18 subjects (3hrs sleep) 10 days observations on each

Note. Comments in English

1 Data Management

Load and view data structure

data(sleepstudy, package="lme4")  
dim(sleepstudy)  
## [1] 180   3
head(sleepstudy)
##   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
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

2 Visdualization 01

Display the relationship between Reaction and Sleep Deprivation per subject

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

3 Analysis

library(lme4)
## Loading required package: Matrix
sleeplmList <- lmList(Reaction ~ Days |Subject, data = sleepstudy)

Compute the mean of coefficients for variale 1 & 2 (Reaction time and Days)

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

mean(coef(sleeplmList)[,1]) 
## [1] 251.4051
mean(coef(sleeplmList)[,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

Compute the variance of coefficients for variale 1 & 2 (Reaction time and Days)

# - SSR/SST
var(coef(sleeplmList)[,1]) 
## [1] 838.3423
var(coef(sleeplmList)[,2]) 
## [1] 43.01034

Compute the quantile of coefficients for variale 1 & 2 (Reaction time and Days)

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

Compute stem & leaf plot for the coefficients (variable 2= Days )

#     
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

Compute correlations for the coefficients of variable 1 & 2 (Reaction, Days)

cor(coef(sleeplmList)[,1], coef(sleeplmList)[,2])  
## [1] -0.1375534

Comparing models with different “Random effects”, REML True

# first 2 do REML, second pair match 
# Bates p.67 in doing MLE (REML FALSE)

sleeplmer <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
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

Comparing models with different “Fixed effects”, REML False

sleeplmer2 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, 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

4 Visdualization 02

Plot the relationship between Deprivation and Reaction Time per subject

library(ggplot2)
 ggplot(sleepstudy, aes(Days, Reaction)) +
 geom_smooth(method = "lm", se = F, lwd = .6) +
 geom_point(size = rel(0.5)) +
 facet_wrap(~ Subject, nrow = 1) +
 labs(x = "Deprivation(in days)", y="Reaction time")+
 theme(axis.text.x = element_text(size = 6))
## `geom_smooth()` using formula 'y ~ x'

Plot the relationship between Deprivation and Reaction Time, using quardratic x to predict y (specifying the distribution points of Reaction time) per day

ggplot(sleepstudy, aes(Days, Reaction))+
 geom_point()+
 stat_smooth(method = "lm", formula = y ~ poly(x, 2)) +
 labs(x = "Deprivation (in days)", y = "Reaction time")+
 theme(legend.position = c(.9, .2))

The end ```