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[, dstrat2 := as.integer(baseline_isi_high == "ISI <= 7" & baseline_stage_high == "Stage >= 3")]
d.all[, dstrat4 := as.integer(baseline_isi_high == "ISI >= 8" & baseline_stage_high == "Stage >= 3")]
## select only initially high stage (>= 3) cancer at screening
## note that this means we cannot use dstrat1 or dstrat3
## as those only occur in women with stage 1 or 2
d.all <- d.all[baseline_stage_high == "Stage >= 3"]
## 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 Symptoms
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 = "06_AdvancedStage_isi.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 06_AdvancedStage_isi.inp
## Wrote data to: 06_AdvancedStage_isi.dat
## model information, coefficients, and effect sizes
summary(isifit)
## Primary ISI Analysis
##
## Estimated using ML
## Number of obs: 42, number of (free) parameters: 13
##
## Model: Chi2(df = 9) = 12.581, p = 0.1825
## Baseline model: Chi2(df = 14) = 88.083, p = 0
##
## Fit Indices:
##
## CFI = 0.952, TLI = 0.925, SRMR = 0.17
## RMSEA = 0.097, 90% CI [0, 0.213], p < .05 = 0.25
## AIC = 775.677, BIC = 798.267
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 5.695 1.864 0.002 2.041 9.349
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND -1.515 1.820 0.405 -5.083 2.052
## 16 SLOPE2<-DCOND 0.847 1.601 0.597 -2.290 3.984
## 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 8.383 1.642 0.000 5.164 11.601
## 25 SLOPE<-Intercepts -2.520 1.416 0.075 -5.295 0.255
## 26 SLOPE2<-Intercepts -2.836 1.197 0.018 -5.182 -0.491
## 34 ADJREF 12.857 0.741 0.000 11.404 14.310
## 35 MTAU1 12.857 0.741 0.000 11.404 14.310
## 36 MTAU2 11.597 0.888 0.000 9.856 13.339
## 37 MTAU3 10.337 1.425 0.000 7.544 13.131
## 38 MTAU4 7.501 1.591 0.000 4.383 10.619
## 39 MCBT1 12.857 0.741 0.000 11.404 14.310
## 40 MCBT2 10.840 0.820 0.000 9.232 12.447
## 41 MCBT3 8.822 1.254 0.000 6.365 11.280
## 42 MCBT4 6.833 1.411 0.000 4.067 9.599
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 -0.758 0.910 0.405 -2.542 1.026
## 45 MDIF3 -1.515 1.820 0.405 -5.083 2.052
## 46 MDIF4 -0.668 2.044 0.744 -4.675 3.339
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 -0.155 0.187 0.407 -0.522 0.212
## 49 BD3 -0.310 0.374 0.407 -1.044 0.424
## 50 BD4 -0.137 0.419 0.744 -0.958 0.684
## 51 S1CBT -4.035 1.246 0.001 -6.477 -1.593
## 52 S2CBT -1.989 1.057 0.060 -4.062 0.083
## 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 = "06_AdvancedStage_sri.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 06_AdvancedStage_sri.inp
## Wrote data to: 06_AdvancedStage_sri.dat
## model information, coefficients, and effect sizes
summary(srifit)
## Primary SRI Analysis
##
## Estimated using ML
## Number of obs: 42, number of (free) parameters: 13
##
## Model: Chi2(df = 9) = 17.198, p = 0.0457
## Baseline model: Chi2(df = 14) = 83.608, p = 0
##
## Fit Indices:
##
## CFI = 0.882, TLI = 0.817, SRMR = 0.136
## RMSEA = 0.147, 90% CI [0.02, 0.252], p < .05 = 0.076
## AIC = 897.837, BIC = 920.427
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 6.904 2.696 0.010 1.620 12.187
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND -3.458 2.577 0.180 -8.508 1.592
## 16 SLOPE2<-DCOND -0.895 3.099 0.773 -6.970 5.180
## 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 51.193 2.410 0.000 46.469 55.917
## 25 SLOPE<-Intercepts -0.588 2.031 0.772 -4.568 3.393
## 26 SLOPE2<-Intercepts -2.582 2.320 0.266 -7.128 1.965
## 34 ADJREF 56.618 1.148 0.000 54.367 58.868
## 35 MTAU1 56.618 1.148 0.000 54.367 58.868
## 36 MTAU2 56.324 1.291 0.000 53.793 58.854
## 37 MTAU3 56.030 2.019 0.000 52.072 59.988
## 38 MTAU4 53.448 2.594 0.000 48.364 58.533
## 39 MCBT1 56.618 1.148 0.000 54.367 58.868
## 40 MCBT2 54.595 1.201 0.000 52.240 56.949
## 41 MCBT3 52.572 1.788 0.000 49.068 56.076
## 42 MCBT4 49.095 2.413 0.000 44.365 53.825
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 -1.729 1.288 0.180 -4.254 0.796
## 45 MDIF3 -3.458 2.577 0.180 -8.508 1.592
## 46 MDIF4 -4.353 3.249 0.180 -10.722 2.016
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 -0.227 0.171 0.184 -0.562 0.108
## 49 BD3 -0.454 0.342 0.184 -1.125 0.216
## 50 BD4 -0.572 0.431 0.185 -1.417 0.273
## 51 S1CBT -4.046 1.805 0.025 -7.583 -0.508
## 52 S2CBT -3.477 2.147 0.105 -7.685 0.731
## 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 = "06_AdvancedStage_sd.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 06_AdvancedStage_sd.inp
## Wrote data to: 06_AdvancedStage_sd.dat
## model information, coefficients, and effect sizes
summary(sdfit)
## Primary SD Analysis
##
## Estimated using ML
## Number of obs: 42, number of (free) parameters: 12
##
## Model: Chi2(df = 10) = 29.439, p = 0.0011
## Baseline model: Chi2(df = 14) = 72.64, p = 0
##
## Fit Indices:
##
## CFI = 0.669, TLI = 0.536, SRMR = 0.489
## RMSEA = 0.215, 90% CI [0.128, 0.307], p < .05 = 0.003
## AIC = 916.915, BIC = 937.767
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 9.138 2.842 0.001 3.568 14.708
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND -1.234 2.850 0.665 -6.819 4.352
## 16 SLOPE2<-DCOND -0.286 3.083 0.926 -6.330 5.757
## 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 47.949 2.488 0.000 43.074 52.825
## 25 SLOPE<-Intercepts -3.312 2.196 0.132 -7.617 0.993
## 26 SLOPE2<-Intercepts -3.185 2.306 0.167 -7.706 1.335
## 34 ADJREF 55.129 1.095 0.000 52.983 57.275
## 35 MTAU1 55.129 1.095 0.000 52.983 57.275
## 36 MTAU2 53.473 1.381 0.000 50.766 56.180
## 37 MTAU3 51.817 2.242 0.000 47.422 56.212
## 38 MTAU4 48.632 2.588 0.000 43.559 53.705
## 39 MCBT1 55.129 1.095 0.000 52.983 57.275
## 40 MCBT2 52.856 1.272 0.000 50.363 55.349
## 41 MCBT3 50.583 1.973 0.000 46.717 54.450
## 42 MCBT4 47.112 2.392 0.000 42.424 51.799
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 -0.617 1.425 0.665 -3.410 2.176
## 45 MDIF3 -1.234 2.850 0.665 -6.819 4.352
## 46 MDIF4 -1.520 3.216 0.637 -7.824 4.784
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 -0.083 0.192 0.665 -0.460 0.294
## 49 BD3 -0.166 0.384 0.665 -0.920 0.587
## 50 BD4 -0.205 0.434 0.637 -1.055 0.646
## 51 S1CBT -4.546 1.927 0.018 -8.322 -0.770
## 52 S2CBT -3.471 2.121 0.102 -7.628 0.685
## 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 = "06_AdvancedStage_fa.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 06_AdvancedStage_fa.inp
## Wrote data to: 06_AdvancedStage_fa.dat
## model information, coefficients, and effect sizes
summary(fafit)
## Primary FA Analysis
##
## Estimated using ML
## Fit Indices:
##
## CFI = NA, TLI = NA, SRMR = NA
## RMSEA = NA, 90% CI [NA, NA], p < .05 = NA
## AIC = NA, BIC = NA
cbind(
coef(fafit, param = c("reg", "exp", "new")),
confint(fafit, param = c("reg", "exp", "new"))[,-1])
## Error in `[.data.frame`(out, , c("Label", estimate, se, pvalue)): object 'se' not found
## 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))
})
## Error in `[.data.frame`(out, , c("Label", estimate, se, pvalue)): object 'se' not found
## 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]))
## Error in ggplot(est, aes(Timepoint + ifelse(Condition == "CBT-I+Light", : object 'est' not found
p.fa <- set_palette(p.fa, "jco")
## Error in set_palette(p.fa, "jco"): object 'p.fa' not found
print(p.fa)
## Error in print(p.fa): object 'p.fa' not found
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 = "06_AdvancedStage_dep.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 06_AdvancedStage_dep.inp
## Wrote data to: 06_AdvancedStage_dep.dat
## model information, coefficients, and effect sizes
summary(depfit)
## Primary DEP Analysis
##
## Estimated using ML
## Number of obs: 42, number of (free) parameters: 13
##
## Model: Chi2(df = 9) = 4.806, p = 0.8509
## Baseline model: Chi2(df = 14) = 74.515, p = 0
##
## Fit Indices:
##
## CFI = 1, TLI = 1, SRMR = 0.094
## RMSEA = 0, 90% CI [0, 0.097], p < .05 = 0.888
## AIC = 897.647, BIC = 920.237
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 3.513 3.066 0.252 -2.496 9.522
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND 2.694 2.780 0.333 -2.755 8.143
## 16 SLOPE2<-DCOND -4.952 2.572 0.054 -9.992 0.088
## 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 47.741 2.733 0.000 42.384 53.099
## 25 SLOPE<-Intercepts -3.562 2.153 0.098 -7.782 0.657
## 26 SLOPE2<-Intercepts 3.520 1.894 0.063 -0.191 7.232
## 34 ADJREF 50.502 1.290 0.000 47.973 53.031
## 35 MTAU1 50.502 1.290 0.000 47.973 53.031
## 36 MTAU2 48.721 1.478 0.000 45.824 51.617
## 37 MTAU3 46.940 2.241 0.000 42.548 51.331
## 38 MTAU4 50.460 2.128 0.000 46.288 54.632
## 39 MCBT1 50.502 1.290 0.000 47.973 53.031
## 40 MCBT2 50.068 1.385 0.000 47.354 52.782
## 41 MCBT3 49.633 1.993 0.000 45.727 53.540
## 42 MCBT4 48.202 2.012 0.000 44.259 52.144
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 1.347 1.390 0.333 -1.377 4.071
## 45 MDIF3 2.694 2.780 0.333 -2.755 8.143
## 46 MDIF4 -2.258 2.561 0.378 -7.278 2.762
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 0.158 0.164 0.335 -0.164 0.480
## 49 BD3 0.317 0.328 0.335 -0.327 0.960
## 50 BD4 -0.265 0.302 0.380 -0.858 0.327
## 51 S1CBT -0.868 1.899 0.647 -4.591 2.854
## 52 S2CBT -1.432 1.762 0.416 -4.884 2.021
## 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 = "06_AdvancedStage_anx.dat", run = TRUE, writeData = "always", hashfilename = FALSE)
## Wrote model to: 06_AdvancedStage_anx.inp
## Wrote data to: 06_AdvancedStage_anx.dat
## model information, coefficients, and effect sizes
summary(anxfit)
## Primary ANX Analysis
##
## Estimated using ML
## Number of obs: 42, number of (free) parameters: 13
##
## Model: Chi2(df = 9) = 17.378, p = 0.0431
## Baseline model: Chi2(df = 14) = 88.59, p = 0
##
## Fit Indices:
##
## CFI = 0.888, TLI = 0.825, SRMR = 0.148
## RMSEA = 0.149, 90% CI [0.026, 0.253], p < .05 = 0.072
## AIC = 924.448, BIC = 947.038
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 2.176 3.560 0.541 -4.801 9.152
## 14 INT<-DCOND 0.000 0.000 999.000 0.000 0.000
## 15 SLOPE<-DCOND 2.076 2.867 0.469 -3.544 7.696
## 16 SLOPE2<-DCOND -4.772 3.662 0.193 -11.949 2.405
## 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 50.487 3.156 0.000 44.300 56.674
## 25 SLOPE<-Intercepts -2.927 2.177 0.179 -7.193 1.340
## 26 SLOPE2<-Intercepts 2.597 2.706 0.337 -2.707 7.900
## 34 ADJREF 52.196 1.461 0.000 49.332 55.061
## 35 MTAU1 52.196 1.461 0.000 49.332 55.061
## 36 MTAU2 50.733 1.686 0.000 47.429 54.037
## 37 MTAU3 49.270 2.432 0.000 44.503 54.037
## 38 MTAU4 51.866 2.277 0.000 47.403 56.330
## 39 MCBT1 52.196 1.461 0.000 49.332 55.061
## 40 MCBT2 51.771 1.582 0.000 48.670 54.872
## 41 MCBT3 51.346 2.140 0.000 47.151 55.540
## 42 MCBT4 49.170 2.099 0.000 45.056 53.285
## 43 MDIF1 0.000 0.000 1.000 0.000 0.000
## 44 MDIF2 1.038 1.434 0.469 -1.772 3.848
## 45 MDIF3 2.076 2.867 0.469 -3.544 7.696
## 46 MDIF4 -2.696 2.773 0.331 -8.131 2.739
## 47 BD1 0.000 0.000 1.000 0.000 0.000
## 48 BD2 0.108 0.149 0.470 -0.185 0.400
## 49 BD3 0.216 0.299 0.470 -0.370 0.801
## 50 BD4 -0.280 0.290 0.334 -0.847 0.288
## 51 S1CBT -0.851 1.849 0.645 -4.475 2.774
## 52 S2CBT -2.175 2.371 0.359 -6.823 2.472
## 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", "E", "F"),
ncol=2, nrow=3,
legend = "bottom",
common.legend = TRUE)
print(p.all)
pdf(file = "06_advancedstage_figure_pros.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
## advanced stage. Results showed no significant differences between conditions
## except in education
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)
## Relaxation M (SD)/N (%) CBT+ M (SD)/N (%)
## 1: Age 52.65 (10.46) 47.36 (10.77)
## 2: CHILD2
## 3: >= 1 Child 13 (68.4) 13 (56.5)
## 4: No Children 6 (31.6) 10 (43.5)
## 5: CHILDLV18
## 6: Child < 18y 4 (21.1) 7 (30.4)
## 7: No Child < 18y 15 (78.9) 16 (69.6)
## 8: ETH2
## 9: Non-White 1 (5.3) 3 (13.0)
## 10: White 18 (94.7) 20 (87.0)
## 11: MARI2
## 12: Married/living as 12 (63.2) 15 (65.2)
## 13: Other 7 (36.8) 8 (34.8)
## 14: EMPL3
## 15: Employed 12 (63.2) 13 (56.5)
## 16: Retired 4 (21.1) 4 (17.4)
## 17: Not employed 3 (15.8) 6 (26.1)
## 18: EDU3
## 19: < Bachelor 11 (57.9) 4 (17.4)
## 20: Bachelor 1 (5.3) 11 (47.8)
## 21: > Bachelor 7 (36.8) 8 (34.8)
## 22: INCOM3
## 23: <= $50,000 4 (22.2) 4 (18.2)
## 24: $50,001 - $100,000 8 (44.4) 9 (40.9)
## 25: > $100,000 6 (33.3) 9 (40.9)
## 26: SES_1 6.47 (1.93) 6.45 (1.50)
## 27: SMO
## 28: No 18 (94.7) 23 (100.0)
## 29: Yes 1 (5.3) 0 (0.0)
## 30: CANSTG
## 31: 1 1 (5.3) 0 (0.0)
## 32: 2 2 (10.5) 0 (0.0)
## 33: 3 5 (26.3) 10 (45.5)
## 34: 4 11 (57.9) 12 (54.5)
## 35: RADIATION
## 36: Current 3 (15.8) 3 (14.3)
## 37: No 7 (36.8) 6 (28.6)
## 38: Planned 6 (31.6) 9 (42.9)
## 39: UnknownPlan 3 (15.8) 3 (14.3)
## 40: BIOLOGICAL
## 41: Current 8 (47.1) 6 (30.0)
## 42: No 3 (17.6) 8 (40.0)
## 43: Planned 2 (11.8) 1 (5.0)
## 44: UnknownPlan 4 (23.5) 5 (25.0)
## 45: HORMONAL
## 46: Current 5 (26.3) 3 (14.3)
## 47: No 7 (36.8) 7 (33.3)
## 48: Planned 2 (10.5) 5 (23.8)
## 49: UnknownPlan 5 (26.3) 6 (28.6)
## 50: MENOPAUSE
## 51: Current 2 (10.5) 6 (30.0)
## 52: No 8 (42.1) 8 (40.0)
## 53: Past 9 (47.4) 6 (30.0)
## 54: SURG
## 55: No 10 (52.6) 7 (31.8)
## 56: Yes 9 (47.4) 15 (68.2)
## 57: CHEM
## 58: No 0 (0.0) 1 (4.5)
## 59: Yes 19 (100.0) 21 (95.5)
## 60: MHTX
## 61: No 10 (55.6) 6 (27.3)
## 62: Yes 8 (44.4) 16 (72.7)
## 63: SleepTX
## 64: No 13 (68.4) 10 (45.5)
## 65: Yes 6 (31.6) 12 (54.5)
## 66: ISI_Total 12.79 (4.50) 12.96 (6.66)
## Relaxation M (SD)/N (%) CBT+ M (SD)/N (%)
## Test
## 1: t(df=40) = 1.60, p = .117, d = 0.50
## 2: Chi-square = 0.62, simulated, p = .530, Phi = 0.12
## 3:
## 4:
## 5: Chi-square = 0.47, simulated, p = .727, Phi = 0.11
## 6:
## 7:
## 8: Chi-square = 0.73, simulated, p = .613, Phi = 0.13
## 9:
## 10:
## 11: Chi-square = 0.02, simulated, p = 1.000, Phi = 0.02
## 12:
## 13:
## 14: Chi-square = 0.67, simulated, p = .760, Cramer's V = 0.13
## 15:
## 16:
## 17:
## 18: Chi-square = 11.39, simulated, p = .003, Cramer's V = 0.52
## 19:
## 20:
## 21:
## 22: Chi-square = 0.26, simulated, p = .845, Cramer's V = 0.08
## 23:
## 24:
## 25:
## 26: t(df=39) = 0.04, p = .972, d = 0.01
## 27: Chi-square = 1.24, simulated, p = .452, Phi = 0.17
## 28:
## 29:
## 30: Chi-square = 4.51, simulated, p = .176, Cramer's V = 0.33
## 31:
## 32:
## 33:
## 34:
## 35: Chi-square = 0.58, simulated, p = .967, Cramer's V = 0.12
## 36:
## 37:
## 38:
## 39:
## 40: Chi-square = 2.78, simulated, p = .511, Cramer's V = 0.27
## 41:
## 42:
## 43:
## 44:
## 45: Chi-square = 1.78, simulated, p = .658, Cramer's V = 0.21
## 46:
## 47:
## 48:
## 49:
## 50: Chi-square = 2.58, simulated, p = .325, Cramer's V = 0.26
## 51:
## 52:
## 53:
## 54: Chi-square = 1.82, simulated, p = .216, Phi = 0.21
## 55:
## 56:
## 57: Chi-square = 0.89, simulated, p = 1.000, Phi = 0.15
## 58:
## 59:
## 60: Chi-square = 3.30, simulated, p = .106, Phi = 0.29
## 61:
## 62:
## 63: Chi-square = 2.18, simulated, p = .208, Phi = 0.23
## 64:
## 65:
## 66: t(df=40) = -0.09, p = .926, d = 0.03
## 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