Setup, Data, Packages
if (Sys.info()[["user"]] == "jwile") {
loc.data <- "g:/My Drive/CBTI+Light project/Recruitment and data collection/"
} else {
loc.data <- "/Volumes/GoogleDrive/My Drive/CBTI+Light project/Recruitment and data collection/"
}
## graphics packages
library(ggplot2)
library(ggpubr)
library(scales)
## tools and miscellaneous packages
library(JWileymisc)
## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car
## data related packages
library(reshape2)
library(data.table)
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:reshape2':
##
## dcast, melt
## modelling packages
library(MplusAutomation)
## Version: 0.8
## We work hard to write this free software. Please help us get credit by citing:
##
## Hallquist, M. N. & Wiley, J. F. (2018). MplusAutomation: An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus. Structural Equation Modeling, 25, 621-638. doi: 10.1080/10705511.2017.1402334.
##
## -- see citation("MplusAutomation").
##
## Attaching package: 'MplusAutomation'
## The following object is masked from 'package:JWileymisc':
##
## cd
theme_set(theme_pubr())
Sys.setenv(TZ = "Australia/Melbourne")
SEback <- function(x) 100 - ((100 - x)^2)
d.all <- readRDS(file.path(loc.data, "Data", "SleepWell_Surveys.RDS"))
## create stratification factors
d.all[, dcond := as.integer(condition == "CBT+")]
d.all[, dstrat3 := as.integer(baseline_isi_high == "ISI >= 8" & baseline_stage_high == "Stage <= 2")]
d.all[, dstrat4 := as.integer(baseline_isi_high == "ISI >= 8" & baseline_stage_high == "Stage >= 3")]
## select only initially higher ISI at screening
## note that this means we cannot use dstrat1 or dstrat2
## as those only occur in women with low initial ISI levels
d.all <- d.all[baseline_isi_high == "ISI >= 8"]
## sleep data
d.sleep <- readRDS(file.path(loc.data, "Data", "SleepWell_Sleep.RDS"))
## merge survey and sleep data
## only take women in the second dataset, those with high screening ISI
d.all2 <- merge(
d.sleep,
d.all[Survey == "Baseline"],
by = "ID",
all.y=TRUE)
## weighted.mean(c(1, 4, 2, 5), w = c(.3, .2, .25, .25))
## check that weight the differences relative to first by weights matches
## 1 + 3 * .2 + 1 * .25 + 4 * .25
mcsetup <- "
NEW(
adjref
mtau1 mtau2 mtau3 mtau4
mcbt1 mcbt2 mcbt3 mcbt4
mdif1 mdif2 mdif3 mdif4
bd1 bd2 bd3 bd4
s1cbt s2cbt);
adjref = rm1 + i3 * %0.4f;
mtau1 = adjref + rm2 * 0 + rm3 * 0;
mtau2 = adjref + rm2 * 0.5 + rm3 * 0;
mtau3 = adjref + rm2 * 1 + rm3 * 0;
mtau4 = adjref + rm2 * 1 + rm3 * 1;
mcbt1 = (adjref) + (rm2 + s4) * 0 + (rm3 + s24) * 0;
mcbt2 = (adjref) + (rm2 + s4) * 0.5 + (rm3 + s24) * 0;
mcbt3 = (adjref) + (rm2 + s4) * 1 + (rm3 + s24) * 0;
mcbt4 = (adjref) + (rm2 + s4) * 1 + (rm3 + s24) * 1;
mdif1 = mcbt1 - mtau1;
mdif2 = mcbt2 - mtau2;
mdif3 = mcbt3 - mtau3;
mdif4 = mcbt4 - mtau4;
bd1 = mdif1 / sqrt(rv1 + resvar);
bd2 = mdif2 / sqrt(rv1 + resvar);
bd3 = mdif3 / sqrt(rv1 + resvar);
bd4 = mdif4 / sqrt(rv1 + resvar);
s1cbt = rm2 + s4;
s2cbt = rm3 + s24;
rv1 > .001;
rv2 > .001;
rv3 > .001;
"
Insomnia Symptons
isidat <- copy(reshape(d.all[, .(
ID, dcond, dstrat4, ISI_Total, Survey)],
v.names = "ISI_Total",
timevar = "Survey",
idvar = c("ID", "dcond", "dstrat4"),
direction = "wide"))
setnames(isidat, names(isidat),
c("id", "dcond", "dstrat4",
"ISI0", "ISI1", "ISI2", "ISI3", "ISI4"))
isidat <- isidat[!is.na(ISI1) | !is.na(ISI2) | !is.na(ISI3) | !is.na(ISI4)]
isifit <- mplusModeler(mplusObject(
TITLE = "Primary ISI Analysis;",
ANALYSIS = "ESTIMATOR = ML;",
MODEL = "
[ISI1 - ISI4@0]; ! fix residual intercepts to 0
ISI1* ISI2* ISI3* ISI4* (resvar); ! homogenous residual variance
int BY ISI1@1 ISI2@1 ISI3@1 ISI4@1;
slope BY ISI1@0 ISI2@0.5 ISI3@1 ISI4@1;
slope2 BY ISI1@0 ISI2@0 ISI3@0 ISI4@1;
int ON dstrat4 dcond@0 (i3 - i4);
slope ON dcond (s4);
slope2 ON dcond (s24);
[int* slope* slope2*] (rm1 - rm3);
int* slope* slope2*0 (rv1-rv3);
int WITH slope*0;
slope2 WITH int*0 slope*0;
",
OUTPUT = "STDYX; CINTERVAL;",
MODELCONSTRAINT = sprintf(
mcsetup,
mean(isidat$dstrat4)),
usevariables = c("ISI1", "ISI2", "ISI3", "ISI4",
"dstrat4", "dcond"),
rdata = as.data.frame(isidat)),
dataout = "03_HighISI_isi.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 03_HighISI_isi.inp
## Wrote data to: 03_HighISI_isi.dat
## model information, coefficients, and effect sizes
summary(isifit)
## Primary ISI Analysis
##
## Estimated using ML
## Number of obs: 67, number of (free) parameters: 13
##
## Model: Chi2(df = 9) = 17.337, p = 0.0437
## Baseline model: Chi2(df = 14) = 83.518, p = 0
##
## Fit Indices:
##
## CFI = 0.88, TLI = 0.813, SRMR = 0.109
## RMSEA = 0.118, 90% CI [0.019, 0.2], p < .05 = 0.093
## AIC = 1223.323, BIC = 1251.984
cbind(
coef(isifit, param = c("reg", "exp", "new")),
confint(isifit, param = c("reg", "exp", "new"))[,-1])
## Label est se pval LowerCI UpperCI
## 13 INT<-DSTRAT4 -1.172 1.116 0.294 -3.358 1.015
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND -3.173 1.417 0.025 -5.950 -0.395
## 16 SLOPE2<-DCOND 4.030 2.010 0.045 0.092 7.969
## 20 ISI1<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 21 ISI2<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 22 ISI3<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 23 ISI4<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 24 INT<-Intercepts 15.411 0.805 0.000 13.834 16.988
## 25 SLOPE<-Intercepts -2.944 1.062 0.006 -5.026 -0.862
## 26 SLOPE2<-Intercepts -3.493 1.527 0.022 -6.487 -0.500
## 34 ADJREF 14.834 0.587 0.000 13.684 15.984
## 35 MTAU1 14.834 0.587 0.000 13.684 15.984
## 36 MTAU2 13.362 0.695 0.000 12.000 14.724
## 37 MTAU3 11.890 1.089 0.000 9.756 14.024
## 38 MTAU4 8.397 1.346 0.000 5.758 11.035
## 39 MCBT1 14.834 0.587 0.000 13.684 15.984
## 40 MCBT2 11.776 0.669 0.000 10.465 13.086
## 41 MCBT3 8.717 1.025 0.000 6.708 10.727
## 42 MCBT4 9.254 1.127 0.000 7.045 11.464
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 -1.586 0.709 0.025 -2.975 -0.197
## 45 MDIF3 -3.173 1.417 0.025 -5.950 -0.395
## 46 MDIF4 0.858 1.688 0.611 -2.450 4.166
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 -0.323 0.147 0.028 -0.611 -0.035
## 49 BD3 -0.646 0.294 0.028 -1.222 -0.070
## 50 BD4 0.175 0.344 0.612 -0.500 0.849
## 51 S1CBT -6.117 1.002 0.000 -8.080 -4.154
## 52 S2CBT 0.537 1.313 0.683 -2.036 3.110
## create all the data for the plot, including means, effect sizes, Ns
## confidence intervals and significance
isi.plotdat <- list()
isi.plotdat <- within(isi.plotdat, {
est <- as.data.table(cbind(coef(isifit), confint(isifit)[, -1]))[35:42, ]
est[, Condition := rep(c("TAU+", "CBT-I+Light"), each = 4)]
est[, Timepoint := c(0:3, 0:3)]
freq <- xtabs(~ condition + Survey, data = d.all[!is.na(ISI_Total)])
ES <- c(sprintf("ES = %s",
gsub("^0", "",
format(round(abs(subset(coef(isifit, param = c("new")),
grepl("BD", Label))$est), 2),
digits = 2, nsmall = 2))))[-1]
labels <- c(
sprintf("Baseline\nCBT-I+Light N: %d \nTAU+ N: %d ",
freq["CBT+", "Baseline"],
freq["Relaxation", "Baseline"]),
sprintf("Midpoint\n%d\n%d",
freq["CBT+", "Midpoint"],
freq["Relaxation", "Midpoint"]),
sprintf("Post\n%d\n%d",
freq["CBT+", "Post"],
freq["Relaxation", "Post"]),
sprintf("Follow-up\n%d\n%d",
freq["CBT+", "Follow-Up"],
freq["Relaxation", "Follow-Up"]))
sig <- symnum(subset(coef(isifit), grepl("MDIF", Label))$pval,
legend = FALSE, na = "",
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "n.s."))[-1]
yval.sig <- with(est, max(UpperCI) + (abs(diff(range(c(LowerCI, UpperCI)))) * .02))
yval.es <- with(est, min(LowerCI) - (abs(diff(range(c(LowerCI, UpperCI)))) * .03))
})
## make the plot
p.isi <- with(isi.plotdat,
ggplot(est, aes(Timepoint + ifelse(Condition=="CBT-I+Light", .03, -.03),
est, ymin = LowerCI, ymax = UpperCI,
colour = Condition, shape = Condition)) +
geom_hline(yintercept = 8, linetype = 3, colour = "grey50") +
geom_line(size = .1, alpha = .5) + geom_pointrange(size=1) +
scale_x_continuous("", breaks = 0:3, labels = labels) +
scale_y_continuous("Insomnia Symptoms (Adjusted)",
breaks = seq(5, 18, 1)) +
annotate("text", x = 1, y = yval.sig, label = sig[1]) +
annotate("text", x = 2, y = yval.sig, label = sig[2]) +
annotate("text", x = 3, y = yval.sig, label = sig[3]) +
annotate("text", x = 1, y = yval.es, label = ES[1]) +
annotate("text", x = 2, y = yval.es, label = ES[2]) +
annotate("text", x = 3, y = yval.es, label = ES[3]))
p.isi <- set_palette(p.isi, "jco")
print(p.isi)
Sleep Related Impairment
sridat <- copy(reshape(d.all[, .(
ID, dcond, dstrat4, SRI_TScore, Survey)],
v.names = "SRI_TScore",
timevar = "Survey",
idvar = c("ID", "dcond", "dstrat4"),
direction = "wide"))
setnames(sridat, names(sridat),
c("id", "dcond", "dstrat4",
"SRI0", "SRI1", "SRI2", "SRI3", "SRI4"))
sridat <- sridat[!is.na(SRI1) | !is.na(SRI2) | !is.na(SRI3) | !is.na(SRI4)]
srifit <- mplusModeler(mplusObject(
TITLE = "Primary SRI Analysis;",
ANALYSIS = "ESTIMATOR = ML;",
MODEL = "
[SRI1 - SRI4@0]; ! fix residual intercepts to 0
SRI1* SRI2* SRI3* SRI4* (resvar); ! homogenous residual variance
int BY SRI1@1 SRI2@1 SRI3@1 SRI4@1;
slope BY SRI1@0 SRI2@0.5 SRI3@1 SRI4@1;
slope2 BY SRI1@0 SRI2@0 SRI3@0 SRI4@1;
int ON dstrat4 dcond@0 (i3 - i4);
slope ON dcond (s4);
slope2 ON dcond (s24);
[int* slope* slope2*] (rm1 - rm3);
int* slope* slope2*0 (rv1-rv3);
int WITH slope*0;
slope2 WITH int*0 slope*0;
",
OUTPUT = "STDYX; CINTERVAL;",
MODELCONSTRAINT = sprintf(
mcsetup,
mean(sridat$dstrat4)),
usevariables = c("SRI1", "SRI2", "SRI3", "SRI4",
"dstrat4", "dcond"),
rdata = as.data.frame(sridat)),
dataout = "03_HighISI_sri.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 03_HighISI_sri.inp
## Wrote data to: 03_HighISI_sri.dat
## model information, coefficients, and effect sizes
summary(srifit)
## Primary SRI Analysis
##
## Estimated using ML
## Number of obs: 67, number of (free) parameters: 13
##
## Model: Chi2(df = 9) = 10.574, p = 0.3061
## Baseline model: Chi2(df = 14) = 78.059, p = 0
##
## Fit Indices:
##
## CFI = 0.975, TLI = 0.962, SRMR = 0.07
## RMSEA = 0.051, 90% CI [0, 0.152], p < .05 = 0.433
## AIC = 1380.965, BIC = 1409.626
cbind(
coef(srifit, param = c("reg", "exp", "new")),
confint(srifit, param = c("reg", "exp", "new"))[,-1])
## Label est se pval LowerCI UpperCI
## 13 INT<-DSTRAT4 0.033 1.673 0.984 -3.246 3.311
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND -3.602 2.012 0.073 -7.545 0.340
## 16 SLOPE2<-DCOND 0.016 2.624 0.995 -5.128 5.159
## 20 SRI1<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 21 SRI2<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 22 SRI3<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 23 SRI4<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 24 INT<-Intercepts 58.429 1.209 0.000 56.061 60.798
## 25 SLOPE<-Intercepts -2.380 1.522 0.118 -5.362 0.602
## 26 SLOPE2<-Intercepts -2.694 1.980 0.174 -6.575 1.187
## 34 ADJREF 58.445 0.875 0.000 56.730 60.161
## 35 MTAU1 58.445 0.875 0.000 56.730 60.161
## 36 MTAU2 57.255 1.021 0.000 55.254 59.257
## 37 MTAU3 56.065 1.574 0.000 52.981 59.150
## 38 MTAU4 53.371 1.940 0.000 49.570 57.173
## 39 MCBT1 58.445 0.875 0.000 56.730 60.161
## 40 MCBT2 55.454 0.975 0.000 53.544 57.365
## 41 MCBT3 52.463 1.457 0.000 49.608 55.319
## 42 MCBT4 49.785 1.702 0.000 46.449 53.121
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 -1.801 1.006 0.073 -3.772 0.170
## 45 MDIF3 -3.602 2.012 0.073 -7.545 0.340
## 46 MDIF4 -3.586 2.416 0.138 -8.321 1.148
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 -0.244 0.138 0.077 -0.515 0.027
## 49 BD3 -0.488 0.276 0.077 -1.029 0.053
## 50 BD4 -0.486 0.330 0.141 -1.133 0.161
## 51 S1CBT -5.982 1.406 0.000 -8.738 -3.226
## 52 S2CBT -2.678 1.727 0.121 -6.063 0.706
## create all the data for the plot, including means, effect sizes, Ns
## confidence intervals and significance
sri.plotdat <- list()
sri.plotdat <- within(sri.plotdat, {
est <- as.data.table(cbind(coef(srifit), confint(srifit)[, -1]))[35:42, ]
est[, Condition := rep(c("TAU+", "CBT-I+Light"), each = 4)]
est[, Timepoint := c(0:3, 0:3)]
freq <- xtabs(~ condition + Survey, data = d.all[!is.na(SRI_Total)])
ES <- c(sprintf("ES = %s",
gsub("^0", "",
format(round(abs(subset(coef(srifit, param = c("new")),
grepl("BD", Label))$est), 2),
digits = 2, nsmall = 2))))[-1]
labels <- c(
sprintf("Baseline\nCBT-I+Light N: %d \nTAU+ N: %d ",
freq["CBT+", "Baseline"],
freq["Relaxation", "Baseline"]),
sprintf("Midpoint\n%d\n%d",
freq["CBT+", "Midpoint"],
freq["Relaxation", "Midpoint"]),
sprintf("Post\n%d\n%d",
freq["CBT+", "Post"],
freq["Relaxation", "Post"]),
sprintf("Follow-up\n%d\n%d",
freq["CBT+", "Follow-Up"],
freq["Relaxation", "Follow-Up"]))
sig <- symnum(subset(coef(srifit), grepl("MDIF", Label))$pval,
legend = FALSE, na = "",
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "n.s."))[-1]
yval.sig <- with(est, max(UpperCI) + (abs(diff(range(c(LowerCI, UpperCI)))) * .02))
yval.es <- with(est, min(LowerCI) - (abs(diff(range(c(LowerCI, UpperCI)))) * .03))
})
## make the plot
p.sri <- with(sri.plotdat,
ggplot(est, aes(Timepoint + ifelse(Condition=="CBT-I+Light", .03, -.03),
est, ymin = LowerCI, ymax = UpperCI,
colour = Condition, shape = Condition)) +
geom_hline(yintercept = 50, linetype = 3, colour = "grey50") +
geom_line(size = .1, alpha = .5) + geom_pointrange(size=1) +
scale_x_continuous("", breaks = 0:3, labels = labels) +
scale_y_continuous("Sleep Related Impairment (Adjusted)",
breaks = seq(40, 60, 2)) +
annotate("text", x = 1, y = yval.sig, label = sig[1]) +
annotate("text", x = 2, y = yval.sig, label = sig[2]) +
annotate("text", x = 3, y = yval.sig, label = sig[3]) +
annotate("text", x = 1, y = yval.es, label = ES[1]) +
annotate("text", x = 2, y = yval.es, label = ES[2]) +
annotate("text", x = 3, y = yval.es, label = ES[3]))
p.sri <- set_palette(p.sri, "jco")
print(p.sri)
Sleep Disturbance Symptoms
sddat <- copy(reshape(d.all[, .(
ID, dcond, dstrat4, SD_TScore, Survey)],
v.names = "SD_TScore",
timevar = "Survey",
idvar = c("ID", "dcond", "dstrat4"),
direction = "wide"))
setnames(sddat, names(sddat),
c("id", "dcond", "dstrat4",
"SD0", "SD1", "SD2", "SD3", "SD4"))
sddat <- sddat[!is.na(SD1) | !is.na(SD2) | !is.na(SD3) | !is.na(SD4)]
sdfit <- mplusModeler(mplusObject(
TITLE = "Primary SD Analysis;",
ANALYSIS = "ESTIMATOR = ML;",
MODEL = "
[SD1 - SD4@0]; ! fix residual intercepts to 0
SD1* SD2* SD3* SD4* (resvar); ! homogenous residual variance
int BY SD1@1 SD2@1 SD3@1 SD4@1;
slope BY SD1@0 SD2@0.5 SD3@1 SD4@1;
slope2 BY SD1@0 SD2@0 SD3@0 SD4@1;
int ON dstrat4 dcond@0 (i3 - i4);
slope ON dcond (s4);
slope2 ON dcond (s24);
[int* slope* slope2*] (rm1 - rm3);
int* slope* slope2*0 (rv1-rv3);
int WITH slope*0;
slope2 WITH int*0 slope@0;
",
OUTPUT = "STDYX; CINTERVAL;",
MODELCONSTRAINT = sprintf(
mcsetup,
mean(sddat$dstrat4)),
usevariables = c("SD1", "SD2", "SD3", "SD4",
"dstrat4", "dcond"),
rdata = as.data.frame(sddat)),
dataout = "03_HighISI_sd.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 03_HighISI_sd.inp
## Wrote data to: 03_HighISI_sd.dat
## model information, coefficients, and effect sizes
summary(sdfit)
## Primary SD Analysis
##
## Estimated using ML
## Number of obs: 67, number of (free) parameters: 12
##
## Model: Chi2(df = 10) = 20.642, p = 0.0237
## Baseline model: Chi2(df = 14) = 71.29, p = 0
##
## Fit Indices:
##
## CFI = 0.814, TLI = 0.74, SRMR = 0.283
## RMSEA = 0.126, 90% CI [0.044, 0.203], p < .05 = 0.059
## AIC = 1397.562, BIC = 1424.018
cbind(
coef(sdfit, param = c("reg", "exp", "new")),
confint(sdfit, param = c("reg", "exp", "new"))[,-1])
## Label est se pval LowerCI UpperCI
## 13 INT<-DSTRAT4 -1.101 1.756 0.531 -4.542 2.341
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND -2.034 1.943 0.295 -5.843 1.775
## 16 SLOPE2<-DCOND 4.871 2.751 0.077 -0.521 10.264
## 20 SD1<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 21 SD2<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 22 SD3<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 23 SD4<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 24 INT<-Intercepts 57.985 1.290 0.000 55.457 60.514
## 25 SLOPE<-Intercepts -4.663 1.472 0.002 -7.547 -1.779
## 26 SLOPE2<-Intercepts -3.973 2.077 0.056 -8.045 0.099
## 34 ADJREF 57.443 0.951 0.000 55.578 59.308
## 35 MTAU1 57.443 0.951 0.000 55.578 59.308
## 36 MTAU2 55.112 1.041 0.000 53.071 57.153
## 37 MTAU3 52.781 1.531 0.000 49.779 55.782
## 38 MTAU4 48.808 1.940 0.000 45.005 52.611
## 39 MCBT1 57.443 0.951 0.000 55.578 59.308
## 40 MCBT2 54.095 1.005 0.000 52.126 56.064
## 41 MCBT3 50.747 1.439 0.000 47.926 53.568
## 42 MCBT4 51.645 1.677 0.000 48.357 54.933
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 -1.017 0.972 0.295 -2.921 0.888
## 45 MDIF3 -2.034 1.943 0.295 -5.843 1.775
## 46 MDIF4 2.838 2.494 0.255 -2.051 7.726
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 -0.127 0.122 0.298 -0.367 0.112
## 49 BD3 -0.255 0.245 0.298 -0.734 0.225
## 50 BD4 0.355 0.314 0.257 -0.259 0.970
## 51 S1CBT -6.696 1.384 0.000 -9.409 -3.984
## 52 S2CBT 0.898 1.812 0.620 -2.653 4.449
## create all the data for the plot, including means, effect sizes, Ns
## confidence intervals and significance
sd.plotdat <- list()
sd.plotdat <- within(sd.plotdat, {
est <- as.data.table(cbind(coef(sdfit), confint(sdfit)[, -1]))[35:42, ]
est[, Condition := rep(c("TAU+", "CBT-I+Light"), each = 4)]
est[, Timepoint := c(0:3, 0:3)]
freq <- xtabs(~ condition + Survey, data = d.all[!is.na(SD_Total)])
ES <- c(sprintf("ES = %s",
gsub("^0", "",
format(round(abs(subset(coef(sdfit, param = c("new")),
grepl("BD", Label))$est), 2),
digits = 2, nsmall = 2))))[-1]
labels <- c(
sprintf("Baseline\nCBT-I+Light N: %d \nTAU+ N: %d ",
freq["CBT+", "Baseline"],
freq["Relaxation", "Baseline"]),
sprintf("Midpoint\n%d\n%d",
freq["CBT+", "Midpoint"],
freq["Relaxation", "Midpoint"]),
sprintf("Post\n%d\n%d",
freq["CBT+", "Post"],
freq["Relaxation", "Post"]),
sprintf("Follow-up\n%d\n%d",
freq["CBT+", "Follow-Up"],
freq["Relaxation", "Follow-Up"]))
sig <- symnum(subset(coef(sdfit), grepl("MDIF", Label))$pval,
legend = FALSE, na = "",
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "n.s."))[-1]
yval.sig <- with(est, max(UpperCI) + (abs(diff(range(c(LowerCI, UpperCI)))) * .02))
yval.es <- with(est, min(LowerCI) - (abs(diff(range(c(LowerCI, UpperCI)))) * .03))
})
## make the plot
p.sd <- with(sd.plotdat,
ggplot(est, aes(Timepoint + ifelse(Condition=="CBT-I+Light", .03, -.03),
est, ymin = LowerCI, ymax = UpperCI,
colour = Condition, shape = Condition)) +
geom_hline(yintercept = 50, linetype = 3, colour = "grey50") +
geom_line(size = .1, alpha = .5) + geom_pointrange(size=1) +
scale_x_continuous("", breaks = 0:3, labels = labels) +
scale_y_continuous("Sleep Disturbance (Adjusted)",
breaks = seq(40, 58, 2)) +
annotate("text", x = 1, y = yval.sig, label = sig[1]) +
annotate("text", x = 2, y = yval.sig, label = sig[2]) +
annotate("text", x = 3, y = yval.sig, label = sig[3]) +
annotate("text", x = 1, y = yval.es, label = ES[1]) +
annotate("text", x = 2, y = yval.es, label = ES[2]) +
annotate("text", x = 3, y = yval.es, label = ES[3]))
p.sd <- set_palette(p.sd, "jco")
print(p.sd)
Fatigue Symptoms
fadat <- copy(reshape(d.all[, .(
ID, dcond, dstrat4, FA_TScore, Survey)],
v.names = "FA_TScore",
timevar = "Survey",
idvar = c("ID", "dcond", "dstrat4"),
direction = "wide"))
setnames(fadat, names(fadat),
c("id", "dcond", "dstrat4",
"FA0", "FA1", "FA2", "FA3", "FA4"))
fadat <- fadat[!is.na(FA1) | !is.na(FA2) | !is.na(FA3) | !is.na(FA4)]
fafit <- mplusModeler(mplusObject(
TITLE = "Primary FA Analysis;",
ANALYSIS = "ESTIMATOR = ML;",
MODEL = "
[FA1 - FA4@0]; ! fix residual intercepts to 0
FA1* FA2* FA3* FA4* (resvar); ! homogenous residual variance
int BY FA1@1 FA2@1 FA3@1 FA4@1;
slope BY FA1@0 FA2@0.5 FA3@1 FA4@1;
slope2 BY FA1@0 FA2@0 FA3@0 FA4@1;
int ON dstrat4 dcond@0 (i3 - i4);
slope ON dcond (s4);
slope2 ON dcond (s24);
[int* slope* slope2*] (rm1 - rm3);
int* slope* slope2*0 (rv1-rv3);
int WITH slope*0;
slope2 WITH int*0 slope*0;
",
OUTPUT = "STDYX; CINTERVAL;",
MODELCONSTRAINT = sprintf(
mcsetup,
mean(fadat$dstrat4)),
usevariables = c("FA1", "FA2", "FA3", "FA4",
"dstrat4", "dcond"),
rdata = as.data.frame(fadat)),
dataout = "03_HighISI_fa.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 03_HighISI_fa.inp
## Wrote data to: 03_HighISI_fa.dat
## model information, coefficients, and effect sizes
summary(fafit)
## Primary FA Analysis
##
## Estimated using ML
## Number of obs: 67, number of (free) parameters: 13
##
## Model: Chi2(df = 9) = 12.438, p = 0.1897
## Baseline model: Chi2(df = 14) = 90.029, p = 0
##
## Fit Indices:
##
## CFI = 0.955, TLI = 0.93, SRMR = 0.169
## RMSEA = 0.076, 90% CI [0, 0.167], p < .05 = 0.3
## AIC = 1373.881, BIC = 1402.542
cbind(
coef(fafit, param = c("reg", "exp", "new")),
confint(fafit, param = c("reg", "exp", "new"))[,-1])
## Label est se pval LowerCI UpperCI
## 13 INT<-DSTRAT4 -0.753 1.545 0.626 -3.781 2.276
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND -5.403 1.946 0.005 -9.216 -1.590
## 16 SLOPE2<-DCOND -0.671 2.730 0.806 -6.022 4.680
## 20 FA1<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 21 FA2<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 22 FA3<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 23 FA4<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 24 INT<-Intercepts 60.876 1.143 0.000 58.636 63.117
## 25 SLOPE<-Intercepts -0.906 1.441 0.530 -3.730 1.918
## 26 SLOPE2<-Intercepts -1.451 2.133 0.496 -5.632 2.729
## 34 ADJREF 60.505 0.847 0.000 58.846 62.165
## 35 MTAU1 60.505 0.847 0.000 58.846 62.165
## 36 MTAU2 60.052 1.066 0.000 57.964 62.141
## 37 MTAU3 59.600 1.610 0.000 56.444 62.756
## 38 MTAU4 58.148 2.174 0.000 53.886 62.410
## 39 MCBT1 60.505 0.847 0.000 58.846 62.165
## 40 MCBT2 57.351 1.023 0.000 55.346 59.356
## 41 MCBT3 54.197 1.500 0.000 51.257 57.136
## 42 MCBT4 52.074 1.905 0.000 48.340 55.809
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 -2.702 0.973 0.005 -4.608 -0.795
## 45 MDIF3 -5.403 1.946 0.005 -9.216 -1.590
## 46 MDIF4 -6.074 2.882 0.035 -11.723 -0.424
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 -0.380 0.141 0.007 -0.656 -0.104
## 49 BD3 -0.759 0.282 0.007 -1.311 -0.207
## 50 BD4 -0.854 0.412 0.038 -1.661 -0.046
## 51 S1CBT -6.309 1.322 0.000 -8.900 -3.717
## 52 S2CBT -2.122 1.877 0.258 -5.801 1.557
## create all the data for the plot, including means, effect sizes, Ns
## confidence intervals and significance
fa.plotdat <- list()
fa.plotdat <- within(fa.plotdat, {
est <- as.data.table(cbind(coef(fafit), confint(fafit)[, -1]))[35:42, ]
est[, Condition := rep(c("TAU+", "CBT-I+Light"), each = 4)]
est[, Timepoint := c(0:3, 0:3)]
freq <- xtabs(~ condition + Survey, data = d.all[!is.na(FA_Total)])
ES <- c(sprintf("ES = %s",
gsub("^0", "",
format(round(abs(subset(coef(fafit, param = c("new")),
grepl("BD", Label))$est), 2),
digits = 2, nsmall = 2))))[-1]
labels <- c(
sprintf("Baseline\nCBT-I+Light N: %d \nTAU+ N: %d ",
freq["CBT+", "Baseline"],
freq["Relaxation", "Baseline"]),
sprintf("Midpoint\n%d\n%d",
freq["CBT+", "Midpoint"],
freq["Relaxation", "Midpoint"]),
sprintf("Post\n%d\n%d",
freq["CBT+", "Post"],
freq["Relaxation", "Post"]),
sprintf("Follow-up\n%d\n%d",
freq["CBT+", "Follow-Up"],
freq["Relaxation", "Follow-Up"]))
sig <- symnum(subset(coef(fafit), grepl("MDIF", Label))$pval,
legend = FALSE, na = "",
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "n.s."))[-1]
yval.sig <- with(est, max(UpperCI) + (abs(diff(range(c(LowerCI, UpperCI)))) * .02))
yval.es <- with(est, min(LowerCI) - (abs(diff(range(c(LowerCI, UpperCI)))) * .03))
})
## make the plot
p.fa <- with(fa.plotdat,
ggplot(est, aes(Timepoint + ifelse(Condition=="CBT-I+Light", .03, -.03),
est, ymin = LowerCI, ymax = UpperCI,
colour = Condition, shape = Condition)) +
geom_hline(yintercept = 50, linetype = 3, colour = "grey50") +
geom_line(size = .1, alpha = .5) + geom_pointrange(size=1) +
scale_x_continuous("", breaks = 0:3, labels = labels) +
scale_y_continuous("Fatigue Symptoms (Adjusted)",
breaks = seq(46, 62, 2)) +
annotate("text", x = 1, y = yval.sig, label = sig[1]) +
annotate("text", x = 2, y = yval.sig, label = sig[2]) +
annotate("text", x = 3, y = yval.sig, label = sig[3]) +
annotate("text", x = 1, y = yval.es, label = ES[1]) +
annotate("text", x = 2, y = yval.es, label = ES[2]) +
annotate("text", x = 3, y = yval.es, label = ES[3]))
p.fa <- set_palette(p.fa, "jco")
print(p.fa)
Depression Symptoms
depdat <- copy(reshape(d.all[, .(
ID, dcond, dstrat4, DEP_TScore, Survey)],
v.names = "DEP_TScore",
timevar = "Survey",
idvar = c("ID", "dcond", "dstrat4"),
direction = "wide"))
setnames(depdat, names(depdat),
c("id", "dcond", "dstrat4",
"DEP0", "DEP1", "DEP2", "DEP3", "DEP4"))
depdat <- depdat[!is.na(DEP1) | !is.na(DEP2) | !is.na(DEP3) | !is.na(DEP4)]
depfit <- mplusModeler(mplusObject(
TITLE = "Primary DEP Analysis;",
ANALYSIS = "ESTIMATOR = ML;",
MODEL = "
[DEP1 - DEP4@0]; ! fix residual intercepts to 0
DEP1* DEP2* DEP3* DEP4* (resvar); ! homogenous residual variance
int BY DEP1@1 DEP2@1 DEP3@1 DEP4@1;
slope BY DEP1@0 DEP2@0.5 DEP3@1 DEP4@1;
slope2 BY DEP1@0 DEP2@0 DEP3@0 DEP4@1;
int ON dstrat4 dcond@0 (i3 - i4);
slope ON dcond (s4);
slope2 ON dcond (s24);
[int* slope* slope2*] (rm1 - rm3);
int* slope* slope2*0 (rv1-rv3);
int WITH slope*0;
slope2 WITH int*0 slope*0;
",
OUTPUT = "STDYX; CINTERVAL;",
MODELCONSTRAINT = sprintf(
mcsetup,
mean(depdat$dstrat4)),
usevariables = c("DEP1", "DEP2", "DEP3", "DEP4",
"dstrat4", "dcond"),
rdata = as.data.frame(depdat)),
dataout = "03_HighISI_dep.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 03_HighISI_dep.inp
## Wrote data to: 03_HighISI_dep.dat
## model information, coefficients, and effect sizes
summary(depfit)
## Primary DEP Analysis
##
## Estimated using ML
## Number of obs: 67, number of (free) parameters: 13
##
## Model: Chi2(df = 9) = 10.469, p = 0.3138
## Baseline model: Chi2(df = 14) = 130.223, p = 0
##
## Fit Indices:
##
## CFI = 0.987, TLI = 0.98, SRMR = 0.089
## RMSEA = 0.049, 90% CI [0, 0.151], p < .05 = 0.442
## AIC = 1380.871, BIC = 1409.532
cbind(
coef(depfit, param = c("reg", "exp", "new")),
confint(depfit, param = c("reg", "exp", "new"))[,-1])
## Label est se pval LowerCI UpperCI
## 13 INT<-DSTRAT4 -2.761 2.143 0.198 -6.962 1.440
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND 2.107 1.891 0.265 -1.600 5.814
## 16 SLOPE2<-DCOND -2.437 2.261 0.281 -6.867 1.994
## 20 DEP1<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 21 DEP2<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 22 DEP3<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 23 DEP4<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 24 INT<-Intercepts 53.972 1.528 0.000 50.976 56.968
## 25 SLOPE<-Intercepts -2.659 1.419 0.061 -5.440 0.122
## 26 SLOPE2<-Intercepts -0.198 1.675 0.906 -3.482 3.086
## 34 ADJREF 52.612 1.102 0.000 50.453 54.772
## 35 MTAU1 52.612 1.102 0.000 50.453 54.772
## 36 MTAU2 51.283 1.174 0.000 48.982 53.584
## 37 MTAU3 49.953 1.596 0.000 46.825 53.082
## 38 MTAU4 49.755 1.844 0.000 46.142 53.369
## 39 MCBT1 52.612 1.102 0.000 50.453 54.772
## 40 MCBT2 52.336 1.141 0.000 50.100 54.572
## 41 MCBT3 52.060 1.504 0.000 49.112 55.008
## 42 MCBT4 49.425 1.682 0.000 46.128 52.722
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 1.053 0.946 0.265 -0.800 2.907
## 45 MDIF3 2.107 1.891 0.265 -1.600 5.814
## 46 MDIF4 -0.330 2.193 0.880 -4.628 3.968
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 0.115 0.104 0.267 -0.088 0.318
## 49 BD3 0.230 0.207 0.267 -0.176 0.637
## 50 BD4 -0.036 0.240 0.880 -0.506 0.434
## 51 S1CBT -0.552 1.322 0.676 -3.143 2.038
## 52 S2CBT -2.635 1.493 0.078 -5.561 0.291
## create all the data for the plot, including means, effect sizes, Ns
## confidence intervals and significance
dep.plotdat <- list()
dep.plotdat <- within(dep.plotdat, {
est <- as.data.table(cbind(coef(depfit), confint(depfit)[, -1]))[35:42, ]
est[, Condition := rep(c("TAU+", "CBT-I+Light"), each = 4)]
est[, Timepoint := c(0:3, 0:3)]
freq <- xtabs(~ condition + Survey, data = d.all[!is.na(DEP_Total)])
ES <- c(sprintf("ES = %s",
gsub("^0", "",
format(round(abs(subset(coef(depfit, param = c("new")),
grepl("BD", Label))$est), 2),
digits = 2, nsmall = 2))))[-1]
labels <- c(
sprintf("Baseline\nCBT-I+Light N: %d \nTAU+ N: %d ",
freq["CBT+", "Baseline"],
freq["Relaxation", "Baseline"]),
sprintf("Midpoint\n%d\n%d",
freq["CBT+", "Midpoint"],
freq["Relaxation", "Midpoint"]),
sprintf("Post\n%d\n%d",
freq["CBT+", "Post"],
freq["Relaxation", "Post"]),
sprintf("Follow-up\n%d\n%d",
freq["CBT+", "Follow-Up"],
freq["Relaxation", "Follow-Up"]))
sig <- symnum(subset(coef(depfit), grepl("MDIF", Label))$pval,
legend = FALSE, na = "",
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "n.s."))[-1]
yval.sig <- with(est, max(UpperCI) + (abs(diff(range(c(LowerCI, UpperCI)))) * .02))
yval.es <- with(est, min(LowerCI) - (abs(diff(range(c(LowerCI, UpperCI)))) * .03))
})
## make the plot
p.dep <- with(dep.plotdat,
ggplot(est, aes(Timepoint + ifelse(Condition=="CBT-I+Light", .03, -.03),
est, ymin = LowerCI, ymax = UpperCI,
colour = Condition, shape = Condition)) +
geom_hline(yintercept = 50, linetype = 3, colour = "grey50") +
geom_line(size = .1, alpha = .5) + geom_pointrange(size=1) +
scale_x_continuous("", breaks = 0:3, labels = labels) +
scale_y_continuous("Depressive Symptoms (Adjusted)",
breaks = seq(44, 54, 2)) +
annotate("text", x = 1, y = yval.sig, label = sig[1]) +
annotate("text", x = 2, y = yval.sig, label = sig[2]) +
annotate("text", x = 3, y = yval.sig, label = sig[3]) +
annotate("text", x = 1, y = yval.es, label = ES[1]) +
annotate("text", x = 2, y = yval.es, label = ES[2]) +
annotate("text", x = 3, y = yval.es, label = ES[3]))
p.dep <- set_palette(p.dep, "jco")
print(p.dep)
Anxiety Symptoms
anxdat <- copy(reshape(d.all[, .(
ID, dcond, dstrat4, ANX_TScore, Survey)],
v.names = "ANX_TScore",
timevar = "Survey",
idvar = c("ID", "dcond", "dstrat4"),
direction = "wide"))
setnames(anxdat, names(anxdat),
c("id", "dcond", "dstrat4",
"ANX0", "ANX1", "ANX2", "ANX3", "ANX4"))
anxdat <- anxdat[!is.na(ANX1) | !is.na(ANX2) | !is.na(ANX3) | !is.na(ANX4)]
anxfit <- mplusModeler(mplusObject(
TITLE = "Primary ANX Analysis;",
ANALYSIS = "ESTIMATOR = ML;",
MODEL = "
[ANX1 - ANX4@0]; ! fix residual intercepts to 0
ANX1* ANX2* ANX3* ANX4* (resvar); ! homogenous residual variance
int BY ANX1@1 ANX2@1 ANX3@1 ANX4@1;
slope BY ANX1@0 ANX2@0.5 ANX3@1 ANX4@1;
slope2 BY ANX1@0 ANX2@0 ANX3@0 ANX4@1;
int ON dstrat4 dcond@0 (i3 - i4);
slope ON dcond (s4);
slope2 ON dcond (s24);
[int* slope* slope2*] (rm1 - rm3);
int* slope* slope2*0 (rv1-rv3);
int WITH slope*0;
slope2 WITH int*0 slope*0;
",
OUTPUT = "STDYX; CINTERVAL;",
MODELCONSTRAINT = sprintf(
mcsetup,
mean(anxdat$dstrat4)),
usevariables = c("ANX1", "ANX2", "ANX3", "ANX4",
"dstrat4", "dcond"),
rdata = as.data.frame(anxdat)),
dataout = "03_HighISI_anx.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 03_HighISI_anx.inp
## Wrote data to: 03_HighISI_anx.dat
## model information, coefficients, and effect sizes
summary(anxfit)
## Primary ANX Analysis
##
## Estimated using ML
## Number of obs: 67, number of (free) parameters: 13
##
## Model: Chi2(df = 9) = 12.898, p = 0.1673
## Baseline model: Chi2(df = 14) = 99.676, p = 0
##
## Fit Indices:
##
## CFI = 0.955, TLI = 0.929, SRMR = 0.148
## RMSEA = 0.08, 90% CI [0, 0.171], p < .05 = 0.272
## AIC = 1429.206, BIC = 1457.867
cbind(
coef(anxfit, param = c("reg", "exp", "new")),
confint(anxfit, param = c("reg", "exp", "new"))[,-1])
## Label est se pval LowerCI UpperCI
## 13 INT<-DSTRAT4 -1.295 2.102 0.538 -5.416 2.826
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND 2.485 2.266 0.273 -1.957 6.927
## 16 SLOPE2<-DCOND -3.704 2.760 0.180 -9.113 1.706
## 20 ANX1<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 21 ANX2<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 22 ANX3<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 23 ANX4<-Intercepts 0.000 0.000 999.000 0.000 0.000
## 24 INT<-Intercepts 54.330 1.542 0.000 51.308 57.352
## 25 SLOPE<-Intercepts -2.704 1.731 0.118 -6.096 0.689
## 26 SLOPE2<-Intercepts -1.011 2.098 0.630 -5.122 3.101
## 34 ADJREF 53.692 1.138 0.000 51.463 55.922
## 35 MTAU1 53.692 1.138 0.000 51.463 55.922
## 36 MTAU2 52.340 1.241 0.000 49.909 54.772
## 37 MTAU3 50.989 1.812 0.000 47.437 54.540
## 38 MTAU4 49.978 2.059 0.000 45.943 54.013
## 39 MCBT1 53.692 1.138 0.000 51.463 55.922
## 40 MCBT2 53.583 1.188 0.000 51.255 55.911
## 41 MCBT3 53.473 1.670 0.000 50.200 56.747
## 42 MCBT4 48.759 1.770 0.000 45.290 52.228
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 1.242 1.133 0.273 -0.979 3.464
## 45 MDIF3 2.485 2.266 0.273 -1.957 6.927
## 46 MDIF4 -1.219 2.534 0.630 -6.186 3.748
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 0.131 0.120 0.275 -0.104 0.367
## 49 BD3 0.262 0.240 0.275 -0.208 0.733
## 50 BD4 -0.129 0.268 0.631 -0.653 0.396
## 51 S1CBT -0.219 1.590 0.890 -3.335 2.897
## 52 S2CBT -4.715 1.788 0.008 -8.218 -1.211
## create all the data for the plot, including means, effect sizes, Ns
## confidence intervals and significance
anx.plotdat <- list()
anx.plotdat <- within(anx.plotdat, {
est <- as.data.table(cbind(coef(anxfit), confint(anxfit)[, -1]))[35:42, ]
est[, Condition := rep(c("TAU+", "CBT-I+Light"), each = 4)]
est[, Timepoint := c(0:3, 0:3)]
freq <- xtabs(~ condition + Survey, data = d.all[!is.na(ANX_Total)])
ES <- c(sprintf("ES = %s",
gsub("^0", "",
format(round(abs(subset(coef(anxfit, param = c("new")),
grepl("BD", Label))$est), 2),
digits = 2, nsmall = 2))))[-1]
labels <- c(
sprintf("Baseline\nCBT-I+Light N: %d \nTAU+ N: %d ",
freq["CBT+", "Baseline"],
freq["Relaxation", "Baseline"]),
sprintf("Midpoint\n%d\n%d",
freq["CBT+", "Midpoint"],
freq["Relaxation", "Midpoint"]),
sprintf("Post\n%d\n%d",
freq["CBT+", "Post"],
freq["Relaxation", "Post"]),
sprintf("Follow-up\n%d\n%d",
freq["CBT+", "Follow-Up"],
freq["Relaxation", "Follow-Up"]))
sig <- symnum(subset(coef(anxfit), grepl("MDIF", Label))$pval,
legend = FALSE, na = "",
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "n.s."))[-1]
yval.sig <- with(est, max(UpperCI) + (abs(diff(range(c(LowerCI, UpperCI)))) * .02))
yval.es <- with(est, min(LowerCI) - (abs(diff(range(c(LowerCI, UpperCI)))) * .03))
})
## make the plot
p.anx <- with(anx.plotdat,
ggplot(est, aes(Timepoint + ifelse(Condition=="CBT-I+Light", .03, -.03),
est, ymin = LowerCI, ymax = UpperCI,
colour = Condition, shape = Condition)) +
geom_hline(yintercept = 50, linetype = 3, colour = "grey50") +
geom_line(size = .1, alpha = .5) + geom_pointrange(size=1) +
scale_x_continuous("", breaks = 0:3, labels = labels) +
scale_y_continuous("Anxiety Symptoms (Adjusted)",
breaks = seq(40, 56, 2)) +
annotate("text", x = 1, y = yval.sig, label = sig[1]) +
annotate("text", x = 2, y = yval.sig, label = sig[2]) +
annotate("text", x = 3, y = yval.sig, label = sig[3]) +
annotate("text", x = 1, y = yval.es, label = ES[1]) +
annotate("text", x = 2, y = yval.es, label = ES[2]) +
annotate("text", x = 3, y = yval.es, label = ES[3]))
p.anx <- set_palette(p.anx, "jco")
print(p.anx)
Final Panel Plot
p.all <- ggarrange(
p.isi + theme(plot.margin = margin(2, 14, 0, 32, "pt")) + coord_cartesian(clip="off"),
p.sri + theme(plot.margin = margin(2, 14, 0, 32, "pt")) + coord_cartesian(clip="off"),
p.sd + theme(plot.margin = margin(2, 14, 0, 32, "pt")) + coord_cartesian(clip="off"),
p.fa + theme(plot.margin = margin(2, 14, 0, 32, "pt")) + coord_cartesian(clip="off"),
p.dep + theme(plot.margin = margin(2, 14, 0, 32, "pt")) + coord_cartesian(clip="off"),
p.anx + theme(plot.margin = margin(2, 14, 0, 32, "pt")) + coord_cartesian(clip="off"),
labels = c("A", "B", "C", "D", "E", "F"),
ncol=2, nrow=3,
legend = "bottom",
common.legend = TRUE)
print(p.all)
pdf(file = "figure_pros_highisi.pdf", width = 9, height = 12)
print(p.all)
dev.off()
## png
## 2
# Demographics Table
## used to examine whether conditions differ on any demographics in that subset with
## high initial ISI. Results showed no significant differences between conditions
## in those with high initial ISI.
egltable(c(
"Age", "CHILD2", "CHILDLV18", "ETH2", "MARI2",
"EMPL3", "EDU3", "INCOM3", "SES_1", "SMO",
"CANSTG", "RADIATION", "BIOLOGICAL", "HORMONAL", "MENOPAUSE", "SURG", "CHEM",
"MHTX", "SleepTX", "ISI_Total"),
g = "condition",
data = d.all[Survey=="Baseline"], simChisq=TRUE)
## Warning in chisq.test(tabs, correct = FALSE, simulate.p.value = simChisq, :
## cannot compute simulated p-value with zero marginals
## Warning in chisq.test(tabs, correct = FALSE, simulate.p.value = simChisq, : Chi-
## squared approximation may be incorrect
## Relaxation M (SD)/N (%) CBT+ M (SD)/N (%)
## 1: Age 51.67 (11.64) 48.36 (10.08)
## 2: CHILD2
## 3: >= 1 Child 21 (63.6) 21 (61.8)
## 4: No Children 12 (36.4) 13 (38.2)
## 5: CHILDLV18
## 6: Child < 18y 10 (30.3) 12 (35.3)
## 7: No Child < 18y 23 (69.7) 22 (64.7)
## 8: ETH2
## 9: Non-White 2 (6.1) 4 (11.8)
## 10: White 31 (93.9) 30 (88.2)
## 11: MARI2
## 12: Married/living as 17 (51.5) 23 (67.6)
## 13: Other 16 (48.5) 11 (32.4)
## 14: EMPL3
## 15: Employed 21 (63.6) 21 (61.8)
## 16: Retired 7 (21.2) 6 (17.6)
## 17: Not employed 5 (15.2) 7 (20.6)
## 18: EDU3
## 19: < Bachelor 14 (42.4) 10 (29.4)
## 20: Bachelor 9 (27.3) 10 (29.4)
## 21: > Bachelor 10 (30.3) 14 (41.2)
## 22: INCOM3
## 23: <= $50,000 6 (20.0) 7 (22.6)
## 24: $50,001 - $100,000 16 (53.3) 11 (35.5)
## 25: > $100,000 8 (26.7) 13 (41.9)
## 26: SES_1 6.05 (1.86) 6.41 (1.41)
## 27: SMO
## 28: No 32 (100.0) 34 (100.0)
## 29: Yes 0 (0.0) 0 (0.0)
## 30: CANSTG
## 31: 1 4 (12.5) 4 (12.1)
## 32: 2 11 (34.4) 10 (30.3)
## 33: 3 8 (25.0) 9 (27.3)
## 34: 4 9 (28.1) 10 (30.3)
## 35: RADIATION
## 36: Current 3 (9.4) 4 (12.1)
## 37: No 9 (28.1) 4 (12.1)
## 38: Planned 16 (50.0) 19 (57.6)
## 39: UnknownPlan 4 (12.5) 6 (18.2)
## 40: BIOLOGICAL
## 41: Current 11 (37.9) 7 (25.0)
## 42: No 6 (20.7) 12 (42.9)
## 43: Planned 4 (13.8) 4 (14.3)
## 44: UnknownPlan 8 (27.6) 5 (17.9)
## 45: HORMONAL
## 46: Current 7 (22.6) 4 (12.9)
## 47: No 10 (32.3) 8 (25.8)
## 48: Planned 8 (25.8) 7 (22.6)
## 49: UnknownPlan 6 (19.4) 12 (38.7)
## 50: MENOPAUSE
## 51: Current 8 (25.0) 8 (26.7)
## 52: No 12 (37.5) 14 (46.7)
## 53: Past 12 (37.5) 8 (26.7)
## 54: SURG
## 55: No 11 (34.4) 10 (30.3)
## 56: Yes 21 (65.6) 23 (69.7)
## 57: CHEM
## 58: No 0 (0.0) 1 (3.0)
## 59: Yes 32 (100.0) 32 (97.0)
## 60: MHTX
## 61: No 14 (45.2) 11 (33.3)
## 62: Yes 17 (54.8) 22 (66.7)
## 63: SleepTX
## 64: No 15 (46.9) 14 (42.4)
## 65: Yes 17 (53.1) 19 (57.6)
## 66: ISI_Total 14.97 (4.69) 15.09 (5.10)
## Relaxation M (SD)/N (%) CBT+ M (SD)/N (%)
## Test
## 1: t(df=65) = 1.25, p = .218, d = 0.30
## 2: Chi-square = 0.03, simulated, p = 1.000, Phi = 0.02
## 3:
## 4:
## 5: Chi-square = 0.19, simulated, p = .796, Phi = 0.05
## 6:
## 7:
## 8: Chi-square = 0.67, simulated, p = .672, Phi = 0.10
## 9:
## 10:
## 11: Chi-square = 1.81, simulated, p = .218, Phi = 0.16
## 12:
## 13:
## 14: Chi-square = 0.40, simulated, p = .883, Cramer's V = 0.08
## 15:
## 16:
## 17:
## 18: Chi-square = 1.37, simulated, p = .508, Cramer's V = 0.14
## 19:
## 20:
## 21:
## 22: Chi-square = 2.18, simulated, p = .384, Cramer's V = 0.19
## 23:
## 24:
## 25:
## 26: t(df=62) = -0.86, p = .396, d = 0.21
## 27: Chi-square = NaN, simulated, NA, Phi = NaN
## 28:
## 29:
## 30: Chi-square = 0.14, simulated, p = 1.000, Cramer's V = 0.05
## 31:
## 32:
## 33:
## 34:
## 35: Chi-square = 2.71, simulated, p = .430, Cramer's V = 0.20
## 36:
## 37:
## 38:
## 39:
## 40: Chi-square = 3.56, simulated, p = .307, Cramer's V = 0.25
## 41:
## 42:
## 43:
## 44:
## 45: Chi-square = 3.11, simulated, p = .392, Cramer's V = 0.22
## 46:
## 47:
## 48:
## 49:
## 50: Chi-square = 0.89, simulated, p = .654, Cramer's V = 0.12
## 51:
## 52:
## 53:
## 54: Chi-square = 0.12, simulated, p = .795, Phi = 0.04
## 55:
## 56:
## 57: Chi-square = 0.98, simulated, p = 1.000, Phi = 0.12
## 58:
## 59:
## 60: Chi-square = 0.94, simulated, p = .443, Phi = 0.12
## 61:
## 62:
## 63: Chi-square = 0.13, simulated, p = .805, Phi = 0.04
## 64:
## 65:
## 66: t(df=65) = -0.10, p = .923, d = 0.02
## Test
Session Info for Documentation
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] MplusAutomation_0.8 data.table_1.14.0 reshape2_1.4.4
## [4] JWileymisc_1.2.1 scales_1.1.1 ggpubr_0.4.0
## [7] ggplot2_3.3.3 rmarkdown_2.8 knitr_1.33
##
## loaded via a namespace (and not attached):
## [1] TH.data_1.0-10 VGAM_1.1-5 minqa_1.2.4
## [4] colorspace_2.0-1 ggsignif_0.6.1 ellipsis_0.3.2
## [7] rio_0.5.26 estimability_1.3 htmlTable_2.2.1
## [10] base64enc_0.1-3 rstudioapi_0.13 mice_3.13.0
## [13] farver_2.1.0 lavaan_0.6-8 MatrixModels_0.5-0
## [16] fansi_0.5.0 mvtnorm_1.1-1 codetools_0.2-18
## [19] splines_4.1.0 mnormt_2.0.2 robustbase_0.93-7
## [22] texreg_1.37.5 Formula_1.2-4 jsonlite_1.7.2
## [25] nloptr_1.2.2.2 broom_0.7.6 cluster_2.1.2
## [28] png_0.1-7 httr_1.4.2 compiler_4.1.0
## [31] emmeans_1.6.0 backports_1.2.1 assertthat_0.2.1
## [34] Matrix_1.3-3 htmltools_0.5.1.1 quantreg_5.85
## [37] tools_4.1.0 coda_0.19-4 gtable_0.3.0
## [40] glue_1.4.2 dplyr_1.0.6 ggthemes_4.2.4
## [43] Rcpp_1.0.6 carData_3.0-4 cellranger_1.1.0
## [46] jquerylib_0.1.4 vctrs_0.3.8 nlme_3.1-152
## [49] conquer_1.0.2 psych_2.1.3 xfun_0.23
## [52] stringr_1.4.0 proto_1.0.0 openxlsx_4.2.3
## [55] lme4_1.1-27 lifecycle_1.0.0 extraoperators_0.1.1
## [58] rstatix_0.7.0 polspline_1.1.19 DEoptimR_1.0-9
## [61] MASS_7.3-54 zoo_1.8-9 hms_1.1.0
## [64] parallel_4.1.0 sandwich_3.0-1 SparseM_1.81
## [67] RColorBrewer_1.1-2 yaml_2.2.1 curl_4.3.1
## [70] gridExtra_2.3 pander_0.6.3 sass_0.4.0
## [73] rms_6.2-0 rpart_4.1-15 latticeExtra_0.6-29
## [76] stringi_1.6.2 highr_0.9 checkmate_2.0.0
## [79] boot_1.3-28 zip_2.1.1 rlang_0.4.11
## [82] pkgconfig_2.0.3 matrixStats_0.58.0 evaluate_0.14
## [85] lattice_0.20-44 purrr_0.3.4 labeling_0.4.2
## [88] htmlwidgets_1.5.3 cowplot_1.1.1 tidyselect_1.1.1
## [91] ggsci_2.9 plyr_1.8.6 magrittr_2.0.1
## [94] R6_2.5.0 generics_0.1.0 multcompView_0.1-8
## [97] Hmisc_4.5-0 multcomp_1.4-17 DBI_1.1.1
## [100] gsubfn_0.7 pillar_1.6.1 haven_2.4.1
## [103] foreign_0.8-81 withr_2.4.2 mgcv_1.8-35
## [106] survival_3.2-11 abind_1.4-5 nnet_7.3-16
## [109] tibble_3.1.2 crayon_1.4.1 car_3.0-10
## [112] utf8_1.2.1 tmvnsim_1.0-2 jpeg_0.1-8.1
## [115] grid_4.1.0 readxl_1.3.1 pbivnorm_0.6.0
## [118] forcats_0.5.1 digest_0.6.27 xtable_1.8-4
## [121] tidyr_1.1.3 stats4_4.1.0 munsell_0.5.0
## [124] bslib_0.2.5.1