rm(list=ls())


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

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


# Data location for reading datasets and storing data ---------------------------

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

# Read Data ---------------------------

raw_dt_ncigs_smkd <- read.csv("/Volumes/caas/MERITS/RESOURCES/Data Processing/Data Analysis/Outcome Analyses 2016/Manuscript Revisions_2020/Aditya Khanna Analyses 2022/NCIGSSMDK.csv")


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


# Create covariate data objects ---------------------------

covariates_dt <- 
  cbind.data.frame(
    cAGE = raw_dt_ncigs_smkd$AGE-mean(raw_dt_ncigs_smkd$AGE),
    MI = as.factor(raw_dt_ncigs_smkd$MI),
    CBT = as.factor(raw_dt_ncigs_smkd$CBT),
    TIME=as.factor(raw_dt_ncigs_smkd$TIME),
    #CBT=as.factor(raw_dt_ncigs_smkd$CBT),
    M3ENVDAYS = raw_dt_ncigs_smkd$M3ENVDAYS,
    M6ENVDAYS = raw_dt_ncigs_smkd$M6ENVDAYS ,
    MIBYCBT = as.factor(raw_dt_ncigs_smkd$MIBYCBT),
    MIBYTIME = as.factor(raw_dt_ncigs_smkd$MIBYTIME),
    CBTBYTIME = as.factor(raw_dt_ncigs_smkd$CBTBYTIME),
    MIBYCBTBYTIME = as.factor(raw_dt_ncigs_smkd$MIBYCBTBYTIME),
    id = raw_dt_ncigs_smkd$id,
    log_unctrldays = raw_dt_ncigs_smkd$log_unctrldays
  )


# Create dataframe ---------------------------
ncigs_dt <- as.data.frame(
  cbind.data.frame(NCIGS=raw_dt_ncigs_smkd$NCIGS,
                   covariates_dt)
)

class(ncigs_dt)
## [1] "data.frame"
dim(ncigs_dt)
## [1] 729  13
summary(ncigs_dt)  
##      NCIGS             cAGE         MI      CBT     TIME      M3ENVDAYS    
##  Min.   :   0.0   Min.   :-2.8828   0:366   0:327   0:243   Min.   : 0.00  
##  1st Qu.:  85.0   1st Qu.:-0.8253   1:363   1:402   1:243   1st Qu.: 0.00  
##  Median : 213.0   Median : 0.1144                   2:243   Median : 1.00  
##  Mean   : 329.4   Mean   : 0.0000                           Mean   :17.71  
##  3rd Qu.: 460.0   3rd Qu.: 0.8103                           3rd Qu.:35.00  
##  Max.   :2400.0   Max.   : 2.3363                           Max.   :58.00  
##    M6ENVDAYS     MIBYCBT MIBYTIME CBTBYTIME MIBYCBTBYTIME       id       
##  Min.   : 0.00   0:525   0:487    0:461     0:593         Min.   :  1.0  
##  1st Qu.: 0.00   1:204   1:121    1:134     1: 68         1st Qu.: 86.0  
##  Median : 0.00           2:121    2:134     2: 68         Median :163.0  
##  Mean   :15.23                                            Mean   :162.8  
##  3rd Qu.:34.00                                            3rd Qu.:240.0  
##  Max.   :58.00                                            Max.   :316.0  
##  log_unctrldays 
##  Min.   :1.099  
##  1st Qu.:3.332  
##  Median :3.434  
##  Mean   :3.560  
##  3rd Qu.:4.111  
##  Max.   :4.111
colnames(ncigs_dt)
##  [1] "NCIGS"          "cAGE"           "MI"             "CBT"           
##  [5] "TIME"           "M3ENVDAYS"      "M6ENVDAYS"      "MIBYCBT"       
##  [9] "MIBYTIME"       "CBTBYTIME"      "MIBYCBTBYTIME"  "id"            
## [13] "log_unctrldays"
## restrict to complete cases
ncigs_dt_na.omit <- na.omit(ncigs_dt)  
dim(ncigs_dt_na.omit)  
## [1] 729  13
summary(ncigs_dt_na.omit)
##      NCIGS             cAGE         MI      CBT     TIME      M3ENVDAYS    
##  Min.   :   0.0   Min.   :-2.8828   0:366   0:327   0:243   Min.   : 0.00  
##  1st Qu.:  85.0   1st Qu.:-0.8253   1:363   1:402   1:243   1st Qu.: 0.00  
##  Median : 213.0   Median : 0.1144                   2:243   Median : 1.00  
##  Mean   : 329.4   Mean   : 0.0000                           Mean   :17.71  
##  3rd Qu.: 460.0   3rd Qu.: 0.8103                           3rd Qu.:35.00  
##  Max.   :2400.0   Max.   : 2.3363                           Max.   :58.00  
##    M6ENVDAYS     MIBYCBT MIBYTIME CBTBYTIME MIBYCBTBYTIME       id       
##  Min.   : 0.00   0:525   0:487    0:461     0:593         Min.   :  1.0  
##  1st Qu.: 0.00   1:204   1:121    1:134     1: 68         1st Qu.: 86.0  
##  Median : 0.00           2:121    2:134     2: 68         Median :163.0  
##  Mean   :15.23                                            Mean   :162.8  
##  3rd Qu.:34.00                                            3rd Qu.:240.0  
##  Max.   :58.00                                            Max.   :316.0  
##  log_unctrldays 
##  Min.   :1.099  
##  1st Qu.:3.332  
##  Median :3.434  
##  Mean   :3.560  
##  3rd Qu.:4.111  
##  Max.   :4.111
colnames(ncigs_dt_na.omit)
##  [1] "NCIGS"          "cAGE"           "MI"             "CBT"           
##  [5] "TIME"           "M3ENVDAYS"      "M6ENVDAYS"      "MIBYCBT"       
##  [9] "MIBYTIME"       "CBTBYTIME"      "MIBYCBTBYTIME"  "id"            
## [13] "log_unctrldays"
# DV: NCIGS ---------------------------

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

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

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


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

# modeling
## fit main model
ncigs_main <- geem(formula =   
                     NCIGS~
                     cAGE +
                     MI +
                     CBT +
                     TIME + 
                     offset( log_unctrldays),
                   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.369000  0.10060   0.08185 28.94000 0.0000
## cAGE         0.104500  0.04729   0.04654  2.24500 0.0248
## MI1          0.082880  0.10640   0.09894  0.83770 0.4022
## CBT1         0.003578  0.10680   0.10200  0.03509 0.9720
## TIME1       -0.540100  0.05169   0.06700 -8.06100 0.0000
## TIME2       -0.448700  0.06663   0.06756 -6.64200 0.0000
## 
##  Estimated Correlation Parameter:  0.6627 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.958 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
## fit 3-way interaction model
ncigs_interaction <- geem(formula =   
                            NCIGS~
                            cAGE +
                            MI +
                            CBT +
                            TIME + 
                            MI*CBT+
                            MI*TIME+
                            CBT*TIME+
                            MI*CBT*TIME+
                            offset( log_unctrldays),
                          #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.29100  0.13100   0.10190 22.48000 0.000000
## cAGE             0.10340  0.04746   0.04657  2.22000 0.026430
## MI1              0.19110  0.18790   0.14440  1.32300 0.185700
## CBT1             0.19190  0.17800   0.15030  1.27700 0.201500
## TIME1           -0.49560  0.10610   0.15700 -3.15800 0.001591
## TIME2           -0.40790  0.13710   0.14550 -2.80400 0.005054
## MI1:CBT1        -0.29790  0.25270   0.20160 -1.47800 0.139400
## MI1:TIME1       -0.06360  0.15210   0.20490 -0.31030 0.756300
## MI1:TIME2        0.01873  0.19650   0.19490  0.09607 0.923500
## CBT1:TIME1      -0.16070  0.14420   0.20590 -0.78070 0.435000
## CBT1:TIME2      -0.16360  0.18630   0.20460 -0.79970 0.423900
## MI1:CBT1:TIME1   0.28210  0.20480   0.26930  1.04700 0.295000
## MI1:CBT1:TIME2   0.15070  0.26460   0.26780  0.56290 0.573500
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_TIMEref1) #ref: MI=0, CBT=0, time=1
##                     Estimates Model SE Robust SE     wald        p
## (Intercept)          1.796000  0.13120   0.15090 11.90000 0.000000
## cAGE                 0.103400  0.04746   0.04657  2.22000 0.026430
## MI1                  0.127500  0.18820   0.22660  0.56270 0.573600
## CBT1                 0.031210  0.17830   0.20120  0.15510 0.876700
## TIME_ref10           0.495600  0.10610   0.15700  3.15800 0.001591
## TIME_ref12           0.087670  0.10610   0.06703  1.30800 0.190900
## MI1:CBT1            -0.015880  0.25300   0.28630 -0.05547 0.955800
## MI1:TIME_ref10       0.063600  0.15210   0.20490  0.31030 0.756300
## MI1:TIME_ref12       0.082330  0.15220   0.12680  0.64950 0.516000
## CBT1:TIME_ref10      0.160700  0.14420   0.20590  0.78070 0.435000
## CBT1:TIME_ref12     -0.002899  0.14430   0.09798 -0.02959 0.976400
## MI1:CBT1:TIME_ref10 -0.282100  0.20480   0.26930 -1.04700 0.295000
## MI1:CBT1:TIME_ref12 -0.131300  0.20490   0.16750 -0.78410 0.433000
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_TIMEref2) #ref: MI=0, CBT=0, time=2
##                     Estimates Model SE Robust SE     wald        p
## (Intercept)          1.883000  0.13110   0.13190 14.28000 0.000000
## cAGE                 0.103400  0.04746   0.04657  2.22000 0.026430
## MI1                  0.209800  0.18800   0.20490  1.02400 0.305700
## CBT1                 0.028310  0.17810   0.18950  0.14940 0.881200
## TIME_ref20           0.407900  0.13710   0.14550  2.80400 0.005054
## TIME_ref21          -0.087670  0.10610   0.06703 -1.30800 0.190900
## MI1:CBT1            -0.147200  0.25280   0.26910 -0.54710 0.584300
## MI1:TIME_ref20      -0.018730  0.19650   0.19490 -0.09607 0.923500
## MI1:TIME_ref21      -0.082330  0.15220   0.12680 -0.64950 0.516000
## CBT1:TIME_ref20      0.163600  0.18630   0.20460  0.79970 0.423900
## CBT1:TIME_ref21      0.002899  0.14430   0.09798  0.02959 0.976400
## MI1:CBT1:TIME_ref20 -0.150700  0.26460   0.26780 -0.56290 0.573500
## MI1:CBT1:TIME_ref21  0.131300  0.20490   0.16750  0.78410 0.433000
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_MIref1) #ref: MI=1, CBT=0, time=0
##                     Estimates Model SE Robust SE     wald         p
## (Intercept)           2.48200  0.13430   0.10250 24.21000 0.000e+00
## cAGE                  0.10340  0.04746   0.04657  2.22000 2.643e-02
## MI_ref10             -0.19110  0.18790   0.14440 -1.32300 1.857e-01
## CBT1                 -0.10600  0.17910   0.13400 -0.79110 4.289e-01
## TIME1                -0.55920  0.10900   0.13180 -4.24400 2.193e-05
## TIME2                -0.38920  0.14080   0.12970 -3.00000 2.699e-03
## MI_ref10:CBT1         0.29790  0.25270   0.20160  1.47800 1.394e-01
## MI_ref10:TIME1        0.06360  0.15210   0.20490  0.31030 7.563e-01
## MI_ref10:TIME2       -0.01873  0.19650   0.19490 -0.09607 9.235e-01
## CBT1:TIME1            0.12130  0.14540   0.17370  0.69860 4.848e-01
## CBT1:TIME2           -0.01291  0.18790   0.17270 -0.07474 9.404e-01
## MI_ref10:CBT1:TIME1  -0.28210  0.20480   0.26930 -1.04700 2.950e-01
## MI_ref10:CBT1:TIME2  -0.15070  0.26460   0.26780 -0.56290 5.735e-01
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, time=1
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)                1.92300  0.13450   0.16940 11.36000 0.000e+00
## cAGE                       0.10340  0.04746   0.04657  2.22000 2.643e-02
## MI_ref10                  -0.12750  0.18820   0.22660 -0.56270 5.736e-01
## CBT1                       0.01533  0.17930   0.20380  0.07525 9.400e-01
## TIME_ref10                 0.55920  0.10900   0.13180  4.24400 2.193e-05
## TIME_ref12                 0.17000  0.10900   0.10760  1.58000 1.140e-01
## MI_ref10:CBT1              0.01588  0.25300   0.28630  0.05547 9.558e-01
## MI_ref10:TIME_ref10       -0.06360  0.15210   0.20490 -0.31030 7.563e-01
## MI_ref10:TIME_ref12       -0.08233  0.15220   0.12680 -0.64950 5.160e-01
## CBT1:TIME_ref10           -0.12130  0.14540   0.17370 -0.69860 4.848e-01
## CBT1:TIME_ref12           -0.13420  0.14550   0.13590 -0.98810 3.231e-01
## MI_ref10:CBT1:TIME_ref10   0.28210  0.20480   0.26930  1.04700 2.950e-01
## MI_ref10:CBT1:TIME_ref12   0.13130  0.20490   0.16750  0.78410 4.330e-01
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_MIref1_TIMEref2) #ref: MI=1, CBT=0, time=2
##                          Estimates Model SE Robust SE     wald        p
## (Intercept)                2.09300  0.13440   0.15600 13.42000 0.000000
## cAGE                       0.10340  0.04746   0.04657  2.22000 0.026430
## MI_ref10                  -0.20980  0.18800   0.20490 -1.02400 0.305700
## CBT1                      -0.11890  0.17920   0.18980 -0.62660 0.530900
## TIME_ref20                 0.38920  0.14080   0.12970  3.00000 0.002699
## TIME_ref21                -0.17000  0.10900   0.10760 -1.58000 0.114000
## MI_ref10:CBT1              0.14720  0.25280   0.26910  0.54710 0.584300
## MI_ref10:TIME_ref20        0.01873  0.19650   0.19490  0.09607 0.923500
## MI_ref10:TIME_ref21        0.08233  0.15220   0.12680  0.64950 0.516000
## CBT1:TIME_ref20            0.01291  0.18790   0.17270  0.07474 0.940400
## CBT1:TIME_ref21            0.13420  0.14550   0.13590  0.98810 0.323100
## MI_ref10:CBT1:TIME_ref20   0.15070  0.26460   0.26780  0.56290 0.573500
## MI_ref10:CBT1:TIME_ref21  -0.13130  0.20490   0.16750 -0.78410 0.433000
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_CBTref1) #ref: MI=0, CBT=1, time=0
##                     Estimates Model SE Robust SE    wald         p
## (Intercept)            2.4830  0.12030   0.10810 22.9800 0.000e+00
## cAGE                   0.1034  0.04746   0.04657  2.2200 2.643e-02
## MI1                   -0.1068  0.16880   0.13770 -0.7758 4.379e-01
## CBT_ref10             -0.1919  0.17800   0.15030 -1.2770 2.015e-01
## TIME1                 -0.6563  0.09768   0.13320 -4.9260 8.400e-07
## TIME2                 -0.5716  0.12620   0.14380 -3.9740 7.080e-05
## MI1:CBT_ref10          0.2979  0.25270   0.20160  1.4780 1.394e-01
## MI1:TIME1              0.2185  0.13710   0.17480  1.2500 2.114e-01
## MI1:TIME2              0.1694  0.17720   0.18350  0.9232 3.559e-01
## CBT_ref10:TIME1        0.1607  0.14420   0.20590  0.7807 4.350e-01
## CBT_ref10:TIME2        0.1636  0.18630   0.20460  0.7997 4.239e-01
## MI1:CBT_ref10:TIME1   -0.2821  0.20480   0.26930 -1.0470 2.950e-01
## MI1:CBT_ref10:TIME2   -0.1507  0.26460   0.26780 -0.5629 5.735e-01
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, time=1
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)               1.827000  0.12050   0.13170 13.87000 0.000e+00
## cAGE                      0.103400  0.04746   0.04657  2.22000 2.643e-02
## MI1                       0.111600  0.16900   0.17340  0.64380 5.197e-01
## CBT_ref10                -0.031210  0.17830   0.20120 -0.15510 8.767e-01
## TIME_ref10                0.656300  0.09768   0.13320  4.92600 8.400e-07
## TIME_ref12                0.084770  0.09772   0.07145  1.18600 2.355e-01
## MI1:CBT_ref10             0.015880  0.25300   0.28630  0.05547 9.558e-01
## MI1:TIME_ref10           -0.218500  0.13710   0.17480 -1.25000 2.114e-01
## MI1:TIME_ref12           -0.049010  0.13720   0.10950 -0.44760 6.544e-01
## CBT_ref10:TIME_ref10     -0.160700  0.14420   0.20590 -0.78070 4.350e-01
## CBT_ref10:TIME_ref12      0.002899  0.14430   0.09798  0.02959 9.764e-01
## MI1:CBT_ref10:TIME_ref10  0.282100  0.20480   0.26930  1.04700 2.950e-01
## MI1:CBT_ref10:TIME_ref12  0.131300  0.20490   0.16750  0.78410 4.330e-01
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_CBTref1_TIMEref2) #ref: MI=0, CBT=1, time=2
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)               1.912000  0.12040   0.13430 14.23000 0.0000000
## cAGE                      0.103400  0.04746   0.04657  2.22000 0.0264300
## MI1                       0.062610  0.16900   0.17210  0.36380 0.7160000
## CBT_ref10                -0.028310  0.17810   0.18950 -0.14940 0.8812000
## TIME_ref20                0.571600  0.12620   0.14380  3.97400 0.0000708
## TIME_ref21               -0.084770  0.09772   0.07145 -1.18600 0.2355000
## MI1:CBT_ref10             0.147200  0.25280   0.26910  0.54710 0.5843000
## MI1:TIME_ref20           -0.169400  0.17720   0.18350 -0.92320 0.3559000
## MI1:TIME_ref21            0.049010  0.13720   0.10950  0.44760 0.6544000
## CBT_ref10:TIME_ref20     -0.163600  0.18630   0.20460 -0.79970 0.4239000
## CBT_ref10:TIME_ref21     -0.002899  0.14430   0.09798 -0.02959 0.9764000
## MI1:CBT_ref10:TIME_ref20  0.150700  0.26460   0.26780  0.56290 0.5735000
## MI1:CBT_ref10:TIME_ref21 -0.131300  0.20490   0.16750 -0.78410 0.4330000
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, time=0
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)                2.37600  0.11850   0.08570 27.73000 0.0000000
## cAGE                       0.10340  0.04746   0.04657  2.22000 0.0264300
## MI_ref10                   0.10680  0.16880   0.13770  0.77580 0.4379000
## CBT_ref10                  0.10600  0.17910   0.13400  0.79110 0.4289000
## TIME1                     -0.43790  0.09622   0.11310 -3.87000 0.0001089
## TIME2                     -0.40210  0.12440   0.11400 -3.52700 0.0004204
## MI_ref10:CBT_ref10        -0.29790  0.25270   0.20160 -1.47800 0.1394000
## MI_ref10:TIME1            -0.21850  0.13710   0.17480 -1.25000 0.2114000
## MI_ref10:TIME2            -0.16940  0.17720   0.18350 -0.92320 0.3559000
## CBT_ref10:TIME1           -0.12130  0.14540   0.17370 -0.69860 0.4848000
## CBT_ref10:TIME2            0.01291  0.18790   0.17270  0.07474 0.9404000
## MI_ref10:CBT_ref10:TIME1   0.28210  0.20480   0.26930  1.04700 0.2950000
## MI_ref10:CBT_ref10:TIME2   0.15070  0.26460   0.26780  0.56290 0.5735000
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, time=1
##                               Estimates Model SE Robust SE     wald         p
## (Intercept)                     1.93900  0.11860   0.11300 17.16000 0.0000000
## cAGE                            0.10340  0.04746   0.04657  2.22000 0.0264300
## MI_ref10                       -0.11160  0.16900   0.17340 -0.64380 0.5197000
## CBT_ref10                      -0.01533  0.17930   0.20380 -0.07525 0.9400000
## TIME_ref10                      0.43790  0.09622   0.11310  3.87000 0.0001089
## TIME_ref12                      0.03575  0.09626   0.08297  0.43090 0.6665000
## MI_ref10:CBT_ref10             -0.01588  0.25300   0.28630 -0.05547 0.9558000
## MI_ref10:TIME_ref10             0.21850  0.13710   0.17480  1.25000 0.2114000
## MI_ref10:TIME_ref12             0.04901  0.13720   0.10950  0.44760 0.6544000
## CBT_ref10:TIME_ref10            0.12130  0.14540   0.17370  0.69860 0.4848000
## CBT_ref10:TIME_ref12            0.13420  0.14550   0.13590  0.98810 0.3231000
## MI_ref10:CBT_ref10:TIME_ref10  -0.28210  0.20480   0.26930 -1.04700 0.2950000
## MI_ref10:CBT_ref10:TIME_ref12  -0.13130  0.20490   0.16750 -0.78410 0.4330000
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(ncigs_interaction_MIref1_CBTref1_TIMEref2) #ref: MI=1, CBT=1, time=2
##                               Estimates Model SE Robust SE     wald         p
## (Intercept)                     1.97400  0.11860   0.10790 18.31000 0.0000000
## cAGE                            0.10340  0.04746   0.04657  2.22000 0.0264300
## MI_ref10                       -0.06261  0.16900   0.17210 -0.36380 0.7160000
## CBT_ref10                       0.11890  0.17920   0.18980  0.62660 0.5309000
## TIME_ref20                      0.40210  0.12440   0.11400  3.52700 0.0004204
## TIME_ref21                     -0.03575  0.09626   0.08297 -0.43090 0.6665000
## MI_ref10:CBT_ref10             -0.14720  0.25280   0.26910 -0.54710 0.5843000
## MI_ref10:TIME_ref20             0.16940  0.17720   0.18350  0.92320 0.3559000
## MI_ref10:TIME_ref21            -0.04901  0.13720   0.10950 -0.44760 0.6544000
## CBT_ref10:TIME_ref20           -0.01291  0.18790   0.17270 -0.07474 0.9404000
## CBT_ref10:TIME_ref21           -0.13420  0.14550   0.13590 -0.98810 0.3231000
## MI_ref10:CBT_ref10:TIME_ref20  -0.15070  0.26460   0.26780 -0.56290 0.5735000
## MI_ref10:CBT_ref10:TIME_ref21   0.13130  0.20490   0.16750  0.78410 0.4330000
## 
##  Estimated Correlation Parameter:  0.6707 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.9519 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729