rm(list=ls())


# Load libraries ---------------------------

library(haven)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(geeM)
## Loading required package: Matrix
library(magrittr)


# Set wd for generated output ---------------------------

setwd("/Volumes/caas/MERITS/RESOURCES/Data Processing/Data Analysis/Outcome Analyses 2016/Manuscript Revisions_2020/Aditya Khanna Analyses 2022/RESOURCES-AK-R")


# Load data ---------------------------

load(file="gee-formulae.RData")


# Specify corstr = ar1 ---------------------------

corstr = "ar1"


# DV = AVCIG ---------------------------

#relevel MI and CBT
AVCIG_nobl_dt$MI_ref1 <- AVCIG_nobl_dt$MI %>% relevel("1")
AVCIG_nobl_dt$CBT_ref1 <- AVCIG_nobl_dt$CBT %>% relevel("1")
AVCIG_nobl_dt$TIME_ref1 <- AVCIG_nobl_dt$TIME %>% relevel("1")

# ignore missing values
## create an indicator column for any missing values 
AVCIG_nobl_dt$any_na <- AVCIG_nobl_dt %>% apply(1, function(x){any(is.na(x))}) 

## left join the dataset by the id column using the group_by() function
AVCIG_nobl_dt %<>% left_join(AVCIG_nobl_dt %>% group_by(id) %>% 
                               summarise(any_na2 = any(any_na)), by="id")


## filter out rows with any NAs
AVCIG_nobl_dt %<>% filter(any_na2 != T)

# modeling
## main model
AVCIG_nobl_main <- geem(formula = AVCIG.nobl.main,
                       data = AVCIG_nobl_dt,
                       id=id,
                       family = MASS::negative.binomial(1),
                       corstr = corstr
)


## print main model summary
summary(AVCIG_nobl_main)
##             Estimates Model SE Robust SE    wald        p
## (Intercept)  1.732000 0.108400  0.117200 14.7800 0.000000
## cAGE        -0.027910 0.053040  0.060770 -0.4592 0.646100
## cBLCGSMD     0.049730 0.006412  0.005501  9.0400 0.000000
## MI1          0.081820 0.117900  0.116900  0.6998 0.484100
## CBT1        -0.001529 0.118500  0.116700 -0.0131 0.989500
## TIME1        0.123500 0.045150  0.045430  2.7170 0.006579
## 
##  Estimated Correlation Parameter:  0.7434 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8339 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  246    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  492
## interaction model
AVCIG_nobl_interaction <- geem(formula = AVCIG.nobl.interaction,
                               #mi=0, cbt=0, time=0
                        data = AVCIG_nobl_dt,
                        id=id,
                        family = MASS::negative.binomial(1),
                        corstr = corstr
)

AVCIG_nobl_interaction_TIMEref1 <- geem(formula = 
                                 AVCIG ~ cAGE + cBLCGSMD + MI + CBT + TIME_ref1 + MI * CBT + MI * TIME_ref1 + 
                                 CBT * TIME_ref1 + MI * CBT * TIME_ref1,
                               #mi=0, cbt=0, time=1
                               data = AVCIG_nobl_dt,
                               id=id,
                               family = MASS::negative.binomial(1),
                               corstr = corstr
)


AVCIG_nobl_interaction_MIref1 <- 
  geem(formula = AVCIG ~ cAGE + cBLCGSMD + MI_ref1 + CBT + TIME + MI_ref1 * CBT + MI_ref1 * TIME + 
         CBT * TIME + MI_ref1 * CBT * TIME,
       #ref: Mi=1, Cbt=0, Time=0
       data = AVCIG_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

AVCIG_nobl_interaction_MIref1_TIMEref1 <- 
  geem(formula = AVCIG ~ cAGE + cBLCGSMD + MI_ref1 + CBT + TIME_ref1 + 
         MI_ref1 * CBT + MI_ref1 * TIME_ref1 + 
         CBT * TIME_ref1 + MI_ref1 * CBT * TIME_ref1,
       #ref: Mi=1, Cbt=0, Time=1
       data = AVCIG_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

AVCIG_nobl_interaction_CBTref1 <- 
  geem(formula = AVCIG ~ cAGE + cBLCGSMD + MI + CBT_ref1 + TIME + MI * CBT_ref1 + MI * 
         TIME + CBT_ref1 * TIME + MI * CBT_ref1 * TIME,
       #ref: Mi=0, Cbt=1, Time = 0
       data = AVCIG_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

AVCIG_nobl_interaction_CBTref1_TIMEref1 <- 
  geem(formula = AVCIG ~ cAGE + cBLCGSMD + MI + CBT_ref1 + TIME_ref1 + 
         MI * CBT_ref1 + MI * TIME_ref1 + CBT_ref1 * TIME_ref1 + MI * CBT_ref1 * TIME_ref1,
       #ref: Mi=0, Cbt=1, Time = 1
       data = AVCIG_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

AVCIG_nobl_interaction_MIref1_CBTref1 <- 
  geem(formula = AVCIG ~ cAGE + cBLCGSMD + MI_ref1 + CBT_ref1 + TIME + MI_ref1 * CBT_ref1 + MI_ref1 * 
         TIME + CBT_ref1 * TIME + MI_ref1 * CBT_ref1 * TIME,
       #ref: Mi=1, Cbt=1, Time=0
       data = AVCIG_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

AVCIG_nobl_interaction_MIref1_CBTref1_TIMEref1 <- 
  geem(formula = AVCIG ~ cAGE + cBLCGSMD + MI_ref1 + CBT_ref1 + TIME_ref1 + 
         MI_ref1 * CBT_ref1 + MI_ref1 * TIME_ref1 + 
         CBT_ref1 * TIME_ref1 + MI_ref1 * CBT_ref1 * TIME_ref1,
       #ref: Mi=1, Cbt=1, Time=1
       data = AVCIG_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )


## print interaction model summary
summary(AVCIG_nobl_interaction) #base: MI=0, CBT=0, Time=0 
##                 Estimates Model SE Robust SE      wald      p
## (Intercept)     1.8020000 0.133200  0.147800 12.190000 0.0000
## cAGE           -0.0264300 0.053150  0.061120 -0.432500 0.6654
## cBLCGSMD        0.0509600 0.006433  0.005465  9.326000 0.0000
## MI1            -0.1181000 0.190900  0.194100 -0.608600 0.5428
## CBT1           -0.1225000 0.180400  0.185800 -0.659400 0.5096
## TIME1           0.1129000 0.093470  0.075660  1.493000 0.1356
## MI1:CBT1        0.3341000 0.256000  0.251800  1.327000 0.1845
## MI1:TIME1       0.1265000 0.133400  0.135700  0.932200 0.3512
## CBT1:TIME1      0.0009478 0.126400  0.109500  0.008652 0.9931
## MI1:CBT1:TIME1 -0.1793000 0.179200  0.182100 -0.984600 0.3248
## 
##  Estimated Correlation Parameter:  0.7493 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8299 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  246    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  492
summary(AVCIG_nobl_interaction_TIMEref1) #base: MI=0, CBT=0, Time=1
##                      Estimates Model SE Robust SE      wald      p
## (Intercept)          1.9150000 0.132100  0.123400 15.520000 0.0000
## cAGE                -0.0264300 0.053150  0.061120 -0.432500 0.6654
## cBLCGSMD             0.0509600 0.006433  0.005465  9.326000 0.0000
## MI1                  0.0083580 0.188500  0.190100  0.043960 0.9649
## CBT1                -0.1216000 0.178900  0.170100 -0.714700 0.4748
## TIME_ref10          -0.1129000 0.093470  0.075660 -1.493000 0.1356
## MI1:CBT1             0.1548000 0.253500  0.249900  0.619600 0.5355
## MI1:TIME_ref10      -0.1265000 0.133400  0.135700 -0.932200 0.3512
## CBT1:TIME_ref10     -0.0009478 0.126400  0.109500 -0.008652 0.9931
## MI1:CBT1:TIME_ref10  0.1793000 0.179200  0.182100  0.984600 0.3248
## 
##  Estimated Correlation Parameter:  0.7493 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8299 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  246    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  492
summary(AVCIG_nobl_interaction_MIref1) #ref: MI=1, CBT=0, Time=0
##                     Estimates Model SE Robust SE    wald       p
## (Intercept)           1.68400 0.135900  0.124700 13.5000 0.00000
## cAGE                 -0.02643 0.053150  0.061120 -0.4325 0.66540
## cBLCGSMD              0.05096 0.006433  0.005465  9.3260 0.00000
## MI_ref10              0.11810 0.190900  0.194100  0.6086 0.54280
## CBT1                  0.21160 0.180600  0.170100  1.2440 0.21340
## TIME1                 0.23940 0.095130  0.112600  2.1260 0.03352
## MI_ref10:CBT1        -0.33410 0.256000  0.251800 -1.3270 0.18450
## MI_ref10:TIME1       -0.12650 0.133400  0.135700 -0.9322 0.35120
## CBT1:TIME1           -0.17830 0.127000  0.145500 -1.2260 0.22020
## MI_ref10:CBT1:TIME1   0.17930 0.179200  0.182100  0.9846 0.32480
## 
##  Estimated Correlation Parameter:  0.7493 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8299 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  246    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  492
summary(AVCIG_nobl_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, Time=1
##                          Estimates Model SE Robust SE     wald       p
## (Intercept)               1.923000 0.133500  0.141500 13.59000 0.00000
## cAGE                     -0.026430 0.053150  0.061120 -0.43250 0.66540
## cBLCGSMD                  0.050960 0.006433  0.005465  9.32600 0.00000
## MI_ref10                 -0.008358 0.188500  0.190100 -0.04396 0.96490
## CBT1                      0.033240 0.178600  0.180800  0.18380 0.85420
## TIME_ref10               -0.239400 0.095130  0.112600 -2.12600 0.03352
## MI_ref10:CBT1            -0.154800 0.253500  0.249900 -0.61960 0.53550
## MI_ref10:TIME_ref10       0.126500 0.133400  0.135700  0.93220 0.35120
## CBT1:TIME_ref10           0.178300 0.127000  0.145500  1.22600 0.22020
## MI_ref10:CBT1:TIME_ref10 -0.179300 0.179200  0.182100 -0.98460 0.32480
## 
##  Estimated Correlation Parameter:  0.7493 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8299 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  246    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  492
summary(AVCIG_nobl_interaction_CBTref1) #ref: MI=0, CBT=1, Time=0
##                      Estimates Model SE Robust SE      wald      p
## (Intercept)          1.6790000 0.121000  0.114800 14.630000 0.0000
## cAGE                -0.0264300 0.053150  0.061120 -0.432500 0.6654
## cBLCGSMD             0.0509600 0.006433  0.005465  9.326000 0.0000
## MI1                  0.2160000 0.169800  0.163100  1.324000 0.1854
## CBT_ref10            0.1225000 0.180400  0.185800  0.659400 0.5096
## TIME1                0.1139000 0.085110  0.079230  1.437000 0.1506
## MI1:CBT_ref10       -0.3341000 0.256000  0.251800 -1.327000 0.1845
## MI1:TIME1           -0.0528100 0.119700  0.121500 -0.434800 0.6637
## CBT_ref10:TIME1     -0.0009478 0.126400  0.109500 -0.008652 0.9931
## MI1:CBT_ref10:TIME1  0.1793000 0.179200  0.182100  0.984600 0.3248
## 
##  Estimated Correlation Parameter:  0.7493 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8299 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  246    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  492
summary(AVCIG_nobl_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, Time=1
##                           Estimates Model SE Robust SE      wald      p
## (Intercept)               1.7930000 0.120000  0.117700 15.230000 0.0000
## cAGE                     -0.0264300 0.053150  0.061120 -0.432500 0.6654
## cBLCGSMD                  0.0509600 0.006433  0.005465  9.326000 0.0000
## MI1                       0.1632000 0.168700  0.163200  0.999900 0.3173
## CBT_ref10                 0.1216000 0.178900  0.170100  0.714700 0.4748
## TIME_ref10               -0.1139000 0.085110  0.079230 -1.437000 0.1506
## MI1:CBT_ref10            -0.1548000 0.253500  0.249900 -0.619600 0.5355
## MI1:TIME_ref10            0.0528100 0.119700  0.121500  0.434800 0.6637
## CBT_ref10:TIME_ref10      0.0009478 0.126400  0.109500  0.008652 0.9931
## MI1:CBT_ref10:TIME_ref10 -0.1793000 0.179200  0.182100 -0.984600 0.3248
## 
##  Estimated Correlation Parameter:  0.7493 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8299 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  246    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  492
summary(AVCIG_nobl_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, Time=0
##                          Estimates Model SE Robust SE    wald      p
## (Intercept)                1.89500 0.119200  0.116600 16.2600 0.0000
## cAGE                      -0.02643 0.053150  0.061120 -0.4325 0.6654
## cBLCGSMD                   0.05096 0.006433  0.005465  9.3260 0.0000
## MI_ref10                  -0.21600 0.169800  0.163100 -1.3240 0.1854
## CBT_ref10                 -0.21160 0.180600  0.170100 -1.2440 0.2134
## TIME1                      0.06106 0.084170  0.092070  0.6632 0.5072
## MI_ref10:CBT_ref10         0.33410 0.256000  0.251800  1.3270 0.1845
## MI_ref10:TIME1             0.05281 0.119700  0.121500  0.4348 0.6637
## CBT_ref10:TIME1            0.17830 0.127000  0.145500  1.2260 0.2202
## MI_ref10:CBT_ref10:TIME1  -0.17930 0.179200  0.182100 -0.9846 0.3248
## 
##  Estimated Correlation Parameter:  0.7493 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8299 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  246    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  492
summary(AVCIG_nobl_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, Time=1
##                               Estimates Model SE Robust SE    wald      p
## (Intercept)                     1.95600 0.118700  0.114400 17.0900 0.0000
## cAGE                           -0.02643 0.053150  0.061120 -0.4325 0.6654
## cBLCGSMD                        0.05096 0.006433  0.005465  9.3260 0.0000
## MI_ref10                       -0.16320 0.168700  0.163200 -0.9999 0.3173
## CBT_ref10                      -0.03324 0.178600  0.180800 -0.1838 0.8542
## TIME_ref10                     -0.06106 0.084170  0.092070 -0.6632 0.5072
## MI_ref10:CBT_ref10              0.15480 0.253500  0.249900  0.6196 0.5355
## MI_ref10:TIME_ref10            -0.05281 0.119700  0.121500 -0.4348 0.6637
## CBT_ref10:TIME_ref10           -0.17830 0.127000  0.145500 -1.2260 0.2202
## MI_ref10:CBT_ref10:TIME_ref10   0.17930 0.179200  0.182100  0.9846 0.3248
## 
##  Estimated Correlation Parameter:  0.7493 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8299 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  246    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  492
# DV = Cig/SmkD ---------------------------

#relevel MI and CBT
CGSMD_nobl_dt$MI_ref1 <- CGSMD_nobl_dt$MI %>% relevel("1")
CGSMD_nobl_dt$CBT_ref1 <- CGSMD_nobl_dt$CBT %>% relevel("1")
CGSMD_nobl_dt$TIME_ref1 <- CGSMD_nobl_dt$TIME %>% relevel("1")

# ignore missing values
## create an indicator column for any missing values 
CGSMD_nobl_dt$any_na <- CGSMD_nobl_dt %>% apply(1, function(x){any(is.na(x))}) 

# ## left join the dataset by the id column using the group_by() function
CGSMD_nobl_dt %<>% left_join(CGSMD_nobl_dt %>% group_by(id) %>% 
                               summarise(any_na2 = any(any_na)), by="id")


# ## filter out rows with any NAs
CGSMD_nobl_dt %<>% filter(any_na2 != T)

## main model
CGSMD_nobl_main <- geem(formula = CGSMD.nobl.main,
                        data = CGSMD_nobl_dt,
                        id=id,
                        family = MASS::negative.binomial(1),
                        corstr = corstr
)


## print main model summary
summary(CGSMD_nobl_main)
##             Estimates Model SE Robust SE    wald        p
## (Intercept)   1.96100 0.089830  0.095270 20.5800 0.000000
## cAGE         -0.02582 0.044690  0.050800 -0.5083 0.611200
## cBLCGSMD      0.04462 0.005711  0.004886  9.1330 0.000000
## MI1          -0.02427 0.098520  0.097120 -0.2499 0.802700
## CBT1          0.06931 0.098740  0.097130  0.7137 0.475400
## TIME1         0.09890 0.037940  0.038150  2.5920 0.009532
## 
##  Estimated Correlation Parameter:  0.7412 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.528 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  216    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  432
## interaction model
CGSMD_nobl_interaction <- geem(formula = CGSMD.nobl.interaction,
                               #mi=0, cbt=0, time=0
                               data = CGSMD_nobl_dt,
                               id=id,
                               family = MASS::negative.binomial(1),
                               corstr = corstr
)

CGSMD_nobl_interaction_TIMEref1 <- geem(formula = 
                                 CGSMD ~ cAGE + cBLCGSMD + MI + CBT + TIME_ref1 + MI * CBT + MI * TIME_ref1 + 
                                 CBT * TIME_ref1 + MI * CBT * TIME_ref1,
                               #mi=0, cbt=0, time=1
                               data = CGSMD_nobl_dt,
                               id=id,
                               family = MASS::negative.binomial(1),
                               corstr = corstr
)


CGSMD_nobl_interaction_MIref1 <- 
  geem(formula = CGSMD ~ cAGE + cBLCGSMD + MI_ref1 + CBT + TIME + MI_ref1 * CBT + MI_ref1 * TIME + 
         CBT * TIME + MI_ref1 * CBT * TIME,
       #ref: Mi=1, Cbt=0, time=0
       data = CGSMD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

CGSMD_nobl_interaction_MIref1_TIMEref1 <- 
  geem(formula = CGSMD ~ cAGE + cBLCGSMD + MI_ref1 + CBT + TIME_ref1 + 
         MI_ref1 * CBT + MI_ref1 * TIME_ref1 + 
         CBT * TIME_ref1 + MI_ref1 * CBT * TIME_ref1,
       #ref: Mi=1, Cbt=0, time=1
       data = CGSMD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

CGSMD_nobl_interaction_CBTref1 <- 
  geem(formula = CGSMD ~ cAGE + cBLCGSMD + MI + CBT_ref1 + TIME + MI * CBT_ref1 + MI * 
         TIME + CBT_ref1 * TIME + MI * CBT_ref1 * TIME,
       #ref: Mi=0, Cbt=1, time=0
       data = CGSMD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

CGSMD_nobl_interaction_CBTref1_TIMEref1 <- 
  geem(formula = CGSMD ~ cAGE + cBLCGSMD + MI + CBT_ref1 + TIME_ref1 + MI * CBT_ref1 + MI * 
         TIME_ref1 + CBT_ref1 * TIME_ref1 + MI * CBT_ref1 * TIME_ref1,
       #ref: Mi=0, Cbt=1, time=1
       data = CGSMD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

CGSMD_nobl_interaction_MIref1_CBTref1 <- 
  geem(formula = CGSMD ~ cAGE + cBLCGSMD + MI_ref1 + CBT_ref1 + TIME + MI_ref1 * CBT_ref1 + MI_ref1 * 
         TIME + CBT_ref1 * TIME + MI_ref1 * CBT_ref1 * TIME,
       #ref: Mi=0, Cbt=1, time=0
       data = CGSMD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

CGSMD_nobl_interaction_MIref1_CBTref1_TIMEref1 <- 
  geem(formula = CGSMD ~ cAGE + cBLCGSMD + MI_ref1 + CBT_ref1 + TIME_ref1 + MI_ref1 * CBT_ref1 + MI_ref1 * 
         TIME_ref1 + CBT_ref1 * TIME_ref1 + MI_ref1 * CBT_ref1 * TIME_ref1,
       #ref: Mi=0, Cbt=1, time=1
       data = CGSMD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )


## print interaction model summary
summary(CGSMD_nobl_interaction) #base: MI=0, CBT=0, time=0
##                Estimates Model SE Robust SE    wald      p
## (Intercept)      2.01000 0.110100  0.120000 16.7500 0.0000
## cAGE            -0.02455 0.044860  0.051060 -0.4808 0.6307
## cBLCGSMD         0.04529 0.005742  0.004835  9.3670 0.0000
## MI1             -0.15170 0.155900  0.160100 -0.9479 0.3432
## CBT1            -0.04097 0.153100  0.149300 -0.2744 0.7837
## TIME1            0.07267 0.077410  0.062440  1.1640 0.2445
## MI1:CBT1         0.25960 0.213500  0.207800  1.2490 0.2116
## MI1:TIME1        0.10730 0.109100  0.119300  0.8989 0.3687
## CBT1:TIME1       0.07835 0.107200  0.085290  0.9187 0.3583
## MI1:CBT1:TIME1  -0.24370 0.149300  0.151800 -1.6050 0.1084
## 
##  Estimated Correlation Parameter:  0.7507 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5267 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  216    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  432
summary(CGSMD_nobl_interaction_TIMEref1) #base: MI=0, CBT=0, time=1
##                     Estimates Model SE Robust SE     wald      p
## (Intercept)           2.08300 0.109700  0.104200 19.98000 0.0000
## cAGE                 -0.02455 0.044860  0.051060 -0.48080 0.6307
## cBLCGSMD              0.04529 0.005742  0.004835  9.36700 0.0000
## MI1                  -0.04446 0.154600  0.162400 -0.27370 0.7843
## CBT1                  0.03738 0.152200  0.135800  0.27520 0.7832
## TIME_ref10           -0.07267 0.077410  0.062440 -1.16400 0.2445
## MI1:CBT1              0.01595 0.212100  0.209700  0.07605 0.9394
## MI1:TIME_ref10       -0.10730 0.109100  0.119300 -0.89890 0.3687
## CBT1:TIME_ref10      -0.07835 0.107200  0.085290 -0.91870 0.3583
## MI1:CBT1:TIME_ref10   0.24370 0.149300  0.151800  1.60500 0.1084
## 
##  Estimated Correlation Parameter:  0.7507 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5267 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  216    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  432
summary(CGSMD_nobl_interaction_MIref1) #ref: MI=1, CBT=0, time=0
##                     Estimates Model SE Robust SE    wald       p
## (Intercept)           1.85800 0.109800  0.106200 17.5000 0.00000
## cAGE                 -0.02455 0.044860  0.051060 -0.4808 0.63070
## cBLCGSMD              0.04529 0.005742  0.004835  9.3670 0.00000
## MI_ref10              0.15170 0.155900  0.160100  0.9479 0.34320
## CBT1                  0.21860 0.147900  0.145100  1.5070 0.13190
## TIME1                 0.17990 0.076820  0.101700  1.7700 0.07681
## MI_ref10:CBT1        -0.25960 0.213500  0.207800 -1.2490 0.21160
## MI_ref10:TIME1       -0.10730 0.109100  0.119300 -0.8989 0.36870
## CBT1:TIME1           -0.16530 0.103800  0.125600 -1.3160 0.18800
## MI_ref10:CBT1:TIME1   0.24370 0.149300  0.151800  1.6050 0.10840
## 
##  Estimated Correlation Parameter:  0.7507 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5267 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  216    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  432
summary(CGSMD_nobl_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, time=1
##                          Estimates Model SE Robust SE     wald       p
## (Intercept)                2.03800 0.108500  0.123000 16.57000 0.00000
## cAGE                      -0.02455 0.044860  0.051060 -0.48080 0.63070
## cBLCGSMD                   0.04529 0.005742  0.004835  9.36700 0.00000
## MI_ref10                   0.04446 0.154600  0.162400  0.27370 0.78430
## CBT1                       0.05333 0.146900  0.157900  0.33770 0.73560
## TIME_ref10                -0.17990 0.076820  0.101700 -1.77000 0.07681
## MI_ref10:CBT1             -0.01595 0.212100  0.209700 -0.07605 0.93940
## MI_ref10:TIME_ref10        0.10730 0.109100  0.119300  0.89890 0.36870
## CBT1:TIME_ref10            0.16530 0.103800  0.125600  1.31600 0.18800
## MI_ref10:CBT1:TIME_ref10  -0.24370 0.149300  0.151800 -1.60500 0.10840
## 
##  Estimated Correlation Parameter:  0.7507 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5267 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  216    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  432
summary(CGSMD_nobl_interaction_CBTref1) #ref: MI=0, CBT=1, time=0
##                     Estimates Model SE Robust SE    wald       p
## (Intercept)           1.96900 0.105800  0.090330 21.8000 0.00000
## cAGE                 -0.02455 0.044860  0.051060 -0.4808 0.63070
## cBLCGSMD              0.04529 0.005742  0.004835  9.3670 0.00000
## MI1                   0.10790 0.145000  0.134300  0.8032 0.42190
## CBT_ref10             0.04097 0.153100  0.149300  0.2744 0.78370
## TIME1                 0.15100 0.074220  0.058110  2.5990 0.00935
## MI1:CBT_ref10        -0.25960 0.213500  0.207800 -1.2490 0.21160
## MI1:TIME1            -0.13640 0.101900  0.093840 -1.4540 0.14600
## CBT_ref10:TIME1      -0.07835 0.107200  0.085290 -0.9187 0.35830
## MI1:CBT_ref10:TIME1   0.24370 0.149300  0.151800  1.6050 0.10840
## 
##  Estimated Correlation Parameter:  0.7507 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5267 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  216    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  432
summary(CGSMD_nobl_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, time=1
##                          Estimates Model SE Robust SE     wald       p
## (Intercept)                2.12000 0.104900  0.087900 24.12000 0.00000
## cAGE                      -0.02455 0.044860  0.051060 -0.48080 0.63070
## cBLCGSMD                   0.04529 0.005742  0.004835  9.36700 0.00000
## MI1                       -0.02851 0.144200  0.132600 -0.21510 0.82970
## CBT_ref10                 -0.03738 0.152200  0.135800 -0.27520 0.78320
## TIME_ref10                -0.15100 0.074220  0.058110 -2.59900 0.00935
## MI1:CBT_ref10             -0.01595 0.212100  0.209700 -0.07605 0.93940
## MI1:TIME_ref10             0.13640 0.101900  0.093840  1.45400 0.14600
## CBT_ref10:TIME_ref10       0.07835 0.107200  0.085290  0.91870 0.35830
## MI1:CBT_ref10:TIME_ref10  -0.24370 0.149300  0.151800 -1.60500 0.10840
## 
##  Estimated Correlation Parameter:  0.7507 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5267 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  216    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  432
summary(CGSMD_nobl_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, time=0
##                          Estimates Model SE Robust SE    wald      p
## (Intercept)                2.07700 0.098990  0.099590 20.8500 0.0000
## cAGE                      -0.02455 0.044860  0.051060 -0.4808 0.6307
## cBLCGSMD                   0.04529 0.005742  0.004835  9.3670 0.0000
## MI_ref10                  -0.10790 0.145000  0.134300 -0.8032 0.4219
## CBT_ref10                 -0.21860 0.147900  0.145100 -1.5070 0.1319
## TIME1                      0.01461 0.069870  0.073680  0.1983 0.8428
## MI_ref10:CBT_ref10         0.25960 0.213500  0.207800  1.2490 0.2116
## MI_ref10:TIME1             0.13640 0.101900  0.093840  1.4540 0.1460
## CBT_ref10:TIME1            0.16530 0.103800  0.125600  1.3160 0.1880
## MI_ref10:CBT_ref10:TIME1  -0.24370 0.149300  0.151800 -1.6050 0.1084
## 
##  Estimated Correlation Parameter:  0.7507 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5267 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  216    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  432
summary(CGSMD_nobl_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, time=1
##                               Estimates Model SE Robust SE     wald      p
## (Intercept)                     2.09200 0.098910  0.099930 20.93000 0.0000
## cAGE                           -0.02455 0.044860  0.051060 -0.48080 0.6307
## cBLCGSMD                        0.04529 0.005742  0.004835  9.36700 0.0000
## MI_ref10                        0.02851 0.144200  0.132600  0.21510 0.8297
## CBT_ref10                      -0.05333 0.146900  0.157900 -0.33770 0.7356
## TIME_ref10                     -0.01461 0.069870  0.073680 -0.19830 0.8428
## MI_ref10:CBT_ref10              0.01595 0.212100  0.209700  0.07605 0.9394
## MI_ref10:TIME_ref10            -0.13640 0.101900  0.093840 -1.45400 0.1460
## CBT_ref10:TIME_ref10           -0.16530 0.103800  0.125600 -1.31600 0.1880
## MI_ref10:CBT_ref10:TIME_ref10   0.24370 0.149300  0.151800  1.60500 0.1084
## 
##  Estimated Correlation Parameter:  0.7507 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5267 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  216    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  432
# DV = PCSMKD ---------------------------

#relevel MI and CBT
PCSMKD_nobl_dt$MI_ref1 <- PCSMKD_nobl_dt$MI %>% relevel("1")
PCSMKD_nobl_dt$CBT_ref1 <- PCSMKD_nobl_dt$CBT %>% relevel("1")
PCSMKD_nobl_dt$TIME_ref1 <- PCSMKD_nobl_dt$TIME %>% relevel("1")

# check contrasts
contrasts(PCSMKD_nobl_dt$MI) #ref: MI=0
##   1
## 0 0
## 1 1
contrasts(PCSMKD_nobl_dt$MI_ref1) #ref: MI=1
##   0
## 1 0
## 0 1
contrasts(PCSMKD_nobl_dt$CBT) #ref: CBT=0
##   1
## 0 0
## 1 1
contrasts(PCSMKD_nobl_dt$CBT_ref1) #ref: CBT=1
##   0
## 1 0
## 0 1
# ignore missing values
## create an indicator column for any missing values 
PCSMKD_nobl_dt$any_na <- PCSMKD_nobl_dt %>% apply(1, function(x){any(is.na(x))}) 

# ## left join the dataset by the id column using the group_by() function
PCSMKD_nobl_dt %<>% left_join(PCSMKD_nobl_dt %>% group_by(id) %>% 
                               summarise(any_na2 = any(any_na)), by="id")


# ## filter out rows with any NAs
PCSMKD_nobl_dt %<>% filter(any_na2 != T)

## fit main model
PCSMKD_nobl_main <- geem(formula = PCSMKD.nobl.main,
                       data = PCSMKD_nobl_dt,
                       id=id,
                       family = MASS::negative.binomial(1),
                       corstr = corstr
)


## print main model summary
summary(PCSMKD_nobl_main) #ref: Mi=0, Cbt=0
##             Estimates Model SE Robust SE    wald        p
## (Intercept)  4.333000 0.045790  0.046850 92.4800 0.000000
## cAGE         0.013610 0.022390  0.023650  0.5755 0.564900
## cBLCGSMD     0.008767 0.002774  0.003216  2.7260 0.006414
## MI1          0.102300 0.049550  0.048600  2.1050 0.035310
## CBT1        -0.055010 0.049770  0.048560 -1.1330 0.257300
## TIME1        0.027570 0.024220  0.024300  1.1340 0.256600
## 
##  Estimated Correlation Parameter:  0.6129 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1818 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
## fit MI and CBT 2-way interactions
PCSMKD_MI_interactions <- geem(formula = PCSMDK.MI.interaction,
                               #ref: Mi=0, Cbt=0
                               data = PCSMKD_nobl_dt, 
                               id=id,
                               family = MASS::negative.binomial(1),
                               corstr = corstr
)

PCSMKD_MI_interactions_MIref1_CBT_ref0 <- 
  geem(PCSMKD ~ cAGE + cBLCGSMD + MI_ref1*TIME + CBT, 
       #ref: Mi=1, Cbt=0
       data = PCSMKD_nobl_dt ,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )


PCSMKD_CBT_interactions <- geem(formula = PCSMDK.CBT.interaction,
                                #ref: Mi=0, Cbt=0
                                data = PCSMKD_nobl_dt,
                                id=id,
                                family = MASS::negative.binomial(1),
                                corstr = corstr
)


## print main model summaries
summary(PCSMKD_MI_interactions)
##             Estimates Model SE Robust SE    wald        p
## (Intercept)  4.350000  0.04744  0.048590 89.5300 0.000000
## cAGE         0.013420  0.02245  0.023670  0.5671 0.570700
## cBLCGSMD     0.008770  0.00278  0.003221  2.7230 0.006469
## MI1          0.068350  0.05524  0.055790  1.2250 0.220500
## CBT1        -0.055710  0.04988  0.048600 -1.1460 0.251700
## TIME1       -0.006276  0.03415  0.036440 -0.1722 0.863300
## MI1:TIME1    0.067900  0.04837  0.048550  1.3980 0.162000
## 
##  Estimated Correlation Parameter:  0.6152 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1824 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
summary(PCSMKD_MI_interactions_MIref1_CBT_ref0)
##                Estimates Model SE Robust SE     wald        p
## (Intercept)      4.41800  0.04806  0.043320 102.0000 0.000000
## cAGE             0.01342  0.02245  0.023670   0.5671 0.570700
## cBLCGSMD         0.00877  0.00278  0.003221   2.7230 0.006469
## MI_ref10        -0.06835  0.05524  0.055790  -1.2250 0.220500
## TIME1            0.06162  0.03426  0.032090   1.9200 0.054810
## CBT1            -0.05571  0.04988  0.048600  -1.1460 0.251700
## MI_ref10:TIME1  -0.06790  0.04837  0.048550  -1.3980 0.162000
## 
##  Estimated Correlation Parameter:  0.6152 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1824 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
summary(PCSMKD_CBT_interactions)
##             Estimates Model SE Robust SE    wald        p
## (Intercept)   4.31300  0.04780   0.05075 84.9800 0.000000
## cAGE          0.01357  0.02245   0.02365  0.5735 0.566300
## cBLCGSMD      0.00879  0.00278   0.00322  2.7300 0.006327
## MI1           0.10300  0.04966   0.04864  2.1180 0.034160
## CBT1         -0.01928  0.05549   0.05630 -0.3424 0.732000
## TIME1         0.06690  0.03607   0.03507  1.9070 0.056460
## CBT1:TIME1   -0.07137  0.04858   0.04850 -1.4720 0.141200
## 
##  Estimated Correlation Parameter:  0.6159 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1823 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
## fit 3-way interaction model
PCSMKD_nobl_interaction <- 
  geem(formula = PCSMKD.nobl.interaction,
       #ref: Mi=0, Cbt=0, time=0
       data = PCSMKD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

PCSMKD_nobl_interaction_TIMEref1 <- 
  geem(formula = PCSMKD ~ cAGE + cBLCGSMD + MI + CBT + TIME_ref1 + MI * CBT + MI * 
         TIME_ref1 + CBT * TIME_ref1 + MI * CBT * TIME_ref1,
       #ref: Mi=0, Cbt=0, time=1
       data = PCSMKD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

PCSMKD_nobl_interaction_MIref1 <- 
  geem(formula = PCSMKD ~ cAGE + cBLCGSMD + MI_ref1 + CBT + TIME + MI_ref1 * CBT + MI_ref1 * 
         TIME + CBT * TIME + MI_ref1 * CBT * TIME,
       #ref: Mi=1, Cbt=0, time=0
       data = PCSMKD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

PCSMKD_nobl_interaction_MIref1_TIMEref1 <- 
  geem(formula = PCSMKD ~ cAGE + cBLCGSMD + MI_ref1 + CBT + TIME_ref1 + MI_ref1 * CBT + MI_ref1 * 
         TIME_ref1 + CBT * TIME_ref1 + MI_ref1 * CBT * TIME_ref1,
       #ref: Mi=1, Cbt=0, time=1
       data = PCSMKD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

PCSMKD_nobl_interaction_CBTref1 <- 
  geem(formula = PCSMKD ~ cAGE + cBLCGSMD + MI + CBT_ref1 + TIME + MI * CBT_ref1 + MI * 
         TIME + CBT_ref1 * TIME + MI * CBT_ref1 * TIME,
       #ref: Mi=0, Cbt=1, time=0
       data = PCSMKD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

PCSMKD_nobl_interaction_CBTref1_TIMEref1 <- 
  geem(formula = PCSMKD ~ cAGE + cBLCGSMD + MI + CBT_ref1 + TIME_ref1 + MI * CBT_ref1 + MI * 
         TIME_ref1 + CBT_ref1 * TIME_ref1 + MI * CBT_ref1 * TIME_ref1,
       #ref: Mi=0, Cbt=1, time=1
       data = PCSMKD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

PCSMKD_nobl_interaction_MIref1_CBTref1 <- 
  geem(formula = PCSMKD ~ cAGE + cBLCGSMD + MI_ref1 + CBT_ref1 + TIME + MI_ref1 * CBT_ref1 + MI_ref1 * 
         TIME + CBT_ref1 * TIME + MI_ref1 * CBT_ref1 * TIME,
       #ref: Mi=1, Cbt=1, time=0
       data = PCSMKD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

PCSMKD_nobl_interaction_MIref1_CBTref1_TIMEref1 <- 
  geem(formula = PCSMKD ~ cAGE + cBLCGSMD + MI_ref1 + CBT_ref1 + TIME_ref1 + MI_ref1 * CBT_ref1 + MI_ref1 * 
         TIME_ref1 + CBT_ref1 * TIME_ref1 + MI_ref1 * CBT_ref1 * TIME_ref1,
       #ref: Mi=1, Cbt=1, time=1
       data = PCSMKD_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

## print interaction model summary
summary(PCSMKD_nobl_interaction) #base: MI=0, CBT=0, time=0
##                 Estimates Model SE Robust SE      wald        p
## (Intercept)     4.3300000  0.05797  0.064030 67.620000 0.000000
## cAGE            0.0146700  0.02263  0.023800  0.616200 0.537800
## cBLCGSMD        0.0090460  0.00281  0.003239  2.793000 0.005224
## MI1             0.0676200  0.08316  0.081390  0.830700 0.406100
## CBT1           -0.0186100  0.07880  0.085690 -0.217200 0.828100
## TIME1           0.0789700  0.04993  0.055900  1.413000 0.157800
## MI1:CBT1       -0.0001742  0.11190  0.112500 -0.001549 0.998800
## MI1:TIME1      -0.0253200  0.07158  0.069030 -0.366800 0.713800
## CBT1:TIME1     -0.1607000  0.06790  0.073140 -2.197000 0.027990
## MI1:CBT1:TIME1  0.1746000  0.09640  0.096000  1.819000 0.068990
## 
##  Estimated Correlation Parameter:  0.6249 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1837 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
summary(PCSMKD_nobl_interaction_TIMEref1) #base: MI=0, CBT=0, time=1
##                     Estimates Model SE Robust SE    wald        p
## (Intercept)          4.409000  0.05794  0.054990 80.1700 0.000000
## cAGE                 0.014670  0.02263  0.023800  0.6162 0.537800
## cBLCGSMD             0.009046  0.00281  0.003239  2.7930 0.005224
## MI1                  0.042290  0.08313  0.069030  0.6126 0.540100
## CBT1                -0.179300  0.07880  0.086900 -2.0640 0.039060
## TIME_ref10          -0.078970  0.04993  0.055900 -1.4130 0.157800
## MI1:CBT1             0.174400  0.11190  0.105600  1.6520 0.098580
## MI1:TIME_ref10       0.025320  0.07158  0.069030  0.3668 0.713800
## CBT1:TIME_ref10      0.160700  0.06790  0.073140  2.1970 0.027990
## MI1:CBT1:TIME_ref10 -0.174600  0.09640  0.096000 -1.8190 0.068990
## 
##  Estimated Correlation Parameter:  0.6249 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1837 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
summary(PCSMKD_nobl_interaction_MIref1) #ref: MI=1, CBT=0, time=0
##                      Estimates Model SE Robust SE      wald        p
## (Intercept)          4.3970000  0.05932  0.049960 88.010000 0.000000
## cAGE                 0.0146700  0.02263  0.023800  0.616200 0.537800
## cBLCGSMD             0.0090460  0.00281  0.003239  2.793000 0.005224
## MI_ref10            -0.0676200  0.08316  0.081390 -0.830700 0.406100
## CBT1                -0.0187800  0.07911  0.072500 -0.259100 0.795500
## TIME1                0.0536500  0.05130  0.040500  1.325000 0.185300
## MI_ref10:CBT1        0.0001742  0.11190  0.112500  0.001549 0.998800
## MI_ref10:TIME1       0.0253200  0.07158  0.069030  0.366800 0.713800
## CBT1:TIME1           0.0138500  0.06843  0.062170  0.222800 0.823700
## MI_ref10:CBT1:TIME1 -0.1746000  0.09640  0.096000 -1.819000 0.068990
## 
##  Estimated Correlation Parameter:  0.6249 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1837 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
summary(PCSMKD_nobl_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, time=1
##                          Estimates Model SE Robust SE      wald        p
## (Intercept)               4.451000  0.05930  0.040970 108.70000 0.000000
## cAGE                      0.014670  0.02263  0.023800   0.61620 0.537800
## cBLCGSMD                  0.009046  0.00281  0.003239   2.79300 0.005224
## MI_ref10                 -0.042290  0.08313  0.069030  -0.61260 0.540100
## CBT1                     -0.004933  0.07908  0.059260  -0.08323 0.933700
## TIME_ref10               -0.053650  0.05130  0.040500  -1.32500 0.185300
## MI_ref10:CBT1            -0.174400  0.11190  0.105600  -1.65200 0.098580
## MI_ref10:TIME_ref10      -0.025320  0.07158  0.069030  -0.36680 0.713800
## CBT1:TIME_ref10          -0.013850  0.06843  0.062170  -0.22280 0.823700
## MI_ref10:CBT1:TIME_ref10  0.174600  0.09640  0.096000   1.81900 0.068990
## 
##  Estimated Correlation Parameter:  0.6249 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1837 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
summary(PCSMKD_nobl_interaction_CBTref1) #ref: MI=0, CBT=1, time=0
##                      Estimates Model SE Robust SE      wald        p
## (Intercept)          4.3110000  0.05316  0.056170 76.750000 0.000000
## cAGE                 0.0146700  0.02263  0.023800  0.616200 0.537800
## cBLCGSMD             0.0090460  0.00281  0.003239  2.793000 0.005224
## MI1                  0.0674400  0.07462  0.076760  0.878600 0.379600
## CBT_ref10            0.0186100  0.07880  0.085690  0.217200 0.828100
## TIME1               -0.0817500  0.04601  0.047170 -1.733000 0.083050
## MI1:CBT_ref10        0.0001742  0.11190  0.112500  0.001549 0.998800
## MI1:TIME1            0.1492000  0.06457  0.066710  2.237000 0.025260
## CBT_ref10:TIME1      0.1607000  0.06790  0.073140  2.197000 0.027990
## MI1:CBT_ref10:TIME1 -0.1746000  0.09640  0.096000 -1.819000 0.068990
## 
##  Estimated Correlation Parameter:  0.6249 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1837 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
summary(PCSMKD_nobl_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, time=1
##                          Estimates Model SE Robust SE    wald        p
## (Intercept)               4.229000  0.05319  0.066420 63.6700 0.000000
## cAGE                      0.014670  0.02263  0.023800  0.6162 0.537800
## cBLCGSMD                  0.009046  0.00281  0.003239  2.7930 0.005224
## MI1                       0.216700  0.07462  0.078930  2.7450 0.006043
## CBT_ref10                 0.179300  0.07880  0.086900  2.0640 0.039060
## TIME_ref10                0.081750  0.04601  0.047170  1.7330 0.083050
## MI1:CBT_ref10            -0.174400  0.11190  0.105600 -1.6520 0.098580
## MI1:TIME_ref10           -0.149200  0.06457  0.066710 -2.2370 0.025260
## CBT_ref10:TIME_ref10     -0.160700  0.06790  0.073140 -2.1970 0.027990
## MI1:CBT_ref10:TIME_ref10  0.174600  0.09640  0.096000  1.8190 0.068990
## 
##  Estimated Correlation Parameter:  0.6249 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1837 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
summary(PCSMKD_nobl_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, time=0
##                           Estimates Model SE Robust SE      wald        p
## (Intercept)               4.3790000  0.05233  0.052640 83.180000 0.000000
## cAGE                      0.0146700  0.02263  0.023800  0.616200 0.537800
## cBLCGSMD                  0.0090460  0.00281  0.003239  2.793000 0.005224
## MI_ref10                 -0.0674400  0.07462  0.076760 -0.878600 0.379600
## CBT_ref10                 0.0187800  0.07911  0.072500  0.259100 0.795500
## TIME1                     0.0675000  0.04529  0.047170  1.431000 0.152500
## MI_ref10:CBT_ref10       -0.0001742  0.11190  0.112500 -0.001549 0.998800
## MI_ref10:TIME1           -0.1492000  0.06457  0.066710 -2.237000 0.025260
## CBT_ref10:TIME1          -0.0138500  0.06843  0.062170 -0.222800 0.823700
## MI_ref10:CBT_ref10:TIME1  0.1746000  0.09640  0.096000  1.819000 0.068990
## 
##  Estimated Correlation Parameter:  0.6249 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1837 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
summary(PCSMKD_nobl_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, time=1
##                               Estimates Model SE Robust SE      wald        p
## (Intercept)                    4.446000  0.05231  0.043060 103.30000 0.000000
## cAGE                           0.014670  0.02263  0.023800   0.61620 0.537800
## cBLCGSMD                       0.009046  0.00281  0.003239   2.79300 0.005224
## MI_ref10                      -0.216700  0.07462  0.078930  -2.74500 0.006043
## CBT_ref10                      0.004933  0.07908  0.059260   0.08323 0.933700
## TIME_ref10                    -0.067500  0.04529  0.047170  -1.43100 0.152500
## MI_ref10:CBT_ref10             0.174400  0.11190  0.105600   1.65200 0.098580
## MI_ref10:TIME_ref10            0.149200  0.06457  0.066710   2.23700 0.025260
## CBT_ref10:TIME_ref10           0.013850  0.06843  0.062170   0.22280 0.823700
## MI_ref10:CBT_ref10:TIME_ref10 -0.174600  0.09640  0.096000  -1.81900 0.068990
## 
##  Estimated Correlation Parameter:  0.6249 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1837 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  243    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  486
# DV = PPVABS ---------------------------

#relevel MI and CBT
PPVABS_nobl_dt$MI_ref1 <- PPVABS_nobl_dt$MI %>% relevel("1")
PPVABS_nobl_dt$CBT_ref1 <- PPVABS_nobl_dt$CBT %>% relevel("1")
PPVABS_nobl_dt$TIME_ref1 <- PPVABS_nobl_dt$TIME %>% relevel("1")

# ignore missing values
## create an indicator column for any missing values 
PPVABS_nobl_dt$any_na <- PPVABS_nobl_dt %>% apply(1, function(x){any(is.na(x))}) 

# ## left join the dataset by the id column using the group_by() function
PPVABS_nobl_dt %<>% left_join(PPVABS_nobl_dt %>% group_by(id) %>% 
                               summarise(any_na2 = any(any_na)), by="id")


# ## filter out rows with any NAs
PPVABS_nobl_dt %<>% filter(any_na2 != T)


## fit main model
PPVABS_nobl_main <- geem(formula = PPVABS.nobl.main,
                        data = PPVABS_nobl_dt,
                        id=id,
                        family = binomial,
                        corstr = corstr
)
## Warning in geem(formula = PPVABS.nobl.main, data = PPVABS_nobl_dt, id = id, :
## Did not converge
## print model summary
summary(PPVABS_nobl_main)
##              Estimates Model SE Robust SE       wald         p
## (Intercept) -1.338e+00  0.24400   0.24570 -5.445e+00 5.000e-08
## cAGE        -3.475e-01  0.12080   0.10210 -3.403e+00 6.665e-04
## cBLCGSMD    -2.310e-02  0.01568   0.01911 -1.209e+00 2.267e-01
## MI1         -1.101e-01  0.26170   0.25810 -4.267e-01 6.696e-01
## CBT1         1.659e-01  0.26240   0.26370  6.293e-01 5.291e-01
## TIME1        1.840e-17  0.13640   0.13510  1.361e-16 1.000e+00
## 
##  Estimated Correlation Parameter:  0.569 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.01 
## 
##  Number of GEE iterations: 20 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
## fit 3-way interaction model
PPVABS_nobl_interaction <- 
  geem(formula = PPVABS.nobl.interaction,
       #ref: Mi=0, Cbt=0, time=0
       data = PPVABS_nobl_dt,
       id=id,
       family = binomial,
       corstr = corstr
  )

PPVABS_nobl_interaction_TIMEref1 <- 
  geem(formula = PPVABS ~ cAGE + cBLCGSMD + MI + CBT + TIME_ref1 + MI * CBT + MI * 
         TIME_ref1 + CBT * TIME_ref1 + MI * CBT * TIME_ref1,
       #ref: Mi=0, Cbt=0, time=1
       data = PPVABS_nobl_dt,
       id=id,
       family = binomial,
       corstr = corstr
  )

PPVABS_nobl_interaction_MIref1 <- 
  geem(formula = PPVABS ~ cAGE + cBLCGSMD + MI_ref1 + CBT + TIME + MI_ref1 * CBT + MI_ref1 * 
         TIME + CBT * TIME + MI_ref1 * CBT * TIME,
       #ref: Mi=1, Cbt=0, time=0
       data = PPVABS_nobl_dt,
       id=id,
       family = binomial,
       corstr = corstr
  )

PPVABS_nobl_interaction_MIref1_TIMEref1 <- 
  geem(formula = PPVABS ~ cAGE + cBLCGSMD + MI_ref1 + CBT + TIME_ref1 + MI_ref1 * CBT + MI_ref1 * 
         TIME_ref1 + CBT * TIME_ref1 + MI_ref1 * CBT * TIME_ref1,
       #ref: Mi=1, Cbt=0, time=1
       data = PPVABS_nobl_dt,
       id=id,
       family = binomial,
       corstr = corstr
  )

PPVABS_nobl_interaction_CBTref1 <- 
  geem(formula = PPVABS ~ cAGE + cBLCGSMD + MI + CBT_ref1 + TIME + MI * CBT_ref1 + MI * 
         TIME + CBT_ref1 * TIME + MI * CBT_ref1 * TIME,
       #ref: Mi=0, Cbt=1, time=0
       data = PPVABS_nobl_dt,
       id=id,
       family = binomial,
       corstr = corstr
  )

PPVABS_nobl_interaction_CBTref1_TIMEref1 <- 
  geem(formula = PPVABS ~ cAGE + cBLCGSMD + MI + CBT_ref1 + TIME_ref1 + MI * CBT_ref1 + MI * 
         TIME_ref1 + CBT_ref1 * TIME_ref1 + MI * CBT_ref1 * TIME_ref1,
       #ref: Mi=0, Cbt=1, time=1
       data = PPVABS_nobl_dt,
       id=id,
       family = binomial,
       corstr = corstr
  )

PPVABS_nobl_interaction_MIref1_CBTref1 <- 
  geem(formula = PPVABS ~ cAGE + cBLCGSMD + MI_ref1 + CBT_ref1 + TIME + MI_ref1 * CBT_ref1 + MI_ref1 * 
         TIME + CBT_ref1 * TIME + MI_ref1 * CBT_ref1 * TIME,
       #ref: Mi=0, Cbt=1, time=0
       data = PPVABS_nobl_dt,
       id=id,
       family = binomial,
       corstr = corstr
  )

PPVABS_nobl_interaction_MIref1_CBTref1_TIMEref1 <- 
  geem(formula = PPVABS ~ cAGE + cBLCGSMD + MI_ref1 + CBT_ref1 + TIME_ref1 + MI_ref1 * CBT_ref1 + MI_ref1 * 
         TIME_ref1 + CBT_ref1 * TIME_ref1 + MI_ref1 * CBT_ref1 * TIME_ref1,
       #ref: Mi=0, Cbt=1, time=1
       data = PPVABS_nobl_dt,
       id=id,
       family = binomial,
       corstr = corstr
  )

## print interaction model summary
summary(PPVABS_nobl_interaction) #base: MI=0, CBT=0, TIME=0
##                Estimates Model SE Robust SE     wald         p
## (Intercept)     -1.29600  0.30370    0.3033 -4.27400 1.917e-05
## cAGE            -0.35640  0.12210    0.1029 -3.46200 5.355e-04
## cBLCGSMD        -0.02131  0.01562    0.0189 -1.12700 2.597e-01
## MI1             -0.49950  0.47360    0.4754 -1.05100 2.934e-01
## CBT1             0.16200  0.40480    0.4072  0.39790 6.907e-01
## TIME1           -0.08773  0.27320    0.2647 -0.33150 7.403e-01
## MI1:CBT1         0.50510  0.61170    0.6096  0.82870 4.073e-01
## MI1:TIME1        0.72270  0.41420    0.4566  1.58300 1.134e-01
## CBT1:TIME1       0.01319  0.37110    0.3621  0.03643 9.709e-01
## MI1:CBT1:TIME1  -0.97240  0.55010    0.5651 -1.72100 8.529e-02
## 
##  Estimated Correlation Parameter:  0.5833 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.014 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(PPVABS_nobl_interaction_TIMEref1) #base: MI=0, CBT=0, TIME=1
##                     Estimates Model SE Robust SE     wald         p
## (Intercept)          -1.38400  0.31010    0.3149 -4.39500 1.107e-05
## cAGE                 -0.35640  0.12210    0.1029 -3.46200 5.355e-04
## cBLCGSMD             -0.02131  0.01562    0.0189 -1.12700 2.597e-01
## MI1                   0.22320  0.43250    0.4362  0.51170 6.089e-01
## CBT1                  0.17520  0.41270    0.4168  0.42040 6.742e-01
## TIME_ref10            0.08773  0.27320    0.2647  0.33150 7.403e-01
## MI1:CBT1             -0.46720  0.59440    0.5933 -0.78750 4.310e-01
## MI1:TIME_ref10       -0.72270  0.41420    0.4566 -1.58300 1.134e-01
## CBT1:TIME_ref10      -0.01319  0.37110    0.3621 -0.03643 9.709e-01
## MI1:CBT1:TIME_ref10   0.97240  0.55010    0.5651  1.72100 8.529e-02
## 
##  Estimated Correlation Parameter:  0.5833 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.014 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(PPVABS_nobl_interaction_MIref1) #ref: MI=1, CBT=0, TIME=0
##                     Estimates Model SE Robust SE    wald         p
## (Intercept)          -1.79600  0.36710    0.3754 -4.7830 1.720e-06
## cAGE                 -0.35640  0.12210    0.1029 -3.4620 5.355e-04
## cBLCGSMD             -0.02131  0.01562    0.0189 -1.1270 2.597e-01
## MI_ref10              0.49950  0.47360    0.4754  1.0510 2.934e-01
## CBT1                  0.66720  0.45720    0.4552  1.4660 1.427e-01
## TIME1                 0.63500  0.31130    0.3717  1.7080 8.760e-02
## MI_ref10:CBT1        -0.50510  0.61170    0.6096 -0.8287 4.073e-01
## MI_ref10:TIME1       -0.72270  0.41420    0.4566 -1.5830 1.134e-01
## CBT1:TIME1           -0.95920  0.40620    0.4334 -2.2130 2.690e-02
## MI_ref10:CBT1:TIME1   0.97240  0.55010    0.5651  1.7210 8.529e-02
## 
##  Estimated Correlation Parameter:  0.5833 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.014 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(PPVABS_nobl_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, TIME=1
##                          Estimates Model SE Robust SE    wald         p
## (Intercept)               -1.16100  0.30430    0.3122 -3.7180 0.0002009
## cAGE                      -0.35640  0.12210    0.1029 -3.4620 0.0005355
## cBLCGSMD                  -0.02131  0.01562    0.0189 -1.1270 0.2597000
## MI_ref10                  -0.22320  0.43250    0.4362 -0.5117 0.6089000
## CBT1                      -0.29200  0.42520    0.4226 -0.6910 0.4895000
## TIME_ref10                -0.63500  0.31130    0.3717 -1.7080 0.0876000
## MI_ref10:CBT1              0.46720  0.59440    0.5933  0.7875 0.4310000
## MI_ref10:TIME_ref10        0.72270  0.41420    0.4566  1.5830 0.1134000
## CBT1:TIME_ref10            0.95920  0.40620    0.4334  2.2130 0.0269000
## MI_ref10:CBT1:TIME_ref10  -0.97240  0.55010    0.5651 -1.7210 0.0852900
## 
##  Estimated Correlation Parameter:  0.5833 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.014 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(PPVABS_nobl_interaction_CBTref1) #ref: MI=0, CBT=1, TIME=0
##                     Estimates Model SE Robust SE     wald         p
## (Intercept)         -1.134000  0.27390    0.2730 -4.15400 3.262e-05
## cAGE                -0.356400  0.12210    0.1029 -3.46200 5.355e-04
## cBLCGSMD            -0.021310  0.01562    0.0189 -1.12700 2.597e-01
## MI1                  0.005599  0.38640    0.3780  0.01481 9.882e-01
## CBT_ref10           -0.162000  0.40480    0.4072 -0.39790 6.907e-01
## TIME1               -0.074540  0.25110    0.2471 -0.30170 7.629e-01
## MI1:CBT_ref10       -0.505100  0.61170    0.6096 -0.82870 4.073e-01
## MI1:TIME1           -0.249600  0.36200    0.3327 -0.75030 4.531e-01
## CBT_ref10:TIME1     -0.013190  0.37110    0.3621 -0.03643 9.709e-01
## MI1:CBT_ref10:TIME1  0.972400  0.55010    0.5651  1.72100 8.529e-02
## 
##  Estimated Correlation Parameter:  0.5833 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.014 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(PPVABS_nobl_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, TIME=1
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)               -1.20900  0.27890    0.2741 -4.41000 1.033e-05
## cAGE                      -0.35640  0.12210    0.1029 -3.46200 5.355e-04
## cBLCGSMD                  -0.02131  0.01562    0.0189 -1.12700 2.597e-01
## MI1                       -0.24400  0.40660    0.3979 -0.61330 5.397e-01
## CBT_ref10                 -0.17520  0.41270    0.4168 -0.42040 6.742e-01
## TIME_ref10                 0.07454  0.25110    0.2471  0.30170 7.629e-01
## MI1:CBT_ref10              0.46720  0.59440    0.5933  0.78750 4.310e-01
## MI1:TIME_ref10             0.24960  0.36200    0.3327  0.75030 4.531e-01
## CBT_ref10:TIME_ref10       0.01319  0.37110    0.3621  0.03643 9.709e-01
## MI1:CBT_ref10:TIME_ref10  -0.97240  0.55010    0.5651 -1.72100 8.529e-02
## 
##  Estimated Correlation Parameter:  0.5833 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.014 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(PPVABS_nobl_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, TIME=0
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)              -1.129000  0.27520    0.2611 -4.32300 1.538e-05
## cAGE                     -0.356400  0.12210    0.1029 -3.46200 5.355e-04
## cBLCGSMD                 -0.021310  0.01562    0.0189 -1.12700 2.597e-01
## MI_ref10                 -0.005599  0.38640    0.3780 -0.01481 9.882e-01
## CBT_ref10                -0.667200  0.45720    0.4552 -1.46600 1.427e-01
## TIME1                    -0.324200  0.26080    0.2227 -1.45500 1.456e-01
## MI_ref10:CBT_ref10        0.505100  0.61170    0.6096  0.82870 4.073e-01
## MI_ref10:TIME1            0.249600  0.36200    0.3327  0.75030 4.531e-01
## CBT_ref10:TIME1           0.959200  0.40620    0.4334  2.21300 2.690e-02
## MI_ref10:CBT_ref10:TIME1 -0.972400  0.55010    0.5651 -1.72100 8.529e-02
## 
##  Estimated Correlation Parameter:  0.5833 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.014 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(PPVABS_nobl_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, TIME=1
##                               Estimates Model SE Robust SE    wald         p
## (Intercept)                    -1.45300  0.29890    0.2879 -5.0460 4.500e-07
## cAGE                           -0.35640  0.12210    0.1029 -3.4620 5.355e-04
## cBLCGSMD                       -0.02131  0.01562    0.0189 -1.1270 2.597e-01
## MI_ref10                        0.24400  0.40660    0.3979  0.6133 5.397e-01
## CBT_ref10                       0.29200  0.42520    0.4226  0.6910 4.895e-01
## TIME_ref10                      0.32420  0.26080    0.2227  1.4550 1.456e-01
## MI_ref10:CBT_ref10             -0.46720  0.59440    0.5933 -0.7875 4.310e-01
## MI_ref10:TIME_ref10            -0.24960  0.36200    0.3327 -0.7503 4.531e-01
## CBT_ref10:TIME_ref10           -0.95920  0.40620    0.4334 -2.2130 2.690e-02
## MI_ref10:CBT_ref10:TIME_ref10   0.97240  0.55010    0.5651  1.7210 8.529e-02
## 
##  Estimated Correlation Parameter:  0.5833 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.014 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
# DV: NCIGS ---------------------------

#relevel MI and CBT
ncigs_dt_na.omit$MI_ref1 <- ncigs_dt_na.omit$MI %>% relevel("1")
ncigs_dt_na.omit$CBT_ref1 <- ncigs_dt_na.omit$CBT %>% relevel("1")
ncigs_dt_na.omit$TIME_ref1 <- ncigs_dt_na.omit$TIME %>% relevel("1")
ncigs_dt_na.omit$TIME_ref2 <- ncigs_dt_na.omit$TIME %>% relevel("2")

# ignore missing values
## create an indicator column for any missing values 
ncigs_dt_na.omit$any_na <- ncigs_dt_na.omit %>% apply(1, function(x){any(is.na(x))}) 

# ## left join the dataset by the id column using the group_by() function
ncigs_dt_na.omit %<>% left_join(ncigs_dt_na.omit %>% group_by(id) %>% 
                               summarise(any_na2 = any(any_na)), by="id")


# ## filter out rows with any NAs
ncigs_dt_na.omit %<>% filter(any_na2 != T)

## fit main model
ncigs_main <- geem(formula = ncigs.main,
                  data = ncigs_dt_na.omit,
                  id=id,
                  family = MASS::negative.binomial(1),
                  corstr = corstr
)


## print model summary
summary(ncigs_main)
##             Estimates Model SE Robust SE    wald        p
## (Intercept)   2.40600  0.09502   0.08162 29.4800 0.000000
## cAGE          0.11440  0.04536   0.04437  2.5770 0.009962
## MI1           0.03442  0.10080   0.09391  0.3665 0.714000
## CBT1         -0.01055  0.10090   0.09603 -0.1099 0.912500
## TIME1        -0.55700  0.05347   0.06819 -8.1680 0.000000
## TIME2        -0.47850  0.06816   0.07091 -6.7480 0.000000
## 
##  Estimated Correlation Parameter:  0.6187 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.002 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
## fit 3-way interaction model
ncigs_interaction <- geem(formula = ncigs.interaction,
                          #ref: Mi=0, Cbt=0, time=0
                   data = ncigs_dt_na.omit,
                   id=id,
                   family = MASS::negative.binomial(1),
                   corstr = corstr
)

ncigs_interaction_TIMEref1 <- geem(formula = 
                            NCIGS ~ cAGE + MI + CBT + TIME_ref1 + MI * CBT + MI * TIME_ref1 + CBT * 
                            TIME_ref1 + MI * CBT * TIME_ref1 + offset(log_unctrldays),
                          #ref: Mi=0, Cbt=0, time=1
                          data = ncigs_dt_na.omit,
                          id=id,
                          family = MASS::negative.binomial(1),
                          corstr = corstr
)

ncigs_interaction_TIMEref2 <- geem(formula = 
                                     NCIGS ~ cAGE + MI + CBT + TIME_ref2 + MI * CBT + MI * TIME_ref2 + CBT * 
                                     TIME_ref2 + MI * CBT * TIME_ref2 + offset(log_unctrldays),
                                   #ref: Mi=0, Cbt=0, time=1
                                   data = ncigs_dt_na.omit,
                                   id=id,
                                   family = MASS::negative.binomial(1),
                                   corstr = corstr
)

ncigs_interaction_MIref1 <- 
  geem(formula = NCIGS ~ cAGE + MI_ref1 + CBT + TIME + MI_ref1 * CBT + MI_ref1 * TIME + CBT * 
         TIME + MI_ref1 * CBT * TIME + offset(log_unctrldays),
       #ref: Mi=1, Cbt=0, time=0
       data = ncigs_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

ncigs_interaction_MIref1_TIMEref1 <- 
  geem(formula = NCIGS ~ cAGE + MI_ref1 + CBT + TIME_ref1 + MI_ref1 * CBT + MI_ref1 * TIME_ref1 + CBT * 
         TIME_ref1 + MI_ref1 * CBT * TIME_ref1 + offset(log_unctrldays),
       #ref: Mi=1, Cbt=0, time=1
       data = ncigs_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

ncigs_interaction_MIref1_TIMEref2 <- 
  geem(formula = NCIGS ~ cAGE + MI_ref1 + CBT + TIME_ref2 + MI_ref1 * CBT + MI_ref1 * TIME_ref2 + CBT * 
         TIME_ref2 + MI_ref1 * CBT * TIME_ref2 + offset(log_unctrldays),
       #ref: Mi=1, Cbt=0, time=2
       data = ncigs_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

ncigs_interaction_CBTref1 <- 
  geem(formula = NCIGS ~ cAGE + MI + CBT_ref1 + TIME + MI * CBT_ref1 + MI * TIME + CBT_ref1 * 
         TIME + MI * CBT_ref1 * TIME + offset(log_unctrldays),
       #ref: Mi=0, Cbt=1, time=0
       data = ncigs_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

ncigs_interaction_CBTref1_TIMEref1 <- 
  geem(formula = NCIGS ~ cAGE + MI + CBT_ref1 + TIME_ref1 + MI * CBT_ref1 + MI * TIME_ref1 + CBT_ref1 * 
         TIME_ref1 + MI * CBT_ref1 * TIME_ref1 + offset(log_unctrldays),
       #ref: Mi=0, Cbt=1, time=1
       data = ncigs_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

ncigs_interaction_CBTref1_TIMEref2 <- 
  geem(formula = NCIGS ~ cAGE + MI + CBT_ref1 + TIME_ref2 + MI * CBT_ref1 + MI * TIME_ref2 + CBT_ref1 * 
         TIME_ref2 + MI * CBT_ref1 * TIME_ref2 + offset(log_unctrldays),
       #ref: Mi=0, Cbt=1, time=2
       data = ncigs_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )


ncigs_interaction_MIref1_CBTref1 <- 
  geem(formula = NCIGS ~ cAGE + MI_ref1 + CBT_ref1 + TIME + MI_ref1 * CBT_ref1 + MI_ref1 * TIME + CBT_ref1 * 
         TIME + MI_ref1 * CBT_ref1 * TIME + offset(log_unctrldays),
       #ref: MI_ref1=1, Cbt=1, time=0
       data = ncigs_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

ncigs_interaction_MIref1_CBTref1_TIMEref1 <- 
  geem(formula = NCIGS ~ cAGE + MI_ref1 + CBT_ref1 + TIME_ref1 + MI_ref1 * CBT_ref1 + MI_ref1 * TIME_ref1 + CBT_ref1 * 
         TIME_ref1 + MI_ref1 * CBT_ref1 * TIME_ref1 + offset(log_unctrldays),
       #ref: MI_ref1=1, Cbt=1, time=1
       data = ncigs_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

ncigs_interaction_MIref1_CBTref1_TIMEref2 <- 
  geem(formula = NCIGS ~ cAGE + MI_ref1 + CBT_ref1 + TIME_ref2 + MI_ref1 * CBT_ref1 + MI_ref1 * TIME_ref2 + CBT_ref1 * 
         TIME_ref2 + MI_ref1 * CBT_ref1 * TIME_ref2 + offset(log_unctrldays),
       #ref: MI_ref1=1, Cbt=1, time=2
       data = ncigs_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )


## print model summary
summary(ncigs_interaction) #ref: MI=0, CBT=0, time=0
##                Estimates Model SE Robust SE    wald        p
## (Intercept)      2.32800  0.12470   0.11870 19.6000 0.000000
## cAGE             0.11290  0.04564   0.04463  2.5300 0.011400
## MI1              0.11850  0.17780   0.15330  0.7729 0.439600
## CBT1             0.18630  0.16950   0.15820  1.1780 0.239000
## TIME1           -0.48660  0.10920   0.17520 -2.7770 0.005493
## TIME2           -0.42420  0.14040   0.15860 -2.6760 0.007461
## MI1:CBT1        -0.27090  0.24160   0.20400 -1.3280 0.184200
## MI1:TIME1       -0.08472  0.15680   0.21690 -0.3905 0.696200
## MI1:TIME2        0.04657  0.20090   0.20570  0.2264 0.820900
## CBT1:TIME1      -0.19570  0.14940   0.21670 -0.9032 0.366400
## CBT1:TIME2      -0.21650  0.19140   0.21980 -0.9850 0.324600
## MI1:CBT1:TIME1   0.29440  0.21310   0.27600  1.0670 0.286100
## MI1:CBT1:TIME2   0.15550  0.27220   0.28010  0.5552 0.578800
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_TIMEref1) #ref: MI=0, CBT=0, time=1
##                     Estimates Model SE Robust SE     wald        p
## (Intercept)          1.841000  0.12670   0.14300 12.87000 0.000000
## cAGE                 0.112900  0.04564   0.04463  2.53000 0.011400
## MI1                  0.033810  0.18170   0.21440  0.15770 0.874700
## CBT1                -0.009413  0.17280   0.19130 -0.04920 0.960800
## TIME_ref10           0.486600  0.10920   0.17520  2.77700 0.005493
## TIME_ref12           0.062350  0.11390   0.06659  0.93630 0.349100
## MI1:CBT1             0.023550  0.24640   0.27490  0.08567 0.931700
## MI1:TIME_ref10       0.084720  0.15680   0.21690  0.39050 0.696200
## MI1:TIME_ref12       0.131300  0.16320   0.12910  1.01700 0.309300
## CBT1:TIME_ref10      0.195700  0.14940   0.21670  0.90320 0.366400
## CBT1:TIME_ref12     -0.020810  0.15500   0.10510 -0.19810 0.843000
## MI1:CBT1:TIME_ref10 -0.294400  0.21310   0.27600 -1.06700 0.286100
## MI1:CBT1:TIME_ref12 -0.138900  0.22030   0.17550 -0.79140 0.428700
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_TIMEref2) #ref: MI=0, CBT=0, time=2
##                     Estimates Model SE Robust SE    wald        p
## (Intercept)           1.90300  0.12940   0.12340 15.4200 0.000000
## cAGE                  0.11290  0.04564   0.04463  2.5300 0.011400
## MI1                   0.16510  0.18550   0.19410  0.8505 0.395100
## CBT1                 -0.03022  0.17620   0.18520 -0.1632 0.870300
## TIME_ref20            0.42420  0.14040   0.15860  2.6760 0.007461
## TIME_ref21           -0.06235  0.11390   0.06659 -0.9363 0.349100
## MI1:CBT1             -0.11540  0.25060   0.26350 -0.4377 0.661600
## MI1:TIME_ref20       -0.04657  0.20090   0.20570 -0.2264 0.820900
## MI1:TIME_ref21       -0.13130  0.16320   0.12910 -1.0170 0.309300
## CBT1:TIME_ref20       0.21650  0.19140   0.21980  0.9850 0.324600
## CBT1:TIME_ref21       0.02081  0.15500   0.10510  0.1981 0.843000
## MI1:CBT1:TIME_ref20  -0.15550  0.27220   0.28010 -0.5552 0.578800
## MI1:CBT1:TIME_ref21   0.13890  0.22030   0.17550  0.7914 0.428700
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_MIref1) #ref: MI=1, CBT=0, time=0
##                     Estimates Model SE Robust SE    wald         p
## (Intercept)           2.44600  0.12620   0.09717 25.1700 0.000e+00
## cAGE                  0.11290  0.04564   0.04463  2.5300 1.140e-02
## MI_ref10             -0.11850  0.17780   0.15330 -0.7729 4.396e-01
## CBT1                 -0.08460  0.17170   0.12750 -0.6634 5.071e-01
## TIME1                -0.57130  0.11250   0.12790 -4.4680 7.890e-06
## TIME2                -0.37770  0.14370   0.13110 -2.8800 3.978e-03
## MI_ref10:CBT1         0.27090  0.24160   0.20400  1.3280 1.842e-01
## MI_ref10:TIME1        0.08472  0.15680   0.21690  0.3905 6.962e-01
## MI_ref10:TIME2       -0.04657  0.20090   0.20570 -0.2264 8.209e-01
## CBT1:TIME1            0.09874  0.15200   0.17100  0.5774 5.636e-01
## CBT1:TIME2           -0.06098  0.19350   0.17370 -0.3510 7.256e-01
## MI_ref10:CBT1:TIME1  -0.29440  0.21310   0.27600 -1.0670 2.861e-01
## MI_ref10:CBT1:TIME2  -0.15550  0.27220   0.28010 -0.5552 5.788e-01
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, time=1
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)                1.87500  0.12970   0.16000 11.72000 0.000e+00
## cAGE                       0.11290  0.04564   0.04463  2.53000 1.140e-02
## MI_ref10                  -0.03381  0.18170   0.21440 -0.15770 8.747e-01
## CBT1                       0.01414  0.17530   0.19710  0.07175 9.428e-01
## TIME_ref10                 0.57130  0.11250   0.12790  4.46800 7.890e-06
## TIME_ref12                 0.19360  0.11690   0.11070  1.74900 8.037e-02
## MI_ref10:CBT1             -0.02355  0.24640   0.27490 -0.08567 9.317e-01
## MI_ref10:TIME_ref10       -0.08472  0.15680   0.21690 -0.39050 6.962e-01
## MI_ref10:TIME_ref12       -0.13130  0.16320   0.12910 -1.01700 3.093e-01
## CBT1:TIME_ref10           -0.09874  0.15200   0.17100 -0.57740 5.636e-01
## CBT1:TIME_ref12           -0.15970  0.15650   0.14060 -1.13600 2.561e-01
## MI_ref10:CBT1:TIME_ref10   0.29440  0.21310   0.27600  1.06700 2.861e-01
## MI_ref10:CBT1:TIME_ref12   0.13890  0.22030   0.17550  0.79140 4.287e-01
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_MIref1_TIMEref2) #ref: MI=1, CBT=0, time=2
##                          Estimates Model SE Robust SE    wald        p
## (Intercept)                2.06800  0.13240   0.14870 13.9100 0.000000
## cAGE                       0.11290  0.04564   0.04463  2.5300 0.011400
## MI_ref10                  -0.16510  0.18550   0.19410 -0.8505 0.395100
## CBT1                      -0.14560  0.17780   0.18550 -0.7847 0.432600
## TIME_ref20                 0.37770  0.14370   0.13110  2.8800 0.003978
## TIME_ref21                -0.19360  0.11690   0.11070 -1.7490 0.080370
## MI_ref10:CBT1              0.11540  0.25060   0.26350  0.4377 0.661600
## MI_ref10:TIME_ref20        0.04657  0.20090   0.20570  0.2264 0.820900
## MI_ref10:TIME_ref21        0.13130  0.16320   0.12910  1.0170 0.309300
## CBT1:TIME_ref20            0.06098  0.19350   0.17370  0.3510 0.725600
## CBT1:TIME_ref21            0.15970  0.15650   0.14060  1.1360 0.256100
## MI_ref10:CBT1:TIME_ref20   0.15550  0.27220   0.28010  0.5552 0.578800
## MI_ref10:CBT1:TIME_ref21  -0.13890  0.22030   0.17550 -0.7914 0.428700
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_CBTref1) #ref: MI=0, CBT=1, time=0
##                     Estimates Model SE Robust SE    wald         p
## (Intercept)            2.5140  0.11480   0.10260 24.5000 0.000e+00
## cAGE                   0.1129  0.04564   0.04463  2.5300 1.140e-02
## MI1                   -0.1524  0.16340   0.13150 -1.1580 2.467e-01
## CBT_ref10             -0.1863  0.16950   0.15820 -1.1780 2.390e-01
## TIME1                 -0.6823  0.10190   0.12750 -5.3530 9.000e-08
## TIME2                 -0.6407  0.13010   0.15230 -4.2080 2.575e-05
## MI1:CBT_ref10          0.2709  0.24160   0.20400  1.3280 1.842e-01
## MI1:TIME1              0.2097  0.14430   0.17070  1.2290 2.192e-01
## MI1:TIME2              0.2021  0.18370   0.19020  1.0630 2.880e-01
## CBT_ref10:TIME1        0.1957  0.14940   0.21670  0.9032 3.664e-01
## CBT_ref10:TIME2        0.2165  0.19140   0.21980  0.9850 3.246e-01
## MI1:CBT_ref10:TIME1   -0.2944  0.21310   0.27600 -1.0670 2.861e-01
## MI1:CBT_ref10:TIME2   -0.1555  0.27220   0.28010 -0.5552 5.788e-01
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, time=1
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)               1.832000  0.11740   0.12570 14.57000 0.000e+00
## cAGE                      0.112900  0.04564   0.04463  2.53000 1.140e-02
## MI1                       0.057360  0.16640   0.17030  0.33680 7.363e-01
## CBT_ref10                 0.009413  0.17280   0.19130  0.04920 9.608e-01
## TIME_ref10                0.682300  0.10190   0.12750  5.35300 9.000e-08
## TIME_ref12                0.041540  0.10520   0.08130  0.51100 6.094e-01
## MI1:CBT_ref10            -0.023550  0.24640   0.27490 -0.08567 9.317e-01
## MI1:TIME_ref10           -0.209700  0.14430   0.17070 -1.22900 2.192e-01
## MI1:TIME_ref12           -0.007612  0.14790   0.11890 -0.06404 9.489e-01
## CBT_ref10:TIME_ref10     -0.195700  0.14940   0.21670 -0.90320 3.664e-01
## CBT_ref10:TIME_ref12      0.020810  0.15500   0.10510  0.19810 8.430e-01
## MI1:CBT_ref10:TIME_ref10  0.294400  0.21310   0.27600  1.06700 2.861e-01
## MI1:CBT_ref10:TIME_ref12  0.138900  0.22030   0.17550  0.79140 4.287e-01
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_CBTref1_TIMEref2) #ref: MI=0, CBT=1, time=2
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)               1.873000  0.11950   0.13640 13.74000 0.000e+00
## cAGE                      0.112900  0.04564   0.04463  2.53000 1.140e-02
## MI1                       0.049750  0.16840   0.17560  0.28330 7.770e-01
## CBT_ref10                 0.030220  0.17620   0.18520  0.16320 8.703e-01
## TIME_ref20                0.640700  0.13010   0.15230  4.20800 2.575e-05
## TIME_ref21               -0.041540  0.10520   0.08130 -0.51100 6.094e-01
## MI1:CBT_ref10             0.115400  0.25060   0.26350  0.43770 6.616e-01
## MI1:TIME_ref20           -0.202100  0.18370   0.19020 -1.06300 2.880e-01
## MI1:TIME_ref21            0.007612  0.14790   0.11890  0.06404 9.489e-01
## CBT_ref10:TIME_ref20     -0.216500  0.19140   0.21980 -0.98500 3.246e-01
## CBT_ref10:TIME_ref21     -0.020810  0.15500   0.10510 -0.19810 8.430e-01
## MI1:CBT_ref10:TIME_ref20  0.155500  0.27220   0.28010  0.55520 5.788e-01
## MI1:CBT_ref10:TIME_ref21 -0.138900  0.22030   0.17550 -0.79140 4.287e-01
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, time=0
##                          Estimates Model SE Robust SE    wald         p
## (Intercept)                2.36200  0.11630   0.08222 28.7200 0.000e+00
## cAGE                       0.11290  0.04564   0.04463  2.5300 1.140e-02
## MI_ref10                   0.15240  0.16340   0.13150  1.1580 2.467e-01
## CBT_ref10                  0.08460  0.17170   0.12750  0.6634 5.071e-01
## TIME1                     -0.47260  0.10210   0.11350 -4.1630 3.145e-05
## TIME2                     -0.43860  0.12970   0.11400 -3.8470 1.194e-04
## MI_ref10:CBT_ref10        -0.27090  0.24160   0.20400 -1.3280 1.842e-01
## MI_ref10:TIME1            -0.20970  0.14430   0.17070 -1.2290 2.192e-01
## MI_ref10:TIME2            -0.20210  0.18370   0.19020 -1.0630 2.880e-01
## CBT_ref10:TIME1           -0.09874  0.15200   0.17100 -0.5774 5.636e-01
## CBT_ref10:TIME2            0.06098  0.19350   0.17370  0.3510 7.256e-01
## MI_ref10:CBT_ref10:TIME1   0.29440  0.21310   0.27600  1.0670 2.861e-01
## MI_ref10:CBT_ref10:TIME2   0.15550  0.27220   0.28010  0.5552 5.788e-01
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, time=1
##                               Estimates Model SE Robust SE     wald         p
## (Intercept)                    1.889000  0.11800   0.11490 16.43000 0.000e+00
## cAGE                           0.112900  0.04564   0.04463  2.53000 1.140e-02
## MI_ref10                      -0.057360  0.16640   0.17030 -0.33680 7.363e-01
## CBT_ref10                     -0.014140  0.17530   0.19710 -0.07175 9.428e-01
## TIME_ref10                     0.472600  0.10210   0.11350  4.16300 3.145e-05
## TIME_ref12                     0.033930  0.10400   0.08674  0.39120 6.956e-01
## MI_ref10:CBT_ref10             0.023550  0.24640   0.27490  0.08567 9.317e-01
## MI_ref10:TIME_ref10            0.209700  0.14430   0.17070  1.22900 2.192e-01
## MI_ref10:TIME_ref12            0.007612  0.14790   0.11890  0.06404 9.489e-01
## CBT_ref10:TIME_ref10           0.098740  0.15200   0.17100  0.57740 5.636e-01
## CBT_ref10:TIME_ref12           0.159700  0.15650   0.14060  1.13600 2.561e-01
## MI_ref10:CBT_ref10:TIME_ref10 -0.294400  0.21310   0.27600 -1.06700 2.861e-01
## MI_ref10:CBT_ref10:TIME_ref12 -0.138900  0.22030   0.17550 -0.79140 4.287e-01
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(ncigs_interaction_MIref1_CBTref1_TIMEref2) #ref: MI=1, CBT=1, time=2
##                               Estimates Model SE Robust SE     wald         p
## (Intercept)                    1.923000  0.11870   0.11070 17.37000 0.0000000
## cAGE                           0.112900  0.04564   0.04463  2.53000 0.0114000
## MI_ref10                      -0.049750  0.16840   0.17560 -0.28330 0.7770000
## CBT_ref10                      0.145600  0.17780   0.18550  0.78470 0.4326000
## TIME_ref20                     0.438600  0.12970   0.11400  3.84700 0.0001194
## TIME_ref21                    -0.033930  0.10400   0.08674 -0.39120 0.6956000
## MI_ref10:CBT_ref10            -0.115400  0.25060   0.26350 -0.43770 0.6616000
## MI_ref10:TIME_ref20            0.202100  0.18370   0.19020  1.06300 0.2880000
## MI_ref10:TIME_ref21           -0.007612  0.14790   0.11890 -0.06404 0.9489000
## CBT_ref10:TIME_ref20          -0.060980  0.19350   0.17370 -0.35100 0.7256000
## CBT_ref10:TIME_ref21          -0.159700  0.15650   0.14060 -1.13600 0.2561000
## MI_ref10:CBT_ref10:TIME_ref20 -0.155500  0.27220   0.28010 -0.55520 0.5788000
## MI_ref10:CBT_ref10:TIME_ref21  0.138900  0.22030   0.17550  0.79140 0.4287000
## 
##  Estimated Correlation Parameter:  0.624 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.001 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
# DV: SMDK  ---------------------------

#relevel MI and CBT
SMDK_dt_na.omit$MI_ref1 <- SMDK_dt_na.omit$MI %>% relevel("1")
SMDK_dt_na.omit$CBT_ref1 <- SMDK_dt_na.omit$CBT %>% relevel("1")
SMDK_dt_na.omit$TIME_ref1 <- SMDK_dt_na.omit$TIME %>% relevel("1")
SMDK_dt_na.omit$TIME_ref2 <- SMDK_dt_na.omit$TIME %>% relevel("2")

# ignore missing values
## create an indicator column for any missing values 
SMDK_dt_na.omit$any_na <- SMDK_dt_na.omit %>% apply(1, function(x){any(is.na(x))}) 

## left join the dataset by the id column using the group_by() function
SMDK_dt_na.omit %<>% left_join(SMDK_dt_na.omit %>% group_by(id) %>% 
                               summarise(any_na2 = any(any_na)), by="id")


## filter out rows with any NAs
SMDK_dt_na.omit %<>% filter(any_na2 != T)



## fit main model
SMDK_main <- geem(formula = SMDK.main,
                 data = SMDK_dt_na.omit,
                 id=id,
                 family = MASS::negative.binomial(1),
                 corstr = corstr
)

## print main model summary
summary(SMDK_main)
##             Estimates Model SE Robust SE   wald         p
## (Intercept)  -0.08032  0.03220   0.02432 -3.303 0.0009558
## cAGE          0.02530  0.01492   0.01339  1.890 0.0588100
## MI1           0.05554  0.03317   0.02938  1.890 0.0587300
## CBT1         -0.03473  0.03321   0.02935 -1.183 0.2367000
## TIME1        -0.21080  0.02223   0.02664 -7.912 0.0000000
## TIME2        -0.19200  0.02703   0.02752 -6.977 0.0000000
## 
##  Estimated Correlation Parameter:  0.4749 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1209 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
## fit interaction model
SMDK_interaction <- geem(formula = SMDK.interaction,
                  data = SMDK_dt_na.omit,
                  id=id,
                  family = MASS::negative.binomial(1),
                  corstr = corstr
)

SMDK_interaction_TIMEref1 <- geem(formula = 
                                     SMDK ~ cAGE + MI + CBT + TIME_ref1 + MI * CBT + MI * TIME_ref1 + CBT * 
                                     TIME_ref1 + MI * CBT * TIME_ref1 + offset(log_unctrldays),
                                   #ref: Mi=0, Cbt=0, time=1
                                   data = SMDK_dt_na.omit,
                                   id=id,
                                   family = MASS::negative.binomial(1),
                                  corstr = corstr
)



SMDK_interaction_TIMEref2 <- geem(formula = 
                                     SMDK ~ cAGE + MI + CBT + TIME_ref2 + MI * CBT + MI * TIME_ref2 + CBT * 
                                     TIME_ref2 + MI * CBT * TIME_ref2 + offset(log_unctrldays),
                                   #ref: Mi=0, Cbt=0, time=1
                                   data = SMDK_dt_na.omit,
                                   id=id,
                                   family = MASS::negative.binomial(1),
                                  corstr = corstr
)

SMDK_interaction_MIref1 <- 
  geem(formula = SMDK ~ cAGE + MI_ref1 + CBT + TIME + MI_ref1 * CBT + MI_ref1 * TIME + CBT * 
         TIME + MI_ref1 * CBT * TIME + offset(log_unctrldays),
       #ref: Mi=1, Cbt=0, time=0
       data = SMDK_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

SMDK_interaction_MIref1_TIMEref1 <- 
  geem(formula = SMDK ~ cAGE + MI_ref1 + CBT + TIME_ref1 + MI_ref1 * CBT + MI_ref1 * TIME_ref1 + CBT * 
         TIME_ref1 + MI_ref1 * CBT * TIME_ref1 + offset(log_unctrldays),
       #ref: Mi=1, Cbt=0, time=1
       data = SMDK_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

SMDK_interaction_MIref1_TIMEref2 <- 
  geem(formula = SMDK ~ cAGE + MI_ref1 + CBT + TIME_ref2 + MI_ref1 * CBT + MI_ref1 * TIME_ref2 + CBT * 
         TIME_ref2 + MI_ref1 * CBT * TIME_ref2 + offset(log_unctrldays),
       #ref: Mi=1, Cbt=0, time=2
       data = SMDK_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

SMDK_interaction_CBTref1 <- 
  geem(formula = SMDK ~ cAGE + MI + CBT_ref1 + TIME + MI * CBT_ref1 + MI * TIME + CBT_ref1 * 
         TIME + MI * CBT_ref1 * TIME + offset(log_unctrldays),
       #ref: Mi=0, Cbt=1, time=0
       data = SMDK_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

SMDK_interaction_CBTref1_TIMEref1 <- 
  geem(formula = SMDK ~ cAGE + MI + CBT_ref1 + TIME_ref1 + MI * CBT_ref1 + MI * TIME_ref1 + CBT_ref1 * 
         TIME_ref1 + MI * CBT_ref1 * TIME_ref1 + offset(log_unctrldays),
       #ref: Mi=0, Cbt=1, time=1
       data = SMDK_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

SMDK_interaction_CBTref1_TIMEref2 <- 
  geem(formula = SMDK ~ cAGE + MI + CBT_ref1 + TIME_ref2 + MI * CBT_ref1 + MI * TIME_ref2 + CBT_ref1 * 
         TIME_ref2 + MI * CBT_ref1 * TIME_ref2 + offset(log_unctrldays),
       #ref: Mi=0, Cbt=1, time=2
       data = SMDK_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )


SMDK_interaction_MIref1_CBTref1 <- 
  geem(formula = SMDK ~ cAGE + MI_ref1 + CBT_ref1 + TIME + MI_ref1 * CBT_ref1 + MI_ref1 * TIME + CBT_ref1 * 
         TIME + MI_ref1 * CBT_ref1 * TIME + offset(log_unctrldays),
       #ref: MI_ref1=1, Cbt=1, time=0
       data = SMDK_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

SMDK_interaction_MIref1_CBTref1_TIMEref1 <- 
  geem(formula = SMDK ~ cAGE + MI_ref1 + CBT_ref1 + TIME_ref1 + MI_ref1 * CBT_ref1 + MI_ref1 * TIME_ref1 + CBT_ref1 * 
         TIME_ref1 + MI_ref1 * CBT_ref1 * TIME_ref1 + offset(log_unctrldays),
       #ref: MI_ref1=1, Cbt=1, time=1
       data = SMDK_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

SMDK_interaction_MIref1_CBTref1_TIMEref2 <- 
  geem(formula = SMDK ~ cAGE + MI_ref1 + CBT_ref1 + TIME_ref2 + MI_ref1 * CBT_ref1 + MI_ref1 * TIME_ref2 + CBT_ref1 * 
         TIME_ref2 + MI_ref1 * CBT_ref1 * TIME_ref2 + offset(log_unctrldays),
       #ref: MI_ref1=1, Cbt=1, time=2
       data = SMDK_dt_na.omit,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

## print interaction model summary
summary(SMDK_interaction) #ref: MI=0, CBT=0, time=0
##                 Estimates Model SE Robust SE      wald         p
## (Intercept)    -0.0916800  0.04436   0.02252 -4.071000 4.688e-05
## cAGE            0.0265900  0.01515   0.01375  1.933000 5.321e-02
## MI1             0.0211400  0.06328   0.02733  0.773500 4.392e-01
## CBT1            0.0298700  0.06032   0.02805  1.065000 2.870e-01
## TIME1          -0.2095000  0.04589   0.05909 -3.546000 3.917e-04
## TIME2          -0.1404000  0.05625   0.05370 -2.614000 8.961e-03
## MI1:CBT1       -0.0254400  0.08595   0.03711 -0.685600 4.930e-01
## MI1:TIME1       0.0311200  0.06584   0.07422  0.419300 6.750e-01
## MI1:TIME2       0.0313600  0.08037   0.06463  0.485200 6.275e-01
## CBT1:TIME1     -0.0304000  0.06265   0.07843 -0.387600 6.983e-01
## CBT1:TIME2     -0.2042000  0.07665   0.08605 -2.373000 1.765e-02
## MI1:CBT1:TIME1 -0.0008054  0.08935   0.10550 -0.007637 9.939e-01
## MI1:CBT1:TIME2  0.1546000  0.10890   0.10650  1.451000 1.468e-01
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_TIMEref1) #ref: MI=0, CBT=0, time=1
##                      Estimates Model SE Robust SE      wald         p
## (Intercept)         -0.3012000  0.04556   0.05909 -5.098000 3.400e-07
## cAGE                 0.0265900  0.01515   0.01375  1.933000 5.321e-02
## MI1                  0.0522500  0.06541   0.07627  0.685100 4.933e-01
## CBT1                -0.0005307  0.06209   0.07857 -0.006755 9.946e-01
## TIME_ref10           0.2095000  0.04589   0.05909  3.546000 3.917e-04
## TIME_ref12           0.0691600  0.04770   0.05043  1.371000 1.702e-01
## MI1:CBT1            -0.0262500  0.08856   0.10580 -0.248100 8.041e-01
## MI1:TIME_ref10      -0.0311200  0.06584   0.07422 -0.419300 6.750e-01
## MI1:TIME_ref12       0.0002409  0.06835   0.06348  0.003795 9.970e-01
## CBT1:TIME_ref10      0.0304000  0.06265   0.07843  0.387600 6.983e-01
## CBT1:TIME_ref12     -0.1738000  0.06493   0.06856 -2.535000 1.124e-02
## MI1:CBT1:TIME_ref10  0.0008054  0.08935   0.10550  0.007637 9.939e-01
## MI1:CBT1:TIME_ref12  0.1554000  0.09225   0.09248  1.680000 9.296e-02
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_TIMEref2) #ref: MI=0, CBT=0, time=2
##                      Estimates Model SE Robust SE      wald         p
## (Intercept)         -0.2320000  0.04656   0.05286 -4.390000 1.133e-05
## cAGE                 0.0265900  0.01515   0.01375  1.933000 5.321e-02
## MI1                  0.0524900  0.06671   0.06560  0.800200 4.236e-01
## CBT1                -0.1743000  0.06342   0.08520 -2.046000 4.076e-02
## TIME_ref20           0.1404000  0.05625   0.05370  2.614000 8.961e-03
## TIME_ref21          -0.0691600  0.04770   0.05043 -1.371000 1.702e-01
## MI1:CBT1             0.1291000  0.09004   0.10470  1.233000 2.176e-01
## MI1:TIME_ref20      -0.0313600  0.08037   0.06463 -0.485200 6.275e-01
## MI1:TIME_ref21      -0.0002409  0.06835   0.06348 -0.003795 9.970e-01
## CBT1:TIME_ref20      0.2042000  0.07665   0.08605  2.373000 1.765e-02
## CBT1:TIME_ref21      0.1738000  0.06493   0.06856  2.535000 1.124e-02
## MI1:CBT1:TIME_ref20 -0.1546000  0.10890   0.10650 -1.451000 1.468e-01
## MI1:CBT1:TIME_ref21 -0.1554000  0.09225   0.09248 -1.680000 9.296e-02
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_MIref1) #ref: MI=1, CBT=0, time=0
##                      Estimates Model SE Robust SE      wald         p
## (Intercept)         -0.0705500  0.04494   0.01407 -5.014000 5.300e-07
## cAGE                 0.0265900  0.01515   0.01375  1.933000 5.321e-02
## MI_ref10            -0.0211400  0.06328   0.02733 -0.773500 4.392e-01
## CBT1                 0.0044220  0.06109   0.02328  0.190000 8.493e-01
## TIME1               -0.1784000  0.04721   0.04483 -3.980000 6.905e-05
## TIME2               -0.1090000  0.05741   0.03597 -3.030000 2.446e-03
## MI_ref10:CBT1        0.0254400  0.08595   0.03711  0.685600 4.930e-01
## MI_ref10:TIME1      -0.0311200  0.06584   0.07422 -0.419300 6.750e-01
## MI_ref10:TIME2      -0.0313600  0.08037   0.06463 -0.485200 6.275e-01
## CBT1:TIME1          -0.0312000  0.06371   0.07048 -0.442700 6.580e-01
## CBT1:TIME2          -0.0496200  0.07738   0.06280 -0.790300 4.294e-01
## MI_ref10:CBT1:TIME1  0.0008054  0.08935   0.10550  0.007637 9.939e-01
## MI_ref10:CBT1:TIME2 -0.1546000  0.10890   0.10650 -1.451000 1.468e-01
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, time=1
##                           Estimates Model SE Robust SE      wald         p
## (Intercept)              -0.2490000  0.04680   0.04807 -5.179000 2.200e-07
## cAGE                      0.0265900  0.01515   0.01375  1.933000 5.321e-02
## MI_ref10                 -0.0522500  0.06541   0.07627 -0.685100 4.933e-01
## CBT1                     -0.0267800  0.06305   0.07063 -0.379200 7.046e-01
## TIME_ref10                0.1784000  0.04721   0.04483  3.980000 6.905e-05
## TIME_ref12                0.0694000  0.04896   0.03854  1.801000 7.170e-02
## MI_ref10:CBT1             0.0262500  0.08856   0.10580  0.248100 8.041e-01
## MI_ref10:TIME_ref10       0.0311200  0.06584   0.07422  0.419300 6.750e-01
## MI_ref10:TIME_ref12      -0.0002409  0.06835   0.06348 -0.003795 9.970e-01
## CBT1:TIME_ref10           0.0312000  0.06371   0.07048  0.442700 6.580e-01
## CBT1:TIME_ref12          -0.0184200  0.06553   0.06205 -0.296900 7.666e-01
## MI_ref10:CBT1:TIME_ref10 -0.0008054  0.08935   0.10550 -0.007637 9.939e-01
## MI_ref10:CBT1:TIME_ref12 -0.1554000  0.09225   0.09248 -1.680000 9.296e-02
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_MIref1_TIMEref2) #ref: MI=1, CBT=0, time=2
##                           Estimates Model SE Robust SE      wald         p
## (Intercept)              -0.1795000  0.04761   0.03848 -4.665000 3.080e-06
## cAGE                      0.0265900  0.01515   0.01375  1.933000 5.321e-02
## MI_ref10                 -0.0524900  0.06671   0.06560 -0.800200 4.236e-01
## CBT1                     -0.0452000  0.06381   0.06037 -0.748800 4.540e-01
## TIME_ref20                0.1090000  0.05741   0.03597  3.030000 2.446e-03
## TIME_ref21               -0.0694000  0.04896   0.03854 -1.801000 7.170e-02
## MI_ref10:CBT1            -0.1291000  0.09004   0.10470 -1.233000 2.176e-01
## MI_ref10:TIME_ref20       0.0313600  0.08037   0.06463  0.485200 6.275e-01
## MI_ref10:TIME_ref21       0.0002409  0.06835   0.06348  0.003795 9.970e-01
## CBT1:TIME_ref20           0.0496200  0.07738   0.06280  0.790300 4.294e-01
## CBT1:TIME_ref21           0.0184200  0.06553   0.06205  0.296900 7.666e-01
## MI_ref10:CBT1:TIME_ref20  0.1546000  0.10890   0.10650  1.451000 1.468e-01
## MI_ref10:CBT1:TIME_ref21  0.1554000  0.09225   0.09248  1.680000 9.296e-02
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_CBTref1) #ref: MI=0, CBT=1, time=0
##                      Estimates Model SE Robust SE      wald         p
## (Intercept)         -0.0618200  0.04087   0.01628 -3.797000 1.464e-04
## cAGE                 0.0265900  0.01515   0.01375  1.933000 5.321e-02
## MI1                 -0.0043090  0.05816   0.02469 -0.174500 8.615e-01
## CBT_ref10           -0.0298700  0.06032   0.02805 -1.065000 2.870e-01
## TIME1               -0.2399000  0.04265   0.05164 -4.646000 3.390e-06
## TIME2               -0.3445000  0.05207   0.06732 -5.118000 3.100e-07
## MI1:CBT_ref10        0.0254400  0.08595   0.03711  0.685600 4.930e-01
## MI1:TIME1            0.0303100  0.06041   0.07497  0.404300 6.860e-01
## MI1:TIME2            0.1859000  0.07350   0.08471  2.195000 2.818e-02
## CBT_ref10:TIME1      0.0304000  0.06265   0.07843  0.387600 6.983e-01
## CBT_ref10:TIME2      0.2042000  0.07665   0.08605  2.373000 1.765e-02
## MI1:CBT_ref10:TIME1  0.0008054  0.08935   0.10550  0.007637 9.939e-01
## MI1:CBT_ref10:TIME2 -0.1546000  0.10890   0.10650 -1.451000 1.468e-01
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, time=1
##                           Estimates Model SE Robust SE      wald         p
## (Intercept)              -0.3017000  0.04215   0.05122 -5.891000 0.000e+00
## cAGE                      0.0265900  0.01515   0.01375  1.933000 5.321e-02
## MI1                       0.0260000  0.05969   0.07283  0.357000 7.211e-01
## CBT_ref10                 0.0005307  0.06209   0.07857  0.006755 9.946e-01
## TIME_ref10                0.2399000  0.04265   0.05164  4.646000 3.390e-06
## TIME_ref12               -0.1046000  0.04405   0.04644 -2.253000 2.425e-02
## MI1:CBT_ref10             0.0262500  0.08856   0.10580  0.248100 8.041e-01
## MI1:TIME_ref10           -0.0303100  0.06041   0.07497 -0.404300 6.860e-01
## MI1:TIME_ref12            0.1556000  0.06195   0.06724  2.314000 2.066e-02
## CBT_ref10:TIME_ref10     -0.0304000  0.06265   0.07843 -0.387600 6.983e-01
## CBT_ref10:TIME_ref12      0.1738000  0.06493   0.06856  2.535000 1.124e-02
## MI1:CBT_ref10:TIME_ref10 -0.0008054  0.08935   0.10550 -0.007637 9.939e-01
## MI1:CBT_ref10:TIME_ref12 -0.1554000  0.09225   0.09248 -1.680000 9.296e-02
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_CBTref1_TIMEref2) #ref: MI=0, CBT=1, time=2
##                          Estimates Model SE Robust SE   wald         p
## (Intercept)               -0.40640  0.04301   0.06625 -6.134 0.000e+00
## cAGE                       0.02659  0.01515   0.01375  1.933 5.321e-02
## MI1                        0.18160  0.06046   0.08095  2.243 2.487e-02
## CBT_ref10                  0.17430  0.06342   0.08520  2.046 4.076e-02
## TIME_ref20                 0.34450  0.05207   0.06732  5.118 3.100e-07
## TIME_ref21                 0.10460  0.04405   0.04644  2.253 2.425e-02
## MI1:CBT_ref10             -0.12910  0.09004   0.10470 -1.233 2.176e-01
## MI1:TIME_ref20            -0.18590  0.07350   0.08471 -2.195 2.818e-02
## MI1:TIME_ref21            -0.15560  0.06195   0.06724 -2.314 2.066e-02
## CBT_ref10:TIME_ref20      -0.20420  0.07665   0.08605 -2.373 1.765e-02
## CBT_ref10:TIME_ref21      -0.17380  0.06493   0.06856 -2.535 1.124e-02
## MI1:CBT_ref10:TIME_ref20   0.15460  0.10890   0.10650  1.451 1.468e-01
## MI1:CBT_ref10:TIME_ref21   0.15540  0.09225   0.09248  1.680 9.296e-02
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, time=0
##                           Estimates Model SE Robust SE      wald         p
## (Intercept)              -0.0661300  0.04138   0.01856 -3.563000 0.0003670
## cAGE                      0.0265900  0.01515   0.01375  1.933000 0.0532100
## MI_ref10                  0.0043090  0.05816   0.02469  0.174500 0.8615000
## CBT_ref10                -0.0044220  0.06109   0.02328 -0.190000 0.8493000
## TIME1                    -0.2096000  0.04278   0.05438 -3.855000 0.0001159
## TIME2                    -0.1586000  0.05188   0.05149 -3.080000 0.0020670
## MI_ref10:CBT_ref10       -0.0254400  0.08595   0.03711 -0.685600 0.4930000
## MI_ref10:TIME1           -0.0303100  0.06041   0.07497 -0.404300 0.6860000
## MI_ref10:TIME2           -0.1859000  0.07350   0.08471 -2.195000 0.0281800
## CBT_ref10:TIME1           0.0312000  0.06371   0.07048  0.442700 0.6580000
## CBT_ref10:TIME2           0.0496200  0.07738   0.06280  0.790300 0.4294000
## MI_ref10:CBT_ref10:TIME1 -0.0008054  0.08935   0.10550 -0.007637 0.9939000
## MI_ref10:CBT_ref10:TIME2  0.1546000  0.10890   0.10650  1.451000 0.1468000
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, time=1
##                                Estimates Model SE Robust SE      wald         p
## (Intercept)                   -0.2757000  0.04227   0.05182 -5.321000 0.0000001
## cAGE                           0.0265900  0.01515   0.01375  1.933000 0.0532100
## MI_ref10                      -0.0260000  0.05969   0.07283 -0.357000 0.7211000
## CBT_ref10                      0.0267800  0.06305   0.07063  0.379200 0.7046000
## TIME_ref10                     0.2096000  0.04278   0.05438  3.855000 0.0001159
## TIME_ref12                     0.0509800  0.04356   0.04864  1.048000 0.2946000
## MI_ref10:CBT_ref10            -0.0262500  0.08856   0.10580 -0.248100 0.8041000
## MI_ref10:TIME_ref10            0.0303100  0.06041   0.07497  0.404300 0.6860000
## MI_ref10:TIME_ref12           -0.1556000  0.06195   0.06724 -2.314000 0.0206600
## CBT_ref10:TIME_ref10          -0.0312000  0.06371   0.07048 -0.442700 0.6580000
## CBT_ref10:TIME_ref12           0.0184200  0.06553   0.06205  0.296900 0.7666000
## MI_ref10:CBT_ref10:TIME_ref10  0.0008054  0.08935   0.10550  0.007637 0.9939000
## MI_ref10:CBT_ref10:TIME_ref12  0.1554000  0.09225   0.09248  1.680000 0.0929600
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
summary(SMDK_interaction_MIref1_CBTref1_TIMEref2) #ref: MI=1, CBT=1, time=2
##                               Estimates Model SE Robust SE    wald         p
## (Intercept)                    -0.22470  0.04250   0.04663 -4.8200 1.430e-06
## cAGE                            0.02659  0.01515   0.01375  1.9330 5.321e-02
## MI_ref10                       -0.18160  0.06046   0.08095 -2.2430 2.487e-02
## CBT_ref10                       0.04520  0.06381   0.06037  0.7488 4.540e-01
## TIME_ref20                      0.15860  0.05188   0.05149  3.0800 2.067e-03
## TIME_ref21                     -0.05098  0.04356   0.04864 -1.0480 2.946e-01
## MI_ref10:CBT_ref10              0.12910  0.09004   0.10470  1.2330 2.176e-01
## MI_ref10:TIME_ref20             0.18590  0.07350   0.08471  2.1950 2.818e-02
## MI_ref10:TIME_ref21             0.15560  0.06195   0.06724  2.3140 2.066e-02
## CBT_ref10:TIME_ref20           -0.04962  0.07738   0.06280 -0.7903 4.294e-01
## CBT_ref10:TIME_ref21           -0.01842  0.06553   0.06205 -0.2969 7.666e-01
## MI_ref10:CBT_ref10:TIME_ref20  -0.15460  0.10890   0.10650 -1.4510 1.468e-01
## MI_ref10:CBT_ref10:TIME_ref21  -0.15540  0.09225   0.09248 -1.6800 9.296e-02
## 
##  Estimated Correlation Parameter:  0.4817 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1228 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  791
# DV = LONGABS ---------------------------

## fit main model

#relevel MI and CBT
LONGABS_nobl_dt$MI_ref1 <- LONGABS_nobl_dt$MI %>% relevel("1")
LONGABS_nobl_dt$CBT_ref1 <- LONGABS_nobl_dt$CBT %>% relevel("1")
LONGABS_nobl_dt$TIME_ref1 <- LONGABS_nobl_dt$TIME %>% relevel("1")

CLONGABS_nobl_main <- geem(formula = 
                             CLONGABS ~ cAGE + cBLCGSMD + MI + CBT + TIME,
                          data = LONGABS_nobl_dt,
                          id=id,
                          family = MASS::negative.binomial(1),
                          corstr = corstr
)
  

## print main model summary
summary(CLONGABS_nobl_main) #unlogged
##             Estimates Model SE Robust SE     wald        p
## (Intercept)  1.725000  0.23570   0.22090  7.80600 0.000000
## cAGE        -0.128900  0.11610   0.11420 -1.12900 0.258900
## cBLCGSMD    -0.001244  0.01455   0.01597 -0.07790 0.937900
## MI1         -0.747300  0.25370   0.24840 -3.00800 0.002628
## CBT1         0.666200  0.25430   0.24460  2.72400 0.006453
## TIME1        0.008899  0.13590   0.13410  0.06635 0.947100
## 
##  Estimated Correlation Parameter:  0.5516 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.821 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
## fit interaction models

CLONGABS_nobl_interaction <- geem(formula = CLONGABS.nobl.interaction,
                           data = LONGABS_nobl_dt,
                           id=id,
                           family = MASS::negative.binomial(1),
                           corstr = corstr
)



## print interaction models' summary
summary(CLONGABS_nobl_interaction) #unlogged
##                Estimates Model SE Robust SE    wald       p
## (Intercept)     -1.93700  0.24060   0.22920 -8.4500 0.00000
## cAGE            -0.17350  0.09237   0.10340 -1.6780 0.09328
## cBLCGSMD        -0.01200  0.01158   0.01789 -0.6710 0.50220
## MI1             -0.54560  0.35440   0.36740 -1.4850 0.13750
## CBT1             0.31990  0.32360   0.29100  1.0990 0.27160
## TIME1           -0.07626  0.24310   0.23590 -0.3233 0.74650
## MI1:CBT1         0.32140  0.46990   0.47210  0.6809 0.49600
## MI1:TIME1       -0.47490  0.37020   0.40590 -1.1700 0.24200
## CBT1:TIME1       0.32360  0.32630   0.27670  1.1690 0.24220
## MI1:CBT1:TIME1  -0.06697  0.48430   0.50950 -0.1315 0.89540
## 
##  Estimated Correlation Parameter:  0.5029 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  2.94 
## 
##  Number of GEE iterations: 17 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
# DV = LONGABS ---------------------------

# #relevel MI and CBT
LONGABS_nobl_dt$MI_ref1 <- LONGABS_nobl_dt$MI %>% relevel("1")
LONGABS_nobl_dt$CBT_ref1 <- LONGABS_nobl_dt$CBT %>% relevel("1")
LONGABS_nobl_dt$TIME_ref1 <- LONGABS_nobl_dt$TIME %>% relevel("1")

# ignore missing values
## create an indicator column for any missing values 
LONGABS_nobl_dt$any_na <- LONGABS_nobl_dt %>% apply(1, function(x){any(is.na(x))}) 

## left join the dataset by the id column using the group_by() function
LONGABS_nobl_dt %<>% left_join(LONGABS_nobl_dt %>% group_by(id) %>% 
                               summarise(any_na2 = any(any_na)), by="id")


## filter out rows with any NAs
LONGABS_nobl_dt %<>% filter(any_na2 != T)

## main model
LONGABS_nobl_main <- geem(formula = CLONGABS.nobl.main,
                        data = LONGABS_nobl_dt,
                        id=id,
                        family = MASS::negative.binomial(1),
                        corstr = corstr
)
## Warning in geem(formula = CLONGABS.nobl.main, data = LONGABS_nobl_dt, id = id,
## : Did not converge
## print main model summary
summary(LONGABS_nobl_main)
##             Estimates Model SE Robust SE    wald        p
## (Intercept)  -1.97100  0.18990   0.20800 -9.4740 0.000000
## cAGE         -0.14630  0.09258   0.10710 -1.3660 0.171900
## cBLCGSMD     -0.01598  0.01151   0.02026 -0.7886 0.430400
## MI1          -0.59970  0.20440   0.21530 -2.7860 0.005342
## CBT1          0.60390  0.20600   0.22100  2.7330 0.006283
## TIME1        -0.12680  0.11960   0.11920 -1.0640 0.287500
## 
##  Estimated Correlation Parameter:  0.4962 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  2.955 
## 
##  Number of GEE iterations: 20 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
## interaction model
LONGABS_nobl_interaction <- geem(formula = CLONGABS ~ cAGE + cBLCGSMD + MI + CBT + TIME + MI * CBT + MI * TIME +
                                   CBT * TIME + MI * CBT * TIME,
                               #mi=0, cbt=0, time=0
                               data = LONGABS_nobl_dt,
                               id=id,
                               family = MASS::negative.binomial(1),
                               corstr = corstr
)

LONGABS_nobl_interaction_TIMEref1 <- geem(formula =
                                          CLONGABS ~ cAGE + cBLCGSMD + MI + CBT + TIME_ref1 + MI * CBT + MI * TIME_ref1 +
                                          CBT * TIME_ref1 + MI * CBT * TIME_ref1,
                                        #mi=0, cbt=0, time=1
                                        data = LONGABS_nobl_dt,
                                        id=id,
                                        family = MASS::negative.binomial(1),
                                        corstr = corstr
)


LONGABS_nobl_interaction_MIref1 <-
  geem(formula = CLONGABS ~ cAGE + cBLCGSMD + MI_ref1 + CBT + TIME + MI_ref1 * CBT + MI_ref1 * TIME +
         CBT * TIME + MI_ref1 * CBT * TIME,
       #ref: Mi=1, Cbt=0, Time=0
       data = LONGABS_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

LONGABS_nobl_interaction_MIref1_TIMEref1 <-
  geem(formula = CLONGABS ~ cAGE + cBLCGSMD + MI_ref1 + CBT + TIME_ref1 +
         MI_ref1 * CBT + MI_ref1 * TIME_ref1 +
         CBT * TIME_ref1 + MI_ref1 * CBT * TIME_ref1,
       #ref: Mi=1, Cbt=0, Time=1
       data = LONGABS_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

LONGABS_nobl_interaction_CBTref1 <-
  geem(formula = CLONGABS ~ cAGE + cBLCGSMD + MI + CBT_ref1 + TIME + MI * CBT_ref1 + MI *
         TIME + CBT_ref1 * TIME + MI * CBT_ref1 * TIME,
       #ref: Mi=0, Cbt=1, Time = 0
       data = LONGABS_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

LONGABS_nobl_interaction_CBTref1_TIMEref1 <-
  geem(formula = CLONGABS ~ cAGE + cBLCGSMD + MI + CBT_ref1 + TIME_ref1 +
         MI * CBT_ref1 + MI * TIME_ref1 + CBT_ref1 * TIME_ref1 + 
         MI * CBT_ref1 * TIME_ref1,
       #ref: Mi=0, Cbt=1, Time = 1
       data = LONGABS_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

LONGABS_nobl_interaction_MIref1_CBTref1 <-
  geem(formula = CLONGABS ~ cAGE + cBLCGSMD + MI_ref1 + CBT_ref1 + TIME + MI_ref1 * CBT_ref1 + MI_ref1 *
         TIME + CBT_ref1 * TIME + MI_ref1 * CBT_ref1 * TIME,
       #ref: Mi=1, Cbt=1, Time=0
       data = LONGABS_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )

LONGABS_nobl_interaction_MIref1_CBTref1_TIMEref1 <-
  geem(formula = CLONGABS ~ cAGE + cBLCGSMD + MI_ref1 + CBT_ref1 + TIME_ref1 +
         MI_ref1 * CBT_ref1 + MI_ref1 * TIME_ref1 +
         CBT_ref1 * TIME_ref1 + MI_ref1 * CBT_ref1 * TIME_ref1,
       #ref: Mi=1, Cbt=1, Time=1
       data = LONGABS_nobl_dt,
       id=id,
       family = MASS::negative.binomial(1),
       corstr = corstr
  )


## print interaction model summary
summary(LONGABS_nobl_interaction) #base: MI=0, CBT=0, Time=0
##                 Estimates Model SE Robust SE     wald       p
## (Intercept)     1.6550000  0.29820   0.25250  6.55500 0.00000
## cAGE           -0.1379000  0.11640   0.11550 -1.19500 0.23230
## cBLCGSMD        0.0007003  0.01457   0.01521  0.04605 0.96330
## MI1            -0.5604000  0.43700   0.37890 -1.47900 0.13920
## CBT1            0.5545000  0.40080   0.31870  1.74000 0.08184
## TIME1           0.2209000  0.28000   0.26050  0.84810 0.39640
## MI1:CBT1        0.2036000  0.58100   0.51850  0.39270 0.69460
## MI1:TIME1      -0.6803000  0.42240   0.47810 -1.42300 0.15470
## CBT1:TIME1      0.0351800  0.37770   0.30920  0.11380 0.90940
## MI1:CBT1:TIME1  0.2018000  0.55760   0.58310  0.34600 0.72930
## 
##  Estimated Correlation Parameter:  0.545 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.796 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(LONGABS_nobl_interaction_TIMEref1) #base: MI=0, CBT=0, Time=1
##                      Estimates Model SE Robust SE     wald       p
## (Intercept)          1.8760000  0.29340   0.30820  6.08800 0.00000
## cAGE                -0.1379000  0.11640   0.11550 -1.19500 0.23230
## cBLCGSMD             0.0007003  0.01457   0.01521  0.04605 0.96330
## MI1                 -1.2410000  0.45210   0.53510 -2.31900 0.02041
## CBT1                 0.5897000  0.39510   0.35560  1.65800 0.09723
## TIME_ref10          -0.2209000  0.28000   0.26050 -0.84810 0.39640
## MI1:CBT1             0.4054000  0.59250   0.66560  0.60900 0.54250
## MI1:TIME_ref10       0.6803000  0.42240   0.47810  1.42300 0.15470
## CBT1:TIME_ref10     -0.0351800  0.37770   0.30920 -0.11380 0.90940
## MI1:CBT1:TIME_ref10 -0.2018000  0.55760   0.58310 -0.34600 0.72930
## 
##  Estimated Correlation Parameter:  0.545 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.796 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(LONGABS_nobl_interaction_MIref1) #ref: MI=1, CBT=0, Time=0
##                      Estimates Model SE Robust SE     wald         p
## (Intercept)          1.0950000  0.31880   0.27900  3.92400 8.727e-05
## cAGE                -0.1379000  0.11640   0.11550 -1.19500 2.323e-01
## cBLCGSMD             0.0007003  0.01457   0.01521  0.04605 9.633e-01
## MI_ref10             0.5604000  0.43700   0.37890  1.47900 1.392e-01
## CBT1                 0.7581000  0.41910   0.40670  1.86400 6.228e-02
## TIME1               -0.4594000  0.31630   0.40090 -1.14600 2.518e-01
## MI_ref10:CBT1       -0.2036000  0.58100   0.51850 -0.39270 6.946e-01
## MI_ref10:TIME1       0.6803000  0.42240   0.47810  1.42300 1.547e-01
## CBT1:TIME1           0.2369000  0.41020   0.49430  0.47930 6.317e-01
## MI_ref10:CBT1:TIME1 -0.2018000  0.55760   0.58310 -0.34600 7.293e-01
## 
##  Estimated Correlation Parameter:  0.545 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.796 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(LONGABS_nobl_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, Time=1
##                           Estimates Model SE Robust SE     wald       p
## (Intercept)               0.6352000  0.34350   0.42150  1.50700 0.13180
## cAGE                     -0.1379000  0.11640   0.11550 -1.19500 0.23230
## cBLCGSMD                  0.0007003  0.01457   0.01521  0.04605 0.96330
## MI_ref10                  1.2410000  0.45210   0.53510  2.31900 0.02041
## CBT1                      0.9951000  0.44020   0.55310  1.79900 0.07202
## TIME_ref10                0.4594000  0.31630   0.40090  1.14600 0.25180
## MI_ref10:CBT1            -0.4054000  0.59250   0.66560 -0.60900 0.54250
## MI_ref10:TIME_ref10      -0.6803000  0.42240   0.47810 -1.42300 0.15470
## CBT1:TIME_ref10          -0.2369000  0.41020   0.49430 -0.47930 0.63170
## MI_ref10:CBT1:TIME_ref10  0.2018000  0.55760   0.58310  0.34600 0.72930
## 
##  Estimated Correlation Parameter:  0.545 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.796 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(LONGABS_nobl_interaction_CBTref1) #ref: MI=0, CBT=1, Time=0
##                      Estimates Model SE Robust SE     wald       p
## (Intercept)          2.2100000  0.26730   0.19520 11.32000 0.00000
## cAGE                -0.1379000  0.11640   0.11550 -1.19500 0.23230
## cBLCGSMD             0.0007003  0.01457   0.01521  0.04605 0.96330
## MI1                 -0.3568000  0.38190   0.34760 -1.02700 0.30470
## CBT_ref10           -0.5545000  0.40080   0.31870 -1.74000 0.08184
## TIME1                0.2561000  0.25350   0.16670  1.53600 0.12450
## MI1:CBT_ref10       -0.2036000  0.58100   0.51850 -0.39270 0.69460
## MI1:TIME1           -0.4786000  0.36400   0.33300 -1.43700 0.15070
## CBT_ref10:TIME1     -0.0351800  0.37770   0.30920 -0.11380 0.90940
## MI1:CBT_ref10:TIME1 -0.2018000  0.55760   0.58310 -0.34600 0.72930
## 
##  Estimated Correlation Parameter:  0.545 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.796 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(LONGABS_nobl_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, Time=1
##                           Estimates Model SE Robust SE     wald       p
## (Intercept)               2.4660000  0.26420   0.17860 13.81000 0.00000
## cAGE                     -0.1379000  0.11640   0.11550 -1.19500 0.23230
## cBLCGSMD                  0.0007003  0.01457   0.01521  0.04605 0.96330
## MI1                      -0.8354000  0.38220   0.38850 -2.15000 0.03153
## CBT_ref10                -0.5897000  0.39510   0.35560 -1.65800 0.09723
## TIME_ref10               -0.2561000  0.25350   0.16670 -1.53600 0.12450
## MI1:CBT_ref10            -0.4054000  0.59250   0.66560 -0.60900 0.54250
## MI1:TIME_ref10            0.4786000  0.36400   0.33300  1.43700 0.15070
## CBT_ref10:TIME_ref10      0.0351800  0.37770   0.30920  0.11380 0.90940
## MI1:CBT_ref10:TIME_ref10  0.2018000  0.55760   0.58310  0.34600 0.72930
## 
##  Estimated Correlation Parameter:  0.545 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.796 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(LONGABS_nobl_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, Time=0
##                           Estimates Model SE Robust SE     wald       p
## (Intercept)               1.8530000  0.27250   0.28890  6.41200 0.00000
## cAGE                     -0.1379000  0.11640   0.11550 -1.19500 0.23230
## cBLCGSMD                  0.0007003  0.01457   0.01521  0.04605 0.96330
## MI_ref10                  0.3568000  0.38190   0.34760  1.02700 0.30470
## CBT_ref10                -0.7581000  0.41910   0.40670 -1.86400 0.06228
## TIME1                    -0.2225000  0.26110   0.28850 -0.77130 0.44050
## MI_ref10:CBT_ref10        0.2036000  0.58100   0.51850  0.39270 0.69460
## MI_ref10:TIME1            0.4786000  0.36400   0.33300  1.43700 0.15070
## CBT_ref10:TIME1          -0.2369000  0.41020   0.49430 -0.47930 0.63170
## MI_ref10:CBT_ref10:TIME1  0.2018000  0.55760   0.58310  0.34600 0.72930
## 
##  Estimated Correlation Parameter:  0.545 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.796 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
summary(LONGABS_nobl_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, Time=1
##                                Estimates Model SE Robust SE     wald         p
## (Intercept)                    1.6300000  0.27590   0.34740  4.69300 0.0000027
## cAGE                          -0.1379000  0.11640   0.11550 -1.19500 0.2323000
## cBLCGSMD                       0.0007003  0.01457   0.01521  0.04605 0.9633000
## MI_ref10                       0.8354000  0.38220   0.38850  2.15000 0.0315300
## CBT_ref10                     -0.9951000  0.44020   0.55310 -1.79900 0.0720200
## TIME_ref10                     0.2225000  0.26110   0.28850  0.77130 0.4405000
## MI_ref10:CBT_ref10             0.4054000  0.59250   0.66560  0.60900 0.5425000
## MI_ref10:TIME_ref10           -0.4786000  0.36400   0.33300 -1.43700 0.1507000
## CBT_ref10:TIME_ref10           0.2369000  0.41020   0.49430  0.47930 0.6317000
## MI_ref10:CBT_ref10:TIME_ref10 -0.2018000  0.55760   0.58310 -0.34600 0.7293000
## 
##  Estimated Correlation Parameter:  0.545 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.796 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  278    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  556
# Save ---------------------------

save.image(file="gee-main-interaction-models.RData")