Practical significance t test

dmdv.stan

data {
  real<lower = 0> sd_m;
  real<lower = 0> sd_m_diff;
  real<lower = 0> sd_st;
  real<lower = 0> sd_st_r;
  int<lower = 0, upper = 1> nu_choice;
  int<lower = 0> N;
  vector<lower = 0, upper = 1>[N] x;
  vector[N] y;
}
transformed data {
  real n1 = sum(x);
  real n0 = N - n1;
}
parameters {
  real m0;
  real m_diff;
  real ln_st0;
  real ln_st_ratio;
  real<lower = 0> nu;
}
transformed parameters {
  vector[N] mu = m0 + x * m_diff;
  vector[N] sigma = exp(ln_st0 + x * ln_st_ratio);
}
model {
  m0 ~ cauchy(0, sd_m);
  m_diff ~ normal(0, sd_m_diff);
  ln_st0 ~ student_t(3, 0, sd_st);
  ln_st_ratio ~ normal(0, sd_st_r);
  if (nu_choice == 0) nu ~ gamma(2, 0.1);
  else if (nu_choice == 1) nu ~ exponential(1.0 / 29);

  y ~ student_t(nu + nu_choice, mu, sigma);
}
generated quantities {
  vector[N] log_lik;
  vector[N] yhat;
  real m1 = m0 + m_diff;
  real <lower = 0> st0 = exp(ln_st0);
  real <lower = 0> st1 = exp(ln_st0 + ln_st_ratio);
  real<lower = 0> st_ratio = exp(ln_st_ratio);

  for (i in 1:N) {
    yhat[i] = student_t_rng(nu + nu_choice, mu[i], sigma[i]);
    log_lik[i] = student_t_lpdf(y[i] | nu + nu_choice, mu[i], sigma[i]);
  }
}

file_0.R

library(rstan)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

# Uncomment next line and create stan_scripts folder to recreate this
# saveRDS(stan_model(stanc_ret = stanc(file = "stan_scripts/dmdv.stan")), "stan_scripts/dmdv.rds")
dmdv <- readRDS("stan_scripts/dmdv.rds")

dat <- data.frame(
  y = c(
    c(101, 100, 102, 104, 102, 97, 105, 105, 98, 101, 100, 123, 105, 103, 100, 95, 102, 106,
      109, 102, 82, 102, 100, 102, 102, 101, 102, 102, 103, 103, 97, 97, 103, 101, 97, 104,
      96, 103, 124, 101, 101, 100, 101, 101, 104, 100, 101),
    c(99, 101, 100, 101, 102, 100, 97, 101, 104, 101, 102, 102, 100, 105, 88, 101, 100,
      104, 100, 100, 100, 101, 102, 103, 97, 101, 101, 100, 101, 99, 101, 100, 100,
      101, 100, 99, 101, 100, 102, 99, 100, 99)
  ), x = c(rep(1, 47), rep(0, 42))
)

library(dplyr)
library(ggplot2)
library(scales)
theme_set(theme_classic())

image1 <- dat %>%
  mutate(x = factor(x)) %>%
  ggplot(aes(x, y, col = x)) +
  scale_colour_manual(values = c("#00BFC4", "#F8766D")) +
  geom_jitter(height = 0, shape = 1, width = .3) +
  labs(x = "Group assignment", y = "IQ scores", tag = "A") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        axis.title.y = element_blank()) +
  scale_y_continuous(breaks = c(100, 105, sort(dat$y)[c(1, 2, 3, 87, 89)])) +
  guides(color = FALSE) + coord_flip()
image1

stan.dat <- list(
  sd_m = 5, sd_m_diff = 15 / qnorm(.975), sd_st = 1, sd_st_r = log(3) / qnorm(.975),
  nu_choice = 0, N = length(dat$y), y = dat$y, x = dat$x)

fit.dmdv <- sampling(dmdv, data = stan.dat, seed = 123)
print(fit.dmdv, pars = c("m0", "m1", "m_diff", "st0", "st1", "st_ratio", "nu"),
      probs = c(.025, .1, .5, .9, .975), digits_summary = 3)
mdiff.df <- as.data.frame(fit.dmdv, "m_diff")
effect.tab <- cbind(level = seq(0, 5, .01),
                    efficacy = sapply(seq(0, 5, .01), function (x) mean(mdiff.df$m_diff > x))) %>%
  as.data.frame()
effect.tab %>% ggplot(aes(level, efficacy)) + geom_line() +
  scale_y_continuous(labels = percent_format()) + theme_bw() +
  labs(x = "Practically significant treatment effect", y = "",
       subtitle = "Probability that treatment effect is practically significant")