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
# 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")


# 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")

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


## print main model summary
summary(AVCIG_nobl_main)
##             Estimates Model SE Robust SE    wald       p
## (Intercept)   1.77700 0.089150   0.11250 15.7900 0.00000
## cAGE         -0.02810 0.039800   0.05864 -0.4791 0.63180
## cBLCGSMD      0.04724 0.004623   0.00517  9.1360 0.00000
## MI1           0.04773 0.087870   0.11300  0.4224 0.67270
## CBT1         -0.03620 0.088180   0.11260 -0.3215 0.74790
## TIME1         0.09405 0.087700   0.04884  1.9260 0.05413
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.8529 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  274    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  520
## 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)
)

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)
)


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)
  )

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)
  )

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)
  )

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)
  )

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)
  )

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)
  )


## print interaction model summary
summary(AVCIG_nobl_interaction) #base: MI=0, CBT=0, Time=0
##                Estimates Model SE Robust SE    wald      p
## (Intercept)      1.87900 0.126500  0.141500 13.2800 0.0000
## cAGE            -0.02613 0.039720  0.058870 -0.4440 0.6571
## cBLCGSMD         0.04854 0.004616  0.005141  9.4420 0.0000
## MI1             -0.19320 0.181500  0.184500 -1.0470 0.2949
## CBT1            -0.18940 0.172800  0.177700 -1.0660 0.2866
## TIME1            0.01963 0.181500  0.096760  0.2029 0.8392
## MI1:CBT1         0.35650 0.246100  0.242700  1.4690 0.1417
## MI1:TIME1        0.20920 0.259800  0.150000  1.3940 0.1632
## CBT1:TIME1       0.05884 0.246800  0.127500  0.4614 0.6445
## MI1:CBT1:TIME1  -0.21420 0.350900  0.196100 -1.0920 0.2748
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.844 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  274    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  520
summary(AVCIG_nobl_interaction_TIMEref1) #base: MI=0, CBT=0, Time=1
##                     Estimates Model SE Robust SE     wald      p
## (Intercept)           1.89900 0.130800  0.120700 15.72000 0.0000
## cAGE                 -0.02613 0.039720  0.058870 -0.44400 0.6571
## cBLCGSMD              0.04854 0.004616  0.005141  9.44200 0.0000
## MI1                   0.01600 0.187100  0.186300  0.08587 0.9316
## CBT1                 -0.13050 0.177300  0.169300 -0.77090 0.4408
## TIME_ref10           -0.01963 0.181500  0.096760 -0.20290 0.8392
## MI1:CBT1              0.14230 0.251800  0.247900  0.57410 0.5659
## MI1:TIME_ref10       -0.20920 0.259800  0.150000 -1.39400 0.1632
## CBT1:TIME_ref10      -0.05884 0.246800  0.127500 -0.46140 0.6445
## MI1:CBT1:TIME_ref10   0.21420 0.350900  0.196100  1.09200 0.2748
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.844 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  274    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  520
summary(AVCIG_nobl_interaction_MIref1) #ref: MI=1, CBT=0, Time=0
##                     Estimates Model SE Robust SE   wald       p
## (Intercept)           1.68600 0.129800  0.117300 14.370 0.00000
## cAGE                 -0.02613 0.039720  0.058870 -0.444 0.65710
## cBLCGSMD              0.04854 0.004616  0.005141  9.442 0.00000
## MI_ref10              0.19320 0.181500  0.184500  1.047 0.29490
## CBT1                  0.16720 0.174800  0.164700  1.015 0.31010
## TIME1                 0.22890 0.185800  0.113900  2.009 0.04456
## MI_ref10:CBT1        -0.35650 0.246100  0.242700 -1.469 0.14170
## MI_ref10:TIME1       -0.20920 0.259800  0.150000 -1.394 0.16320
## CBT1:TIME1           -0.15540 0.249400  0.148300 -1.048 0.29460
## MI_ref10:CBT1:TIME1   0.21420 0.350900  0.196100  1.092 0.27480
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.844 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  274    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  520
summary(AVCIG_nobl_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, Time=1
##                          Estimates Model SE Robust SE     wald       p
## (Intercept)                1.91500 0.133300  0.139000 13.77000 0.00000
## cAGE                      -0.02613 0.039720  0.058870 -0.44400 0.65710
## cBLCGSMD                   0.04854 0.004616  0.005141  9.44200 0.00000
## MI_ref10                  -0.01600 0.187100  0.186300 -0.08587 0.93160
## CBT1                       0.01177 0.178100  0.179500  0.06557 0.94770
## TIME_ref10                -0.22890 0.185800  0.113900 -2.00900 0.04456
## MI_ref10:CBT1             -0.14230 0.251800  0.247900 -0.57410 0.56590
## MI_ref10:TIME_ref10        0.20920 0.259800  0.150000  1.39400 0.16320
## CBT1:TIME_ref10            0.15540 0.249400  0.148300  1.04800 0.29460
## MI_ref10:CBT1:TIME_ref10  -0.21420 0.350900  0.196100 -1.09200 0.27480
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.844 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  274    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  520
summary(AVCIG_nobl_interaction_CBTref1) #ref: MI=0, CBT=1, Time=0
##                     Estimates Model SE Robust SE     wald      p
## (Intercept)          1.690000 0.117500  0.109600 15.42000 0.0000
## cAGE                -0.026130 0.039720  0.058870 -0.44400 0.6571
## cBLCGSMD             0.048540 0.004616  0.005141  9.44200 0.0000
## MI1                  0.163300 0.165800  0.159100  1.02600 0.3048
## CBT_ref10            0.189400 0.172800  0.177700  1.06600 0.2866
## TIME1                0.078470 0.167100  0.083100  0.94440 0.3450
## MI1:CBT_ref10       -0.356500 0.246100  0.242700 -1.46900 0.1417
## MI1:TIME1           -0.004995 0.235800  0.126200 -0.03958 0.9684
## CBT_ref10:TIME1     -0.058840 0.246800  0.127500 -0.46140 0.6445
## MI1:CBT_ref10:TIME1  0.214200 0.350900  0.196100  1.09200 0.2748
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.844 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  274    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  520
summary(AVCIG_nobl_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, Time=1
##                          Estimates Model SE Robust SE     wald      p
## (Intercept)               1.768000 0.119300  0.119500 14.79000 0.0000
## cAGE                     -0.026130 0.039720  0.058870 -0.44400 0.6571
## cBLCGSMD                  0.048540 0.004616  0.005141  9.44200 0.0000
## MI1                       0.158300 0.167900  0.164700  0.96110 0.3365
## CBT_ref10                 0.130500 0.177300  0.169300  0.77090 0.4408
## TIME_ref10               -0.078470 0.167100  0.083100 -0.94440 0.3450
## MI1:CBT_ref10            -0.142300 0.251800  0.247900 -0.57410 0.5659
## MI1:TIME_ref10            0.004995 0.235800  0.126200  0.03958 0.9684
## CBT_ref10:TIME_ref10      0.058840 0.246800  0.127500  0.46140 0.6445
## MI1:CBT_ref10:TIME_ref10 -0.214200 0.350900  0.196100 -1.09200 0.2748
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.844 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  274    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  520
summary(AVCIG_nobl_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, Time=0
##                          Estimates Model SE Robust SE     wald      p
## (Intercept)               1.853000 0.117000  0.116100 15.96000 0.0000
## cAGE                     -0.026130 0.039720  0.058870 -0.44400 0.6571
## cBLCGSMD                  0.048540 0.004616  0.005141  9.44200 0.0000
## MI_ref10                 -0.163300 0.165800  0.159100 -1.02600 0.3048
## CBT_ref10                -0.167200 0.174800  0.164700 -1.01500 0.3101
## TIME1                     0.073480 0.166300  0.094920  0.77410 0.4389
## MI_ref10:CBT_ref10        0.356500 0.246100  0.242700  1.46900 0.1417
## MI_ref10:TIME1            0.004995 0.235800  0.126200  0.03958 0.9684
## CBT_ref10:TIME1           0.155400 0.249400  0.148300  1.04800 0.2946
## MI_ref10:CBT_ref10:TIME1 -0.214200 0.350900  0.196100 -1.09200 0.2748
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.844 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  274    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  520
summary(AVCIG_nobl_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, Time=1
##                               Estimates Model SE Robust SE     wald      p
## (Intercept)                    1.926000 0.118200  0.114900 16.77000 0.0000
## cAGE                          -0.026130 0.039720  0.058870 -0.44400 0.6571
## cBLCGSMD                       0.048540 0.004616  0.005141  9.44200 0.0000
## MI_ref10                      -0.158300 0.167900  0.164700 -0.96110 0.3365
## CBT_ref10                     -0.011770 0.178100  0.179500 -0.06557 0.9477
## TIME_ref10                    -0.073480 0.166300  0.094920 -0.77410 0.4389
## MI_ref10:CBT_ref10             0.142300 0.251800  0.247900  0.57410 0.5659
## MI_ref10:TIME_ref10           -0.004995 0.235800  0.126200 -0.03958 0.9684
## CBT_ref10:TIME_ref10          -0.155400 0.249400  0.148300 -1.04800 0.2946
## MI_ref10:CBT_ref10:TIME_ref10  0.214200 0.350900  0.196100  1.09200 0.2748
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.844 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  274    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  520
# 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")

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


## print main model summary
summary(CGSMD_nobl_main)
##             Estimates Model SE Robust SE     wald       p
## (Intercept)  1.956000 0.073540  0.091300 21.42000 0.00000
## cAGE        -0.004485 0.033130  0.048740 -0.09201 0.92670
## cBLCGSMD     0.040940 0.003931  0.004366  9.37700 0.00000
## MI1         -0.005918 0.072870  0.092190 -0.06419 0.94880
## CBT1         0.015750 0.073040  0.092660  0.17000 0.86500
## TIME1        0.101100 0.072850  0.042640  2.37000 0.01778
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.562 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  268    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  484
## 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)
)

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)
)


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)
  )

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)
  )

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)
  )

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)
  )

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)
  )

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)
  )


## print interaction model summary
summary(CGSMD_nobl_interaction) #base: MI=0, CBT=0, time=0
##                Estimates Model SE Robust SE     wald       p
## (Intercept)     2.045000 0.104800  0.116600 17.53000 0.00000
## cAGE           -0.002582 0.033080  0.048820 -0.05289 0.95780
## cBLCGSMD        0.041710 0.003925  0.004326  9.64200 0.00000
## MI1            -0.181000 0.149600  0.153900 -1.17600 0.23970
## CBT1           -0.147300 0.142700  0.144400 -1.02000 0.30780
## TIME1          -0.001181 0.150800  0.085310 -0.01385 0.98900
## MI1:CBT1        0.311500 0.203500  0.200900  1.55000 0.12110
## MI1:TIME1       0.190600 0.213400  0.135600  1.40600 0.15970
## CBT1:TIME1      0.181600 0.207000  0.110800  1.63800 0.10140
## MI1:CBT1:TIME1 -0.328200 0.291100  0.172100 -1.90700 0.05653
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  268    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  484
summary(CGSMD_nobl_interaction_TIMEref1) #base: MI=0, CBT=0, time=1
##                     Estimates Model SE Robust SE     wald       p
## (Intercept)          2.044000 0.108700  0.102400 19.96000 0.00000
## cAGE                -0.002582 0.033080  0.048820 -0.05289 0.95780
## cBLCGSMD             0.041710 0.003925  0.004326  9.64200 0.00000
## MI1                  0.009634 0.153000  0.159000  0.06058 0.95170
## CBT1                 0.034240 0.150800  0.135700  0.25230 0.80080
## TIME_ref10           0.001181 0.150800  0.085310  0.01385 0.98900
## MI1:CBT1            -0.016640 0.209300  0.205900 -0.08080 0.93560
## MI1:TIME_ref10      -0.190600 0.213400  0.135600 -1.40600 0.15970
## CBT1:TIME_ref10     -0.181600 0.207000  0.110800 -1.63800 0.10140
## MI1:CBT1:TIME_ref10  0.328200 0.291100  0.172100  1.90700 0.05653
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  268    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  484
summary(CGSMD_nobl_interaction_MIref1) #ref: MI=1, CBT=0, time=0
##                     Estimates Model SE Robust SE     wald       p
## (Intercept)          1.864000 0.106500  0.100300 18.58000 0.00000
## cAGE                -0.002582 0.033080  0.048820 -0.05289 0.95780
## cBLCGSMD             0.041710 0.003925  0.004326  9.64200 0.00000
## MI_ref10             0.181000 0.149600  0.153900  1.17600 0.23970
## CBT1                 0.164200 0.144800  0.140100  1.17200 0.24110
## TIME1                0.189400 0.151000  0.104600  1.81100 0.07020
## MI_ref10:CBT1       -0.311500 0.203500  0.200900 -1.55000 0.12110
## MI_ref10:TIME1      -0.190600 0.213400  0.135600 -1.40600 0.15970
## CBT1:TIME1          -0.146600 0.204600  0.131100 -1.11900 0.26330
## MI_ref10:CBT1:TIME1  0.328200 0.291100  0.172100  1.90700 0.05653
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  268    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  484
summary(CGSMD_nobl_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, time=1
##                          Estimates Model SE Robust SE     wald       p
## (Intercept)               2.053000 0.107300  0.120500 17.04000 0.00000
## cAGE                     -0.002582 0.033080  0.048820 -0.05289 0.95780
## cBLCGSMD                  0.041710 0.003925  0.004326  9.64200 0.00000
## MI_ref10                 -0.009634 0.153000  0.159000 -0.06058 0.95170
## CBT1                      0.017600 0.144700  0.153800  0.11440 0.90890
## TIME_ref10               -0.189400 0.151000  0.104600 -1.81100 0.07020
## MI_ref10:CBT1             0.016640 0.209300  0.205900  0.08080 0.93560
## MI_ref10:TIME_ref10       0.190600 0.213400  0.135600  1.40600 0.15970
## CBT1:TIME_ref10           0.146600 0.204600  0.131100  1.11900 0.26330
## MI_ref10:CBT1:TIME_ref10 -0.328200 0.291100  0.172100 -1.90700 0.05653
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  268    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  484
summary(CGSMD_nobl_interaction_CBTref1) #ref: MI=0, CBT=1, time=0
##                     Estimates Model SE Robust SE     wald       p
## (Intercept)          1.898000 0.096630  0.086580 21.92000 0.00000
## cAGE                -0.002582 0.033080  0.048820 -0.05289 0.95780
## cBLCGSMD             0.041710 0.003925  0.004326  9.64200 0.00000
## MI1                  0.130500 0.137700  0.130500  1.00000 0.31730
## CBT_ref10            0.147300 0.142700  0.144400  1.02000 0.30780
## TIME1                0.180400 0.141900  0.070860  2.54600 0.01091
## MI1:CBT_ref10       -0.311500 0.203500  0.200900 -1.55000 0.12110
## MI1:TIME1           -0.137500 0.197900  0.106100 -1.29600 0.19490
## CBT_ref10:TIME1     -0.181600 0.207000  0.110800 -1.63800 0.10140
## MI1:CBT_ref10:TIME1  0.328200 0.291100  0.172100  1.90700 0.05653
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  268    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  484
summary(CGSMD_nobl_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, time=1
##                          Estimates Model SE Robust SE     wald       p
## (Intercept)               2.078000 0.104200  0.089790 23.14000 0.00000
## cAGE                     -0.002582 0.033080  0.048820 -0.05289 0.95780
## cBLCGSMD                  0.041710 0.003925  0.004326  9.64200 0.00000
## MI1                      -0.007003 0.142500  0.131000 -0.05345 0.95740
## CBT_ref10                -0.034240 0.150800  0.135700 -0.25230 0.80080
## TIME_ref10               -0.180400 0.141900  0.070860 -2.54600 0.01091
## MI1:CBT_ref10             0.016640 0.209300  0.205900  0.08080 0.93560
## MI1:TIME_ref10            0.137500 0.197900  0.106100  1.29600 0.19490
## CBT_ref10:TIME_ref10      0.181600 0.207000  0.110800  1.63800 0.10140
## MI1:CBT_ref10:TIME_ref10 -0.328200 0.291100  0.172100 -1.90700 0.05653
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  268    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  484
summary(CGSMD_nobl_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, time=0
##                          Estimates Model SE Robust SE     wald       p
## (Intercept)               2.028000 0.098080  0.098190 20.65000 0.00000
## cAGE                     -0.002582 0.033080  0.048820 -0.05289 0.95780
## cBLCGSMD                  0.041710 0.003925  0.004326  9.64200 0.00000
## MI_ref10                 -0.130500 0.137700  0.130500 -1.00000 0.31730
## CBT_ref10                -0.164200 0.144800  0.140100 -1.17200 0.24110
## TIME1                     0.042840 0.138000  0.079050  0.54200 0.58790
## MI_ref10:CBT_ref10        0.311500 0.203500  0.200900  1.55000 0.12110
## MI_ref10:TIME1            0.137500 0.197900  0.106100  1.29600 0.19490
## CBT_ref10:TIME1           0.146600 0.204600  0.131100  1.11900 0.26330
## MI_ref10:CBT_ref10:TIME1 -0.328200 0.291100  0.172100 -1.90700 0.05653
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  268    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  484
summary(CGSMD_nobl_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, time=1
##                               Estimates Model SE Robust SE     wald       p
## (Intercept)                    2.071000 0.097140  0.096410 21.48000 0.00000
## cAGE                          -0.002582 0.033080  0.048820 -0.05289 0.95780
## cBLCGSMD                       0.041710 0.003925  0.004326  9.64200 0.00000
## MI_ref10                       0.007003 0.142500  0.131000  0.05345 0.95740
## CBT_ref10                     -0.017600 0.144700  0.153800 -0.11440 0.90890
## TIME_ref10                    -0.042840 0.138000  0.079050 -0.54200 0.58790
## MI_ref10:CBT_ref10            -0.016640 0.209300  0.205900 -0.08080 0.93560
## MI_ref10:TIME_ref10           -0.137500 0.197900  0.106100 -1.29600 0.19490
## CBT_ref10:TIME_ref10          -0.146600 0.204600  0.131100 -1.11900 0.26330
## MI_ref10:CBT_ref10:TIME_ref10  0.328200 0.291100  0.172100  1.90700 0.05653
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  268    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  484
# 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
## fit main model
PCSMKD_nobl_main <- geem(formula = PCSMKD.nobl.main,
                       data = PCSMKD_nobl_dt,
                       id=id,
                       family = MASS::negative.binomial(1)
)


## print main model summary
summary(PCSMKD_nobl_main) #ref: Mi=0, Cbt=0
##             Estimates Model SE Robust SE    wald       p
## (Intercept)  4.355000 0.038350  0.043910 99.2000 0.00000
## cAGE         0.015680 0.017230  0.022700  0.6909 0.48960
## cBLCGSMD     0.006944 0.002055  0.002978  2.3320 0.01972
## MI1          0.088370 0.037970  0.046690  1.8930 0.05838
## CBT1        -0.071120 0.038100  0.046540 -1.5280 0.12650
## TIME1        0.010900 0.037910  0.024880  0.4380 0.66140
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1828 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
## 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)
)

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)
  )


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


## print main model summaries
summary(PCSMKD_MI_interactions)
##             Estimates Model SE Robust SE    wald       p
## (Intercept)  4.376000 0.042560  0.044970 97.3300 0.00000
## cAGE         0.015420 0.017260  0.022740  0.6780 0.49780
## cBLCGSMD     0.006929 0.002059  0.002985  2.3210 0.02028
## MI1          0.046410 0.053010  0.052430  0.8853 0.37600
## CBT1        -0.072190 0.038160  0.046590 -1.5500 0.12120
## TIME1       -0.032010 0.053510  0.036980 -0.8655 0.38680
## MI1:TIME1    0.086340 0.075940  0.049790  1.7340 0.08287
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1835 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
summary(PCSMKD_MI_interactions_MIref1_CBT_ref0)
##                Estimates Model SE Robust SE     wald       p
## (Intercept)     4.423000 0.043090  0.041270 107.2000 0.00000
## cAGE            0.015420 0.017260  0.022740   0.6780 0.49780
## cBLCGSMD        0.006929 0.002059  0.002985   2.3210 0.02028
## MI_ref10       -0.046410 0.053010  0.052430  -0.8853 0.37600
## TIME1           0.054340 0.053890  0.033310   1.6310 0.10280
## CBT1           -0.072190 0.038160  0.046590  -1.5500 0.12120
## MI_ref10:TIME1 -0.086340 0.075940  0.049790  -1.7340 0.08287
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1835 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
summary(PCSMKD_CBT_interactions)
##             Estimates Model SE Robust SE    wald       p
## (Intercept)  4.329000  0.04344  0.047380 91.3600 0.00000
## cAGE         0.015270  0.01727  0.022710  0.6723 0.50140
## cBLCGSMD     0.006982  0.00206  0.002979  2.3440 0.01910
## MI1          0.089610  0.03805  0.046750  1.9170 0.05526
## CBT1        -0.023990  0.05322  0.052940 -0.4532 0.65040
## TIME1        0.063770  0.05641  0.034260  1.8610 0.06270
## CBT1:TIME1  -0.096820  0.07631  0.049260 -1.9660 0.04935
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1836 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
## 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)
  )

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)
  )

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)
  )

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)
  )

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)
  )

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)
  )

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)
  )

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)
  )

## print interaction model summary
summary(PCSMKD_nobl_interaction) #base: MI=0, CBT=0, time=0
##                Estimates Model SE Robust SE     wald       p
## (Intercept)     4.349000 0.055160  0.058660 74.14000 0.00000
## cAGE            0.016360 0.017380  0.022890  0.71470 0.47480
## cBLCGSMD        0.007199 0.002074  0.002995  2.40400 0.01624
## MI1             0.048660 0.078950  0.075900  0.64110 0.52150
## CBT1           -0.021210 0.075230  0.079020 -0.26840 0.78840
## TIME1           0.066670 0.079100  0.053380  1.24900 0.21160
## MI1:CBT1       -0.005168 0.107100  0.105500 -0.04898 0.96090
## MI1:TIME1      -0.006297 0.113300  0.067540 -0.09324 0.92570
## CBT1:TIME1     -0.186700 0.107800  0.073480 -2.54100 0.01106
## MI1:CBT1:TIME1  0.175000 0.153200  0.097850  1.78900 0.07362
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.185 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
summary(PCSMKD_nobl_interaction_TIMEref1) #base: MI=0, CBT=0, time=1
##                     Estimates Model SE Robust SE     wald       p
## (Intercept)          4.415000 0.057020  0.052760 83.68000 0.00000
## cAGE                 0.016360 0.017380  0.022890  0.71470 0.47480
## cBLCGSMD             0.007199 0.002074  0.002995  2.40400 0.01624
## MI1                  0.042360 0.081790  0.066150  0.64040 0.52190
## CBT1                -0.207900 0.077670  0.086540 -2.40200 0.01630
## TIME_ref10          -0.066670 0.079100  0.053380 -1.24900 0.21160
## MI1:CBT1             0.169900 0.110300  0.105900  1.60400 0.10870
## MI1:TIME_ref10       0.006297 0.113300  0.067540  0.09324 0.92570
## CBT1:TIME_ref10      0.186700 0.107800  0.073480  2.54100 0.01106
## MI1:CBT1:TIME_ref10 -0.175000 0.153200  0.097850 -1.78900 0.07362
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.185 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
summary(PCSMKD_nobl_interaction_MIref1) #ref: MI=1, CBT=0, time=0
##                     Estimates Model SE Robust SE     wald       p
## (Intercept)          4.397000 0.056360  0.048050 91.52000 0.00000
## cAGE                 0.016360 0.017380  0.022890  0.71470 0.47480
## cBLCGSMD             0.007199 0.002074  0.002995  2.40400 0.01624
## MI_ref10            -0.048660 0.078950  0.075900 -0.64110 0.52150
## CBT1                -0.026370 0.076030  0.069800 -0.37790 0.70550
## TIME1                0.060380 0.081110  0.041420  1.45800 0.14500
## MI_ref10:CBT1        0.005168 0.107100  0.105500  0.04898 0.96090
## MI_ref10:TIME1       0.006297 0.113300  0.067540  0.09324 0.92570
## CBT1:TIME1          -0.011640 0.108900  0.064650 -0.18010 0.85710
## MI_ref10:CBT1:TIME1 -0.175000 0.153200  0.097850 -1.78900 0.07362
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.185 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
summary(PCSMKD_nobl_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, time=1
##                          Estimates Model SE Robust SE      wald       p
## (Intercept)               4.458000 0.058410  0.039160 113.80000 0.00000
## cAGE                      0.016360 0.017380  0.022890   0.71470 0.47480
## cBLCGSMD                  0.007199 0.002074  0.002995   2.40400 0.01624
## MI_ref10                 -0.042360 0.081790  0.066150  -0.64040 0.52190
## CBT1                     -0.038020 0.078030  0.060480  -0.62860 0.52960
## TIME_ref10               -0.060380 0.081110  0.041420  -1.45800 0.14500
## MI_ref10:CBT1            -0.169900 0.110300  0.105900  -1.60400 0.10870
## MI_ref10:TIME_ref10      -0.006297 0.113300  0.067540  -0.09324 0.92570
## CBT1:TIME_ref10           0.011640 0.108900  0.064650   0.18010 0.85710
## MI_ref10:CBT1:TIME_ref10  0.175000 0.153200  0.097850   1.78900 0.07362
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.185 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
summary(PCSMKD_nobl_interaction_CBTref1) #ref: MI=0, CBT=1, time=0
##                     Estimates Model SE Robust SE     wald       p
## (Intercept)          4.328000 0.051050  0.052200 82.91000 0.00000
## cAGE                 0.016360 0.017380  0.022890  0.71470 0.47480
## cBLCGSMD             0.007199 0.002074  0.002995  2.40400 0.01624
## MI1                  0.043490 0.072190  0.072590  0.59910 0.54910
## CBT_ref10            0.021210 0.075230  0.079020  0.26840 0.78840
## TIME1               -0.120000 0.073250  0.050700 -2.36700 0.01793
## MI1:CBT_ref10        0.005168 0.107100  0.105500  0.04898 0.96090
## MI1:TIME1            0.168700 0.103200  0.070980  2.37800 0.01743
## CBT_ref10:TIME1      0.186700 0.107800  0.073480  2.54100 0.01106
## MI1:CBT_ref10:TIME1 -0.175000 0.153200  0.097850 -1.78900 0.07362
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.185 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
summary(PCSMKD_nobl_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, time=1
##                          Estimates Model SE Robust SE    wald        p
## (Intercept)               4.208000 0.052590  0.067960 61.9100 0.000000
## cAGE                      0.016360 0.017380  0.022890  0.7147 0.474800
## cBLCGSMD                  0.007199 0.002074  0.002995  2.4040 0.016240
## MI1                       0.212200 0.073810  0.082030  2.5870 0.009673
## CBT_ref10                 0.207900 0.077670  0.086540  2.4020 0.016300
## TIME_ref10                0.120000 0.073250  0.050700  2.3670 0.017930
## MI1:CBT_ref10            -0.169900 0.110300  0.105900 -1.6040 0.108700
## MI1:TIME_ref10           -0.168700 0.103200  0.070980 -2.3780 0.017430
## CBT_ref10:TIME_ref10     -0.186700 0.107800  0.073480 -2.5410 0.011060
## MI1:CBT_ref10:TIME_ref10  0.175000 0.153200  0.097850  1.7890 0.073620
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.185 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
summary(PCSMKD_nobl_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, time=0
##                          Estimates Model SE Robust SE     wald       p
## (Intercept)               4.371000 0.051020  0.050660 86.29000 0.00000
## cAGE                      0.016360 0.017380  0.022890  0.71470 0.47480
## cBLCGSMD                  0.007199 0.002074  0.002995  2.40400 0.01624
## MI_ref10                 -0.043490 0.072190  0.072590 -0.59910 0.54910
## CBT_ref10                 0.026370 0.076030  0.069800  0.37790 0.70550
## TIME1                     0.048730 0.072650  0.049700  0.98050 0.32680
## MI_ref10:CBT_ref10       -0.005168 0.107100  0.105500 -0.04898 0.96090
## MI_ref10:TIME1           -0.168700 0.103200  0.070980 -2.37800 0.01743
## CBT_ref10:TIME1           0.011640 0.108900  0.064650  0.18010 0.85710
## MI_ref10:CBT_ref10:TIME1  0.175000 0.153200  0.097850  1.78900 0.07362
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.185 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
summary(PCSMKD_nobl_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, time=1
##                               Estimates Model SE Robust SE    wald        p
## (Intercept)                    4.420000 0.051740  0.046370 95.3100 0.000000
## cAGE                           0.016360 0.017380  0.022890  0.7147 0.474800
## cBLCGSMD                       0.007199 0.002074  0.002995  2.4040 0.016240
## MI_ref10                      -0.212200 0.073810  0.082030 -2.5870 0.009673
## CBT_ref10                      0.038020 0.078030  0.060480  0.6286 0.529600
## TIME_ref10                    -0.048730 0.072650  0.049700 -0.9805 0.326800
## MI_ref10:CBT_ref10             0.169900 0.110300  0.105900  1.6040 0.108700
## MI_ref10:TIME_ref10            0.168700 0.103200  0.070980  2.3780 0.017430
## CBT_ref10:TIME_ref10          -0.011640 0.108900  0.064650 -0.1801 0.857100
## MI_ref10:CBT_ref10:TIME_ref10 -0.175000 0.153200  0.097850 -1.7890 0.073620
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.185 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  273    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  516
# 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")

## fit main model
PPVABS_nobl_main <- geem(formula = PPVABS.nobl.main,
                        data = PPVABS_nobl_dt,
                        id=id,
                        family = binomial
)

## print model summary
summary(PPVABS_nobl_main)
##             Estimates Model SE Robust SE    wald         p
## (Intercept)  -1.37700  0.21090   0.24260 -5.6740 1.000e-08
## cAGE         -0.35500  0.09446   0.09903 -3.5850 3.372e-04
## cBLCGSMD     -0.02391  0.01224   0.01855 -1.2890 1.975e-01
## MI1          -0.09577  0.20580   0.25370 -0.3776 7.057e-01
## CBT1          0.17260  0.20670   0.26010  0.6638 5.068e-01
## TIME1        -0.04156  0.20490   0.13470 -0.3086 7.576e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  1.01 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  295    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  590
## 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
  )

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
  )

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
  )

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
  )

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
  )

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
  )

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
  )

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
  )

## print interaction model summary
summary(PPVABS_nobl_interaction) #base: MI=0, CBT=0, TIME=0
##                Estimates Model SE Robust SE     wald         p
## (Intercept)     -1.36000  0.29890   0.30170 -4.50600 6.590e-06
## cAGE            -0.35670  0.09533   0.09861 -3.61700 2.978e-04
## cBLCGSMD        -0.02402  0.01231   0.01874 -1.28200 1.999e-01
## MI1             -0.40360  0.45830   0.46130 -0.87490 3.816e-01
## CBT1             0.17490  0.40140   0.40420  0.43260 6.653e-01
## TIME1           -0.08703  0.42040   0.26090 -0.33360 7.387e-01
## MI1:CBT1         0.44550  0.59550   0.59550  0.74820 4.544e-01
## MI1:TIME1        0.59840  0.62370   0.45090  1.32700 1.845e-01
## CBT1:TIME1       0.01350  0.57100   0.35700  0.03781 9.698e-01
## MI1:CBT1:TIME1  -0.90760  0.83390   0.56070 -1.61900 1.055e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  1.015 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  295    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  590
summary(PPVABS_nobl_interaction_TIMEref1) #base: MI=0, CBT=0, TIME=1
##                     Estimates Model SE Robust SE     wald         p
## (Intercept)          -1.44700  0.30550   0.31400 -4.60700 4.090e-06
## cAGE                 -0.35670  0.09533   0.09861 -3.61700 2.978e-04
## cBLCGSMD             -0.02402  0.01231   0.01874 -1.28200 1.999e-01
## MI1                   0.19480  0.42710   0.43130  0.45170 6.515e-01
## CBT1                  0.18840  0.40950   0.41440  0.45450 6.495e-01
## TIME_ref10            0.08703  0.42040   0.26090  0.33360 7.387e-01
## MI1:CBT1             -0.46210  0.58820   0.58590 -0.78860 4.303e-01
## MI1:TIME_ref10       -0.59840  0.62370   0.45090 -1.32700 1.845e-01
## CBT1:TIME_ref10      -0.01350  0.57100   0.35700 -0.03781 9.698e-01
## MI1:CBT1:TIME_ref10   0.90760  0.83390   0.56070  1.61900 1.055e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  1.015 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  295    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  590
summary(PPVABS_nobl_interaction_MIref1) #ref: MI=1, CBT=0, TIME=0
##                     Estimates Model SE Robust SE    wald         p
## (Intercept)          -1.76300  0.35010   0.35920 -4.9090 9.200e-07
## cAGE                 -0.35670  0.09533   0.09861 -3.6170 2.978e-04
## cBLCGSMD             -0.02402  0.01231   0.01874 -1.2820 1.999e-01
## MI_ref10              0.40360  0.45830   0.46130  0.8749 3.816e-01
## CBT1                  0.62040  0.43930   0.43990  1.4100 1.585e-01
## TIME1                 0.51130  0.46070   0.36730  1.3920 1.639e-01
## MI_ref10:CBT1        -0.44550  0.59550   0.59550 -0.7482 4.544e-01
## MI_ref10:TIME1       -0.59840  0.62370   0.45090 -1.3270 1.845e-01
## CBT1:TIME1           -0.89410  0.60780   0.43160 -2.0720 3.830e-02
## MI_ref10:CBT1:TIME1   0.90760  0.83390   0.56070  1.6190 1.055e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  1.015 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  295    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  590
summary(PPVABS_nobl_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, TIME=1
##                          Estimates Model SE Robust SE    wald         p
## (Intercept)               -1.25200  0.30030   0.30690 -4.0790 4.531e-05
## cAGE                      -0.35670  0.09533   0.09861 -3.6170 2.978e-04
## cBLCGSMD                  -0.02402  0.01231   0.01874 -1.2820 1.999e-01
## MI_ref10                  -0.19480  0.42710   0.43130 -0.4517 6.515e-01
## CBT1                      -0.27370  0.42040   0.41610 -0.6578 5.107e-01
## TIME_ref10                -0.51130  0.46070   0.36730 -1.3920 1.639e-01
## MI_ref10:CBT1              0.46210  0.58820   0.58590  0.7886 4.303e-01
## MI_ref10:TIME_ref10        0.59840  0.62370   0.45090  1.3270 1.845e-01
## CBT1:TIME_ref10            0.89410  0.60780   0.43160  2.0720 3.830e-02
## MI_ref10:CBT1:TIME_ref10  -0.90760  0.83390   0.56070 -1.6190 1.055e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  1.015 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  295    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  590
summary(PPVABS_nobl_interaction_CBTref1) #ref: MI=0, CBT=1, TIME=0
##                     Estimates Model SE Robust SE     wald         p
## (Intercept)          -1.18500  0.27130   0.26960 -4.39500 1.108e-05
## cAGE                 -0.35670  0.09533   0.09861 -3.61700 2.978e-04
## cBLCGSMD             -0.02402  0.01231   0.01874 -1.28200 1.999e-01
## MI1                   0.04195  0.37980   0.37290  0.11250 9.104e-01
## CBT_ref10            -0.17490  0.40140   0.40420 -0.43260 6.653e-01
## TIME1                -0.07353  0.38640   0.24370 -0.30170 7.629e-01
## MI1:CBT_ref10        -0.44550  0.59550   0.59550 -0.74820 4.544e-01
## MI1:TIME1            -0.30920  0.55350   0.33300 -0.92870 3.531e-01
## CBT_ref10:TIME1      -0.01350  0.57100   0.35700 -0.03781 9.698e-01
## MI1:CBT_ref10:TIME1   0.90760  0.83390   0.56070  1.61900 1.055e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  1.015 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  295    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  590
summary(PPVABS_nobl_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, TIME=1
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)               -1.25800  0.27640   0.27080 -4.64700 3.360e-06
## cAGE                      -0.35670  0.09533   0.09861 -3.61700 2.978e-04
## cBLCGSMD                  -0.02402  0.01231   0.01874 -1.28200 1.999e-01
## MI1                       -0.26730  0.40330   0.39300 -0.68010 4.964e-01
## CBT_ref10                 -0.18840  0.40950   0.41440 -0.45450 6.495e-01
## TIME_ref10                 0.07353  0.38640   0.24370  0.30170 7.629e-01
## MI1:CBT_ref10              0.46210  0.58820   0.58590  0.78860 4.303e-01
## MI1:TIME_ref10             0.30920  0.55350   0.33300  0.92870 3.531e-01
## CBT_ref10:TIME_ref10       0.01350  0.57100   0.35700  0.03781 9.698e-01
## MI1:CBT_ref10:TIME_ref10  -0.90760  0.83390   0.56070 -1.61900 1.055e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  1.015 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  295    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  590
summary(PPVABS_nobl_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, TIME=0
##                          Estimates Model SE Robust SE    wald         p
## (Intercept)               -1.14300  0.26710   0.25720 -4.4440 8.830e-06
## cAGE                      -0.35670  0.09533   0.09861 -3.6170 2.978e-04
## cBLCGSMD                  -0.02402  0.01231   0.01874 -1.2820 1.999e-01
## MI_ref10                  -0.04195  0.37980   0.37290 -0.1125 9.104e-01
## CBT_ref10                 -0.62040  0.43930   0.43990 -1.4100 1.585e-01
## TIME1                     -0.38280  0.39630   0.22670 -1.6890 9.130e-02
## MI_ref10:CBT_ref10         0.44550  0.59550   0.59550  0.7482 4.544e-01
## MI_ref10:TIME1             0.30920  0.55350   0.33300  0.9287 3.531e-01
## CBT_ref10:TIME1            0.89410  0.60780   0.43160  2.0720 3.830e-02
## MI_ref10:CBT_ref10:TIME1  -0.90760  0.83390   0.56070 -1.6190 1.055e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  1.015 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  295    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  590
summary(PPVABS_nobl_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, TIME=1
##                               Estimates Model SE Robust SE    wald         p
## (Intercept)                    -1.52600  0.29550   0.28330 -5.3860 7.000e-08
## cAGE                           -0.35670  0.09533   0.09861 -3.6170 2.978e-04
## cBLCGSMD                       -0.02402  0.01231   0.01874 -1.2820 1.999e-01
## MI_ref10                        0.26730  0.40330   0.39300  0.6801 4.964e-01
## CBT_ref10                       0.27370  0.42040   0.41610  0.6578 5.107e-01
## TIME_ref10                      0.38280  0.39630   0.22670  1.6890 9.130e-02
## MI_ref10:CBT_ref10             -0.46210  0.58820   0.58590 -0.7886 4.303e-01
## MI_ref10:TIME_ref10            -0.30920  0.55350   0.33300 -0.9287 3.531e-01
## CBT_ref10:TIME_ref10           -0.89410  0.60780   0.43160 -2.0720 3.830e-02
## MI_ref10:CBT_ref10:TIME_ref10   0.90760  0.83390   0.56070  1.6190 1.055e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  1.015 
## 
##  Number of GEE iterations: 2 
##  Number of Clusters:  295    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  590
# 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")

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


## print model summary
summary(ncigs_main)
##             Estimates Model SE Robust SE     wald       p
## (Intercept)  2.388000  0.07925   0.08408 28.40000 0.00000
## cAGE         0.114700  0.03195   0.04717  2.43200 0.01503
## MI1          0.049220  0.07124   0.09912  0.49660 0.61950
## CBT1        -0.003666  0.07137   0.10160 -0.03607 0.97120
## TIME1       -0.538000  0.08603   0.06840 -7.86500 0.00000
## TIME2       -0.466700  0.08698   0.07019 -6.64800 0.00000
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9935 
## 
##  Number of GEE iterations: 2 
##  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)
)

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)
)

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)
)

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)
  )

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)
  )

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)
  )

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)
  )

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)
  )

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)
  )


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)
  )

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)
  )

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)
  )


## print model summary
summary(ncigs_interaction) #ref: MI=0, CBT=0, time=0
##                Estimates Model SE Robust SE    wald        p
## (Intercept)       2.3230  0.12390   0.12020 19.3400 0.000000
## cAGE              0.1139  0.03201   0.04747  2.3990 0.016460
## MI1               0.1224  0.17670   0.15320  0.7989 0.424400
## CBT1              0.1758  0.16860   0.15800  1.1120 0.265900
## TIME1            -0.4595  0.17710   0.17510 -2.6230 0.008703
## TIME2            -0.4610  0.18020   0.16150 -2.8550 0.004300
## MI1:CBT1         -0.2682  0.24020   0.20320 -1.3200 0.186900
## MI1:TIME1        -0.1059  0.25370   0.21800 -0.4856 0.627300
## MI1:TIME2         0.1067  0.25750   0.20830  0.5122 0.608500
## CBT1:TIME1       -0.1986  0.24180   0.21650 -0.9173 0.359000
## CBT1:TIME2       -0.1559  0.24530   0.21830 -0.7141 0.475200
## MI1:CBT1:TIME1    0.3147  0.34480   0.27670  1.1370 0.255500
## MI1:CBT1:TIME2    0.1038  0.34880   0.27940  0.3714 0.710400
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.864000  0.12710   0.14060 13.26000 0.000000
## cAGE                 0.113900  0.03201   0.04747  2.39900 0.016460
## MI1                  0.016560  0.18270   0.21460  0.07713 0.938500
## CBT1                -0.022860  0.17380   0.19020 -0.12020 0.904300
## TIME_ref10           0.459500  0.17710   0.17510  2.62300 0.008703
## TIME_ref12          -0.001503  0.18240   0.08536 -0.01760 0.986000
## MI1:CBT1             0.046410  0.24790   0.27500  0.16880 0.866000
## MI1:TIME_ref10       0.105900  0.25370   0.21800  0.48560 0.627300
## MI1:TIME_ref12       0.212600  0.26180   0.13660  1.55600 0.119800
## CBT1:TIME_ref10      0.198600  0.24180   0.21650  0.91730 0.359000
## CBT1:TIME_ref12      0.042750  0.24880   0.11850  0.36090 0.718200
## MI1:CBT1:TIME_ref10 -0.314700  0.34480   0.27670 -1.13700 0.255500
## MI1:CBT1:TIME_ref12 -0.210900  0.35420   0.18190 -1.15900 0.246300
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.862000  0.13120   0.13120 14.2000 0.00000
## cAGE                 0.113900  0.03201   0.04747  2.3990 0.01646
## MI1                  0.229100  0.18810   0.20110  1.1390 0.25460
## CBT1                 0.019890  0.17850   0.19030  0.1045 0.91680
## TIME_ref20           0.461000  0.18020   0.16150  2.8550 0.00430
## TIME_ref21           0.001503  0.18240   0.08536  0.0176 0.98600
## MI1:CBT1            -0.164500  0.25350   0.26820 -0.6133 0.53970
## MI1:TIME_ref20      -0.106700  0.25750   0.20830 -0.5122 0.60850
## MI1:TIME_ref21      -0.212600  0.26180   0.13660 -1.5560 0.11980
## CBT1:TIME_ref20      0.155900  0.24530   0.21830  0.7141 0.47520
## CBT1:TIME_ref21     -0.042750  0.24880   0.11850 -0.3609 0.71820
## MI1:CBT1:TIME_ref20 -0.103800  0.34880   0.27940 -0.3714 0.71040
## MI1:CBT1:TIME_ref21  0.210900  0.35420   0.18190  1.1590 0.24630
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.12560   0.09504 25.7300 0.000e+00
## cAGE                  0.11390  0.03201   0.04747  2.3990 1.646e-02
## MI_ref10             -0.12240  0.17670   0.15320 -0.7989 4.244e-01
## CBT1                 -0.09249  0.17090   0.12650 -0.7310 4.648e-01
## TIME1                -0.56540  0.18160   0.12950 -4.3650 1.272e-05
## TIME2                -0.35430  0.18400   0.13190 -2.6870 7.217e-03
## MI_ref10:CBT1         0.26820  0.24020   0.20320  1.3200 1.869e-01
## MI_ref10:TIME1        0.10590  0.25370   0.21800  0.4856 6.273e-01
## MI_ref10:TIME2       -0.10670  0.25750   0.20830 -0.5122 6.085e-01
## CBT1:TIME1            0.11600  0.24580   0.17210  0.6741 5.002e-01
## CBT1:TIME2           -0.05211  0.24800   0.17440 -0.2988 7.651e-01
## MI_ref10:CBT1:TIME1  -0.31470  0.34480   0.27670 -1.1370 2.555e-01
## MI_ref10:CBT1:TIME2  -0.10380  0.34880   0.27940 -0.3714 7.104e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.88000  0.13120   0.16300 11.54000 0.000e+00
## cAGE                       0.11390  0.03201   0.04747  2.39900 1.646e-02
## MI_ref10                  -0.01656  0.18270   0.21460 -0.07713 9.385e-01
## CBT1                       0.02355  0.17670   0.19900  0.11830 9.058e-01
## TIME_ref10                 0.56540  0.18160   0.12950  4.36500 1.272e-05
## TIME_ref12                 0.21100  0.18780   0.10700  1.97300 4.850e-02
## MI_ref10:CBT1             -0.04641  0.24790   0.27500 -0.16880 8.660e-01
## MI_ref10:TIME_ref10       -0.10590  0.25370   0.21800 -0.48560 6.273e-01
## MI_ref10:TIME_ref12       -0.21260  0.26180   0.13660 -1.55600 1.198e-01
## CBT1:TIME_ref10           -0.11600  0.24580   0.17210 -0.67410 5.002e-01
## CBT1:TIME_ref12           -0.16820  0.25200   0.13820 -1.21700 2.236e-01
## MI_ref10:CBT1:TIME_ref10   0.31470  0.34480   0.27670  1.13700 2.555e-01
## MI_ref10:CBT1:TIME_ref12   0.21090  0.35420   0.18190  1.15900 2.463e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.09200  0.13450   0.15150 13.8000 0.000000
## cAGE                       0.11390  0.03201   0.04747  2.3990 0.016460
## MI_ref10                  -0.22910  0.18810   0.20110 -1.1390 0.254600
## CBT1                      -0.14460  0.17980   0.18720 -0.7724 0.439900
## TIME_ref20                 0.35430  0.18400   0.13190  2.6870 0.007217
## TIME_ref21                -0.21100  0.18780   0.10700 -1.9730 0.048500
## MI_ref10:CBT1              0.16450  0.25350   0.26820  0.6133 0.539700
## MI_ref10:TIME_ref20        0.10670  0.25750   0.20830  0.5122 0.608500
## MI_ref10:TIME_ref21        0.21260  0.26180   0.13660  1.5560 0.119800
## CBT1:TIME_ref20            0.05211  0.24800   0.17440  0.2988 0.765100
## CBT1:TIME_ref21            0.16820  0.25200   0.13820  1.2170 0.223600
## MI_ref10:CBT1:TIME_ref20   0.10380  0.34880   0.27940  0.3714 0.710400
## MI_ref10:CBT1:TIME_ref21  -0.21090  0.35420   0.18190 -1.1590 0.246300
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.4990  0.11430   0.10070 24.8100 0.000e+00
## cAGE                   0.1139  0.03201   0.04747  2.3990 1.646e-02
## MI1                   -0.1458  0.16270   0.13070 -1.1160 2.645e-01
## CBT_ref10             -0.1758  0.16860   0.15800 -1.1120 2.659e-01
## TIME1                 -0.6581  0.16460   0.12740 -5.1660 2.400e-07
## TIME2                 -0.6169  0.16640   0.14700 -4.1960 2.720e-05
## MI1:CBT_ref10          0.2682  0.24020   0.20320  1.3200 1.869e-01
## MI1:TIME1              0.2088  0.23350   0.17050  1.2250 2.206e-01
## MI1:TIME2              0.2104  0.23520   0.18620  1.1300 2.585e-01
## CBT_ref10:TIME1        0.1986  0.24180   0.21650  0.9173 3.590e-01
## CBT_ref10:TIME2        0.1559  0.24530   0.21830  0.7141 4.752e-01
## MI1:CBT_ref10:TIME1   -0.3147  0.34480   0.27670 -1.1370 2.555e-01
## MI1:CBT_ref10:TIME2   -0.1038  0.34880   0.27940 -0.3714 7.104e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.841000  0.11840   0.12660 14.54000 0.000e+00
## cAGE                      0.113900  0.03201   0.04747  2.39900 1.646e-02
## MI1                       0.062970  0.16750   0.17040  0.36940 7.118e-01
## CBT_ref10                 0.022860  0.17380   0.19020  0.12020 9.043e-01
## TIME_ref10                0.658100  0.16460   0.12740  5.16600 2.400e-07
## TIME_ref12                0.041250  0.16930   0.08224  0.50150 6.160e-01
## MI1:CBT_ref10            -0.046410  0.24790   0.27500 -0.16880 8.660e-01
## MI1:TIME_ref10           -0.208800  0.23350   0.17050 -1.22500 2.206e-01
## MI1:TIME_ref12            0.001647  0.23850   0.12000  0.01373 9.890e-01
## CBT_ref10:TIME_ref10     -0.198600  0.24180   0.21650 -0.91730 3.590e-01
## CBT_ref10:TIME_ref12     -0.042750  0.24880   0.11850 -0.36090 7.182e-01
## MI1:CBT_ref10:TIME_ref10  0.314700  0.34480   0.27670  1.13700 2.555e-01
## MI1:CBT_ref10:TIME_ref12  0.210900  0.35420   0.18190  1.15900 2.463e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.882000  0.12100   0.13630 13.81000 0.0000000
## cAGE                      0.113900  0.03201   0.04747  2.39900 0.0164600
## MI1                       0.064610  0.16990   0.17490  0.36940 0.7118000
## CBT_ref10                -0.019890  0.17850   0.19030 -0.10450 0.9168000
## TIME_ref20                0.616900  0.16640   0.14700  4.19600 0.0000272
## TIME_ref21               -0.041250  0.16930   0.08224 -0.50150 0.6160000
## MI1:CBT_ref10             0.164500  0.25350   0.26820  0.61330 0.5397000
## MI1:TIME_ref20           -0.210400  0.23520   0.18620 -1.13000 0.2585000
## MI1:TIME_ref21           -0.001647  0.23850   0.12000 -0.01373 0.9890000
## CBT_ref10:TIME_ref20     -0.155900  0.24530   0.21830 -0.71410 0.4752000
## CBT_ref10:TIME_ref21      0.042750  0.24880   0.11850  0.36090 0.7182000
## MI1:CBT_ref10:TIME_ref20  0.103800  0.34880   0.27940  0.37140 0.7104000
## MI1:CBT_ref10:TIME_ref21 -0.210900  0.35420   0.18190 -1.15900 0.2463000
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.35300  0.11590   0.08323 28.2800 0.000e+00
## cAGE                       0.11390  0.03201   0.04747  2.3990 1.646e-02
## MI_ref10                   0.14580  0.16270   0.13070  1.1160 2.645e-01
## CBT_ref10                  0.09249  0.17090   0.12650  0.7310 4.648e-01
## TIME1                     -0.44930  0.16570   0.11330 -3.9670 7.287e-05
## TIME2                     -0.40640  0.16630   0.11440 -3.5540 3.800e-04
## MI_ref10:CBT_ref10        -0.26820  0.24020   0.20320 -1.3200 1.869e-01
## MI_ref10:TIME1            -0.20880  0.23350   0.17050 -1.2250 2.206e-01
## MI_ref10:TIME2            -0.21040  0.23520   0.18620 -1.1300 2.585e-01
## CBT_ref10:TIME1           -0.11600  0.24580   0.17210 -0.6741 5.002e-01
## CBT_ref10:TIME2            0.05211  0.24800   0.17440  0.2988 7.651e-01
## MI_ref10:CBT_ref10:TIME1   0.31470  0.34480   0.27670  1.1370 2.555e-01
## MI_ref10:CBT_ref10:TIME2   0.10380  0.34880   0.27940  0.3714 7.104e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.904000  0.11840   0.11420 16.68000 0.000e+00
## cAGE                           0.113900  0.03201   0.04747  2.39900 1.646e-02
## MI_ref10                      -0.062970  0.16750   0.17040 -0.36940 7.118e-01
## CBT_ref10                     -0.023550  0.17670   0.19900 -0.11830 9.058e-01
## TIME_ref10                     0.449300  0.16570   0.11330  3.96700 7.287e-05
## TIME_ref12                     0.042900  0.16810   0.08743  0.49060 6.237e-01
## MI_ref10:CBT_ref10             0.046410  0.24790   0.27500  0.16880 8.660e-01
## MI_ref10:TIME_ref10            0.208800  0.23350   0.17050  1.22500 2.206e-01
## MI_ref10:TIME_ref12           -0.001647  0.23850   0.12000 -0.01373 9.890e-01
## CBT_ref10:TIME_ref10           0.116000  0.24580   0.17210  0.67410 5.002e-01
## CBT_ref10:TIME_ref12           0.168200  0.25200   0.13820  1.21700 2.236e-01
## MI_ref10:CBT_ref10:TIME_ref10 -0.314700  0.34480   0.27670 -1.13700 2.555e-01
## MI_ref10:CBT_ref10:TIME_ref12 -0.210900  0.35420   0.18190 -1.15900 2.463e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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.947000  0.11920   0.10970 17.74000 0.00000
## cAGE                           0.113900  0.03201   0.04747  2.39900 0.01646
## MI_ref10                      -0.064610  0.16990   0.17490 -0.36940 0.71180
## CBT_ref10                      0.144600  0.17980   0.18720  0.77240 0.43990
## TIME_ref20                     0.406400  0.16630   0.11440  3.55400 0.00038
## TIME_ref21                    -0.042900  0.16810   0.08743 -0.49060 0.62370
## MI_ref10:CBT_ref10            -0.164500  0.25350   0.26820 -0.61330 0.53970
## MI_ref10:TIME_ref20            0.210400  0.23520   0.18620  1.13000 0.25850
## MI_ref10:TIME_ref21            0.001647  0.23850   0.12000  0.01373 0.98900
## CBT_ref10:TIME_ref20          -0.052110  0.24800   0.17440 -0.29880 0.76510
## CBT_ref10:TIME_ref21          -0.168200  0.25200   0.13820 -1.21700 0.22360
## MI_ref10:CBT_ref10:TIME_ref20 -0.103800  0.34880   0.27940 -0.37140 0.71040
## MI_ref10:CBT_ref10:TIME_ref21  0.210900  0.35420   0.18190  1.15900 0.24630
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.9902 
## 
##  Number of GEE iterations: 2 
##  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")

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

## print main model summary
summary(SMDK_main)
##             Estimates Model SE Robust SE    wald         p
## (Intercept)  -0.08744  0.02814   0.02599 -3.3650 0.0007657
## cAGE          0.02691  0.01136   0.01462  1.8410 0.0655900
## MI1           0.05722  0.02532   0.03208  1.7840 0.0744200
## CBT1         -0.03167  0.02537   0.03213 -0.9857 0.3243000
## TIME1        -0.20410  0.03059   0.02651 -7.7000 0.0000000
## TIME2        -0.18840  0.03088   0.02748 -6.8550 0.0000000
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1209 
## 
##  Number of GEE iterations: 2 
##  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)
)

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)
)



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)
)

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)
  )

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)
  )

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)
  )

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)
  )

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)
  )

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)
  )


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)
  )

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)
  )

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)
  )

## print interaction model summary
summary(SMDK_interaction) #ref: MI=0, CBT=0, time=0
##                Estimates Model SE Robust SE     wald         p
## (Intercept)    -0.092470  0.04440   0.02272 -4.07000 4.696e-05
## cAGE            0.027810  0.01150   0.01494  1.86200 6.259e-02
## MI1             0.021000  0.06331   0.02753  0.76290 4.455e-01
## CBT1            0.022800  0.06040   0.02782  0.81950 4.125e-01
## TIME1          -0.208900  0.06370   0.05969 -3.49900 4.672e-04
## TIME2          -0.148400  0.06464   0.05593 -2.65300 7.978e-03
## MI1:CBT1       -0.026240  0.08605   0.03654 -0.71830 4.726e-01
## MI1:TIME1       0.028330  0.09122   0.07527  0.37630 7.067e-01
## MI1:TIME2       0.044790  0.09232   0.06750  0.66350 5.070e-01
## CBT1:TIME1     -0.018060  0.08686   0.07851 -0.23010 8.180e-01
## CBT1:TIME2     -0.195900  0.08801   0.08757 -2.23700 2.528e-02
## MI1:CBT1:TIME1  0.003268  0.12380   0.10530  0.03104 9.752e-01
## MI1:CBT1:TIME2  0.157400  0.12510   0.10710  1.47000 1.416e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.301300  0.04583   0.05928 -5.08300 3.700e-07
## cAGE                 0.027810  0.01150   0.01494  1.86200 6.259e-02
## MI1                  0.049330  0.06592   0.07684  0.64200 5.209e-01
## CBT1                 0.004735  0.06256   0.07880  0.06008 9.521e-01
## TIME_ref10           0.208900  0.06370   0.05969  3.49900 4.672e-04
## TIME_ref12           0.060480  0.06564   0.05059  1.19500 2.320e-01
## MI1:CBT1            -0.022980  0.08926   0.10610 -0.21650 8.286e-01
## MI1:TIME_ref10      -0.028330  0.09122   0.07527 -0.37630 7.067e-01
## MI1:TIME_ref12       0.016460  0.09417   0.06473  0.25430 7.992e-01
## CBT1:TIME_ref10      0.018060  0.08686   0.07851  0.23010 8.180e-01
## CBT1:TIME_ref12     -0.177800  0.08950   0.06960 -2.55500 1.062e-02
## MI1:CBT1:TIME_ref10 -0.003268  0.12380   0.10530 -0.03104 9.752e-01
## MI1:CBT1:TIME_ref12  0.154100  0.12730   0.09416  1.63700 1.017e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.24080  0.04712   0.05465 -4.4070 1.046e-05
## cAGE                  0.02781  0.01150   0.01494  1.8620 6.259e-02
## MI1                   0.06579  0.06748   0.06710  0.9804 3.269e-01
## CBT1                 -0.17310  0.06414   0.08734 -1.9820 4.749e-02
## TIME_ref20            0.14840  0.06464   0.05593  2.6530 7.978e-03
## TIME_ref21           -0.06048  0.06564   0.05059 -1.1950 2.320e-01
## MI1:CBT1              0.13110  0.09097   0.10620  1.2350 2.169e-01
## MI1:TIME_ref20       -0.04479  0.09232   0.06750 -0.6635 5.070e-01
## MI1:TIME_ref21       -0.01646  0.09417   0.06473 -0.2543 7.992e-01
## CBT1:TIME_ref20       0.19590  0.08801   0.08757  2.2370 2.528e-02
## CBT1:TIME_ref21       0.17780  0.08950   0.06960  2.5550 1.062e-02
## MI1:CBT1:TIME_ref20  -0.15740  0.12510   0.10710 -1.4700 1.416e-01
## MI1:CBT1:TIME_ref21  -0.15410  0.12730   0.09416 -1.6370 1.017e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.071460  0.04502   0.01403 -5.09500 3.500e-07
## cAGE                 0.027810  0.01150   0.01494  1.86200 6.259e-02
## MI_ref10            -0.021000  0.06331   0.02753 -0.76290 4.455e-01
## CBT1                -0.003446  0.06121   0.02253 -0.15290 8.785e-01
## TIME1               -0.180500  0.06530   0.04552 -3.96600 7.312e-05
## TIME2               -0.103600  0.06591   0.03789 -2.73400 6.262e-03
## MI_ref10:CBT1        0.026240  0.08605   0.03654  0.71830 4.726e-01
## MI_ref10:TIME1      -0.028330  0.09122   0.07527 -0.37630 7.067e-01
## MI_ref10:TIME2      -0.044790  0.09232   0.06750 -0.66350 5.070e-01
## CBT1:TIME1          -0.014790  0.08828   0.07007 -0.21110 8.328e-01
## CBT1:TIME2          -0.038500  0.08885   0.06165 -0.62460 5.323e-01
## MI_ref10:CBT1:TIME1 -0.003268  0.12380   0.10530 -0.03104 9.752e-01
## MI_ref10:CBT1:TIME2 -0.157400  0.12510   0.10710 -1.47000 1.416e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.252000  0.04733   0.04879 -5.16500 2.400e-07
## cAGE                      0.027810  0.01150   0.01494  1.86200 6.259e-02
## MI_ref10                 -0.049330  0.06592   0.07684 -0.64200 5.209e-01
## CBT1                     -0.018240  0.06363   0.07089 -0.25730 7.969e-01
## TIME_ref10                0.180500  0.06530   0.04552  3.96600 7.312e-05
## TIME_ref12                0.076940  0.06754   0.04027  1.91100 5.606e-02
## MI_ref10:CBT1             0.022980  0.08926   0.10610  0.21650 8.286e-01
## MI_ref10:TIME_ref10       0.028330  0.09122   0.07527  0.37630 7.067e-01
## MI_ref10:TIME_ref12      -0.016460  0.09417   0.06473 -0.25430 7.992e-01
## CBT1:TIME_ref10           0.014790  0.08828   0.07007  0.21110 8.328e-01
## CBT1:TIME_ref12          -0.023710  0.09055   0.06329 -0.37460 7.080e-01
## MI_ref10:CBT1:TIME_ref10  0.003268  0.12380   0.10530  0.03104 9.752e-01
## MI_ref10:CBT1:TIME_ref12 -0.154100  0.12730   0.09416 -1.63700 1.017e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.17510  0.04821   0.03866 -4.5290 5.940e-06
## cAGE                       0.02781  0.01150   0.01494  1.8620 6.259e-02
## MI_ref10                  -0.06579  0.06748   0.06710 -0.9804 3.269e-01
## CBT1                      -0.04195  0.06444   0.05990 -0.7003 4.837e-01
## TIME_ref20                 0.10360  0.06591   0.03789  2.7340 6.262e-03
## TIME_ref21                -0.07694  0.06754   0.04027 -1.9110 5.606e-02
## MI_ref10:CBT1             -0.13110  0.09097   0.10620 -1.2350 2.169e-01
## MI_ref10:TIME_ref20        0.04479  0.09232   0.06750  0.6635 5.070e-01
## MI_ref10:TIME_ref21        0.01646  0.09417   0.06473  0.2543 7.992e-01
## CBT1:TIME_ref20            0.03850  0.08885   0.06165  0.6246 5.323e-01
## CBT1:TIME_ref21            0.02371  0.09055   0.06329  0.3746 7.080e-01
## MI_ref10:CBT1:TIME_ref20   0.15740  0.12510   0.10710  1.4700 1.416e-01
## MI_ref10:CBT1:TIME_ref21   0.15410  0.12730   0.09416  1.6370 1.017e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.069670  0.04094   0.01557 -4.47500 7.630e-06
## cAGE                 0.027810  0.01150   0.01494  1.86200 6.259e-02
## MI1                 -0.005242  0.05827   0.02355 -0.22250 8.239e-01
## CBT_ref10           -0.022800  0.06040   0.02782 -0.81950 4.125e-01
## TIME1               -0.226900  0.05905   0.05119 -4.43300 9.300e-06
## TIME2               -0.344300  0.05973   0.06756 -5.09600 3.500e-07
## MI1:CBT_ref10        0.026240  0.08605   0.03654  0.71830 4.726e-01
## MI1:TIME1            0.031590  0.08376   0.07383  0.42790 6.687e-01
## MI1:TIME2            0.202200  0.08436   0.08318  2.43100 1.507e-02
## CBT_ref10:TIME1      0.018060  0.08686   0.07851  0.23010 8.180e-01
## CBT_ref10:TIME2      0.195900  0.08801   0.08757  2.23700 2.528e-02
## MI1:CBT_ref10:TIME1 -0.003268  0.12380   0.10530 -0.03104 9.752e-01
## MI1:CBT_ref10:TIME2 -0.157400  0.12510   0.10710 -1.47000 1.416e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.296600  0.04255   0.05138 -5.77300 1.000e-08
## cAGE                      0.027810  0.01150   0.01494  1.86200 6.259e-02
## MI1                       0.026350  0.06017   0.07267  0.36260 7.169e-01
## CBT_ref10                -0.004735  0.06256   0.07880 -0.06008 9.521e-01
## TIME_ref10                0.226900  0.05905   0.05119  4.43300 9.300e-06
## TIME_ref12               -0.117400  0.06084   0.04781 -2.45500 1.410e-02
## MI1:CBT_ref10             0.022980  0.08926   0.10610  0.21650 8.286e-01
## MI1:TIME_ref10           -0.031590  0.08376   0.07383 -0.42790 6.687e-01
## MI1:TIME_ref12            0.170600  0.08568   0.06837  2.49500 1.259e-02
## CBT_ref10:TIME_ref10     -0.018060  0.08686   0.07851 -0.23010 8.180e-01
## CBT_ref10:TIME_ref12      0.177800  0.08950   0.06960  2.55500 1.062e-02
## MI1:CBT_ref10:TIME_ref10  0.003268  0.12380   0.10530  0.03104 9.752e-01
## MI1:CBT_ref10:TIME_ref12 -0.154100  0.12730   0.09416 -1.63700 1.017e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.41390  0.04349   0.06757 -6.126 0.000e+00
## cAGE                       0.02781  0.01150   0.01494  1.862 6.259e-02
## MI1                        0.19690  0.06100   0.08161  2.413 1.582e-02
## CBT_ref10                  0.17310  0.06414   0.08734  1.982 4.749e-02
## TIME_ref20                 0.34430  0.05973   0.06756  5.096 3.500e-07
## TIME_ref21                 0.11740  0.06084   0.04781  2.455 1.410e-02
## MI1:CBT_ref10             -0.13110  0.09097   0.10620 -1.235 2.169e-01
## MI1:TIME_ref20            -0.20220  0.08436   0.08318 -2.431 1.507e-02
## MI1:TIME_ref21            -0.17060  0.08568   0.06837 -2.495 1.259e-02
## CBT_ref10:TIME_ref20      -0.19590  0.08801   0.08757 -2.237 2.528e-02
## CBT_ref10:TIME_ref21      -0.17780  0.08950   0.06960 -2.555 1.062e-02
## MI1:CBT_ref10:TIME_ref20   0.15740  0.12510   0.10710  1.470 1.416e-01
## MI1:CBT_ref10:TIME_ref21   0.15410  0.12730   0.09416  1.637 1.017e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.074910  0.04146   0.01767 -4.24000 2.239e-05
## cAGE                      0.027810  0.01150   0.01494  1.86200 6.259e-02
## MI_ref10                  0.005242  0.05827   0.02355  0.22250 8.239e-01
## CBT_ref10                 0.003446  0.06121   0.02253  0.15290 8.785e-01
## TIME1                    -0.195300  0.05940   0.05326 -3.66700 2.453e-04
## TIME2                    -0.142100  0.05958   0.04868 -2.91900 3.513e-03
## MI_ref10:CBT_ref10       -0.026240  0.08605   0.03654 -0.71830 4.726e-01
## MI_ref10:TIME1           -0.031590  0.08376   0.07383 -0.42790 6.687e-01
## MI_ref10:TIME2           -0.202200  0.08436   0.08318 -2.43100 1.507e-02
## CBT_ref10:TIME1           0.014790  0.08828   0.07007  0.21110 8.328e-01
## CBT_ref10:TIME2           0.038500  0.08885   0.06165  0.62460 5.323e-01
## MI_ref10:CBT_ref10:TIME1  0.003268  0.12380   0.10530  0.03104 9.752e-01
## MI_ref10:CBT_ref10:TIME2  0.157400  0.12510   0.10710  1.47000 1.416e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.270200  0.04254   0.05147 -5.25000 1.500e-07
## cAGE                           0.027810  0.01150   0.01494  1.86200 6.259e-02
## MI_ref10                      -0.026350  0.06017   0.07267 -0.36260 7.169e-01
## CBT_ref10                      0.018240  0.06363   0.07089  0.25730 7.969e-01
## TIME_ref10                     0.195300  0.05940   0.05326  3.66700 2.453e-04
## TIME_ref12                     0.053230  0.06033   0.04888  1.08900 2.762e-01
## MI_ref10:CBT_ref10            -0.022980  0.08926   0.10610 -0.21650 8.286e-01
## MI_ref10:TIME_ref10            0.031590  0.08376   0.07383  0.42790 6.687e-01
## MI_ref10:TIME_ref12           -0.170600  0.08568   0.06837 -2.49500 1.259e-02
## CBT_ref10:TIME_ref10          -0.014790  0.08828   0.07007 -0.21110 8.328e-01
## CBT_ref10:TIME_ref12           0.023710  0.09055   0.06329  0.37460 7.080e-01
## MI_ref10:CBT_ref10:TIME_ref10 -0.003268  0.12380   0.10530 -0.03104 9.752e-01
## MI_ref10:CBT_ref10:TIME_ref12  0.154100  0.12730   0.09416  1.63700 1.017e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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.21700  0.04278   0.04596 -4.7220 2.340e-06
## cAGE                            0.02781  0.01150   0.01494  1.8620 6.259e-02
## MI_ref10                       -0.19690  0.06100   0.08161 -2.4130 1.582e-02
## CBT_ref10                       0.04195  0.06444   0.05990  0.7003 4.837e-01
## TIME_ref20                      0.14210  0.05958   0.04868  2.9190 3.513e-03
## TIME_ref21                     -0.05323  0.06033   0.04888 -1.0890 2.762e-01
## MI_ref10:CBT_ref10              0.13110  0.09097   0.10620  1.2350 2.169e-01
## MI_ref10:TIME_ref20             0.20220  0.08436   0.08318  2.4310 1.507e-02
## MI_ref10:TIME_ref21             0.17060  0.08568   0.06837  2.4950 1.259e-02
## CBT_ref10:TIME_ref20           -0.03850  0.08885   0.06165 -0.6246 5.323e-01
## CBT_ref10:TIME_ref21           -0.02371  0.09055   0.06329 -0.3746 7.080e-01
## MI_ref10:CBT_ref10:TIME_ref20  -0.15740  0.12510   0.10710 -1.4700 1.416e-01
## MI_ref10:CBT_ref10:TIME_ref21  -0.15410  0.12730   0.09416 -1.6370 1.017e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  0.1229 
## 
##  Number of GEE iterations: 2 
##  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)
)
  

## print main model summary
summary(CLONGABS_nobl_main) #unlogged
##             Estimates Model SE Robust SE     wald        p
## (Intercept)  1.710000  0.20750   0.22290  7.67300 0.000000
## cAGE        -0.130200  0.09316   0.11490 -1.13300 0.257200
## cBLCGSMD    -0.001452  0.01168   0.01595 -0.09105 0.927500
## MI1         -0.748200  0.20370   0.24840 -3.01200 0.002599
## CBT1         0.692700  0.20510   0.24750  2.79800 0.005135
## TIME1        0.008417  0.20280   0.13420  0.06270 0.950000
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  4.816 
## 
##  Number of GEE iterations: 3 
##  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)
)



## print interaction models' summary
summary(CLONGABS_nobl_interaction) #unlogged
##                Estimates Model SE Robust SE    wald      p
## (Intercept)     -1.94000  0.26250    0.2233 -8.6880 0.0000
## cAGE            -0.12350  0.08254    0.1451 -0.8514 0.3945
## cBLCGSMD        -0.03319  0.01003    0.0395 -0.8403 0.4007
## MI1             -0.54300  0.39020    0.3564 -1.5240 0.1276
## CBT1             0.25260  0.35260    0.3218  0.7851 0.4324
## TIME1           -0.03922  0.37340    0.2055 -0.1908 0.8486
## MI1:CBT1         0.41830  0.51530    0.4710  0.8881 0.3745
## MI1:TIME1       -0.58830  0.57090    0.4344 -1.3540 0.1756
## CBT1:TIME1       0.29280  0.49970    0.2599  1.1270 0.2598
## MI1:CBT1:TIME1   0.07867  0.74520    0.5274  0.1492 0.8814
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  3.479 
## 
##  Number of GEE iterations: 9 
##  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")

## main model
LONGABS_nobl_main <- geem(formula = CLONGABS.nobl.main,
                        data = LONGABS_nobl_dt,
                        id=id,
                        family = MASS::negative.binomial(1)
)


## print main model summary
summary(LONGABS_nobl_main)
##             Estimates Model SE Robust SE    wald       p
## (Intercept) -1.983000  0.17630   0.21210 -9.3500 0.00000
## cAGE        -0.188300  0.07857   0.10100 -1.8640 0.06238
## cBLCGSMD    -0.009257  0.01004   0.01623 -0.5705 0.56830
## MI1         -0.660400  0.17400   0.21870 -3.0200 0.00253
## CBT1         0.638000  0.17580   0.22000  2.9000 0.00373
## TIME1       -0.115800  0.17330   0.12700 -0.9118 0.36190
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  3.148 
## 
##  Number of GEE iterations: 9 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
## 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)
)

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)
)


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)
  )

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)
  )

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)
  )

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)
  )

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)
  )

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)
  )


## print interaction model summary
summary(LONGABS_nobl_interaction) #base: MI=0, CBT=0, Time=0
##                 Estimates Model SE Robust SE     wald       p
## (Intercept)     1.6520000  0.29880   0.25290  6.53200 0.00000
## cAGE           -0.1498000  0.09382   0.11560 -1.29600 0.19500
## cBLCGSMD        0.0006783  0.01175   0.01518  0.04469 0.96440
## MI1            -0.5704000  0.43700   0.37950 -1.50300 0.13280
## CBT1            0.5529000  0.40100   0.31830  1.73700 0.08243
## TIME1           0.2183000  0.41700   0.26400  0.82670 0.40840
## MI1:CBT1        0.1861000  0.57940   0.52240  0.35630 0.72160
## MI1:TIME1      -0.6656000  0.62790   0.48010 -1.38600 0.16560
## CBT1:TIME1      0.0476600  0.56180   0.31130  0.15310 0.87830
## MI1:CBT1:TIME1  0.2107000  0.82790   0.58470  0.36040 0.71860
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  4.841 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
summary(LONGABS_nobl_interaction_TIMEref1) #base: MI=0, CBT=0, Time=1
##                      Estimates Model SE Robust SE     wald       p
## (Intercept)          1.8700000  0.29400   0.31000  6.03200 0.00000
## cAGE                -0.1498000  0.09382   0.11560 -1.29600 0.19500
## cBLCGSMD             0.0006783  0.01175   0.01518  0.04469 0.96440
## MI1                 -1.2360000  0.45380   0.53310 -2.31900 0.02042
## CBT1                 0.6005000  0.39630   0.35660  1.68400 0.09221
## TIME_ref10          -0.2183000  0.41700   0.26400 -0.82670 0.40840
## MI1:CBT1             0.3968000  0.59460   0.66450  0.59720 0.55040
## MI1:TIME_ref10       0.6656000  0.62790   0.48010  1.38600 0.16560
## CBT1:TIME_ref10     -0.0476600  0.56180   0.31130 -0.15310 0.87830
## MI1:CBT1:TIME_ref10 -0.2107000  0.82790   0.58470 -0.36040 0.71860
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  4.841 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
summary(LONGABS_nobl_interaction_MIref1) #ref: MI=1, CBT=0, Time=0
##                      Estimates Model SE Robust SE     wald         p
## (Intercept)          1.0820000  0.31850   0.27930  3.87300 0.0001076
## cAGE                -0.1498000  0.09382   0.11560 -1.29600 0.1950000
## cBLCGSMD             0.0006783  0.01175   0.01518  0.04469 0.9644000
## MI_ref10             0.5704000  0.43700   0.37950  1.50300 0.1328000
## CBT1                 0.7390000  0.41750   0.41090  1.79900 0.0720700
## TIME1               -0.4474000  0.46940   0.40120 -1.11500 0.2649000
## MI_ref10:CBT1       -0.1861000  0.57940   0.52240 -0.35630 0.7216000
## MI_ref10:TIME1       0.6656000  0.62790   0.48010  1.38600 0.1656000
## CBT1:TIME1           0.2584000  0.60820   0.49600  0.52080 0.6025000
## MI_ref10:CBT1:TIME1 -0.2107000  0.82790   0.58470 -0.36040 0.7186000
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  4.841 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
summary(LONGABS_nobl_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, Time=1
##                           Estimates Model SE Robust SE     wald       p
## (Intercept)               0.6342000  0.34550   0.41690  1.52100 0.12820
## cAGE                     -0.1498000  0.09382   0.11560 -1.29600 0.19500
## cBLCGSMD                  0.0006783  0.01175   0.01518  0.04469 0.96440
## MI_ref10                  1.2360000  0.45380   0.53310  2.31900 0.02042
## CBT1                      0.9974000  0.44250   0.55080  1.81100 0.07019
## TIME_ref10                0.4474000  0.46940   0.40120  1.11500 0.26490
## MI_ref10:CBT1            -0.3968000  0.59460   0.66450 -0.59720 0.55040
## MI_ref10:TIME_ref10      -0.6656000  0.62790   0.48010 -1.38600 0.16560
## CBT1:TIME_ref10          -0.2584000  0.60820   0.49600 -0.52080 0.60250
## MI_ref10:CBT1:TIME_ref10  0.2107000  0.82790   0.58470  0.36040 0.71860
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  4.841 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
summary(LONGABS_nobl_interaction_CBTref1) #ref: MI=0, CBT=1, Time=0
##                      Estimates Model SE Robust SE     wald       p
## (Intercept)          2.2050000  0.26690   0.19380 11.38000 0.00000
## cAGE                -0.1498000  0.09382   0.11560 -1.29600 0.19500
## cBLCGSMD             0.0006783  0.01175   0.01518  0.04469 0.96440
## MI1                 -0.3843000  0.37990   0.35170 -1.09300 0.27460
## CBT_ref10           -0.5529000  0.40100   0.31830 -1.73700 0.08243
## TIME1                0.2659000  0.37640   0.16570  1.60500 0.10860
## MI1:CBT_ref10       -0.1861000  0.57940   0.52240 -0.35630 0.72160
## MI1:TIME1           -0.4549000  0.53960   0.33330 -1.36500 0.17230
## CBT_ref10:TIME1     -0.0476600  0.56180   0.31130 -0.15310 0.87830
## MI1:CBT_ref10:TIME1 -0.2107000  0.82790   0.58470 -0.36040 0.71860
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  4.841 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
summary(LONGABS_nobl_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, Time=1
##                           Estimates Model SE Robust SE     wald       p
## (Intercept)               2.4710000  0.26540   0.17750 13.92000 0.00000
## cAGE                     -0.1498000  0.09382   0.11560 -1.29600 0.19500
## cBLCGSMD                  0.0006783  0.01175   0.01518  0.04469 0.96440
## MI1                      -0.8392000  0.38380   0.38880 -2.15900 0.03088
## CBT_ref10                -0.6005000  0.39630   0.35660 -1.68400 0.09221
## TIME_ref10               -0.2659000  0.37640   0.16570 -1.60500 0.10860
## MI1:CBT_ref10            -0.3968000  0.59460   0.66450 -0.59720 0.55040
## MI1:TIME_ref10            0.4549000  0.53960   0.33330  1.36500 0.17230
## CBT_ref10:TIME_ref10      0.0476600  0.56180   0.31130  0.15310 0.87830
## MI1:CBT_ref10:TIME_ref10  0.2107000  0.82790   0.58470  0.36040 0.71860
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  4.841 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
summary(LONGABS_nobl_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, Time=0
##                           Estimates Model SE Robust SE     wald       p
## (Intercept)               1.8210000  0.27020   0.29380  6.19600 0.00000
## cAGE                     -0.1498000  0.09382   0.11560 -1.29600 0.19500
## cBLCGSMD                  0.0006783  0.01175   0.01518  0.04469 0.96440
## MI_ref10                  0.3843000  0.37990   0.35170  1.09300 0.27460
## CBT_ref10                -0.7390000  0.41750   0.41090 -1.79900 0.07207
## TIME1                    -0.1890000  0.38670   0.28970 -0.65240 0.51410
## MI_ref10:CBT_ref10        0.1861000  0.57940   0.52240  0.35630 0.72160
## MI_ref10:TIME1            0.4549000  0.53960   0.33330  1.36500 0.17230
## CBT_ref10:TIME1          -0.2584000  0.60820   0.49600 -0.52080 0.60250
## MI_ref10:CBT_ref10:TIME1  0.2107000  0.82790   0.58470  0.36040 0.71860
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  4.841 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
summary(LONGABS_nobl_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, Time=1
##                                Estimates Model SE Robust SE     wald         p
## (Intercept)                    1.6320000  0.27690   0.34850  4.68200 2.840e-06
## cAGE                          -0.1498000  0.09382   0.11560 -1.29600 1.950e-01
## cBLCGSMD                       0.0006783  0.01175   0.01518  0.04469 9.644e-01
## MI_ref10                       0.8392000  0.38380   0.38880  2.15900 3.088e-02
## CBT_ref10                     -0.9974000  0.44250   0.55080 -1.81100 7.019e-02
## TIME_ref10                     0.1890000  0.38670   0.28970  0.65240 5.141e-01
## MI_ref10:CBT_ref10             0.3968000  0.59460   0.66450  0.59720 5.504e-01
## MI_ref10:TIME_ref10           -0.4549000  0.53960   0.33330 -1.36500 1.723e-01
## CBT_ref10:TIME_ref10           0.2584000  0.60820   0.49600  0.52080 6.025e-01
## MI_ref10:CBT_ref10:TIME_ref10 -0.2107000  0.82790   0.58470 -0.36040 7.186e-01
## 
##  Estimated Correlation Parameter:  0 
##  Correlation Structure:  independence 
##  Est. Scale Parameter:  4.841 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  282    Maximum Cluster Size:  2 
##  Number of observations with nonzero weight:  560
# Save ---------------------------

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