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


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

corstr = "ar1"


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

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

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


## print main model summary
summary(AVCIG_nobl_main)
##             Estimates Model SE Robust SE    wald       p
## (Intercept)   1.84800 0.105300  0.111600 16.5700 0.00000
## cAGE         -0.02492 0.052540  0.058530 -0.4258 0.67030
## cBLCGSMD      0.04747 0.006111  0.005162  9.1960 0.00000
## MI1           0.05612 0.116100  0.113600  0.4940 0.62130
## CBT1         -0.17070 0.109500  0.118500 -1.4400 0.14980
## TIME1         0.09425 0.045210  0.048540  1.9420 0.05218
## 
##  Estimated Correlation Parameter:  0.7353 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8569 
## 
##  Number of GEE iterations: 6 
##  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),
                        corstr = corstr
)

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


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

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

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

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

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

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


## print interaction model summary
summary(AVCIG_nobl_interaction) #base: MI=0, CBT=0, Time=0
##                Estimates Model SE Robust SE     wald       p
## (Intercept)     1.922000 0.127000  0.143800 13.37000 0.00000
## cAGE           -0.023870 0.052770  0.058810 -0.40590 0.68480
## cBLCGSMD        0.047940 0.006148  0.005109  9.38300 0.00000
## MI1            -0.115900 0.182000  0.184500 -0.62850 0.52970
## CBT1           -0.267300 0.164200  0.188500 -1.41800 0.15610
## TIME1           0.002878 0.091190  0.094910  0.03032 0.97580
## MI1:CBT1        0.224000 0.234800  0.255400  0.87680 0.38060
## MI1:TIME1       0.243200 0.130000  0.138300  1.75900 0.07861
## CBT1:TIME1      0.089640 0.125600  0.124300  0.72100 0.47090
## MI1:CBT1:TIME1 -0.269700 0.178300  0.190000 -1.42000 0.15570
## 
##  Estimated Correlation Parameter:  0.7488 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8521 
## 
##  Number of GEE iterations: 6 
##  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.925000 0.131400  0.130100 14.80000 0.00000
## cAGE                -0.023870 0.052770  0.058810 -0.40590 0.68480
## cBLCGSMD             0.047940 0.006148  0.005109  9.38300 0.00000
## MI1                  0.127300 0.187700  0.188600  0.67490 0.49980
## CBT1                -0.177700 0.168500  0.183600 -0.96750 0.33330
## TIME_ref10          -0.002878 0.091190  0.094910 -0.03032 0.97580
## MI1:CBT1            -0.045700 0.240300  0.261200 -0.17490 0.86110
## MI1:TIME_ref10      -0.243200 0.130000  0.138300 -1.75900 0.07861
## CBT1:TIME_ref10     -0.089640 0.125600  0.124300 -0.72100 0.47090
## MI1:CBT1:TIME_ref10  0.269700 0.178300  0.190000  1.42000 0.15570
## 
##  Estimated Correlation Parameter:  0.7488 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8521 
## 
##  Number of GEE iterations: 7 
##  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.80600 0.129500  0.114500 15.7700 0.00000
## cAGE                 -0.02387 0.052770  0.058810 -0.4059 0.68480
## cBLCGSMD              0.04794 0.006148  0.005109  9.3830 0.00000
## MI_ref10              0.11590 0.182000  0.184500  0.6285 0.52970
## CBT1                 -0.04336 0.166900  0.172600 -0.2512 0.80160
## TIME1                 0.24610 0.092660  0.100200  2.4560 0.01404
## MI_ref10:CBT1        -0.22400 0.234800  0.255400 -0.8768 0.38060
## MI_ref10:TIME1       -0.24320 0.130000  0.138300 -1.7590 0.07861
## CBT1:TIME1           -0.18000 0.126500  0.143200 -1.2570 0.20880
## MI_ref10:CBT1:TIME1   0.26970 0.178300  0.190000  1.4200 0.15570
## 
##  Estimated Correlation Parameter:  0.7488 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8521 
## 
##  Number of GEE iterations: 7 
##  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)                2.05200 0.133000  0.133500 15.3700 0.00000
## cAGE                      -0.02387 0.052770  0.058810 -0.4059 0.68480
## cBLCGSMD                   0.04794 0.006148  0.005109  9.3830 0.00000
## MI_ref10                  -0.12730 0.187700  0.188600 -0.6749 0.49980
## CBT1                      -0.22340 0.170100  0.184500 -1.2110 0.22590
## TIME_ref10                -0.24610 0.092660  0.100200 -2.4560 0.01404
## MI_ref10:CBT1              0.04570 0.240300  0.261200  0.1749 0.86110
## MI_ref10:TIME_ref10        0.24320 0.130000  0.138300  1.7590 0.07861
## CBT1:TIME_ref10            0.18000 0.126500  0.143200  1.2570 0.20880
## MI_ref10:CBT1:TIME_ref10  -0.26970 0.178300  0.190000 -1.4200 0.15570
## 
##  Estimated Correlation Parameter:  0.7488 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.8521 
## 
##  Number of GEE iterations: 7 
##  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.643000 0.118200  0.130600 12.59000 0.00000
## cAGE                -0.024580 0.052480  0.058960 -0.41680 0.67680
## cBLCGSMD             0.049330 0.006084  0.005194  9.49800 0.00000
## MI1                  0.222400 0.166300  0.176100  1.26300 0.20650
## CBT_ref10            0.293000 0.161400  0.230800  1.27000 0.20410
## TIME1                0.051300 0.085740  0.091560  0.56030 0.57530
## MI1:CBT_ref10       -0.488700 0.230500  0.287100 -1.70200 0.08869
## MI1:TIME1            0.020050 0.120600  0.131900  0.15200 0.87920
## CBT_ref10:TIME1      0.008267 0.127700  0.110800  0.07463 0.94050
## MI1:CBT_ref10:TIME1  0.152900 0.181600  0.184900  0.82700 0.40820
## 
##  Estimated Correlation Parameter:  0.7399 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.846 
## 
##  Number of GEE iterations: 5 
##  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.695000 0.120100  0.155000 10.93000 0.0000
## cAGE                     -0.024580 0.052480  0.058960 -0.41680 0.6768
## cBLCGSMD                  0.049330 0.006084  0.005194  9.49800 0.0000
## MI1                       0.242500 0.168500  0.196300  1.23500 0.2167
## CBT_ref10                 0.301300 0.165400  0.226300  1.33200 0.1830
## TIME_ref10               -0.051300 0.085740  0.091560 -0.56030 0.5753
## MI1:CBT_ref10            -0.335800 0.235600  0.293000 -1.14600 0.2517
## MI1:TIME_ref10           -0.020050 0.120600  0.131900 -0.15200 0.8792
## CBT_ref10:TIME_ref10     -0.008267 0.127700  0.110800 -0.07463 0.9405
## MI1:CBT_ref10:TIME_ref10 -0.152900 0.181600  0.184900 -0.82700 0.4082
## 
##  Estimated Correlation Parameter:  0.7399 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.846 
## 
##  Number of GEE iterations: 5 
##  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.86600 0.117100  0.118600 15.7300 0.00000
## cAGE                      -0.02458 0.052480  0.058960 -0.4168 0.67680
## cBLCGSMD                   0.04933 0.006084  0.005194  9.4980 0.00000
## MI_ref10                  -0.22240 0.166300  0.176100 -1.2630 0.20650
## CBT_ref10                 -0.19570 0.164000  0.170000 -1.1510 0.24980
## TIME1                      0.07135 0.084860  0.094900  0.7518 0.45220
## MI_ref10:CBT_ref10         0.48870 0.230500  0.287100  1.7020 0.08869
## MI_ref10:TIME1            -0.02005 0.120600  0.131900 -0.1520 0.87920
## CBT_ref10:TIME1            0.16120 0.129000  0.147800  1.0900 0.27550
## MI_ref10:CBT_ref10:TIME1  -0.15290 0.181600  0.184900 -0.8270 0.40820
## 
##  Estimated Correlation Parameter:  0.7399 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.846 
## 
##  Number of GEE iterations: 5 
##  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.93700 0.118200  0.121500 15.9500 0.0000
## cAGE                           -0.02458 0.052480  0.058960 -0.4168 0.6768
## cBLCGSMD                        0.04933 0.006084  0.005194  9.4980 0.0000
## MI_ref10                       -0.24250 0.168500  0.196300 -1.2350 0.2167
## CBT_ref10                      -0.03453 0.167100  0.183800 -0.1878 0.8510
## TIME_ref10                     -0.07135 0.084860  0.094900 -0.7518 0.4522
## MI_ref10:CBT_ref10              0.33580 0.235600  0.293000  1.1460 0.2517
## MI_ref10:TIME_ref10             0.02005 0.120600  0.131900  0.1520 0.8792
## CBT_ref10:TIME_ref10           -0.16120 0.129000  0.147800 -1.0900 0.2755
## MI_ref10:CBT_ref10:TIME_ref10   0.15290 0.181600  0.184900  0.8270 0.4082
## 
##  Estimated Correlation Parameter:  0.7399 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.846 
## 
##  Number of GEE iterations: 5 
##  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),
                        corstr = corstr
)


## print main model summary
summary(CGSMD_nobl_main)
##             Estimates Model SE Robust SE     wald       p
## (Intercept)  2.085000 0.088400  0.094870 21.98000 0.00000
## cAGE        -0.002086 0.044150  0.048260 -0.04322 0.96550
## cBLCGSMD     0.041330 0.005246  0.004316  9.57400 0.00000
## MI1          0.002836 0.097180  0.093440  0.03035 0.97580
## CBT1        -0.221100 0.083500  0.105900 -2.08800 0.03681
## TIME1        0.101000 0.036460  0.041910  2.40900 0.01599
## 
##  Estimated Correlation Parameter:  0.7525 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5691 
## 
##  Number of GEE iterations: 8 
##  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),
                               corstr = corstr
)

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


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

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

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

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

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

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


## print interaction model summary
summary(CGSMD_nobl_interaction) #base: MI=0, CBT=0, time=0
##                 Estimates Model SE Robust SE     wald       p
## (Intercept)     2.1950000 0.104600   0.11920 18.42000 0.00000
## cAGE            0.0005014 0.044090   0.04842  0.01036 0.99170
## cBLCGSMD        0.0424200 0.005259   0.00428  9.91000 0.00000
## MI1            -0.2359000 0.149500   0.16350 -1.44300 0.14910
## CBT1           -0.4099000 0.122000   0.16020 -2.55900 0.01051
## TIME1           0.0156900 0.073150   0.07804  0.20100 0.84070
## MI1:CBT1        0.4026000 0.179700   0.22660  1.77600 0.07567
## MI1:TIME1       0.1964000 0.103600   0.12340  1.59200 0.11140
## CBT1:TIME1      0.1435000 0.104800   0.10140  1.41600 0.15690
## MI1:CBT1:TIME1 -0.3240000 0.146100   0.16110 -2.01100 0.04427
## 
##  Estimated Correlation Parameter:  0.7632 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5636 
## 
##  Number of GEE iterations: 8 
##  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.2110000 0.108300   0.11770 18.77000 0.00000
## cAGE                 0.0005014 0.044090   0.04842  0.01036 0.99170
## cBLCGSMD             0.0424200 0.005259   0.00428  9.91000 0.00000
## MI1                 -0.0394700 0.152800   0.17190 -0.22960 0.81840
## CBT1                -0.2664000 0.129000   0.16050 -1.66000 0.09695
## TIME_ref10          -0.0156900 0.073150   0.07804 -0.20100 0.84070
## MI1:CBT1             0.0785200 0.184600   0.23410  0.33540 0.73730
## MI1:TIME_ref10      -0.1964000 0.103600   0.12340 -1.59200 0.11140
## CBT1:TIME_ref10     -0.1435000 0.104800   0.10140 -1.41600 0.15690
## MI1:CBT1:TIME_ref10  0.3240000 0.146100   0.16110  2.01100 0.04427
## 
##  Estimated Correlation Parameter:  0.7632 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5636 
## 
##  Number of GEE iterations: 8 
##  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.9590000 0.106600   0.11190 17.51000 0.00000
## cAGE                 0.0005014 0.044090   0.04842  0.01036 0.99170
## cBLCGSMD             0.0424200 0.005259   0.00428  9.91000 0.00000
## MI_ref10             0.2359000 0.149500   0.16350  1.44300 0.14910
## CBT1                -0.0073280 0.131600   0.16100 -0.04551 0.96370
## TIME1                0.2121000 0.073440   0.09512  2.23000 0.02575
## MI_ref10:CBT1       -0.4026000 0.179700   0.22660 -1.77600 0.07567
## MI_ref10:TIME1      -0.1964000 0.103600   0.12340 -1.59200 0.11140
## CBT1:TIME1          -0.1805000 0.101900   0.12470 -1.44800 0.14770
## MI_ref10:CBT1:TIME1  0.3240000 0.146100   0.16110  2.01100 0.04427
## 
##  Estimated Correlation Parameter:  0.7632 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5636 
## 
##  Number of GEE iterations: 9 
##  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.1710000 0.107200   0.12370 17.55000 0.00000
## cAGE                      0.0005014 0.044090   0.04842  0.01036 0.99170
## cBLCGSMD                  0.0424200 0.005259   0.00428  9.91000 0.00000
## MI_ref10                  0.0394700 0.152800   0.17190  0.22960 0.81840
## CBT1                     -0.1879000 0.131400   0.16990 -1.10600 0.26890
## TIME_ref10               -0.2121000 0.073440   0.09512 -2.23000 0.02575
## MI_ref10:CBT1            -0.0785200 0.184600   0.23410 -0.33540 0.73730
## MI_ref10:TIME_ref10       0.1964000 0.103600   0.12340  1.59200 0.11140
## CBT1:TIME_ref10           0.1805000 0.101900   0.12470  1.44800 0.14770
## MI_ref10:CBT1:TIME_ref10 -0.3240000 0.146100   0.16110 -2.01100 0.04427
## 
##  Estimated Correlation Parameter:  0.7632 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5636 
## 
##  Number of GEE iterations: 8 
##  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.909000 0.096660  0.104400 18.29000 0.00000
## cAGE                -0.003254 0.043760  0.048720 -0.06679 0.94670
## cBLCGSMD             0.041540 0.005173  0.004349  9.55200 0.00000
## MI1                  0.109500 0.137700  0.153400  0.71340 0.47560
## CBT_ref10            0.122500 0.128100  0.196100  0.62490 0.53210
## TIME1                0.154700 0.071760  0.076060  2.03300 0.04201
## MI1:CBT_ref10       -0.264900 0.184100  0.255600 -1.03600 0.30000
## MI1:TIME1           -0.107600 0.100200  0.111300 -0.96700 0.33350
## CBT_ref10:TIME1     -0.126300 0.107400  0.093820 -1.34600 0.17830
## MI1:CBT_ref10:TIME1  0.263800 0.150500  0.157500  1.67500 0.09384
## 
##  Estimated Correlation Parameter:  0.7437 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 4 
##  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.064000 0.104300  0.126600 16.310000 0.00000
## cAGE                     -0.003254 0.043760  0.048720 -0.066790 0.94670
## cBLCGSMD                  0.041540 0.005173  0.004349  9.552000 0.00000
## MI1                       0.001852 0.142500  0.165600  0.011180 0.99110
## CBT_ref10                -0.003770 0.135100  0.193600 -0.019480 0.98450
## TIME_ref10               -0.154700 0.071760  0.076060 -2.033000 0.04201
## MI1:CBT_ref10            -0.001065 0.189000  0.261500 -0.004075 0.99670
## MI1:TIME_ref10            0.107600 0.100200  0.111300  0.967000 0.33350
## CBT_ref10:TIME_ref10      0.126300 0.107400  0.093820  1.346000 0.17830
## MI1:CBT_ref10:TIME_ref10 -0.263800 0.150500  0.157500 -1.675000 0.09384
## 
##  Estimated Correlation Parameter:  0.7437 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 5 
##  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.019000 0.098120  0.112700 17.91000 0.00000
## cAGE                     -0.003254 0.043760  0.048720 -0.06679 0.94670
## cBLCGSMD                  0.041540 0.005173  0.004349  9.55200 0.00000
## MI_ref10                 -0.109500 0.137700  0.153400 -0.71340 0.47560
## CBT_ref10                -0.142400 0.131700  0.164400 -0.86590 0.38650
## TIME1                     0.047050 0.069920  0.081220  0.57930 0.56240
## MI_ref10:CBT_ref10        0.264900 0.184100  0.255600  1.03600 0.30000
## MI_ref10:TIME1            0.107600 0.100200  0.111300  0.96700 0.33350
## CBT_ref10:TIME1           0.137600 0.105500  0.126400  1.08800 0.27640
## MI_ref10:CBT_ref10:TIME1 -0.263800 0.150500  0.157500 -1.67500 0.09384
## 
##  Estimated Correlation Parameter:  0.7437 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 4 
##  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.066000 0.097150  0.107500 19.210000 0.00000
## cAGE                          -0.003254 0.043760  0.048720 -0.066790 0.94670
## cBLCGSMD                       0.041540 0.005173  0.004349  9.552000 0.00000
## MI_ref10                      -0.001852 0.142500  0.165600 -0.011180 0.99110
## CBT_ref10                     -0.004835 0.131600  0.174700 -0.027680 0.97790
## TIME_ref10                    -0.047050 0.069920  0.081220 -0.579300 0.56240
## MI_ref10:CBT_ref10             0.001065 0.189000  0.261500  0.004075 0.99670
## MI_ref10:TIME_ref10           -0.107600 0.100200  0.111300 -0.967000 0.33350
## CBT_ref10:TIME_ref10          -0.137600 0.105500  0.126400 -1.088000 0.27640
## MI_ref10:CBT_ref10:TIME_ref10  0.263800 0.150500  0.157500  1.675000 0.09384
## 
##  Estimated Correlation Parameter:  0.7437 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.5571 
## 
##  Number of GEE iterations: 5 
##  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),
                       corstr = corstr
)


## print main model summary
summary(PCSMKD_nobl_main) #ref: Mi=0, Cbt=0
##             Estimates Model SE Robust SE    wald       p
## (Intercept)  4.371000 0.044350  0.046030 94.9600 0.00000
## cAGE         0.016090 0.021920  0.022860  0.7041 0.48140
## cBLCGSMD     0.006992 0.002615  0.002987  2.3410 0.01923
## MI1          0.089370 0.048310  0.046880  1.9060 0.05660
## CBT1        -0.099800 0.046680  0.052030 -1.9180 0.05507
## TIME1        0.010360 0.023770  0.024890  0.4160 0.67740
## 
##  Estimated Correlation Parameter:  0.609 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1839 
## 
##  Number of GEE iterations: 4 
##  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),
                               corstr = corstr
)

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


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


## print main model summaries
summary(PCSMKD_MI_interactions)
##             Estimates Model SE Robust SE    wald       p
## (Intercept)  4.393000 0.045800  0.046890 93.6800 0.00000
## cAGE         0.015850 0.021980  0.022910  0.6916 0.48920
## cBLCGSMD     0.006978 0.002622  0.002994  2.3300 0.01978
## MI1          0.046730 0.053210  0.052550  0.8893 0.37390
## CBT1        -0.102000 0.046790  0.051990 -1.9610 0.04989
## TIME1       -0.033300 0.033450  0.036840 -0.9039 0.36610
## MI1:TIME1    0.087810 0.047470  0.049750  1.7650 0.07758
## 
##  Estimated Correlation Parameter:  0.6118 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1846 
## 
##  Number of GEE iterations: 4 
##  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.440000 0.046250  0.043540 102.0000 0.00000
## cAGE            0.015850 0.021980  0.022910   0.6916 0.48920
## cBLCGSMD        0.006978 0.002622  0.002994   2.3300 0.01978
## MI_ref10       -0.046730 0.053210  0.052550  -0.8893 0.37390
## TIME1           0.054510 0.033680  0.033440   1.6300 0.10310
## CBT1           -0.102000 0.046790  0.051990  -1.9610 0.04989
## MI_ref10:TIME1 -0.087810 0.047470  0.049750  -1.7650 0.07758
## 
##  Estimated Correlation Parameter:  0.6118 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1846 
## 
##  Number of GEE iterations: 4 
##  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.352000 0.045990  0.048840 89.1000 0.00000
## cAGE         0.015810 0.021980  0.022870  0.6913 0.48930
## cBLCGSMD     0.007024 0.002622  0.002988  2.3510 0.01874
## MI1          0.090370 0.048440  0.046950  1.9250 0.05422
## CBT1        -0.066120 0.051890  0.057360 -1.1530 0.24910
## TIME1        0.051110 0.035230  0.034360  1.4880 0.13680
## CBT1:TIME1  -0.073710 0.048120  0.047900 -1.5390 0.12380
## 
##  Estimated Correlation Parameter:  0.612 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1846 
## 
##  Number of GEE iterations: 4 
##  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),
       corstr = corstr
  )

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

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

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

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

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

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

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

## print interaction model summary
summary(PCSMKD_nobl_interaction) #base: MI=0, CBT=0, time=0
##                Estimates Model SE Robust SE     wald       p
## (Intercept)     4.363000  0.05542  0.060770 71.80000 0.00000
## cAGE            0.016260  0.02218  0.023020  0.70640 0.48000
## cBLCGSMD        0.007137  0.00265  0.002996  2.38200 0.01720
## MI1             0.068330  0.07937  0.079370  0.86080 0.38930
## CBT1           -0.048210  0.07311  0.083860 -0.57490 0.56530
## TIME1           0.049990  0.04884  0.053920  0.92700 0.35390
## MI1:CBT1       -0.039660  0.10450  0.114300 -0.34700 0.72860
## MI1:TIME1       0.002832  0.06995  0.067920  0.04171 0.96670
## CBT1:TIME1     -0.155800  0.06733  0.071500 -2.17900 0.02930
## MI1:CBT1:TIME1  0.158500  0.09555  0.095300  1.66300 0.09628
## 
##  Estimated Correlation Parameter:  0.6209 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1861 
## 
##  Number of GEE iterations: 4 
##  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.413000  0.05729  0.055950 78.88000 0.00000
## cAGE                 0.016260  0.02218  0.023020  0.70640 0.48000
## cBLCGSMD             0.007137  0.00265  0.002996  2.38200 0.01720
## MI1                  0.071160  0.08222  0.070150  1.01400 0.31040
## CBT1                -0.204000  0.07548  0.090500 -2.25500 0.02416
## TIME_ref10          -0.049990  0.04884  0.053920 -0.92700 0.35390
## MI1:CBT1             0.118800  0.10770  0.114900  1.03400 0.30100
## MI1:TIME_ref10      -0.002832  0.06995  0.067920 -0.04171 0.96670
## CBT1:TIME_ref10      0.155800  0.06733  0.071500  2.17900 0.02930
## MI1:CBT1:TIME_ref10 -0.158500  0.09555  0.095300 -1.66300 0.09628
## 
##  Estimated Correlation Parameter:  0.6209 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1861 
## 
##  Number of GEE iterations: 4 
##  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.432000  0.05655  0.051110 86.72000 0.00000
## cAGE                 0.016260  0.02218  0.023020  0.70640 0.48000
## cBLCGSMD             0.007137  0.00265  0.002996  2.38200 0.01720
## MI_ref10            -0.068330  0.07937  0.079370 -0.86080 0.38930
## CBT1                -0.087880  0.07439  0.077790 -1.13000 0.25860
## TIME1                0.052820  0.05008  0.041290  1.27900 0.20080
## MI_ref10:CBT1        0.039660  0.10450  0.114300  0.34700 0.72860
## MI_ref10:TIME1      -0.002832  0.06995  0.067920 -0.04171 0.96670
## CBT1:TIME1           0.002679  0.06780  0.062990  0.04253 0.96610
## MI_ref10:CBT1:TIME1 -0.158500  0.09555  0.095300 -1.66300 0.09628
## 
##  Estimated Correlation Parameter:  0.6209 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1861 
## 
##  Number of GEE iterations: 4 
##  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.485000  0.05861  0.041730 107.50000 0.00000
## cAGE                      0.016260  0.02218  0.023020   0.70640 0.48000
## cBLCGSMD                  0.007137  0.00265  0.002996   2.38200 0.01720
## MI_ref10                 -0.071160  0.08222  0.070150  -1.01400 0.31040
## CBT1                     -0.085200  0.07635  0.070400  -1.21000 0.22620
## TIME_ref10               -0.052820  0.05008  0.041290  -1.27900 0.20080
## MI_ref10:CBT1            -0.118800  0.10770  0.114900  -1.03400 0.30100
## MI_ref10:TIME_ref10       0.002832  0.06995  0.067920   0.04171 0.96670
## CBT1:TIME_ref10          -0.002679  0.06780  0.062990  -0.04253 0.96610
## MI_ref10:CBT1:TIME_ref10  0.158500  0.09555  0.095300   1.66300 0.09628
## 
##  Estimated Correlation Parameter:  0.6209 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1861 
## 
##  Number of GEE iterations: 4 
##  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.31200 0.051220  0.054080 79.7300 0.000000
## cAGE                  0.01720 0.022190  0.023090  0.7447 0.456400
## cBLCGSMD              0.00732 0.002645  0.003005  2.4360 0.014860
## MI1                   0.05612 0.072420  0.074520  0.7530 0.451400
## CBT_ref10             0.05516 0.072720  0.078140  0.7059 0.480200
## TIME1                -0.12420 0.045270  0.051600 -2.4070 0.016080
## MI1:CBT_ref10        -0.02230 0.103400  0.106700 -0.2090 0.834400
## MI1:TIME1             0.17490 0.063760  0.071940  2.4320 0.015020
## CBT_ref10:TIME1       0.19610 0.067500  0.073700  2.6610 0.007781
## MI1:CBT_ref10:TIME1  -0.18910 0.095970  0.097380 -1.9420 0.052160
## 
##  Estimated Correlation Parameter:  0.6204 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1861 
## 
##  Number of GEE iterations: 4 
##  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.18800 0.052750  0.070750 59.1900 0.000000
## cAGE                       0.01720 0.022190  0.023090  0.7447 0.456400
## cBLCGSMD                   0.00732 0.002645  0.003005  2.4360 0.014860
## MI1                        0.23110 0.074020  0.085630  2.6980 0.006966
## CBT_ref10                  0.25130 0.075040  0.086550  2.9040 0.003688
## TIME_ref10                 0.12420 0.045270  0.051600  2.4070 0.016080
## MI1:CBT_ref10             -0.21140 0.106400  0.108000 -1.9580 0.050250
## MI1:TIME_ref10            -0.17490 0.063760  0.071940 -2.4320 0.015020
## CBT_ref10:TIME_ref10      -0.19610 0.067500  0.073700 -2.6610 0.007781
## MI1:CBT_ref10:TIME_ref10   0.18910 0.095970  0.097380  1.9420 0.052160
## 
##  Estimated Correlation Parameter:  0.6204 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1861 
## 
##  Number of GEE iterations: 4 
##  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.36800 0.051170  0.051460 84.8900 0.00000
## cAGE                       0.01720 0.022190  0.023090  0.7447 0.45640
## cBLCGSMD                   0.00732 0.002645  0.003005  2.4360 0.01486
## MI_ref10                  -0.05612 0.072420  0.074520 -0.7530 0.45140
## CBT_ref10                  0.03287 0.073260  0.072440  0.4537 0.65000
## TIME1                      0.05073 0.044890  0.050150  1.0110 0.31180
## MI_ref10:CBT_ref10         0.02230 0.103400  0.106700  0.2090 0.83440
## MI_ref10:TIME1            -0.17490 0.063760  0.071940 -2.4320 0.01502
## CBT_ref10:TIME1            0.00705 0.068210  0.063670  0.1107 0.91180
## MI_ref10:CBT_ref10:TIME1   0.18910 0.095970  0.097380  1.9420 0.05216
## 
##  Estimated Correlation Parameter:  0.6204 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1861 
## 
##  Number of GEE iterations: 4 
##  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.41900 0.051890  0.048630 90.8700 0.000000
## cAGE                            0.01720 0.022190  0.023090  0.7447 0.456400
## cBLCGSMD                        0.00732 0.002645  0.003005  2.4360 0.014860
## MI_ref10                       -0.23110 0.074020  0.085630 -2.6980 0.006966
## CBT_ref10                       0.03992 0.075170  0.064020  0.6235 0.532900
## TIME_ref10                     -0.05073 0.044890  0.050150 -1.0110 0.311800
## MI_ref10:CBT_ref10              0.21140 0.106400  0.108000  1.9580 0.050250
## MI_ref10:TIME_ref10             0.17490 0.063760  0.071940  2.4320 0.015020
## CBT_ref10:TIME_ref10           -0.00705 0.068210  0.063670 -0.1107 0.911800
## MI_ref10:CBT_ref10:TIME_ref10  -0.18910 0.095970  0.097380 -1.9420 0.052160
## 
##  Estimated Correlation Parameter:  0.6204 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1861 
## 
##  Number of GEE iterations: 4 
##  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,
                        corstr = corstr
)

## print model summary
summary(PPVABS_nobl_main)
##             Estimates Model SE Robust SE    wald         p
## (Intercept)  -1.37800  0.24030   0.24270 -5.6780 1.000e-08
## cAGE         -0.35530  0.11790   0.09902 -3.5880 3.328e-04
## cBLCGSMD     -0.02397  0.01528   0.01857 -1.2910 1.968e-01
## MI1          -0.09669  0.25680   0.25360 -0.3812 7.031e-01
## CBT1          0.17600  0.25800   0.26010  0.6767 4.986e-01
## TIME1        -0.04155  0.13650   0.13470 -0.3085 7.577e-01
## 
##  Estimated Correlation Parameter:  0.5563 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.011 
## 
##  Number of GEE iterations: 3 
##  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,
       corstr = corstr
  )

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

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

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

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

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

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

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

## print interaction model summary
summary(PPVABS_nobl_interaction) #base: MI=0, CBT=0, TIME=0
##                Estimates Model SE Robust SE     wald         p
## (Intercept)     -1.35800  0.30110   0.30090 -4.51400 6.370e-06
## cAGE            -0.36400  0.11900   0.09947 -3.65900 2.529e-04
## cBLCGSMD        -0.02180  0.01515   0.01825 -1.19500 2.321e-01
## MI1             -0.40040  0.45840   0.46050 -0.86940 3.847e-01
## CBT1             0.17370  0.40170   0.40410  0.42980 6.674e-01
## TIME1           -0.08657  0.27550   0.26100 -0.33170 7.401e-01
## MI1:CBT1         0.44670  0.59550   0.59450  0.75140 4.524e-01
## MI1:TIME1        0.58750  0.41040   0.45100  1.30300 1.927e-01
## CBT1:TIME1       0.01296  0.37410   0.35690  0.03632 9.710e-01
## MI1:CBT1:TIME1  -0.90330  0.54820   0.56060 -1.61100 1.071e-01
## 
##  Estimated Correlation Parameter:  0.5693 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.012 
## 
##  Number of GEE iterations: 4 
##  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.44500  0.30770   0.31250 -4.62400 3.760e-06
## cAGE                 -0.36400  0.11900   0.09947 -3.65900 2.529e-04
## cBLCGSMD             -0.02180  0.01515   0.01825 -1.19500 2.321e-01
## MI1                   0.18720  0.42790   0.43120  0.43400 6.643e-01
## CBT1                  0.18660  0.40980   0.41380  0.45090 6.520e-01
## TIME_ref10            0.08657  0.27550   0.26100  0.33170 7.401e-01
## MI1:CBT1             -0.45660  0.58880   0.58670 -0.77830 4.364e-01
## MI1:TIME_ref10       -0.58750  0.41040   0.45100 -1.30300 1.927e-01
## CBT1:TIME_ref10      -0.01296  0.37410   0.35690 -0.03632 9.710e-01
## MI1:CBT1:TIME_ref10   0.90330  0.54820   0.56060  1.61100 1.071e-01
## 
##  Estimated Correlation Parameter:  0.5693 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.012 
## 
##  Number of GEE iterations: 4 
##  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.7590  0.34920   0.35790 -4.9140 8.900e-07
## cAGE                  -0.3640  0.11900   0.09947 -3.6590 2.529e-04
## cBLCGSMD              -0.0218  0.01515   0.01825 -1.1950 2.321e-01
## MI_ref10               0.4004  0.45840   0.46050  0.8694 3.847e-01
## CBT1                   0.6203  0.43820   0.43840  1.4150 1.570e-01
## TIME1                  0.5010  0.30420   0.36760  1.3630 1.729e-01
## MI_ref10:CBT1         -0.4467  0.59550   0.59450 -0.7514 4.524e-01
## MI_ref10:TIME1        -0.5875  0.41040   0.45100 -1.3030 1.927e-01
## CBT1:TIME1            -0.8903  0.40070   0.43190 -2.0610 3.926e-02
## MI_ref10:CBT1:TIME1    0.9033  0.54820   0.56060  1.6110 1.071e-01
## 
##  Estimated Correlation Parameter:  0.5693 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.012 
## 
##  Number of GEE iterations: 4 
##  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.2580  0.30030   0.30800 -4.0840 4.421e-05
## cAGE                       -0.3640  0.11900   0.09947 -3.6590 2.529e-04
## cBLCGSMD                   -0.0218  0.01515   0.01825 -1.1950 2.321e-01
## MI_ref10                   -0.1872  0.42790   0.43120 -0.4340 6.643e-01
## CBT1                       -0.2700  0.42040   0.41810 -0.6457 5.185e-01
## TIME_ref10                 -0.5010  0.30420   0.36760 -1.3630 1.729e-01
## MI_ref10:CBT1               0.4566  0.58880   0.58670  0.7783 4.364e-01
## MI_ref10:TIME_ref10         0.5875  0.41040   0.45100  1.3030 1.927e-01
## CBT1:TIME_ref10             0.8903  0.40070   0.43190  2.0610 3.926e-02
## MI_ref10:CBT1:TIME_ref10   -0.9033  0.54820   0.56060 -1.6110 1.071e-01
## 
##  Estimated Correlation Parameter:  0.5693 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.012 
## 
##  Number of GEE iterations: 4 
##  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.27110   0.26960 -4.39500 1.108e-05
## cAGE                 -0.36400  0.11900   0.09947 -3.65900 2.529e-04
## cBLCGSMD             -0.02180  0.01515   0.01825 -1.19500 2.321e-01
## MI1                   0.04630  0.37910   0.37260  0.12430 9.011e-01
## CBT_ref10            -0.17370  0.40170   0.40410 -0.42980 6.674e-01
## TIME1                -0.07361  0.25310   0.24340 -0.30240 7.624e-01
## MI1:CBT_ref10        -0.44670  0.59550   0.59450 -0.75140 4.524e-01
## MI1:TIME1            -0.31580  0.36330   0.33280 -0.94890 3.427e-01
## CBT_ref10:TIME1      -0.01296  0.37410   0.35690 -0.03632 9.710e-01
## MI1:CBT_ref10:TIME1   0.90330  0.54820   0.56060  1.61100 1.071e-01
## 
##  Estimated Correlation Parameter:  0.5693 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.012 
## 
##  Number of GEE iterations: 4 
##  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.27610   0.27090 -4.64600 3.390e-06
## cAGE                      -0.36400  0.11900   0.09947 -3.65900 2.529e-04
## cBLCGSMD                  -0.02180  0.01515   0.01825 -1.19500 2.321e-01
## MI1                       -0.26950  0.40290   0.39420 -0.68360 4.942e-01
## CBT_ref10                 -0.18660  0.40980   0.41380 -0.45090 6.520e-01
## TIME_ref10                 0.07361  0.25310   0.24340  0.30240 7.624e-01
## MI1:CBT_ref10              0.45660  0.58880   0.58670  0.77830 4.364e-01
## MI1:TIME_ref10             0.31580  0.36330   0.33280  0.94890 3.427e-01
## CBT_ref10:TIME_ref10       0.01296  0.37410   0.35690  0.03632 9.710e-01
## MI1:CBT_ref10:TIME_ref10  -0.90330  0.54820   0.56060 -1.61100 1.071e-01
## 
##  Estimated Correlation Parameter:  0.5693 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.012 
## 
##  Number of GEE iterations: 4 
##  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.1380  0.26720   0.25640 -4.4400 9.010e-06
## cAGE                       -0.3640  0.11900   0.09947 -3.6590 2.529e-04
## cBLCGSMD                   -0.0218  0.01515   0.01825 -1.1950 2.321e-01
## MI_ref10                   -0.0463  0.37910   0.37260 -0.1243 9.011e-01
## CBT_ref10                  -0.6203  0.43820   0.43840 -1.4150 1.570e-01
## TIME1                      -0.3894  0.26070   0.22680 -1.7170 8.595e-02
## MI_ref10:CBT_ref10          0.4467  0.59550   0.59450  0.7514 4.524e-01
## MI_ref10:TIME1              0.3158  0.36330   0.33280  0.9489 3.427e-01
## CBT_ref10:TIME1             0.8903  0.40070   0.43190  2.0610 3.926e-02
## MI_ref10:CBT_ref10:TIME1   -0.9033  0.54820   0.56060 -1.6110 1.071e-01
## 
##  Estimated Correlation Parameter:  0.5693 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.012 
## 
##  Number of GEE iterations: 4 
##  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.5280  0.29590   0.28530 -5.3560 9.000e-08
## cAGE                            -0.3640  0.11900   0.09947 -3.6590 2.529e-04
## cBLCGSMD                        -0.0218  0.01515   0.01825 -1.1950 2.321e-01
## MI_ref10                         0.2695  0.40290   0.39420  0.6836 4.942e-01
## CBT_ref10                        0.2700  0.42040   0.41810  0.6457 5.185e-01
## TIME_ref10                       0.3894  0.26070   0.22680  1.7170 8.595e-02
## MI_ref10:CBT_ref10              -0.4566  0.58880   0.58670 -0.7783 4.364e-01
## MI_ref10:TIME_ref10             -0.3158  0.36330   0.33280 -0.9489 3.427e-01
## CBT_ref10:TIME_ref10            -0.8903  0.40070   0.43190 -2.0610 3.926e-02
## MI_ref10:CBT_ref10:TIME_ref10    0.9033  0.54820   0.56060  1.6110 1.071e-01
## 
##  Estimated Correlation Parameter:  0.5693 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  1.012 
## 
##  Number of GEE iterations: 4 
##  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),
                  corstr = corstr
)


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

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

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

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

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

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

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

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

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


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

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

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


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

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

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

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

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



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

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

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

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

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

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

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


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

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

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

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

## fit main model

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

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

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

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



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

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

## main model
LONGABS_nobl_main <- geem(formula = CLONGABS.nobl.main,
                        data = LONGABS_nobl_dt,
                        id=id,
                        family = MASS::negative.binomial(1),
                        corstr = corstr
)
## Warning in geem(formula = CLONGABS.nobl.main, data = LONGABS_nobl_dt, id = id,
## : Did not converge
## print main model summary
summary(LONGABS_nobl_main)
##             Estimates Model SE Robust SE    wald        p
## (Intercept)  -1.97800  0.19070   0.20860 -9.4820 0.000000
## cAGE         -0.14990  0.09251   0.10720 -1.3980 0.162000
## cBLCGSMD     -0.01688  0.01152   0.02085 -0.8097 0.418100
## MI1          -0.60010  0.20470   0.21560 -2.7830 0.005387
## CBT1          0.59740  0.20630   0.22100  2.7030 0.006866
## TIME1        -0.12070  0.12000   0.11970 -1.0080 0.313400
## 
##  Estimated Correlation Parameter:  0.4972 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  2.986 
## 
##  Number of GEE iterations: 20 
##  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),
                               corstr = corstr
)

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


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

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

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

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

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

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


## print interaction model summary
summary(LONGABS_nobl_interaction) #base: MI=0, CBT=0, Time=0
##                Estimates Model SE Robust SE     wald       p
## (Intercept)     1.658000  0.30040    0.2523  6.57100 0.00000
## cAGE           -0.142800  0.11730    0.1151 -1.24000 0.21490
## cBLCGSMD        0.001266  0.01465    0.0151  0.08382 0.93320
## MI1            -0.555400  0.43860    0.3758 -1.47800 0.13940
## CBT1            0.551400  0.40170    0.3158  1.74600 0.08081
## TIME1           0.222800  0.27930    0.2609  0.85390 0.39310
## MI1:CBT1        0.148400  0.57940    0.5181  0.28640 0.77460
## MI1:TIME1      -0.653700  0.42030    0.4683 -1.39600 0.16270
## CBT1:TIME1      0.036570  0.37680    0.3096  0.11810 0.90600
## MI1:CBT1:TIME1  0.193400  0.55550    0.5800  0.33350 0.73880
## 
##  Estimated Correlation Parameter:  0.5538 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.871 
## 
##  Number of GEE iterations: 5 
##  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.881000  0.29560    0.3079  6.10700 0.00000
## cAGE                -0.142800  0.11730    0.1151 -1.24000 0.21490
## cBLCGSMD             0.001266  0.01465    0.0151  0.08382 0.93320
## MI1                 -1.209000  0.45460    0.5246 -2.30500 0.02118
## CBT1                 0.588000  0.39710    0.3539  1.66100 0.09662
## TIME_ref10          -0.222800  0.27930    0.2609 -0.85390 0.39310
## MI1:CBT1             0.341800  0.59390    0.6597  0.51810 0.60440
## MI1:TIME_ref10       0.653700  0.42030    0.4683  1.39600 0.16270
## CBT1:TIME_ref10     -0.036570  0.37680    0.3096 -0.11810 0.90600
## MI1:CBT1:TIME_ref10 -0.193400  0.55550    0.5800 -0.33350 0.73880
## 
##  Estimated Correlation Parameter:  0.5538 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.871 
## 
##  Number of GEE iterations: 5 
##  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.103000  0.31910    0.2746  4.01500 5.946e-05
## cAGE                -0.142800  0.11730    0.1151 -1.24000 2.149e-01
## cBLCGSMD             0.001266  0.01465    0.0151  0.08382 9.332e-01
## MI_ref10             0.555400  0.43860    0.3758  1.47800 1.394e-01
## CBT1                 0.699800  0.41640    0.4074  1.71800 8.583e-02
## TIME1               -0.430900  0.31410    0.3887 -1.10800 2.677e-01
## MI_ref10:CBT1       -0.148400  0.57940    0.5181 -0.28640 7.746e-01
## MI_ref10:TIME1       0.653700  0.42030    0.4683  1.39600 1.627e-01
## CBT1:TIME1           0.230000  0.40810    0.4908  0.46860 6.393e-01
## MI_ref10:CBT1:TIME1 -0.193400  0.55550    0.5800 -0.33350 7.388e-01
## 
##  Estimated Correlation Parameter:  0.5538 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.871 
## 
##  Number of GEE iterations: 5 
##  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.671700  0.34520    0.4084  1.64500 0.10010
## cAGE                     -0.142800  0.11730    0.1151 -1.24000 0.21490
## cBLCGSMD                  0.001266  0.01465    0.0151  0.08382 0.93320
## MI_ref10                  1.209000  0.45460    0.5246  2.30500 0.02118
## CBT1                      0.929800  0.44070    0.5465  1.70100 0.08885
## TIME_ref10                0.430900  0.31410    0.3887  1.10800 0.26770
## MI_ref10:CBT1            -0.341800  0.59390    0.6597 -0.51810 0.60440
## MI_ref10:TIME_ref10      -0.653700  0.42030    0.4683 -1.39600 0.16270
## CBT1:TIME_ref10          -0.230000  0.40810    0.4908 -0.46860 0.63930
## MI_ref10:CBT1:TIME_ref10  0.193400  0.55550    0.5800  0.33350 0.73880
## 
##  Estimated Correlation Parameter:  0.5538 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.871 
## 
##  Number of GEE iterations: 5 
##  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.2080000  0.26660   0.19320 11.43000 0.00000
## cAGE                -0.1429000  0.11660   0.11540 -1.23900 0.21540
## cBLCGSMD             0.0008395  0.01458   0.01514  0.05545 0.95580
## MI1                 -0.3857000  0.37940   0.34980 -1.10200 0.27030
## CBT_ref10           -0.5545000  0.40110   0.31780 -1.74500 0.08101
## TIME1                0.2601000  0.25240   0.16560  1.57000 0.11630
## MI1:CBT_ref10       -0.1882000  0.57840   0.51970 -0.36210 0.71730
## MI1:TIME1           -0.4446000  0.36180   0.33120 -1.34200 0.17950
## CBT_ref10:TIME1     -0.0394600  0.37680   0.30960 -0.12750 0.89860
## MI1:CBT_ref10:TIME1 -0.2289000  0.55630   0.58540 -0.39100 0.69580
## 
##  Estimated Correlation Parameter:  0.5489 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.827 
## 
##  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.4690000  0.26510   0.17800 13.87000 0.00000
## cAGE                     -0.1429000  0.11660   0.11540 -1.23900 0.21540
## cBLCGSMD                  0.0008395  0.01458   0.01514  0.05545 0.95580
## MI1                      -0.8302000  0.38320   0.38710 -2.14500 0.03197
## CBT_ref10                -0.5940000  0.39640   0.35610 -1.66800 0.09529
## TIME_ref10               -0.2601000  0.25240   0.16560 -1.57000 0.11630
## MI1:CBT_ref10            -0.4171000  0.59360   0.66640 -0.62580 0.53140
## MI1:TIME_ref10            0.4446000  0.36180   0.33120  1.34200 0.17950
## CBT_ref10:TIME_ref10      0.0394600  0.37680   0.30960  0.12750 0.89860
## MI1:CBT_ref10:TIME_ref10  0.2289000  0.55630   0.58540  0.39100 0.69580
## 
##  Estimated Correlation Parameter:  0.5489 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.827 
## 
##  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.8230000  0.26990   0.29220  6.23700 0.00000
## cAGE                     -0.1429000  0.11660   0.11540 -1.23900 0.21540
## cBLCGSMD                  0.0008395  0.01458   0.01514  0.05545 0.95580
## MI_ref10                  0.3857000  0.37940   0.34980  1.10200 0.27030
## CBT_ref10                -0.7427000  0.41560   0.40820 -1.82000 0.06882
## TIME1                    -0.1845000  0.25920   0.28710 -0.64260 0.52050
## MI_ref10:CBT_ref10        0.1882000  0.57840   0.51970  0.36210 0.71730
## MI_ref10:TIME1            0.4446000  0.36180   0.33120  1.34200 0.17950
## CBT_ref10:TIME1          -0.2684000  0.40930   0.49750 -0.53940 0.58960
## MI_ref10:CBT_ref10:TIME1  0.2289000  0.55630   0.58540  0.39100 0.69580
## 
##  Estimated Correlation Parameter:  0.5489 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.827 
## 
##  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.6380000  0.27640   0.34600  4.73500 2.190e-06
## cAGE                          -0.1429000  0.11660   0.11540 -1.23900 2.154e-01
## cBLCGSMD                       0.0008395  0.01458   0.01514  0.05545 9.558e-01
## MI_ref10                       0.8302000  0.38320   0.38710  2.14500 3.197e-02
## CBT_ref10                     -1.0110000  0.44070   0.55330 -1.82700 6.765e-02
## TIME_ref10                     0.1845000  0.25920   0.28710  0.64260 5.205e-01
## MI_ref10:CBT_ref10             0.4171000  0.59360   0.66640  0.62580 5.314e-01
## MI_ref10:TIME_ref10           -0.4446000  0.36180   0.33120 -1.34200 1.795e-01
## CBT_ref10:TIME_ref10           0.2684000  0.40930   0.49750  0.53940 5.896e-01
## MI_ref10:CBT_ref10:TIME_ref10 -0.2289000  0.55630   0.58540 -0.39100 6.958e-01
## 
##  Estimated Correlation Parameter:  0.5489 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  4.827 
## 
##  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")