Data preparation

Prepare dataset to proceed with analysis

setwd("C:/Users/hyejung_jung/Documents/PACt-MD")
baseline_complete <- read_excel("PACt-MD_DataReq78_Susmita_Baseline_375.xlsx") 

#ADS total score
summary(baseline_complete$gt_saa_ads) # 94 NA's, i.e., 281  non-missing values
sd(baseline_complete$gt_saa_ads,na.rm=TRUE)
tblFun <- function(x){
    tbl <- table(x)
    res <- cbind(tbl,round(prop.table(tbl)*100,2))
    colnames(res) <- c('Count','Percentage')
    res
}
tblFun(baseline_complete$gt_saa_ads)
#boxplot(baseline_complete$gt_saa_ads)
#ACB total score
summary(baseline_complete$gt_saa_acb) # 253 NA's, i.e., 122 non-missing values
sd(baseline_complete$gt_saa_acb,na.rm=TRUE)
tblFun(baseline_complete$gt_saa_acb)
#boxplot(baseline_complete$gt_saa_acb)

# we need the dataset below to extract complete (375) baseline global cog composite score to merge to the dataset above
cog_composite_long <- read_csv("PACt-MD_LongitudinalData_375.csv")
#table(cog_composite_long$timepoint)
baseline_complete_long <- subset(cog_composite_long, timepoint == 0)

# merge dependent variables to the complete dataset (375 obs)
baseline_complete <- baseline_complete %>% filter(!is.na(gt_saa_acb))
summary(baseline_complete$gt_saa_acb) 
sd(baseline_complete$gt_saa_acb) #Mean(SD)=1.79(1.34)

baseline_complete$ACB <- as.numeric(baseline_complete$gt_saa_acb >= 3)

yn <- c("Yes","No")
dx1 <- c("MCI","MDD-MCI","MDD+MCI")
dx2 <- c("MCI","MDD-MCI","MDD+MCI","Dementia")
dx3 <- c("Normal","MCI","MDD-MCI","MDD+MCI","Dementia")
ethn <- c("Non-hispanic","Hispanic")
race <- c("White", "African descent", "Asian", "Native North American", "Pacific Islander", "Other", "Mixed")
edu <- c("less than 7th grade", "junior high (9th grade)", "partial high school (10 or 11th grade)", "high school graduate", "partial college (at least one year)", "college education", "graduate degree")
edu1 <- c("High school or less", "High school or less", "High school or less", "High school or less", "Some or complete undergraduate", "Some or complete undergraduate", "Graduate degree")
occu <- c("farm/day laborer", "unskilled/service worker", "machine operator/semiskilled worker", "skilled manual worker/craftsman/police and fire services", "clerical/sales; small farm owner", "technicians/semiprofessional/supervisor/office manager", "small business/farm owner; teacher; low level manager; salaried worker", "mid-level manager/professional; mid-size business owner; military officer", "senior manager/professional; owner/CEO of a large business", "retired")
lih <- c("Low", "Intermediate", "High")
highlow <- c("Low", "High")

basel <- baseline_complete %>% 
  mutate(site = factor(site, levels = 1:5, labels = c("SMH","SBH","UHN","CMH","BYC")),
         rand_condition = factor(rand_condition, levels = c("X","B"), labels = c("Sham","Active")),
         randstrat = factor(provisional_dx, levels = 1:3, labels = c("MCI","w/0 MDD","/w MDD")),
         visit_date_t1_class = as.Date(visit_date_t1_class, format = "%m/%d/%Y"),
ACB = factor(ACB, levels = 0:1, labels = c("Low", "High")),
         provisional_dx = factor(provisional_dx, levels = 1:3, labels = dx1),
         demo_gender = factor(demo_gender, levels = 1:2, labels = c("Male","Female")),
         ethnic = factor(demo_ethnicity, levels = 1:2, labels = ethn),
         race = factor(demo_race, levels = c(1,2,3,4,5,7,8), labels = race),
         educ = factor(sdemo_sbjct_highest_ed, levels = 1:7, labels = edu1),
         sdemo_sbjct_highest_ed = factor(sdemo_sbjct_highest_ed, levels = 1:7, labels = edu),
         sdemo_sbjct_highest_occ = factor(sdemo_sbjct_highest_occ, levels = 1:10, labels = occu),
         #gt_apoediplotype = as.factor(gt_apoediplotype),
         analysis_diagnosis = factor(analysis_diagnosis, levels = 1:4, labels = dx2),
         MDD = factor(MDD, levels = 0:1, labels = c("No","Yes")),
         MCI = factor(MCI, levels = 0:1, labels = c("No","Yes")),
         #vs_bmi = as.numeric(vs_bmi),
        # gt_apoe_nonevs.some = factor(gt_apoe_nonevs.some, levels = 0:1, labels = c("None", "One or Two"))
       )

# Read and process follow-up data
ctyp <- paste(c("c", rep("n", 26)), collapse = "")
reasfin <- c("Complete","Dropout","Progressed","Withdrawn")
vismeth <- c("In-person","Phone","Videoconference")

fup <- read_csv("PACt-MD_LongitudinalData_375.csv", col_types = ctyp) %>% 
  filter(!is.na(timepoint)) %>%
  mutate(reason_final_visit = factor(reason_final_visit, levels = 1:4, labels = reasfin),
         np_completed = factor(np_completed, levels = 1:2, labels = yn),
         diabetes = factor(diabetes, levels = 1:2, labels = yn),
         hypertension = factor(hypertension, levels = 1:2, labels = yn),
         analysis_diagnosis1 = factor(analysis_diagnosis, levels = 0:4, dx3),
         assessment_method = factor(assessment_method, levels = 1:3, labels = vismeth),
         #MCI = factor(MCI, levels = 0:1, labels = yn),
         gt_apoe_nonevs.some = factor(gt_apoe_nonevs.some, levels = 0:1, 
                                      labels = c("None", "One or Two")))
fup <- fup %>% filter(!(record_id == "PAC01_CMH_0110095"))
#fup2 <- filter(fup, gt_apoe_nonevs.some == "None")   

# Build the data frame for the primary analysis (H1 and H3).

visit0 <- filter(fup, timepoint == 0) %>% 
  dplyr::select(!c(timepoint,reason_final_visit,days_since_t1,intervention_sessions_attended))
names(visit0) <- c("record_id", paste(names(visit0)[-1], "_0", sep=""))

fuponly <- filter(fup, timepoint != 0) %>% 
  full_join(visit0)

#summary(basel$Frontal_TGC_2B_TC)
#basel_265 <- basel %>% filter(!is.na(Frontal_TGC_2B_TC)) # N = 261

data_complete <- left_join(dplyr::select(basel, record_id, site, rand_condition,
ACB, MDD, MCI, randstrat,
age = age_enrollment_years, educ), fuponly) %>%
  #filter(randstrat == "MCI") %>%
  #filter(gt_apoe_nonevs.some == "One or Two") %>%
  #filter(gt_apoe_nonevs.some == "None") %>%
  #filter(MDD == "Yes") %>%
  #filter(randstrat == "MDD") %>%
  #filter(MCI == "Yes") %>%
  #filter(TGC_Group_Median == "High") %>%
  #filter(TGC_Group_Median == "Low") %>%
  mutate(months = days_since_t1 / 30.4375,
         year = round(days_since_t1 / 365.25))
#length(unique(data_complete$record_id))

#check <- data_complete %>% filter(!is.na(Frontal_TGC_2B_TC))
#length(unique(check$record_id))

Distribution of original ACB total score:

tblFun(baseline_complete$gt_saa_acb)
##   Count Percentage
## 1    77      63.11
## 2    18      14.75
## 3    16      13.11
## 4     6       4.92
## 6     2       1.64
## 7     3       2.46
baseline_complete %>% 
   summarize(n=n(), mean=mean(gt_saa_acb), sd=sd(gt_saa_acb))
## # A tibble: 1 × 3
##       n  mean    sd
##   <int> <dbl> <dbl>
## 1   122  1.79  1.34

ACB total score was dichotomized as low if < 3, and high if >= 3.The distribution of the binary ACB is as follows:

tblFun(baseline_complete$ACB)
##   Count Percentage
## 0    95      77.87
## 1    27      22.13

Complete case analysis

We performed a complete case analysis including N = 122 who are not missing with ACB values. The analysis employed a linear mixed-effect model where the random effect was subject-based and both random intercept and slope were used in the model. The model adjusted for the randomization stratum and included time since randomization, randomized treatment and the time by treatment interaction.

  • Dependent variables: Global cognition, Verbal memory, Visual memory, Speed of processing, Working memory, Language, Executive function

  • A brief interpretation of each result has included for the model with global cognition. The other models can be interpreted similarly.

#Function for regression model 
regression <- function(dependent, dependent_baseline, moderator, data) {
#Two-way interaction model
prim_model1 <- lmer(as.formula(paste(dependent, "~", dependent_baseline, "+", moderator,
                                     "+ age + educ + randstrat + rand_condition * months + (months | record_id)")),data = data)

pars <- getME(prim_model1, c("theta","fixef"))
prim_model1 <- update(prim_model1, start = pars) 

#Null model regression model 
prim_model1_null <- update(prim_model1, ~.-rand_condition:months - rand_condition)

#Three-way interaction model
sg_model1 <- lmer(as.formula(paste(dependent, "~", dependent_baseline, "+", moderator, 
                                   "+ age + educ + randstrat + rand_condition * months *", moderator, "+ (months | record_id)")),data = data)

pars <- getME(sg_model1, c("theta","fixef"))
sg_model1 <- update(sg_model1, start=pars)
pars <- getME(sg_model1, c("theta","fixef"))
sg_model1 <- update(sg_model1, start=pars)

#LRT
lrt1 <- anova(prim_model1,prim_model1_null)
p_value1 <- lrt1$`Pr(>Chisq)`[2]
lrt2 <- anova(prim_model1,sg_model1)
p_value2 <- lrt2$`Pr(>Chisq)`[2]

#Estimated marginal means
prim1_emmG <- ref_grid(prim_model1, at = list(months = 1:6 * 12))
prim1_emm <- emmeans(prim1_emmG, pairwise ~ rand_condition | months)

P.c <- tidy(prim1_emm$contrasts, conf.int = TRUE) %>% 
  select(Year = months, Contrast = contrast, Estimate = estimate,
         Lower = conf.low, Upper = conf.high, P = p.value) %>% 
  mutate(Year = Year/12)


sg1_emmG <- ref_grid(sg_model1, at = list(months = 1:6 * 12))


return(list(
    prim_model1 = prim_model1,
    prim_model1_null = prim_model1_null,
    sg_model1 = sg_model1,
    lrt1 = lrt1,
    lrt2 = lrt2,
    p_value1 = p_value1,
    p_value2 = p_value2,
    sg1_emmG = sg1_emmG,
    P.c = P.c
))}

Baseline ACB on Global cognition

Table 1 summarizes the results of three regression models: the null model, the two-way interaction model, and the three-way interaction model. At the bottom of the table, it is indicated that N = 119 distinct record IDs (participants) were included in the analysis for global cognition.

result1 <- regression(dependent = "globalcog_composite",
                     dependent_baseline = "globalcog_composite_0",
                     moderator = "ACB",
                     data = data_complete )

#Results of two models
tab_model(result1$prim_model1_null,result1$prim_model1,result1$sg_model1,digits = 3, 
          dv.labels = c("Null model", "Model with two-way interaction", "Model with three-way interaction" ), 
          title = "Table 1. Summary of regression models of baseline ACB on global cognition")
Table 1. Summary of regression models of baseline ACB on global cognition
  Null model Model with two-way interaction Model with three-way interaction
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 1.148 0.387 – 1.909 0.003 1.063 0.292 – 1.834 0.007 1.059 0.284 – 1.834 0.007
globalcog composite 0 0.826 0.738 – 0.913 <0.001 0.828 0.740 – 0.915 <0.001 0.832 0.744 – 0.919 <0.001
ACB [High] 0.017 -0.111 – 0.146 0.790 0.016 -0.112 – 0.144 0.802 -0.003 -0.184 – 0.178 0.974
age -0.015 -0.025 – -0.004 0.007 -0.014 -0.025 – -0.003 0.010 -0.014 -0.025 – -0.003 0.011
educ [Some or complete
undergraduate] 		0.024 -0.124 – 0.172 0.746 0.039 -0.110 – 0.188 0.605 0.039 -0.111 – 0.189 0.611
educ [Graduate degree] -0.034 -0.203 – 0.135 0.694 -0.013 -0.184 – 0.158 0.882 -0.013 -0.185 – 0.159 0.882
randstrat [w/0 MDD] 0.082 -0.052 – 0.216 0.229 0.084 -0.050 – 0.218 0.218 0.085 -0.050 – 0.219 0.216
randstrat [/w MDD] -0.028 -0.172 – 0.117 0.706 -0.028 -0.172 – 0.117 0.707 -0.028 -0.173 – 0.117 0.708
months -0.008 -0.010 – -0.006 <0.001 -0.008 -0.011 – -0.006 <0.001 -0.008 -0.011 – -0.005 <0.001
rand condition [Active] 0.063 -0.045 – 0.171 0.252 0.062 -0.061 – 0.184 0.325
rand condition [Active] ×
months 		0.001 -0.003 – 0.005 0.533 -0.000 -0.005 – 0.004 0.919
ACB [High] × rand
condition [Active] 		0.009 -0.248 – 0.265 0.948
ACB [High] × months -0.001 -0.007 – 0.006 0.844
(ACB [High] × rand
condition [Active]) ×
months 		0.006 -0.003 – 0.016 0.204
Random Effects
σ2 0.06 0.06 0.06
τ00 0.05 record_id 0.05 record_id 0.05 record_id
τ11 0.00 record_id.months 0.00 record_id.months 0.00 record_id.months
ρ01 0.14 record_id 0.13 record_id 0.10 record_id
ICC 0.70 0.70 0.70
N 119 record_id 119 record_id 119 record_id
Observations 523 523 523
Marginal R2 / Conditional R2 0.695 / 0.908 0.695 / 0.908 0.700 / 0.909

Comparision of null model vs two-way interaction modelz

We can see that our result shows nonsignificant p-value (p = 0.32), therefore we can conclude that adding the two-way interaction to the model did not provide significant improved fit over the null model.

kable(result1$lrt1)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1_null 13 299.5816 354.9561 -136.7908 273.5816 NA NA NA
prim_model1 15 301.2969 365.1906 -135.6485 271.2969 2.284648 2 0.3190766

Comparision of two-way interaction model vs three-way interaction model

The result (p-value = 0.43) shows a nonsignificant moderator effect of ACB in the model with global cognition.

kable(result1$lrt2)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1 15 301.2969 365.1906 -135.6485 271.2969 NA NA NA
sg_model1 18 304.5637 381.2362 -134.2819 268.5637 2.733187 3 0.4346169

Estimated marginal means

  • Two-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in two-way interaction model.

kable(result1$P.c, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Contrast Estimate Lower Upper P
1 Sham - Active -0.079 -0.188 0.030 0.153
2 Sham - Active -0.095 -0.224 0.035 0.150
3 Sham - Active -0.110 -0.274 0.053 0.184
4 Sham - Active -0.126 -0.331 0.078 0.224
5 Sham - Active -0.142 -0.391 0.107 0.261
6 Sham - Active -0.158 -0.453 0.138 0.292

The following plot shows the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(globalcog_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(globalcog_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

emmip(result1$prim_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12)),
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      theme(plot.title = element_text(size = 12))

  • Three-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in three-way interaction model.

sg1_emm <- emmeans(result1$sg1_emmG, pairwise ~ rand_condition | ACB * months)

sg1.tbl <- tidy(sg1_emm$contrasts, conf.int = TRUE) %>% 
  select(Year = months, 'Sub-Group' = ACB, Contrast = contrast, 
         Estimate = estimate, Lower = conf.low, Upper = conf.high, P = p.value) %>% 
  mutate(Year = Year/12)

kable(sg1.tbl, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Sub-Group Contrast Estimate Lower Upper P
1 Low Sham - Active -0.059 -0.182 0.064 0.347
1 High Sham - Active -0.142 -0.367 0.082 0.211
2 Low Sham - Active -0.056 -0.202 0.091 0.453
2 High Sham - Active -0.215 -0.477 0.048 0.108
3 Low Sham - Active -0.053 -0.238 0.133 0.574
3 High Sham - Active -0.287 -0.616 0.042 0.087
4 Low Sham - Active -0.050 -0.283 0.183 0.672
4 High Sham - Active -0.359 -0.770 0.052 0.086
5 Low Sham - Active -0.047 -0.331 0.237 0.744
5 High Sham - Active -0.432 -0.932 0.069 0.090
6 Low Sham - Active -0.044 -0.381 0.293 0.796
6 High Sham - Active -0.504 -1.098 0.090 0.096

The following plots show the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(globalcog_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(globalcog_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

pldata <- emmip(result1$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "Low"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = Low") +
      theme(plot.title = element_text(size = 12))

pldata2 <- emmip(result1$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "High"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = High") +
      theme(plot.title = element_text(size = 12))

ggarrange(pldata, pldata2, nrow =1, legend = "bottom")

Baseline ACB on Verbal memory

result2 <- regression(dependent = "verbalmemory_composite",
                     dependent_baseline = "verbalmemory_composite_0",
                     moderator = "ACB",
                     data = data_complete )

#Results of two models
tab_model(result2$prim_model1_null,result2$prim_model1,result2$sg_model1,digits = 3, 
          dv.labels = c("Null model", "Model with two-way interaction", "Model with three-way interaction" ), 
          title = "Table 2. Summary of regression models of baseline ACB on verbal memory")
Table 2. Summary of regression models of baseline ACB on verbal memory
  Null model Model with two-way interaction Model with three-way interaction
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 2.965 1.348 – 4.583 <0.001 2.857 1.213 – 4.501 0.001 2.896 1.254 – 4.538 0.001
verbalmemory composite 0 0.667 0.559 – 0.775 <0.001 0.675 0.566 – 0.783 <0.001 0.673 0.564 – 0.782 <0.001
ACB [High] 0.092 -0.184 – 0.369 0.512 0.099 -0.178 – 0.376 0.483 -0.100 -0.552 – 0.353 0.665
age -0.044 -0.066 – -0.022 <0.001 -0.043 -0.065 – -0.021 <0.001 -0.043 -0.065 – -0.021 <0.001
educ [Some or complete
undergraduate] 		-0.244 -0.570 – 0.082 0.142 -0.217 -0.549 – 0.115 0.200 -0.225 -0.557 – 0.107 0.184
educ [Graduate degree] -0.082 -0.444 – 0.280 0.658 -0.053 -0.423 – 0.317 0.780 -0.060 -0.431 – 0.311 0.750
randstrat [w/0 MDD] 0.491 0.210 – 0.773 0.001 0.489 0.207 – 0.771 0.001 0.499 0.217 – 0.780 0.001
randstrat [/w MDD] 0.118 -0.198 – 0.435 0.462 0.115 -0.202 – 0.432 0.476 0.108 -0.209 – 0.424 0.505
months -0.003 -0.006 – 0.000 0.051 -0.005 -0.010 – -0.001 0.014 -0.007 -0.012 – -0.002 0.008
rand condition [Active] 0.010 -0.264 – 0.284 0.942 -0.003 -0.315 – 0.308 0.983
rand condition [Active] ×
months 		0.005 -0.001 – 0.011 0.132 0.004 -0.003 – 0.011 0.300
ACB [High] × rand
condition [Active] 		0.049 -0.594 – 0.691 0.882
ACB [High] × months 0.005 -0.005 – 0.016 0.316
(ACB [High] × rand
condition [Active]) ×
months 		0.005 -0.010 – 0.020 0.535
Random Effects
σ2 0.29 0.29 0.29
τ00 0.35 record_id 0.36 record_id 0.37 record_id
τ11 0.00 record_id.months 0.00 record_id.months 0.00 record_id.months
ρ01 -0.45 record_id -0.46 record_id -0.49 record_id
ICC 0.52 0.53 0.52
N 117 record_id 117 record_id 117 record_id
Observations 520 520 520
Marginal R2 / Conditional R2 0.589 / 0.804 0.591 / 0.806 0.595 / 0.807

Comparision of null model vs two-way interaction model

The results shows that the higher order models (two-way interaction model and three-way interaction model) did not improve the model fit compared to lower order models (null model and two-way interaction model, respectively).

kable(result2$lrt1)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1_null 13 1081.340 1136.64 -527.6701 1055.340 NA NA NA
prim_model1 15 1081.793 1145.60 -525.8965 1051.793 3.547252 2 0.1697165

Comparision of two-way interaction model vs three-way interaction model

kable(result2$lrt2)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1 15 1081.793 1145.60 -525.8965 1051.793 NA NA NA
sg_model1 18 1083.091 1159.66 -523.5454 1047.091 4.702079 3 0.1949581

Estimated marginal means

  • Two-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in two-way interaction model.

kable(result2$P.c, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Contrast Estimate Lower Upper P
1 Sham - Active -0.070 -0.311 0.172 0.570
2 Sham - Active -0.129 -0.359 0.101 0.269
3 Sham - Active -0.188 -0.433 0.056 0.129
4 Sham - Active -0.248 -0.529 0.033 0.083
5 Sham - Active -0.307 -0.640 0.026 0.070
6 Sham - Active -0.366 -0.760 0.027 0.068

The following plot shows the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(verbalmemory_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(verbalmemory_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

emmip(result2$prim_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12)),
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      theme(plot.title = element_text(size = 12))

  • Three-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in three-way interaction model.

sg1_emm <- emmeans(result2$sg1_emmG, pairwise ~ rand_condition | ACB * months)

sg1.tbl <- tidy(sg1_emm$contrasts, conf.int = TRUE) %>% 
  select(Year = months, 'Sub-Group' = ACB, Contrast = contrast, 
         Estimate = estimate, Lower = conf.low, Upper = conf.high, P = p.value) %>% 
  mutate(Year = Year/12)

kable(sg1.tbl, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Sub-Group Contrast Estimate Lower Upper P
1 Low Sham - Active -0.043 -0.318 0.233 0.759
1 High Sham - Active -0.149 -0.644 0.346 0.553
2 Low Sham - Active -0.089 -0.350 0.172 0.502
2 High Sham - Active -0.252 -0.717 0.212 0.283
3 Low Sham - Active -0.135 -0.411 0.141 0.335
3 High Sham - Active -0.356 -0.845 0.132 0.151
4 Low Sham - Active -0.181 -0.498 0.136 0.260
4 High Sham - Active -0.460 -1.021 0.101 0.107
5 Low Sham - Active -0.227 -0.602 0.148 0.232
5 High Sham - Active -0.564 -1.230 0.103 0.096
6 Low Sham - Active -0.273 -0.717 0.171 0.224
6 High Sham - Active -0.667 -1.459 0.124 0.097

The following plots show the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(verbalmemory_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(verbalmemory_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

pldata <- emmip(result2$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "Low"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = Low") +
      theme(plot.title = element_text(size = 12))

pldata2 <- emmip(result2$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "High"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = High") +
      theme(plot.title = element_text(size = 12))

ggarrange(pldata, pldata2, nrow =1, legend = "bottom")

Baseline ACB on Visual memory

result3 <- regression(dependent = "vsmemory_composite",
                     dependent_baseline = "vsmemory_composite_0",
                     moderator = "ACB",
                     data = data_complete )

#Results of two models
tab_model(result3$prim_model1_null,result3$prim_model1,result3$sg_model1,digits = 3, 
          dv.labels = c("Null model", "Model with two-way interaction", "Model with three-way interaction" ), 
          title = "Table 3. Summary of regression models of baseline ACB on visual memory")
Table 3. Summary of regression models of baseline ACB on visual memory
  Null model Model with two-way interaction Model with three-way interaction
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 1.152 -0.444 – 2.748 0.157 0.946 -0.661 – 2.553 0.248 0.958 -0.656 – 2.572 0.244
vsmemory composite 0 0.475 0.371 – 0.578 <0.001 0.477 0.374 – 0.580 <0.001 0.479 0.376 – 0.583 <0.001
ACB [High] 0.035 -0.235 – 0.306 0.798 0.029 -0.240 – 0.298 0.835 -0.137 -0.553 – 0.278 0.516
age -0.019 -0.041 – 0.002 0.080 -0.018 -0.040 – 0.003 0.096 -0.018 -0.040 – 0.004 0.103
educ [Some or complete
undergraduate] 		0.176 -0.139 – 0.490 0.273 0.217 -0.100 – 0.534 0.180 0.221 -0.098 – 0.539 0.175
educ [Graduate degree] 0.271 -0.071 – 0.614 0.120 0.325 -0.023 – 0.673 0.067 0.332 -0.018 – 0.682 0.063
randstrat [w/0 MDD] 0.462 0.191 – 0.733 0.001 0.467 0.197 – 0.737 0.001 0.466 0.195 – 0.737 0.001
randstrat [/w MDD] -0.046 -0.358 – 0.265 0.770 -0.041 -0.351 – 0.269 0.797 -0.046 -0.358 – 0.265 0.770
months -0.008 -0.012 – -0.004 <0.001 -0.007 -0.012 – -0.001 0.012 -0.008 -0.014 – -0.002 0.011
rand condition [Active] 0.201 -0.048 – 0.451 0.113 0.161 -0.122 – 0.443 0.264
rand condition [Active] ×
months 		-0.002 -0.010 – 0.005 0.582 -0.003 -0.011 – 0.006 0.553
ACB [High] × rand
condition [Active] 		0.185 -0.402 – 0.771 0.536
ACB [High] × months 0.004 -0.009 – 0.018 0.506
(ACB [High] × rand
condition [Active]) ×
months 		0.001 -0.017 – 0.020 0.881
Random Effects
σ2 0.48 0.48 0.48
τ00 0.16 record_id 0.15 record_id 0.16 record_id
τ11 0.00 record_id.months 0.00 record_id.months 0.00 record_id.months
ρ01 1.00 record_id 1.00 record_id 1.00 record_id
ICC 0.35    
N 118 record_id 118 record_id 118 record_id
Observations 473 473 473
Marginal R2 / Conditional R2 0.412 / 0.616 0.520 / NA 0.522 / NA

Comparision of null model vs two-way interaction model

The results shows that the higher order models (two-way interaction model and three-way interaction model) did not improve the model fit compared to lower order models (null model and two-way interaction model, respectively).

kable(result3$lrt1)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1_null 13 1134.924 1188.992 -554.4619 1108.924 NA NA NA
prim_model1 15 1136.204 1198.590 -553.1019 1106.204 2.719894 2 0.2566744

Comparision of two-way interaction model vs three-way interaction model

kable(result3$lrt2)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1 15 1136.204 1198.590 -553.1019 1106.204 NA NA NA
sg_model1 18 1140.254 1215.118 -552.1269 1104.254 1.95006 3 0.5828363

Estimated marginal means

  • Two-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in two-way interaction model.

kable(result3$P.c, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Contrast Estimate Lower Upper P
1 Sham - Active -0.176 -0.401 0.049 0.124
2 Sham - Active -0.151 -0.385 0.084 0.205
3 Sham - Active -0.125 -0.402 0.151 0.370
4 Sham - Active -0.100 -0.440 0.240 0.559
5 Sham - Active -0.075 -0.489 0.340 0.720
6 Sham - Active -0.049 -0.545 0.446 0.843

The following plot shows the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(vsmemory_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(vsmemory_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

emmip(result3$prim_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12)),
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      theme(plot.title = element_text(size = 12))

  • Three-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in three-way interaction model.

sg1_emm <- emmeans(result3$sg1_emmG, pairwise ~ rand_condition | ACB * months)

sg1.tbl <- tidy(sg1_emm$contrasts, conf.int = TRUE) %>% 
  select(Year = months, 'Sub-Group' = ACB, Contrast = contrast, 
         Estimate = estimate, Lower = conf.low, Upper = conf.high, P = p.value) %>% 
  mutate(Year = Year/12)

kable(sg1.tbl, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Sub-Group Contrast Estimate Lower Upper P
1 Low Sham - Active -0.130 -0.386 0.126 0.317
1 High Sham - Active -0.332 -0.793 0.130 0.157
2 Low Sham - Active -0.099 -0.366 0.168 0.463
2 High Sham - Active -0.318 -0.798 0.162 0.192
3 Low Sham - Active -0.068 -0.383 0.246 0.667
3 High Sham - Active -0.304 -0.878 0.270 0.296
4 Low Sham - Active -0.038 -0.424 0.349 0.847
4 High Sham - Active -0.290 -1.005 0.425 0.423
5 Low Sham - Active -0.007 -0.477 0.464 0.977
5 High Sham - Active -0.276 -1.157 0.605 0.535
6 Low Sham - Active 0.024 -0.538 0.586 0.933
6 High Sham - Active -0.262 -1.322 0.797 0.624

The following plots show the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(vsmemory_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(vsmemory_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

pldata <- emmip(result3$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "Low"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = Low") +
      theme(plot.title = element_text(size = 12))

pldata2 <- emmip(result3$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "High"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = High") +
      theme(plot.title = element_text(size = 12))

ggarrange(pldata, pldata2, nrow =1, legend = "bottom")

Baseline ACB on Speed of processing

result4 <- regression(dependent = "procspeed_composite",
                     dependent_baseline = "procspeed_composite_0",
                     moderator = "ACB",
                     data = data_complete )

#Results of two models
tab_model(result4$prim_model1_null,result4$prim_model1,result4$sg_model1,digits = 3, 
          dv.labels = c("Null model", "Model with two-way interaction", "Model with three-way interaction" ), 
          title = "Table 4. Summary of regression models of baseline ACB on speed of processing")
Table 4. Summary of regression models of baseline ACB on speed of processing
  Null model Model with two-way interaction Model with three-way interaction
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 2.131 0.743 – 3.519 0.003 2.160 0.738 – 3.581 0.003 2.097 0.686 – 3.507 0.004
procspeed composite 0 0.740 0.643 – 0.836 <0.001 0.739 0.641 – 0.837 <0.001 0.758 0.659 – 0.857 <0.001
ACB [High] 0.173 -0.062 – 0.407 0.149 0.173 -0.063 – 0.409 0.150 -0.012 -0.355 – 0.332 0.947
age -0.034 -0.053 – -0.015 0.001 -0.034 -0.053 – -0.015 0.001 -0.032 -0.052 – -0.013 0.001
educ [Some or complete
undergraduate] 		0.227 -0.047 – 0.500 0.104 0.226 -0.051 – 0.504 0.110 0.231 -0.045 – 0.506 0.100
educ [Graduate degree] 0.165 -0.137 – 0.467 0.284 0.163 -0.145 – 0.472 0.298 0.172 -0.134 – 0.478 0.270
randstrat [w/0 MDD] 0.018 -0.211 – 0.247 0.878 0.017 -0.213 – 0.247 0.883 0.013 -0.215 – 0.241 0.912
randstrat [/w MDD] -0.149 -0.416 – 0.119 0.276 -0.150 -0.418 – 0.119 0.275 -0.152 -0.418 – 0.115 0.264
months -0.010 -0.013 – -0.007 <0.001 -0.011 -0.015 – -0.007 <0.001 -0.011 -0.016 – -0.006 <0.001
rand condition [Active] -0.039 -0.247 – 0.169 0.712 -0.108 -0.341 – 0.124 0.360
rand condition [Active] ×
months 		0.003 -0.004 – 0.009 0.416 0.001 -0.006 – 0.008 0.788
ACB [High] × rand
condition [Active] 		0.325 -0.167 – 0.818 0.195
ACB [High] × months -0.001 -0.011 – 0.010 0.905
(ACB [High] × rand
condition [Active]) ×
months 		0.007 -0.009 – 0.022 0.402
Random Effects
σ2 0.28 0.27 0.27
τ00 0.14 record_id 0.14 record_id 0.14 record_id
τ11 0.00 record_id.months 0.00 record_id.months 0.00 record_id.months
ρ01 0.35 record_id 0.34 record_id 0.27 record_id
ICC 0.51 0.51 0.51
N 118 record_id 118 record_id 118 record_id
Observations 533 533 533
Marginal R2 / Conditional R2 0.578 / 0.794 0.575 / 0.794 0.583 / 0.797

Comparision of null model vs two-way interaction model

The higher order models did not improve the model fit compared to lower order models.

kable(result4$lrt1)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1_null 13 1048.440 1104.061 -511.2199 1022.440 NA NA NA
prim_model1 15 1051.756 1115.934 -510.8781 1021.756 0.683536 2 0.710513

Comparision of two-way interaction model vs three-way interaction model

kable(result4$lrt2)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1 15 1051.756 1115.934 -510.8781 1021.756 NA NA NA
sg_model1 18 1053.445 1130.458 -508.7223 1017.445 4.31163 3 0.2297206

Estimated marginal means

  • Two-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in two-way interaction model.

kable(result4$P.c, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Contrast Estimate Lower Upper P
1 Sham - Active 0.007 -0.191 0.205 0.944
2 Sham - Active -0.025 -0.241 0.191 0.819
3 Sham - Active -0.057 -0.315 0.201 0.662
4 Sham - Active -0.089 -0.403 0.225 0.575
5 Sham - Active -0.121 -0.500 0.258 0.527
6 Sham - Active -0.153 -0.601 0.295 0.499

The following plot shows the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(procspeed_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(procspeed_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

emmip(result4$prim_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12)),
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      theme(plot.title = element_text(size = 12))

  • Three-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in three-way interaction model.

sg1_emm <- emmeans(result4$sg1_emmG, pairwise ~ rand_condition | ACB * months)

sg1.tbl <- tidy(sg1_emm$contrasts, conf.int = TRUE) %>% 
  select(Year = months, 'Sub-Group' = ACB, Contrast = contrast, 
         Estimate = estimate, Lower = conf.low, Upper = conf.high, P = p.value) %>% 
  mutate(Year = Year/12)

kable(sg1.tbl, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Sub-Group Contrast Estimate Lower Upper P
1 Low Sham - Active 0.096 -0.124 0.316 0.388
1 High Sham - Active -0.309 -0.720 0.102 0.139
2 Low Sham - Active 0.084 -0.157 0.325 0.492
2 High Sham - Active -0.401 -0.843 0.042 0.076
3 Low Sham - Active 0.072 -0.219 0.363 0.626
3 High Sham - Active -0.493 -1.020 0.035 0.067
4 Low Sham - Active 0.060 -0.298 0.417 0.742
4 High Sham - Active -0.585 -1.230 0.061 0.075
5 Low Sham - Active 0.047 -0.386 0.480 0.829
5 High Sham - Active -0.676 -1.458 0.105 0.089
6 Low Sham - Active 0.035 -0.479 0.549 0.892
6 High Sham - Active -0.768 -1.695 0.159 0.103

The following plots show the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(procspeed_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(procspeed_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

pldata <- emmip(result4$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "Low"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = Low") +
      theme(plot.title = element_text(size = 12))

pldata2 <- emmip(result4$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "High"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = High") +
      theme(plot.title = element_text(size = 12))

ggarrange(pldata, pldata2, nrow =1, legend = "bottom")

Baseline ACB on Working memory

result5 <- regression(dependent = "wm_composite",
                     dependent_baseline = "wm_composite_0",
                     moderator = "ACB",
                     data = data_complete )

#Results of two models
tab_model(result5$prim_model1_null,result5$prim_model1,result5$sg_model1,digits = 3, 
          dv.labels = c("Null model", "Model with two-way interaction", "Model with three-way interaction" ), 
          title = "Table 5. Summary of regression models of baseline ACB on working memory")
Table 5. Summary of regression models of baseline ACB on working memory
  Null model Model with two-way interaction Model with three-way interaction
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.810 -0.298 – 1.918 0.151 0.849 -0.276 – 1.974 0.139 0.828 -0.302 – 1.957 0.150
wm composite 0 0.686 0.582 – 0.790 <0.001 0.692 0.585 – 0.798 <0.001 0.692 0.585 – 0.800 <0.001
ACB [High] 0.033 -0.157 – 0.223 0.730 0.034 -0.157 – 0.224 0.729 0.099 -0.185 – 0.382 0.495
age -0.009 -0.024 – 0.006 0.248 -0.009 -0.024 – 0.006 0.246 -0.009 -0.024 – 0.006 0.250
educ [Some or complete
undergraduate] 		0.058 -0.159 – 0.275 0.602 0.046 -0.174 – 0.267 0.680 0.049 -0.172 – 0.271 0.662
educ [Graduate degree] -0.023 -0.268 – 0.223 0.856 -0.041 -0.294 – 0.213 0.753 -0.033 -0.288 – 0.221 0.798
randstrat [w/0 MDD] -0.064 -0.257 – 0.129 0.512 -0.070 -0.265 – 0.124 0.477 -0.068 -0.263 – 0.128 0.496
randstrat [/w MDD] 0.124 -0.090 – 0.339 0.255 0.125 -0.091 – 0.340 0.255 0.127 -0.089 – 0.343 0.248
months -0.003 -0.006 – -0.001 0.008 -0.003 -0.006 – 0.001 0.115 -0.001 -0.005 – 0.003 0.605
rand condition [Active] -0.034 -0.206 – 0.138 0.695 -0.032 -0.229 – 0.164 0.745
rand condition [Active] ×
months 		-0.001 -0.006 – 0.004 0.688 -0.002 -0.008 – 0.003 0.392
ACB [High] × rand
condition [Active] 		-0.004 -0.408 – 0.400 0.983
ACB [High] × months -0.007 -0.016 – 0.001 0.074
(ACB [High] × rand
condition [Active]) ×
months 		0.006 -0.006 – 0.018 0.359
Random Effects
σ2 0.15 0.15 0.15
τ00 0.11 record_id 0.11 record_id 0.11 record_id
τ11 0.00 record_id.months 0.00 record_id.months 0.00 record_id.months
ρ01 1.00 record_id 1.00 record_id 1.00 record_id
ICC 0.48    
N 119 record_id 119 record_id 119 record_id
Observations 408 408 408
Marginal R2 / Conditional R2 0.542 / 0.763 0.697 / NA 0.697 / NA

Comparision of null model vs two-way interaction model

The higher order models did not improve the model fit compared to lower order models.

kable(result5$lrt1)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1_null 13 555.2433 607.3898 -264.6216 529.2433 NA NA NA
prim_model1 15 558.7236 618.8926 -264.3618 528.7236 0.5197011 2 0.7711668

Comparision of two-way interaction model vs three-way interaction model

kable(result5$lrt2)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1 15 558.7236 618.8926 -264.3618 528.7236 NA NA NA
sg_model1 18 561.1184 633.3212 -262.5592 525.1184 3.605188 3 0.3073737

Estimated marginal means

  • Two-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in two-way interaction model.

kable(result5$P.c, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Contrast Estimate Lower Upper P
1 Sham - Active 0.046 -0.115 0.207 0.569
2 Sham - Active 0.058 -0.112 0.229 0.498
3 Sham - Active 0.071 -0.129 0.270 0.485
4 Sham - Active 0.083 -0.159 0.324 0.498
5 Sham - Active 0.095 -0.196 0.385 0.518
6 Sham - Active 0.107 -0.237 0.451 0.538

The following plot shows the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(wm_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(wm_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

emmip(result5$prim_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12)),
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      theme(plot.title = element_text(size = 12))

  • Three-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in three-way interaction model.

sg1_emm <- emmeans(result5$sg1_emmG, pairwise ~ rand_condition | ACB * months)

sg1.tbl <- tidy(sg1_emm$contrasts, conf.int = TRUE) %>% 
  select(Year = months, 'Sub-Group' = ACB, Contrast = contrast, 
         Estimate = estimate, Lower = conf.low, Upper = conf.high, P = p.value) %>% 
  mutate(Year = Year/12)

kable(sg1.tbl, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Sub-Group Contrast Estimate Lower Upper P
1 Low Sham - Active 0.061 -0.122 0.245 0.509
1 High Sham - Active -0.002 -0.329 0.326 0.992
2 Low Sham - Active 0.090 -0.103 0.284 0.357
2 High Sham - Active -0.040 -0.385 0.305 0.819
3 Low Sham - Active 0.119 -0.107 0.345 0.297
3 High Sham - Active -0.078 -0.485 0.328 0.703
4 Low Sham - Active 0.148 -0.124 0.420 0.282
4 High Sham - Active -0.117 -0.613 0.379 0.641
5 Low Sham - Active 0.177 -0.150 0.504 0.284
5 High Sham - Active -0.155 -0.756 0.446 0.610
6 Low Sham - Active 0.206 -0.181 0.593 0.291
6 High Sham - Active -0.194 -0.908 0.521 0.592

The following plots show the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(wm_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(wm_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

pldata <- emmip(result5$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "Low"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = Low") +
      theme(plot.title = element_text(size = 12))

pldata2 <- emmip(result5$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "High"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = High") +
      theme(plot.title = element_text(size = 12))

ggarrange(pldata, pldata2, nrow =1, legend = "bottom")

Baseline ACB on Language

result6 <- regression(dependent = "language_composite",
                     dependent_baseline = "language_composite_0",
                     moderator = "ACB",
                     data = data_complete )

#Results of two models
tab_model(result6$prim_model1_null,result6$prim_model1,result6$sg_model1,digits = 3, 
          dv.labels = c("Null model", "Model with two-way interaction", "Model with three-way interaction" ), 
          title = "Table 6. Summary of regression models of baseline ACB on language")
Table 6. Summary of regression models of baseline ACB on language
  Null model Model with two-way interaction Model with three-way interaction
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 0.291 -0.942 – 1.524 0.643 0.311 -0.949 – 1.571 0.628 0.287 -0.976 – 1.549 0.656
language composite 0 0.822 0.728 – 0.917 <0.001 0.823 0.728 – 0.918 <0.001 0.829 0.732 – 0.926 <0.001
ACB [High] 0.022 -0.188 – 0.232 0.835 0.023 -0.188 – 0.234 0.831 0.180 -0.143 – 0.502 0.274
age -0.005 -0.021 – 0.012 0.590 -0.005 -0.022 – 0.012 0.602 -0.004 -0.021 – 0.013 0.609
educ [Some or complete
undergraduate] 		-0.164 -0.409 – 0.082 0.191 -0.161 -0.411 – 0.088 0.205 -0.171 -0.422 – 0.080 0.181
educ [Graduate degree] -0.240 -0.519 – 0.040 0.093 -0.240 -0.526 – 0.046 0.100 -0.255 -0.543 – 0.034 0.083
randstrat [w/0 MDD] 0.252 0.041 – 0.464 0.019 0.250 0.037 – 0.463 0.021 0.245 0.032 – 0.458 0.025
randstrat [/w MDD] -0.043 -0.280 – 0.195 0.724 -0.044 -0.284 – 0.195 0.716 -0.042 -0.282 – 0.198 0.732
months -0.005 -0.007 – -0.002 0.002 -0.006 -0.010 – -0.002 0.002 -0.005 -0.010 – -0.001 0.022
rand condition [Active] -0.059 -0.254 – 0.135 0.549 -0.001 -0.222 – 0.219 0.990
rand condition [Active] ×
months 		0.004 -0.002 – 0.010 0.204 0.003 -0.004 – 0.009 0.435
ACB [High] × rand
condition [Active] 		-0.262 -0.723 – 0.199 0.265
ACB [High] × months -0.004 -0.014 – 0.005 0.361
(ACB [High] × rand
condition [Active]) ×
months 		0.005 -0.009 – 0.019 0.486
Random Effects
σ2 0.14 0.14 0.14
τ00 0.19 record_id 0.19 record_id 0.19 record_id
τ11 0.00 record_id.months 0.00 record_id.months 0.00 record_id.months
ρ01 -0.28 record_id -0.26 record_id -0.26 record_id
ICC 0.63 0.64 0.64
N 119 record_id 119 record_id 119 record_id
Observations 539 539 539
Marginal R2 / Conditional R2 0.659 / 0.875 0.656 / 0.875 0.654 / 0.876

Comparision of null model vs two-way interaction model

The higher order models did not improve the model fit compared to lower order models.

kable(result6$lrt1)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1_null 13 800.7450 856.5113 -387.3725 774.7450 NA NA NA
prim_model1 15 803.1465 867.4923 -386.5733 773.1465 1.598489 2 0.4496685

Comparision of two-way interaction model vs three-way interaction model

kable(result6$lrt2)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1 15 803.1465 867.4923 -386.5733 773.1465 NA NA NA
sg_model1 18 807.4325 884.6474 -385.7163 771.4325 1.714002 3 0.6338253

Estimated marginal means

  • Two-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in two-way interaction model.

kable(result6$P.c, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Contrast Estimate Lower Upper P
1 Sham - Active 0.014 -0.163 0.191 0.873
2 Sham - Active -0.031 -0.215 0.154 0.742
3 Sham - Active -0.076 -0.292 0.141 0.489
4 Sham - Active -0.121 -0.384 0.143 0.365
5 Sham - Active -0.166 -0.485 0.154 0.306
6 Sham - Active -0.211 -0.591 0.170 0.274

The following plot shows the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(language_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(language_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

emmip(result6$prim_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12)),
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      theme(plot.title = element_text(size = 12))

  • Three-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in three-way interaction model.

sg1_emm <- emmeans(result6$sg1_emmG, pairwise ~ rand_condition | ACB * months)

sg1.tbl <- tidy(sg1_emm$contrasts, conf.int = TRUE) %>% 
  select(Year = months, 'Sub-Group' = ACB, Contrast = contrast, 
         Estimate = estimate, Lower = conf.low, Upper = conf.high, P = p.value) %>% 
  mutate(Year = Year/12)

kable(sg1.tbl, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Sub-Group Contrast Estimate Lower Upper P
1 Low Sham - Active -0.031 -0.231 0.170 0.762
1 High Sham - Active 0.172 -0.196 0.541 0.356
2 Low Sham - Active -0.063 -0.274 0.148 0.557
2 High Sham - Active 0.081 -0.300 0.462 0.673
3 Low Sham - Active -0.095 -0.344 0.155 0.452
3 High Sham - Active -0.010 -0.454 0.435 0.966
4 Low Sham - Active -0.127 -0.432 0.179 0.412
4 High Sham - Active -0.101 -0.642 0.441 0.714
5 Low Sham - Active -0.159 -0.530 0.213 0.398
5 High Sham - Active -0.191 -0.849 0.466 0.565
6 Low Sham - Active -0.191 -0.634 0.252 0.394
6 High Sham - Active -0.282 -1.067 0.502 0.477

The following plots show the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(language_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(language_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

pldata <- emmip(result6$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "Low"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = Low") +
      theme(plot.title = element_text(size = 12))

pldata2 <- emmip(result6$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "High"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = High") +
      theme(plot.title = element_text(size = 12))

ggarrange(pldata, pldata2, nrow =1, legend = "bottom")

Baseline ACB on Executive function

result7 <- regression(dependent = "ef_composite",
                     dependent_baseline = "ef_composite_0",
                     moderator = "ACB",
                     data = data_complete )

#Results of two models
tab_model(result7$prim_model1_null,result7$prim_model1,result7$sg_model1,digits = 3, 
          dv.labels = c("Null model", "Model with two-way interaction", "Model with three-way interaction" ), 
          title = "Table 7. Summary of regression models of baseline ACB on executive function")
Table 7. Summary of regression models of baseline ACB on executive function
  Null model Model with two-way interaction Model with three-way interaction
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 1.140 -0.541 – 2.821 0.183 0.807 -0.885 – 2.499 0.349 0.817 -0.883 – 2.516 0.345
ef composite 0 0.615 0.467 – 0.764 <0.001 0.625 0.477 – 0.772 <0.001 0.625 0.476 – 0.773 <0.001
ACB [High] -0.015 -0.298 – 0.267 0.916 -0.032 -0.312 – 0.247 0.820 0.106 -0.311 – 0.523 0.619
age -0.018 -0.041 – 0.005 0.131 -0.016 -0.039 – 0.007 0.176 -0.016 -0.039 – 0.007 0.166
educ [Some or complete
undergraduate] 		-0.064 -0.397 – 0.268 0.704 -0.004 -0.336 – 0.328 0.981 -0.013 -0.347 – 0.321 0.940
educ [Graduate degree] -0.063 -0.428 – 0.303 0.737 0.010 -0.358 – 0.378 0.957 -0.001 -0.371 – 0.369 0.995
randstrat [w/0 MDD] 0.296 0.020 – 0.571 0.035 0.295 0.022 – 0.568 0.034 0.294 0.020 – 0.567 0.036
randstrat [/w MDD] -0.051 -0.373 – 0.272 0.758 -0.034 -0.354 – 0.285 0.832 -0.032 -0.352 – 0.289 0.846
months -0.008 -0.011 – -0.005 <0.001 -0.006 -0.010 – -0.001 0.012 -0.005 -0.010 – -0.000 0.050
rand condition [Active] 0.312 0.059 – 0.565 0.016 0.405 0.119 – 0.690 0.006
rand condition [Active] ×
months 		-0.005 -0.012 – 0.001 0.112 -0.009 -0.016 – -0.001 0.020
ACB [High] × rand
condition [Active] 		-0.422 -1.014 – 0.171 0.163
ACB [High] × months -0.003 -0.014 – 0.007 0.534
(ACB [High] × rand
condition [Active]) ×
months 		0.016 -0.000 – 0.031 0.054
Random Effects
σ2 0.28 0.28 0.28
τ00 0.31 record_id 0.29 record_id 0.29 record_id
τ11 0.00 record_id.months 0.00 record_id.months 0.00 record_id.months
ρ01 -0.15 record_id -0.06 record_id -0.05 record_id
ICC 0.56 0.56 0.57
N 118 record_id 118 record_id 118 record_id
Observations 521 521 521
Marginal R2 / Conditional R2 0.358 / 0.721 0.364 / 0.723 0.369 / 0.729

Comparision of null model vs two-way interaction model

The two-way interaction model shows the significantly improved model fit compared to the null model (p=0.04)

kable(result7$lrt1)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1_null 13 1081.460 1136.785 -527.7302 1055.460 NA NA NA
prim_model1 15 1079.018 1142.854 -524.5088 1049.018 6.442938 2 0.0398964

Comparision of two-way interaction model vs three-way interaction model

The moderator effect of ACB was not significant.

kable(result7$lrt2)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
prim_model1 15 1079.018 1142.854 -524.5088 1049.018 NA NA NA
sg_model1 18 1079.767 1156.371 -521.8837 1043.767 5.250125 3 0.1543717

Estimated marginal means

  • Two-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in two-way interaction model.

kable(result7$P.c, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Contrast Estimate Lower Upper P
1 Sham - Active -0.249 -0.483 -0.014 0.038
2 Sham - Active -0.185 -0.425 0.056 0.131
3 Sham - Active -0.121 -0.392 0.150 0.378
4 Sham - Active -0.057 -0.377 0.262 0.723
5 Sham - Active 0.006 -0.372 0.385 0.973
6 Sham - Active 0.070 -0.375 0.515 0.754

The following plot shows the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(ef_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(ef_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

emmip(result7$prim_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12)),
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      theme(plot.title = element_text(size = 12))

  • Three-way interaction model

The next table estimates the treatment differences at each year (1 to 6) in three-way interaction model.

sg1_emm <- emmeans(result7$sg1_emmG, pairwise ~ rand_condition | ACB * months)

sg1.tbl <- tidy(sg1_emm$contrasts, conf.int = TRUE) %>% 
  select(Year = months, 'Sub-Group' = ACB, Contrast = contrast, 
         Estimate = estimate, Lower = conf.low, Upper = conf.high, P = p.value) %>% 
  mutate(Year = Year/12)

kable(sg1.tbl, booktabs = TRUE, digits = 3) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"))
Year Sub-Group Contrast Estimate Lower Upper P
1 Low Sham - Active -0.299 -0.565 -0.033 0.028
1 High Sham - Active -0.064 -0.547 0.420 0.795
2 Low Sham - Active -0.193 -0.467 0.081 0.165
2 High Sham - Active -0.144 -0.637 0.349 0.564
3 Low Sham - Active -0.087 -0.397 0.223 0.579
3 High Sham - Active -0.224 -0.780 0.332 0.425
4 Low Sham - Active 0.019 -0.346 0.384 0.918
4 High Sham - Active -0.305 -0.963 0.353 0.361
5 Low Sham - Active 0.125 -0.308 0.558 0.568
5 High Sham - Active -0.385 -1.169 0.399 0.332
6 Low Sham - Active 0.231 -0.278 0.740 0.369
6 High Sham - Active -0.465 -1.389 0.458 0.320

The following plots show the estimated response at each year with 95% confidence intervals.

time0 <- filter(fup, timepoint == 0) %>% 
  right_join(select(basel, record_id, site, rand_condition)) %>%
  filter(!is.na(ef_composite)) %>%
  group_by(rand_condition) %>% 
  summarise(yvar = mean(ef_composite)) %>% 
  mutate(months = 0, SE = NA, df = NA, LCL = NA, UCL = NA, tvar = rand_condition, xvar = 0) %>% 
  select(rand_condition, months, yvar, SE, df, LCL, UCL, tvar, xvar) %>% 
  as.data.frame()

pldata <- emmip(result7$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "Low"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + 
      theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = Low") +
      theme(plot.title = element_text(size = 12))

pldata2 <- emmip(result7$sg_model1, rand_condition ~ months, 
      at = list(months = c(2,1:6 * 12), ACB = "High"), 
      type = "response", CIs = TRUE)+ geom_point(aes(x=months,y=yvar,colour=rand_condition), time0, size = 3, shape="circle") + xlab("Year") + 
      scale_x_continuous(breaks = c(0, 1:6 * 12), labels = c(0:6)) + 
      theme_classic() + theme(text = element_text(size = 12), legend.title = element_blank()) + 
      ggtitle("ACB = High") +
      theme(plot.title = element_text(size = 12))

ggarrange(pldata, pldata2, nrow =1, legend = "bottom")