# 下載HSAUR3的package，使用其中的BtheB的data
if (!require("HSAUR3")) install.packages("HSAUR3")
## 載入需要的套件：HSAUR3
## Warning: 套件 'HSAUR3' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：tools
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 NA
# 針對Treated as Usaual與Beat the Blues繪製盒形圖
# 可以發現，隨著月份的增長，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)
library("lme4")
## Warning: 套件 'lme4' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：Matrix
## Warning: 套件 'Matrix' 是用 R 版本 4.1.3 來建造的
BtheB_lmer1 <- lmer(bdi ~ bdi.pre + time + treatment + drug +
length + (1 | subject), data = BtheB_long,
REML = FALSE, na.action = na.omit)
BtheB_lmer2 <- lmer(bdi ~ bdi.pre + time + treatment + drug +
length + (time | subject), data = BtheB_long,
REML = FALSE, na.action = na.omit)
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.7969
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.157
if (!require("multcomp")) install.packages("multcomp")
## 載入需要的套件：multcomp
## Warning: 套件 'multcomp' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：mvtnorm
## 載入需要的套件：survival
## 載入需要的套件：TH.data
## Warning: 套件 'TH.data' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：MASS
##
## 載入套件：'TH.data'
## 下列物件被遮斷自 'package:MASS':
##
## geyser
## 下列物件被遮斷自 'package:HSAUR3':
##
## birds
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)
# 比較隨機效果模型與殘差模型
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)
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))
}
