This plot was initially designed for Lichen Ma’s RCT on VR-ABM for social anxiety - the plots in this document show mock data with similar characteristics.
All code to generate and plot mock data is displayed below.
The plot builds on a function to obtain Jacobson-Truax/Reliable Change Indices which I also wrote for this project. The JTRCI() function can be found here on github: https://github.com/AWKruijt/JT-RCI and as an online dashboard version where you can plot either mock data or use your own data: https://awkruijt.shinyapps.io/JTRCI_dashboard/
# generate data to plot:
# set n, means, and sds for each group as well as correlations between the time points.
# these values are very loosely based on real LSAS data observed in a study by Lichen Ma.
npp <- c(25, 25, 25, 25)
t1means <- c(71, 68, 70, 69)
t1sds <- c(20, 18, 18, 21)
t2means <- c(59, 61, 58, 63)
t2sds <- c(20, 21, 25, 20)
t3means <- c(60, 61, 50, 62)
t3sds <- c(19, 21, 24, 21)
t4means <- c(59, 57, 44, 59)
t4sds <- c(20, 21, 24, 23)
cor12 <- .85
cor13 <- .80
cor14 <- .70
cor23 <- .88
cor24 <- .70
cor34 <- .80
require(MASS)
tdf <- NULL
tdf <- data.frame(NULL)
# generate data separately for each group and bind into tdf:
for(g in 1:length(npp)) {
sigmat <- matrix(c(
t1sds[g]^2, cor12*t1sds[g]*t2sds[g], cor13*t1sds[g]*t3sds[g], cor14*t1sds[g]*t4sds[g],
cor12*t1sds[g]*t2sds[g], t2sds[g]^2, cor23*t2sds[g]*t3sds[g], cor24*t2sds[g]*t4sds[g],
cor13*t1sds[g]*t3sds[g], cor23*t2sds[g]*t3sds[g], t3sds[g]^2, cor34*t3sds[g]*t4sds[g],
cor14*t1sds[g]*t4sds[g], cor24*t2sds[g]*t4sds[g], cor34*t3sds[g]*t4sds[g], t4sds[g]^2),
nrow=4)
gencordat <- mvrnorm(n = npp[g], mu = c(t1means[g], t2means[g], t3means[g], t4means[g]),
Sigma = sigmat, empirical = T)
groupdf <- cbind.data.frame("ppid" = c((nrow(tdf)+1): (nrow(tdf) + npp[g])),
"group" = paste0("Group: ", c(rep(g, npp[g]))),
"t1" = gencordat[,1], "t2" = gencordat[,2], "t3" = gencordat[,3], "t4" = gencordat[,4] )
tdf <-rbind(tdf, groupdf)
}# Obtain Jacobson-Truax classifications:
# download the JTRCI function from github:
devtools::install_github("AWKruijt/JT-RCI")
library(JTRCI)
# create data frame that will hold the JT data for each of the timepoints in long format.
pdf <- cbind.data.frame("ppid" = tdf$ppid, # start plotdf with the ppids
"time" = "t1", # add time var with value BL
"change_abs" = 0, # add change var with value 0 for all at t1/BL
"Sdiff" = 0, # add Sdiff var with value 0 for all at t1/BL
"class_JTRCI" = NA) # add JTRCI classification var with value NA for all at t1
# set the reliabiliyt estimate for the symptom measure:
relEstQ <- .92
# this value is again based on the data by Lichen Ma that I originally created this plot for.
# it was obtained by running psych::alpha() on the item scores.
# for task data, I recommend using Sam Parson's splithalf() package to obtain an estimate of reliability.
# run the JTRCI function on t1-t2 data to get the first set of JT up data (outputed as JTRCIdf by the JTRCI() function)
JTRCI(data = tdf, pre= "t1", post = "t2", reliability = relEstQ, ppid = "ppid", plot = F, table = F)
JTRCIdf$time <- "t2" # add time var with value t2 to the JTRCIdf dataframe
pdf <- rbind(pdf, JTRCIdf [, c("ppid", "change_abs", "Sdiff", "class_JTRCI", "time")]) # rbind JTRCIdf with the t1 data with pdf
# repeat for the two follow-ups:
# t1 to t3
JTRCI(data = tdf, pre= "t1", post = "t3", reliability = relEstQ, ppid = "ppid", plot = F, table = F)
JTRCIdf$time <- "t3"
pdf <- rbind(pdf, JTRCIdf [, c("ppid", "change_abs", "Sdiff", "class_JTRCI", "time")])
# t1 to t4:
JTRCI(data = tdf, pre= "t1", post = "t4", reliability = relEstQ, ppid = "ppid", plot = F, table = F)
JTRCIdf$time <- "t4"
pdf <- rbind(pdf, JTRCIdf [, c("ppid", "change_abs", "Sdiff", "class_JTRCI", "time")])
# merge in group allocation from the tdf dataframe
pdf <- merge (pdf, tdf[, c("ppid", "group")], by = "ppid")# plotterdeplot:
# fix the colors to the classification labels:
pdf$JTclass_RCI <- factor(pdf$class_JTRCI, levels = c("recovered", "non reliably recovered", "improved", "unchanged", "deteriorated"))
# create vector of plot colours linked to specific level names:
cols <- c("recovered" = "orchid3", "non reliably recovered" = "palevioletred3", "improved" = "olivedrab3", "unchanged" = "skyblue3", "deteriorated" = "coral1")
Sdiff <- unique(JTRCIdf$Sdiff) # get the Sdiff value to use for positioning guide lines in the plot.
# a dataframe to define the lines with - this way we can add a legend for the lines to the plot:
linesdf <- data.frame(cc = c(" 1.96 Sdiff", " .84 Sdiff", " no change", " .84 Sdiff", " 1.96 Sdiff"),
ic = c(1.96 * Sdiff, .84 * Sdiff, 0 * Sdiff, -.84 * Sdiff, -1.96 * Sdiff))
linesdf$cc <- factor(linesdf$cc, levels = c (" 1.96 Sdiff", " .84 Sdiff", " no change"))
## if changing anything, notice the added spaces in the factor levels to align the linetype labels (!)
p <- ggplot(data=pdf, aes(x=time, y=change_abs, group = ppid)) +
geom_line(col = "gray80") +
geom_point(aes(col = class_JTRCI), size = 2.5) +
geom_point(data = pdf[!pdf$time == "t1",], size = 2.5, shape = 1, colour = "gray45", alpha = .8) +
facet_grid(~ group) +
expand_limits(y= min(pdf$change_abs)-(.20*(max(pdf$change_abs)-(min(pdf$change_abs))))) +
scale_y_continuous(breaks = scales::pretty_breaks(n = 5)) +
geom_abline(data= linesdf, mapping=aes(slope=0, intercept=ic, linetype= factor(cc)), col = "gray50") +
theme_classic() +
scale_colour_manual(breaks = c("recovered", "non reliably recovered", "improved", "unchanged", "deteriorated"), values = cols) +
scale_linetype_manual(values=c("dashed", "dotted", "1F")) +
theme(strip.background = element_blank()) + #remove facet strip background
labs(title = "longitudinal symptoms Jacobson-Truax plot", x = "time point", y = "symptom change \n(points relative to baseline score)", colour = "Jacobson-Truax \n classification:", linetype = "line guides: \n reliable change") + guides(linetype = guide_legend(order = 2), colour = guide_legend(order = 1, reverse = T))
# note: ggplot will issue a warning about missing values: there are no dots plotted for t1 (baseline) - these will make up the larger part (if not all) of the missing values.
# if there is missing data in the follow-up measures, the available data for that ppid will be plotted as intended and only the missing timepoint will be omitted from the plot.
p