To study the impact of sleep deprivation on cognitive performance, a sample of 18 truck drivers were allowed only 3 hours of sleep during 10 consecutive days and their reaction times were measured several times each day.

Source: Belenky, G., Wesensten, N.J., Thorne, D.R., Thomas, M.L., Sing, H.C., Redmond, D.P., Russo, M.B., & Balkin, T.J. (2003). Patterns of performance degradation and restoration during sleep restriction and subsequent recovery: a sleep does response study. Journal of Sleep Research, 12, 1-12. Reported in Pinheiro, J.C. & Bates, D.M. (2000). Mixed effects models in S and S-PLUS.

Data: sleepstudy{lme4}

Column 1: Mean reaction time in milliseconds Column 2: Days of sleep deprivation Column 3: Truck driver ID

1 Data management

library(lme4) 
## Loading required package: Matrix

2 18 subjects (3hrs sleep) 10 observations on each, p64 of bates

3 the “new” random effects (vs nlme) see Bates book 2010

4 sleep deprivation 3hrs/night truckers

5 

data(sleepstudy, package="lme4") #load data

6 

#nrow( )觀察列的數量、ncol( )觀察行的數量、dim( )可同時看行列的數量

dim(sleepstudy) 
## [1] 180   3

7 

head(sleepstudy,10)
##    Reaction Days Subject
## 1  249.5600    0     308
## 2  258.7047    1     308
## 3  250.8006    2     308
## 4  321.4398    3     308
## 5  356.8519    4     308
## 6  414.6901    5     308
## 7  382.2038    6     308
## 8  290.1486    7     308
## 9  430.5853    8     308
## 10 466.3535    9     308

8 

with(sleepstudy, table(Subject)) #依照受試者表列出睡眠時間有幾筆觀察值
## Subject
## 308 309 310 330 331 332 333 334 335 337 349 350 351 352 369 370 371 372 
##  10  10  10  10  10  10  10  10  10  10  10  10  10  10  10  10  10  10

9 Visdualization

lattice::xyplot(Reaction ~ Days | Subject, data=sleepstudy, 
       xlab="Days of sleep deprivation",
       ylab="Average reaction time (ms)",
       type=c("g","p","r"))

10 Analysis

sleeplmList <- lmList(Reaction ~ Days |Subject, data = sleepstudy) #將data依據每個受試者「剝奪睡眠的天數Days」與「平均反應時間Reaction」之間的關係進行分區(?)

sleeplmList

#mean int and slope match lmer Fixed effects results p.67

mean(coef(sleeplmList)[,1]) #計算變數1的平均數
## [1] 251.4051
mean(coef(sleeplmList)[,2]) #計算變數2的平均數
## [1] 10.46729

#these are too big as they should be, #compare with variance random effects p.67 # mle subtracts off wobble in estimated indiv regressions Y on t

11 - SSR/SST

var(coef(sleeplmList)[,1]) #計算變數1的variance
## [1] 838.3423
var(coef(sleeplmList)[,2]) #計算變數2的variance
## [1] 43.01034

12 

quantile(coef(sleeplmList)[,1]) #計算四分位數
##       0%      25%      50%      75%     100% 
## 203.4842 229.4167 258.0576 273.0255 290.1041
quantile(coef(sleeplmList)[,2]) 
##        0%       25%       50%       75%      100% 
## -2.881034  6.194548 10.432421 13.548395 21.764702

13 

stem(coef(sleeplmList)[,2]) #莖葉圖
## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##   -0 | 3
##   0 | 23
##   0 | 56699
##   1 | 011234
##   1 | 89
##   2 | 02

14 

#計算Reaction和Days的相關係數
cor(coef(sleeplmList)[,1], coef(sleeplmList)[,2])  
## [1] -0.1375534

15 first 2 do REML, second pair match

16 Bates p.67 in doing MLE (REML FALSE)

17 

sleeplmer <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) #每個受試者的Days具有隨機截距項
summary(sleeplmer)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (Days | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1743.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9536 -0.4634  0.0231  0.4634  5.1793 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  Subject  (Intercept) 612.10   24.741       
##           Days         35.07    5.922   0.07
##  Residual             654.94   25.592       
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  251.405      6.825  36.838
## Days          10.467      1.546   6.771
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.138

18 

sleeplmer2 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, REML=FALSE) #每個受試者的Days具有隨機截距項，在比較具有不同“固定效果”的模型時使用REML = FALSE
summary(sleeplmer2)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Reaction ~ Days + (Days | Subject)
##    Data: sleepstudy
## 
##      AIC      BIC   logLik deviance df.resid 
##   1763.9   1783.1   -876.0   1751.9      174 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9416 -0.4656  0.0289  0.4636  5.1793 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  Subject  (Intercept) 565.48   23.780       
##           Days         32.68    5.717   0.08
##  Residual             654.95   25.592       
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  251.405      6.632  37.907
## Days          10.467      1.502   6.968
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.138

19 The end