Two-Way Interaction plots from CBT-TIME Model

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

Model summary for CBT-Time 2-way interaction:

summary(PCSMKD_CBT_interactions)

            Estimates Model SE Robust SE    wald       p
(Intercept)  4.329000  0.04344  0.047380 91.3600 0.00000
cAGE         0.015270  0.01727  0.022710  0.6723 0.50140
cBLCGSMD     0.006982  0.00206  0.002979  2.3440 0.01910
MI1          0.089610  0.03805  0.046750  1.9170 0.05526
CBT1        -0.023990  0.05322  0.052940 -0.4532 0.65040
TIME1        0.063770  0.05641  0.034260  1.8610 0.06270
CBT1:TIME1  -0.096820  0.07631  0.049260 -1.9660 0.04935

 Estimated Correlation Parameter:  0 
 Correlation Structure:  independence 
 Est. Scale Parameter:  0.1836 

 Number of GEE iterations: 2 
 Number of Clusters:  273    Maximum Cluster Size:  2 
 Number of observations with nonzero weight:  516 

Let’s start by extracting all the coefficients

beta_coefs <- (summary(PCSMKD_CBT_interactions))$beta
beta_coefs <- lapply(beta_coefs, function(x){round(exp(x), 5)})
names(beta_coefs) <- summary(PCSMKD_CBT_interactions)$coefnames

beta_coefs

$`(Intercept)`
[1] 75.85613

$cAGE
[1] 1.01538

$cBLCGSMD
[1] 1.00701

$MI1
[1] 1.09375

$CBT1
[1] 0.97629

$TIME1
[1] 1.06584

$`CBT1:TIME1`
[1] 0.90772

Lay out the two-panel interaction plot:

• Panel 1: CBT = 0 (i.e., NicA)

  For group receiving NicA:

  • Plot p1 = % smoking days coefficient reported at time 0 for persons receiving RT (MI=0)

  • Plot p2 = % smoking days coefficient reported at time 1 for persons receiving RT (MI=0)

  • Connect (p1, p2) by a straight line

  • Plot p3 = % smoking days coefficient reported at time 0 for persons receiving MI

  • Plot p4 = % smoking days coefficient reported at time 1 for persons receiving MI

  • Connect (p3, p4) by a straight line

• Panel 2: CBT = 1

  For group receiving CBT:

  • Define points p1, p2, p3, p4 and connect them as described above.

beta_coefs["(Intercept)"]

$`(Intercept)`
[1] 75.85613

#Panel 1: CBT=0
panel1_p1 <- as.numeric(beta_coefs["(Intercept)"])
panel1_p2 <- as.numeric(beta_coefs["(Intercept)"]) + as.numeric(beta_coefs["TIME1"]) 
panel1_p3 <- as.numeric(beta_coefs["(Intercept)"]) + as.numeric(beta_coefs["MI1"]) 
panel1_p4 <- as.numeric(beta_coefs["(Intercept)"]) + as.numeric(beta_coefs["MI1"]) + as.numeric(beta_coefs["TIME1"])
  
#Panel 2: CBT=1
panel2_p1 <- panel1_p1 + as.numeric(beta_coefs["CBT1"])
panel2_p2 <- panel1_p2 + as.numeric(beta_coefs["CBT1"])
panel2_p3 <- panel1_p3 + as.numeric(beta_coefs["CBT1"])
panel2_p4 <- panel1_p4 + as.numeric(beta_coefs["CBT1"])

panel1_slope_l1 <- (panel1_p2-panel1_p1)/(1-0)
panel1_slope_l2 <- (panel1_p4-panel1_p3)/(1-0)

panel1_int_l1 = panel1_p1 - panel1_slope_l1*0
panel1_int_l2 = panel1_p3 - panel1_slope_l2*0
  
#Visualize Panel 1: CBT=0
plot(x=0, y=panel1_p1, type="p", xlab="Time", ylab = "Model Coefficients",
    xaxt="n", yaxt="n", xlim = c(-0.2, 1.2), ylim = c(75, 80), main="CBT=0")
axis(1, at = c(0, 1), labels=c("0", "1"), tick=TRUE)

points(x=1, y=panel1_p2)
points(x=0, y=panel1_p3)

points(x=1, y=panel1_p4)
abline(a=panel1_int_l1, b=panel1_slope_l1)

abline(a=panel1_int_l2, b=panel1_slope_l2)

panel2_slope_l1 <- (panel2_p2-panel2_p1)/(1-0)
panel2_slope_l2 <- (panel2_p4-panel2_p3)/(1-0)

panel2_int_l1 = panel2_p1 - panel2_slope_l1*0
panel2_int_l2 = panel2_p3 - panel2_slope_l2*0

#Visualize Panel 2: CBT=1
plot(x=0, y=panel2_p1, type="p", xlab="Time", ylab = "Model Coefficients",
    xaxt="n", yaxt="n", xlim = c(-0.2, 1.2), ylim = c(75, 80), main="CBT=1")

axis(1, at = c(0, 1), labels=c("0", "1"), tick=TRUE)

points(x=1, y=panel2_p2)
points(x=0, y=panel2_p3)

points(x=1, y=panel2_p4)

abline(a=panel2_int_l1, b=panel2_slope_l1)

abline(a=panel2_int_l2, b=panel2_slope_l2)