Beat the Blues
Torsten Hothorn & Brian S. Everitt
2022-10-12
CHAPTER 13 Analyzing Longitudinal Data I: Computerized Delivery of Cognitive Behavioral Therapy – Beat the Blues
1 Introduction
- drug: 是否服用藥物 (Yes，No).
- length: 病症發作長度 (<6個月，>6個月)
- treatment: 治療方法(TAU，Beat the Blues)
- bdi.pre: 治療前的貝克抑鬱量表II分數
- bdi.2m: 兩個月後的貝克抑鬱量表II分數
- bdi.3m: 三個月後的貝克抑鬱量表II分數
- bdi.5m: 五個月後的貝克抑鬱量表II分數
- bdi.8m: 八個月後的貝克抑鬱量表II分數
2 Data
# 載入資料並更改格式
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
## 2.8m Yes >6m BtheB 32 2 8 20
## 3.8m Yes <6m TAU 25 3 8 NA3 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
library("lme4")
# 個人的截距隨機效果
BtheB_lmer1 <- 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.6975 -0.5026 -0.0638 0.4124 3.8203
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 48.78 6.984
## Residual 25.14 5.014
## Number of obs: 280, groups: subject, 97
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.59239 2.24244 2.494
## bdi.pre 0.63968 0.07789 8.212
## time -0.70476 0.14639 -4.814
## treatmentBtheB -2.32908 1.67036 -1.394
## drugYes -2.82495 1.72684 -1.636
## length>6m 0.19708 1.63832 0.120
##
## 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.157BtheB_lmer2 <- lmer(bdi ~ bdi.pre + time + treatment + drug +
length + (time | 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.6975 -0.5026 -0.0638 0.4124 3.8203
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 48.78 6.984
## Residual 25.14 5.014
## Number of obs: 280, groups: subject, 97
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.59239 2.24244 2.494
## bdi.pre 0.63968 0.07789 8.212
## time -0.70476 0.14639 -4.814
## treatmentBtheB -2.32908 1.67036 -1.394
## drugYes -2.82495 1.72684 -1.636
## length>6m 0.19708 1.63832 0.120
##
## 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.157anova(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 1887.5 1916.6 -935.75 1871.5
## BtheB_lmer2 10 1891.0 1927.4 -935.52 1871.0 0.4542 2 0.7969# 檢驗固定效果
library("multcomp")
cftest(BtheB_lmer1)##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: 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.59239 2.24244 2.494 0.0126 *
## bdi.pre == 0 0.63968 0.07789 8.212 2.22e-16 ***
## time == 0 -0.70476 0.14639 -4.814 1.48e-06 ***
## treatmentBtheB == 0 -2.32908 1.67036 -1.394 0.1632
## drugYes == 0 -2.82495 1.72684 -1.636 0.1019
## length>6m == 0 0.19708 1.63832 0.120 0.9043
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Univariate p values reported)4.2 Random intercepts and Random slopes
# 個人的截距隨機效果與時間隨機斜率效果
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.3118 -0.4989 -0.0696 0.3971 3.7190
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subject (Intercept) 50.7343 7.1228
## time 0.1844 0.4294 -0.18
## Residual 23.9945 4.8984
## Number of obs: 280, groups: subject, 97
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.61391 2.24541 2.500
## bdi.pre 0.64185 0.07784 8.245
## time -0.70214 0.15407 -4.557
## treatmentBtheB -2.37598 1.67129 -1.422
## drugYes -2.86942 1.72766 -1.661
## length>6m 0.14041 1.63854 0.086
##
## 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# 檢驗固定效果
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.61391 2.24541 2.500 0.0124 *
## bdi.pre == 0 0.64185 0.07784 8.245 2.22e-16 ***
## time == 0 -0.70214 0.15407 -4.557 5.18e-06 ***
## treatmentBtheB == 0 -2.37598 1.67129 -1.422 0.1551
## drugYes == 0 -2.86942 1.72766 -1.661 0.0967 .
## length>6m == 0 0.14041 1.63854 0.086 0.9317
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (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 1887.5 1916.6 -935.75 1871.5
## BtheB_lmer2 10 1891.0 1927.4 -935.52 1871.0 0.4542 2 0.79694.3.1 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)
4.4 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
Bates, D. (2005), “Fitting linear mixed models in R,” R News, 5,
27–30, URL http://CRAN.R-project.org/doc/Rnews/.
Bates, D. and Sarkar, D. (2014), lme4: Linear Mixed-Effects Models Using
S4 Classes, URL http://CRAN.R-project.org/package=lme4, R package
version 1.1-5.
Pinheiro, J. C. and Bates, D. M. (2000), Mixed-Effects Models in S and
SPLUS, New York, USA: Springer-Verlag.