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
  )

# DV V: SMKD ---------------------------

## create dataframe
SMDK_dt <- as.data.frame(
  cbind.data.frame(SMDK=raw_dt_ncigs_smkd$SMDK,
                   covariates_dt)
)

class(SMDK_dt)
## [1] "data.frame"
dim(SMDK_dt)
## [1] 729  13
summary(SMDK_dt)  
##       SMDK            cAGE         MI      CBT     TIME      M3ENVDAYS    
##  Min.   : 0.00   Min.   :-2.8828   0:366   0:327   0:243   Min.   : 0.00  
##  1st Qu.:23.00   1st Qu.:-0.8253   1:363   1:402   1:243   1st Qu.: 0.00  
##  Median :30.00   Median : 0.1144                   2:243   Median : 1.00  
##  Mean   :32.21   Mean   : 0.0000                           Mean   :17.71  
##  3rd Qu.:50.00   3rd Qu.: 0.8103                           3rd Qu.:35.00  
##  Max.   :60.00   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(SMDK_dt)
##  [1] "SMDK"           "cAGE"           "MI"             "CBT"           
##  [5] "TIME"           "M3ENVDAYS"      "M6ENVDAYS"      "MIBYCBT"       
##  [9] "MIBYTIME"       "CBTBYTIME"      "MIBYCBTBYTIME"  "id"            
## [13] "log_unctrldays"
## restrict to complete cases
SMDK_dt_na.omit <- na.omit(SMDK_dt)  
dim(SMDK_dt_na.omit)  
## [1] 729  13
summary(SMDK_dt_na.omit)
##       SMDK            cAGE         MI      CBT     TIME      M3ENVDAYS    
##  Min.   : 0.00   Min.   :-2.8828   0:366   0:327   0:243   Min.   : 0.00  
##  1st Qu.:23.00   1st Qu.:-0.8253   1:363   1:402   1:243   1st Qu.: 0.00  
##  Median :30.00   Median : 0.1144                   2:243   Median : 1.00  
##  Mean   :32.21   Mean   : 0.0000                           Mean   :17.71  
##  3rd Qu.:50.00   3rd Qu.: 0.8103                           3rd Qu.:35.00  
##  Max.   :60.00   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(SMDK_dt_na.omit)
##  [1] "SMDK"           "cAGE"           "MI"             "CBT"           
##  [5] "TIME"           "M3ENVDAYS"      "M6ENVDAYS"      "MIBYCBT"       
##  [9] "MIBYTIME"       "CBTBYTIME"      "MIBYCBTBYTIME"  "id"            
## [13] "log_unctrldays"
# DV: SMDK  ---------------------------

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

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

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


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


# modeling
## fit main model
SMDK_main <- geem(formula =   
                    SMDK~
                    cAGE +
                    MI +
                    CBT +
                    TIME + 
                    offset( log_unctrldays),
                  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.09342  0.03467   0.02591 -3.606 0.000311
## cAGE          0.02584  0.01571   0.01412  1.830 0.067220
## MI1           0.07370  0.03531   0.03118  2.363 0.018110
## CBT1         -0.03342  0.03546   0.03130 -1.068 0.285700
## TIME1        -0.20850  0.02321   0.02811 -7.417 0.000000
## TIME2        -0.17740  0.02825   0.02706 -6.559 0.000000
## 
##  Estimated Correlation Parameter:  0.4875 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1225 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
## fit interaction model
SMDK_interaction <- geem(formula =   
                           SMDK~
                           cAGE +
                           MI +
                           CBT +
                           TIME + 
                           MI*CBT+
                           MI*TIME+
                           CBT*TIME+
                           MI*CBT*TIME+
                           offset( log_unctrldays),
                         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.09526  0.04812   0.02427 -3.9250 8.684e-05
## cAGE             0.02696  0.01594   0.01444  1.8670 6.191e-02
## MI1              0.03962  0.06903   0.02776  1.4270 1.536e-01
## CBT1             0.02362  0.06539   0.03073  0.7689 4.420e-01
## TIME1           -0.22940  0.04847   0.06525 -3.5150 4.398e-04
## TIME2           -0.14630  0.05901   0.05720 -2.5580 1.054e-02
## MI1:CBT1        -0.04718  0.09281   0.03847 -1.2260 2.201e-01
## MI1:TIME1        0.05159  0.06945   0.08027  0.6427 5.204e-01
## MI1:TIME2        0.02514  0.08462   0.06868  0.3660 7.143e-01
## CBT1:TIME1      -0.01412  0.06581   0.08549 -0.1652 8.688e-01
## CBT1:TIME2      -0.17110  0.08028   0.08717 -1.9620 4.973e-02
## MI1:CBT1:TIME1   0.01216  0.09340   0.11150  0.1090 9.132e-01
## MI1:CBT1:TIME2   0.17840  0.11400   0.10560  1.6890 9.119e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_TIMEref1) #ref: MI=0, CBT=0, time=1
##                     Estimates Model SE Robust SE    wald         p
## (Intercept)         -0.324600  0.04850   0.06553 -4.9540 7.300e-07
## cAGE                 0.026960  0.01594   0.01444  1.8670 6.191e-02
## MI1                  0.091210  0.06949   0.08152  1.1190 2.632e-01
## CBT1                 0.009504  0.06579   0.08600  0.1105 9.120e-01
## TIME_ref10           0.229400  0.04847   0.06525  3.5150 4.398e-04
## TIME_ref12           0.083070  0.04847   0.05293  1.5690 1.166e-01
## MI1:CBT1            -0.035020  0.09330   0.11180 -0.3133 7.541e-01
## MI1:TIME_ref10      -0.051590  0.06945   0.08027 -0.6427 5.204e-01
## MI1:TIME_ref12      -0.026450  0.06947   0.06539 -0.4045 6.858e-01
## CBT1:TIME_ref10      0.014120  0.06581   0.08549  0.1652 8.688e-01
## CBT1:TIME_ref12     -0.156900  0.06587   0.06913 -2.2700 2.320e-02
## MI1:CBT1:TIME_ref10 -0.012160  0.09340   0.11150 -0.1090 9.132e-01
## MI1:CBT1:TIME_ref12  0.166200  0.09348   0.09164  1.8140 6.974e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_TIMEref2) #ref: MI=0, CBT=0, time=2
##                     Estimates Model SE Robust SE    wald         p
## (Intercept)          -0.24160  0.04813   0.05618 -4.3000 1.711e-05
## cAGE                  0.02696  0.01594   0.01444  1.8670 6.191e-02
## MI1                   0.06476  0.06906   0.06917  0.9361 3.492e-01
## CBT1                 -0.14740  0.06552   0.08690 -1.6960 8.980e-02
## TIME_ref20            0.14630  0.05901   0.05720  2.5580 1.054e-02
## TIME_ref21           -0.08307  0.04847   0.05293 -1.5690 1.166e-01
## MI1:CBT1              0.13120  0.09295   0.10510  1.2490 2.118e-01
## MI1:TIME_ref20       -0.02514  0.08462   0.06868 -0.3660 7.143e-01
## MI1:TIME_ref21        0.02645  0.06947   0.06539  0.4045 6.858e-01
## CBT1:TIME_ref20       0.17110  0.08028   0.08717  1.9620 4.973e-02
## CBT1:TIME_ref21       0.15690  0.06587   0.06913  2.2700 2.320e-02
## MI1:CBT1:TIME_ref20  -0.17840  0.11400   0.10560 -1.6890 9.119e-02
## MI1:CBT1:TIME_ref21  -0.16620  0.09348   0.09164 -1.8140 6.974e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_MIref1) #ref: MI=1, CBT=0, time=0
##                     Estimates Model SE Robust SE     wald         p
## (Intercept)         -0.055640  0.04936   0.01224 -4.54600 5.460e-06
## cAGE                 0.026960  0.01594   0.01444  1.86700 6.191e-02
## MI_ref10            -0.039620  0.06903   0.02776 -1.42700 1.536e-01
## CBT1                -0.023550  0.06579   0.02266 -1.03900 2.987e-01
## TIME1               -0.177800  0.04974   0.04673 -3.80400 1.425e-04
## TIME2               -0.121200  0.06065   0.03802 -3.18700 1.439e-03
## MI_ref10:CBT1        0.047180  0.09281   0.03847  1.22600 2.201e-01
## MI_ref10:TIME1      -0.051590  0.06945   0.08027 -0.64270 5.204e-01
## MI_ref10:TIME2      -0.025140  0.08462   0.06868 -0.36600 7.143e-01
## CBT1:TIME1          -0.001961  0.06628   0.07162 -0.02739 9.782e-01
## CBT1:TIME2           0.007311  0.08089   0.05958  0.12270 9.023e-01
## MI_ref10:CBT1:TIME1 -0.012160  0.09340   0.11150 -0.10900 9.132e-01
## MI_ref10:CBT1:TIME2 -0.178400  0.11400   0.10560 -1.68900 9.119e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_MIref1_TIMEref1) #ref: MI=1, CBT=0, time=1
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)              -0.233400  0.04965   0.04844 -4.81900 1.440e-06
## cAGE                      0.026960  0.01594   0.01444  1.86700 6.191e-02
## MI_ref10                 -0.091210  0.06949   0.08152 -1.11900 2.632e-01
## CBT1                     -0.025510  0.06610   0.07123 -0.35820 7.202e-01
## TIME_ref10                0.177800  0.04974   0.04673  3.80400 1.425e-04
## TIME_ref12                0.056620  0.04976   0.03837  1.47600 1.401e-01
## MI_ref10:CBT1             0.035020  0.09330   0.11180  0.31330 7.541e-01
## MI_ref10:TIME_ref10       0.051590  0.06945   0.08027  0.64270 5.204e-01
## MI_ref10:TIME_ref12       0.026450  0.06947   0.06539  0.40450 6.858e-01
## CBT1:TIME_ref10           0.001961  0.06628   0.07162  0.02739 9.782e-01
## CBT1:TIME_ref12           0.009272  0.06632   0.06015  0.15420 8.775e-01
## MI_ref10:CBT1:TIME_ref10  0.012160  0.09340   0.11150  0.10900 9.132e-01
## MI_ref10:CBT1:TIME_ref12 -0.166200  0.09348   0.09164 -1.81400 6.974e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_MIref1_TIMEref2) #ref: MI=1, CBT=0, time=2
##                          Estimates Model SE Robust SE    wald         p
## (Intercept)              -0.176800  0.04940   0.04008 -4.4110 1.031e-05
## cAGE                      0.026960  0.01594   0.01444  1.8670 6.191e-02
## MI_ref10                 -0.064760  0.06906   0.06917 -0.9361 3.492e-01
## CBT1                     -0.016240  0.06586   0.05858 -0.2773 7.816e-01
## TIME_ref20                0.121200  0.06065   0.03802  3.1870 1.439e-03
## TIME_ref21               -0.056620  0.04976   0.03837 -1.4760 1.401e-01
## MI_ref10:CBT1            -0.131200  0.09295   0.10510 -1.2490 2.118e-01
## MI_ref10:TIME_ref20       0.025140  0.08462   0.06868  0.3660 7.143e-01
## MI_ref10:TIME_ref21      -0.026450  0.06947   0.06539 -0.4045 6.858e-01
## CBT1:TIME_ref20          -0.007311  0.08089   0.05958 -0.1227 9.023e-01
## CBT1:TIME_ref21          -0.009272  0.06632   0.06015 -0.1542 8.775e-01
## MI_ref10:CBT1:TIME_ref20  0.178400  0.11400   0.10560  1.6890 9.119e-02
## MI_ref10:CBT1:TIME_ref21  0.166200  0.09348   0.09164  1.8140 6.974e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_CBTref1) #ref: MI=0, CBT=1, time=0
##                     Estimates Model SE Robust SE    wald         p
## (Intercept)          -0.07164  0.04421   0.01806 -3.9680 7.260e-05
## cAGE                  0.02696  0.01594   0.01444  1.8670 6.191e-02
## MI1                  -0.00756  0.06202   0.02627 -0.2878 7.735e-01
## CBT_ref10            -0.02362  0.06539   0.03073 -0.7689 4.420e-01
## TIME1                -0.24350  0.04452   0.05530 -4.4030 1.069e-05
## TIME2                -0.31730  0.05443   0.06581 -4.8220 1.420e-06
## MI1:CBT_ref10         0.04718  0.09281   0.03847  1.2260 2.201e-01
## MI1:TIME1             0.06375  0.06246   0.07747  0.8229 4.106e-01
## MI1:TIME2             0.20350  0.07634   0.08021  2.5370 1.118e-02
## CBT_ref10:TIME1       0.01412  0.06581   0.08549  0.1652 8.688e-01
## CBT_ref10:TIME2       0.17110  0.08028   0.08717  1.9620 4.973e-02
## MI1:CBT_ref10:TIME1  -0.01216  0.09340   0.11150 -0.1090 9.132e-01
## MI1:CBT_ref10:TIME2  -0.17840  0.11400   0.10560 -1.6890 9.119e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_CBTref1_TIMEref1) #ref: MI=0, CBT=1, time=1
##                          Estimates Model SE Robust SE    wald         p
## (Intercept)              -0.315100  0.04438   0.05522 -5.7060 1.000e-08
## cAGE                      0.026960  0.01594   0.01444  1.8670 6.191e-02
## MI1                       0.056190  0.06224   0.07596  0.7397 4.595e-01
## CBT_ref10                -0.009504  0.06579   0.08600 -0.1105 9.120e-01
## TIME_ref10                0.243500  0.04452   0.05530  4.4030 1.069e-05
## TIME_ref12               -0.073860  0.04460   0.04445 -1.6620 9.658e-02
## MI1:CBT_ref10             0.035020  0.09330   0.11180  0.3133 7.541e-01
## MI1:TIME_ref10           -0.063750  0.06246   0.07747 -0.8229 4.106e-01
## MI1:TIME_ref12            0.139800  0.06255   0.06420  2.1770 2.949e-02
## CBT_ref10:TIME_ref10     -0.014120  0.06581   0.08549 -0.1652 8.688e-01
## CBT_ref10:TIME_ref12      0.156900  0.06587   0.06913  2.2700 2.320e-02
## MI1:CBT_ref10:TIME_ref10  0.012160  0.09340   0.11150  0.1090 9.132e-01
## MI1:CBT_ref10:TIME_ref12 -0.166200  0.09348   0.09164 -1.8140 6.974e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 4 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_CBTref1_TIMEref2) #ref: MI=0, CBT=1, time=2
##                          Estimates Model SE Robust SE   wald         p
## (Intercept)               -0.38900  0.04439   0.06573 -5.918 0.000e+00
## cAGE                       0.02696  0.01594   0.01444  1.867 6.191e-02
## MI1                        0.19590  0.06220   0.07833  2.502 1.236e-02
## CBT_ref10                  0.14740  0.06552   0.08690  1.696 8.980e-02
## TIME_ref20                 0.31730  0.05443   0.06581  4.822 1.420e-06
## TIME_ref21                 0.07386  0.04460   0.04445  1.662 9.658e-02
## MI1:CBT_ref10             -0.13120  0.09295   0.10510 -1.249 2.118e-01
## MI1:TIME_ref20            -0.20350  0.07634   0.08021 -2.537 1.118e-02
## MI1:TIME_ref21            -0.13980  0.06255   0.06420 -2.177 2.949e-02
## CBT_ref10:TIME_ref20      -0.17110  0.08028   0.08717 -1.962 4.973e-02
## CBT_ref10:TIME_ref21      -0.15690  0.06587   0.06913 -2.270 2.320e-02
## MI1:CBT_ref10:TIME_ref20   0.17840  0.11400   0.10560  1.689 9.119e-02
## MI1:CBT_ref10:TIME_ref21   0.16620  0.09348   0.09164  1.814 6.974e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_MIref1_CBTref1) #ref: MI=1, CBT=1, time=0
##                          Estimates Model SE Robust SE     wald         p
## (Intercept)              -0.079200  0.04352   0.01921 -4.12300 3.747e-05
## cAGE                      0.026960  0.01594   0.01444  1.86700 6.191e-02
## MI_ref10                  0.007560  0.06202   0.02627  0.28780 7.735e-01
## CBT_ref10                 0.023550  0.06579   0.02266  1.03900 2.987e-01
## TIME1                    -0.179700  0.04381   0.05428 -3.31100 9.280e-04
## TIME2                    -0.113800  0.05352   0.04589 -2.48100 1.311e-02
## MI_ref10:CBT_ref10       -0.047180  0.09281   0.03847 -1.22600 2.201e-01
## MI_ref10:TIME1           -0.063750  0.06246   0.07747 -0.82290 4.106e-01
## MI_ref10:TIME2           -0.203500  0.07634   0.08021 -2.53700 1.118e-02
## CBT_ref10:TIME1           0.001961  0.06628   0.07162  0.02739 9.782e-01
## CBT_ref10:TIME2          -0.007311  0.08089   0.05958 -0.12270 9.023e-01
## MI_ref10:CBT_ref10:TIME1  0.012160  0.09340   0.11150  0.10900 9.132e-01
## MI_ref10:CBT_ref10:TIME2  0.178400  0.11400   0.10560  1.68900 9.119e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_MIref1_CBTref1_TIMEref1) #ref: MI=1, CBT=1, time=1
##                               Estimates Model SE Robust SE     wald         p
## (Intercept)                   -0.258900  0.04366   0.05235 -4.94600 7.600e-07
## cAGE                           0.026960  0.01594   0.01444  1.86700 6.191e-02
## MI_ref10                      -0.056190  0.06224   0.07596 -0.73970 4.595e-01
## CBT_ref10                      0.025510  0.06610   0.07123  0.35820 7.202e-01
## TIME_ref10                     0.179700  0.04381   0.05428  3.31100 9.280e-04
## TIME_ref12                     0.065890  0.04385   0.04632  1.42200 1.549e-01
## MI_ref10:CBT_ref10            -0.035020  0.09330   0.11180 -0.31330 7.541e-01
## MI_ref10:TIME_ref10            0.063750  0.06246   0.07747  0.82290 4.106e-01
## MI_ref10:TIME_ref12           -0.139800  0.06255   0.06420 -2.17700 2.949e-02
## CBT_ref10:TIME_ref10          -0.001961  0.06628   0.07162 -0.02739 9.782e-01
## CBT_ref10:TIME_ref12          -0.009272  0.06632   0.06015 -0.15420 8.775e-01
## MI_ref10:CBT_ref10:TIME_ref10 -0.012160  0.09340   0.11150 -0.10900 9.132e-01
## MI_ref10:CBT_ref10:TIME_ref12  0.166200  0.09348   0.09164  1.81400 6.974e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729
summary(SMDK_interaction_MIref1_CBTref1_TIMEref2) #ref: MI=1, CBT=1, time=2
##                               Estimates Model SE Robust SE    wald         p
## (Intercept)                   -0.193000  0.04360   0.04285 -4.5050 6.620e-06
## cAGE                           0.026960  0.01594   0.01444  1.8670 6.191e-02
## MI_ref10                      -0.195900  0.06220   0.07833 -2.5020 1.236e-02
## CBT_ref10                      0.016240  0.06586   0.05858  0.2773 7.816e-01
## TIME_ref20                     0.113800  0.05352   0.04589  2.4810 1.311e-02
## TIME_ref21                    -0.065890  0.04385   0.04632 -1.4220 1.549e-01
## MI_ref10:CBT_ref10             0.131200  0.09295   0.10510  1.2490 2.118e-01
## MI_ref10:TIME_ref20            0.203500  0.07634   0.08021  2.5370 1.118e-02
## MI_ref10:TIME_ref21            0.139800  0.06255   0.06420  2.1770 2.949e-02
## CBT_ref10:TIME_ref20           0.007311  0.08089   0.05958  0.1227 9.023e-01
## CBT_ref10:TIME_ref21           0.009272  0.06632   0.06015  0.1542 8.775e-01
## MI_ref10:CBT_ref10:TIME_ref20 -0.178400  0.11400   0.10560 -1.6890 9.119e-02
## MI_ref10:CBT_ref10:TIME_ref21 -0.166200  0.09348   0.09164 -1.8140 6.974e-02
## 
##  Estimated Correlation Parameter:  0.4949 
##  Correlation Structure:  ar1 
##  Est. Scale Parameter:  0.1244 
## 
##  Number of GEE iterations: 3 
##  Number of Clusters:  243    Maximum Cluster Size:  3 
##  Number of observations with nonzero weight:  729