# 下載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
# layout()可以很容易組合多副圖爲一副總括圖形。
# grep與grepl最大的差異，在於它回傳的是「符合pattern」的「位置」或「文字」，而不是一個布林向量。 因此結果的向量長度，會短很多。
# 針對Treated as Usaual&Beat the Blues 畫盒形圖
# 可以看出隨著月份增長，BDI越低
layout(matrix(1:2, nrow = 1))
ylim <- range(BtheB[, grep("bdi", names(BtheB))], na.rm = T)
tau <- subset(BtheB, treatment == "TAU")[,
grep("bdi", names(BtheB))]
boxplot(tau, main = "Treated as Usaual",
ylab = "BDI", xlab = "Time(on 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)
#比較兩種差別，(1 | subject)為subject整體的隨機效果，(time | subject)為subject時間的隨機效果
#結果發現，BtheB_lmer2更加適配
library(lme4)
## Warning: 套件 'lme4' 是用 R 版本 4.1.3 來建造的
## 載入需要的套件：Matrix
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
# A convenience function for univariate testing via z- and t-tests of estimated model coefficients
# 結果發現bdi.pre、time == 0時，模型顯著
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)
# 比繳隨機效果模型與Residuals模型
# 兩者都可看出，其為正相關(?)或常態分佈(?)
layout(matrix(1:2,ncol=2))
qint<-ranef(BtheB_lmer1)$subject[["(Intercept)"]]
qres<-residuals(BtheB_lmer1)
qqnorm(qint,ylab="Estimatedrandomintercepts", xlim=c(-3,3),ylim=c(-20,20), main="Randomintercepts")
qqline(qint)
qqnorm(qres,xlim=c(-3,3),ylim=c(-20,20), ylab="Estimatedresiduals", 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)) }
