Beat the Blues
Ching-Yuan Chang
2022-10-09
Chapter 13: Analyzing Longitudinal Data I: Computerized Delivery of Cognitive Behavioral Therapy – Beat the Blues
1 Introduction
A data frame with 100 observations of 100 patients on the following 8 variables.
- drug: did the patient take anti-depressant drugs (No or Yes).
- length: the length of the current episode of depression, a factor with levels <6m (less than six months) and >6m (more than six months).
- treatment: treatment group, a factor with levels TAU (treatment as usual) and BtheB (Beat the Blues)
- bdi.pre: Beck Depression Inventory II before treatment.
- bdi.2m: Beck Depression Inventory II after two months.
- bdi.3m: Beck Depression Inventory II after one month follow-up.
- bdi.5m: Beck Depression Inventory II after three months follow-up.
- bdi.8m: Beck Depression Inventory II after six months follow-up.
2 Data
# 載入資料並將資料格式改成long-format
data(BtheB, package = "HSAUR3")
BtheB$subject <- factor(rownames(BtheB))
nobs <- nrow(BtheB)
BtheB_long <- reshape(BtheB, idvar = "subject",
varying = c("bdi.2m", "bdi.3m", "bdi.5m", "bdi.8m"),
direction = "long")
BtheB_long$time <- rep(c(2, 3, 5, 8), rep(nobs, 4))
subset(BtheB_long, subject %in% c("1", "2", "3")) drug length treatment bdi.pre subject time bdi
1.2m No >6m TAU 29 1 2 2
2.2m Yes >6m BtheB 32 2 2 16
3.2m Yes <6m TAU 25 3 2 20
1.3m No >6m TAU 29 1 3 2
2.3m Yes >6m BtheB 32 2 3 24
3.3m Yes <6m TAU 25 3 3 NA
1.5m No >6m TAU 29 1 5 NA
2.5m Yes >6m BtheB 32 2 5 17
3.5m Yes <6m TAU 25 3 5 NA
1.8m No >6m TAU 29 1 8 NA
[ reached 'max' / getOption("max.print") -- omitted 2 rows ]The resulting data.frame BtheB_long contains a number of missing values.
3 Visualization
# 對兩個實驗組(治療方法不同)在不同時間點測量的貝克憂鬱量表(BDI)分數作圖
data("BtheB", package = "HSAUR3")
layout(matrix(1:2, nrow = 1))
ylim <- range(BtheB[,grep("bdi", names(BtheB))],
na.rm = TRUE)
tau <- subset(BtheB, treatment == "TAU")[,grep("bdi", names(BtheB))]
boxplot(tau, main = "Treated as Usual", ylab = "BDI",
xlab = "Time (in months)", names = c(0, 2, 3, 5, 8),
ylim = ylim)
btheb <- subset(BtheB, treatment == "BtheB")[,grep("bdi", names(BtheB))]
boxplot(btheb, main = "Beat the Blues", ylab = "BDI",
xlab = "Time (in months)", names = c(0, 2, 3, 5, 8),
ylim = ylim)
4 Models
4.1 Random intercepts only
# 僅包含個人層次的截距隨機效果
BtheB_lmer1 <- lme4::lmer(bdi ~ bdi.pre + time + treatment + drug + length + (1 | subject), data = BtheB_long, REML = FALSE, na.action = na.omit)
summary(BtheB_lmer1)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: bdi ~ bdi.pre + time + treatment + drug + length + (1 | subject)
Data: BtheB_long
AIC BIC logLik deviance df.resid
1887.5 1916.6 -935.7 1871.5 272
Scaled residuals:
Min 1Q Median 3Q Max
-2.697 -0.503 -0.064 0.412 3.820
Random effects:
Groups Name Variance Std.Dev.
subject (Intercept) 48.8 6.98
Residual 25.1 5.01
Number of obs: 280, groups: subject, 97
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.5924 2.2424 2.49
bdi.pre 0.6397 0.0779 8.21
time -0.7048 0.1464 -4.81
treatmentBtheB -2.3291 1.6704 -1.39
drugYes -2.8250 1.7268 -1.64
length>6m 0.1971 1.6383 0.12
Correlation of Fixed Effects:
(Intr) bdi.pr time trtmBB drugYs
bdi.pre -0.682
time -0.238 0.020
tretmntBthB -0.390 0.121 0.018
drugYes -0.073 -0.237 -0.022 -0.323
length>6m -0.243 -0.242 -0.036 0.002 0.157# 檢驗固定效果並計算p值
cftest(BtheB_lmer1)
Simultaneous Tests for General Linear Hypotheses
Fit: lme4::lmer(formula = bdi ~ bdi.pre + time + treatment + drug +
length + (1 | subject), data = BtheB_long, REML = FALSE,
na.action = na.omit)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
(Intercept) == 0 5.5924 2.2424 2.49 0.013
bdi.pre == 0 0.6397 0.0779 8.21 2.2e-16
time == 0 -0.7048 0.1464 -4.81 1.5e-06
treatmentBtheB == 0 -2.3291 1.6704 -1.39 0.163
drugYes == 0 -2.8250 1.7268 -1.64 0.102
length>6m == 0 0.1971 1.6383 0.12 0.904
(Univariate p values reported)4.2 Random slopes and intercepts model
# 包含個人層次的截距隨機效果及時間(以月為單位)隨機斜率效果
BtheB_lmer2 <- lme4::lmer(bdi ~ bdi.pre + time + treatment + drug + length + (time | subject), data = BtheB_long, REML = FALSE, na.action = na.omit)
summary(BtheB_lmer2)Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: bdi ~ bdi.pre + time + treatment + drug + length + (time | subject)
Data: BtheB_long
AIC BIC logLik deviance df.resid
1891.0 1927.4 -935.5 1871.0 270
Scaled residuals:
Min 1Q Median 3Q Max
-2.312 -0.499 -0.070 0.397 3.719
Random effects:
Groups Name Variance Std.Dev. Corr
subject (Intercept) 50.734 7.123
time 0.184 0.429 -0.18
Residual 23.994 4.898
Number of obs: 280, groups: subject, 97
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.6139 2.2454 2.50
bdi.pre 0.6418 0.0778 8.25
time -0.7021 0.1541 -4.56
treatmentBtheB -2.3760 1.6713 -1.42
drugYes -2.8694 1.7277 -1.66
length>6m 0.1404 1.6385 0.09
Correlation of Fixed Effects:
(Intr) bdi.pr time trtmBB drugYs
bdi.pre -0.681
time -0.241 0.019
tretmntBthB -0.389 0.121 0.018
drugYes -0.074 -0.236 -0.022 -0.323
length>6m -0.243 -0.241 -0.035 0.002 0.158# 檢驗固定效果並計算p值
cftest(BtheB_lmer2)
Simultaneous Tests for General Linear Hypotheses
Fit: lme4::lmer(formula = bdi ~ bdi.pre + time + treatment + drug +
length + (time | subject), data = BtheB_long, REML = FALSE,
na.action = na.omit)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
(Intercept) == 0 5.6139 2.2454 2.50 0.012
bdi.pre == 0 0.6418 0.0778 8.25 2.2e-16
time == 0 -0.7021 0.1541 -4.56 5.2e-06
treatmentBtheB == 0 -2.3760 1.6713 -1.42 0.155
drugYes == 0 -2.8694 1.7277 -1.66 0.097
length>6m == 0 0.1404 1.6385 0.09 0.932
(Univariate p values reported)4.3 Compare models
anova(BtheB_lmer1, BtheB_lmer2)Data: BtheB_long
Models:
BtheB_lmer1: bdi ~ bdi.pre + time + treatment + drug + length + (1 | subject)
BtheB_lmer2: bdi ~ bdi.pre + time + treatment + drug + length + (time | subject)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
BtheB_lmer1 8 1888 1917 -935.7 1872
BtheB_lmer2 10 1891 1927 -935.5 1871 0.454 2 0.7974.4 Residuals
layout(matrix(1:2, ncol = 2))
qint <- ranef(BtheB_lmer1)$subject[["(Intercept)"]]
qres <- residuals(BtheB_lmer1)
qqnorm(qint, ylab = "Estimated random intercepts",
xlim = c(-3, 3), ylim = c(-20, 20),
main = "Random intercepts")
qqline(qint)
qqnorm(qres, xlim = c(-3, 3), ylim = c(-20, 20),
ylab = "Estimated residuals",
main = "Residuals")
qqline(qres)
There appear to be no large departures from linearity in either
plot.
4.5 Graph Plotting
bdi <- BtheB[, grep("bdi", names(BtheB))]
plot(1:4, rep(-0.5, 4), type = "n", axes = FALSE,
ylim = c(0, 50), xlab = "Months", ylab = "BDI")
axis(1, at = 1:4, labels = c(0, 2, 3, 5))
axis(2)
for (i in 1:4) {dropout <- is.na(bdi[,i + 1])
points(rep(i, nrow(bdi)) + ifelse(dropout, 0.05, -0.05),
jitter(bdi[,i]), pch = ifelse(dropout, 20, 1))
}
5 References
Hothorn, T., & Everitt, B. S. (2014). A Handbook of Statistical Analyses Using R (3rd Ed.). Chapman and Hall/CRC. https://doi.org/10.1201/b17081. pdf