Generate data
set.seed(123456789)
modi <- data.frame(matrix(NA, nrow = 1000, ncol = 4))
names(modi) <- c("treat", "sex", "outcome", "neuro")
modi$treat <- rbinom(1000, 1, .5)
modi$sex <- rbinom(1000, 1, .5)
modi$neuro <- rnorm(1000) + modi$sex
modi$outcome <- rnorm(1000) + modi$treat*(2 + 2*modi$sex)
Frequentist models
summary(lm(outcome ~ treat*neuro, data = modi))
##
## Call:
## lm(formula = outcome ~ treat * neuro, data = modi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.567 -0.785 0.020 0.763 3.781
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0390 0.0588 0.66 0.51
## treat 2.8469 0.0820 34.71 < 2e-16 ***
## neuro 0.0554 0.0459 1.21 0.23
## treat:neuro 0.3822 0.0666 5.74 1.3e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.2 on 996 degrees of freedom
## Multiple R-squared: 0.633, Adjusted R-squared: 0.632
## F-statistic: 573 on 3 and 996 DF, p-value: <2e-16
summary(lm(outcome ~ treat*neuro + sex, data = modi))
##
## Call:
## lm(formula = outcome ~ treat * neuro + sex, data = modi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.279 -0.757 -0.014 0.753 3.436
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.3601 0.0644 -5.59 3.0e-08 ***
## treat 2.8762 0.0768 37.44 < 2e-16 ***
## neuro -0.0931 0.0447 -2.08 0.038 *
## sex 0.9198 0.0772 11.91 < 2e-16 ***
## treat:neuro 0.3298 0.0625 5.27 1.6e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.1 on 995 degrees of freedom
## Multiple R-squared: 0.679, Adjusted R-squared: 0.678
## F-statistic: 526 on 4 and 995 DF, p-value: <2e-16
summary(lm(outcome ~ treat*neuro + treat*sex, data = modi))
##
## Call:
## lm(formula = outcome ~ treat * neuro + treat * sex, data = modi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9713 -0.6631 -0.0264 0.6996 3.1587
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0158 0.0664 0.24 0.81
## treat 2.1070 0.0929 22.67 <2e-16 ***
## neuro 0.0467 0.0429 1.09 0.28
## sex 0.0536 0.0982 0.55 0.59
## treat:neuro -0.0239 0.0641 -0.37 0.71
## treat:sex 1.8443 0.1433 12.87 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1 on 994 degrees of freedom
## Multiple R-squared: 0.725, Adjusted R-squared: 0.723
## F-statistic: 523 on 5 and 994 DF, p-value: <2e-16
brms models
library(brms)
modif <- brm(outcome ~ treat*neuro, data = modi, file = "models/effect_modification")
summary(modif)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: outcome ~ treat * neuro
## Data: modi (Number of observations: 1000)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.04 0.06 -0.08 0.15 1.00 4156 3016
## treat 2.85 0.08 2.69 3.00 1.00 3924 2893
## neuro 0.06 0.05 -0.04 0.14 1.00 2725 2867
## treat:neuro 0.38 0.07 0.25 0.51 1.00 2710 2873
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.18 0.03 1.13 1.23 1.00 4376 2917
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
modif2 <- brm(outcome ~ treat*neuro + sex, data = modi, file = "models/effect_modification_sex_adj")
summary(modif2)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: outcome ~ treat * neuro + sex
## Data: modi (Number of observations: 1000)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.36 0.06 -0.49 -0.23 1.00 4302 3091
## treat 2.88 0.08 2.73 3.03 1.00 3700 2982
## neuro -0.09 0.05 -0.18 -0.00 1.00 3192 3138
## sex 0.92 0.08 0.77 1.07 1.00 4922 3218
## treat:neuro 0.33 0.06 0.20 0.45 1.00 2919 3143
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.10 0.02 1.06 1.15 1.00 5863 3050
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
modif3 <- brm(outcome ~ treat*neuro + treat*sex, data = modi, file = "models/effect_modification_sex_int_adj")
summary(modif3)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: outcome ~ treat * neuro + treat * sex
## Data: modi (Number of observations: 1000)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.02 0.07 -0.11 0.14 1.00 2880 3070
## treat 2.10 0.09 1.92 2.29 1.00 2566 2539
## neuro 0.05 0.04 -0.04 0.13 1.00 2608 2990
## sex 0.05 0.10 -0.14 0.25 1.00 2536 2727
## treat:neuro -0.02 0.07 -0.15 0.11 1.00 2688 2688
## treat:sex 1.85 0.14 1.58 2.13 1.00 2229 2800
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.02 0.02 0.98 1.07 1.00 3770 2351
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).