Inclass Exercise2

1 load and see

#load package
library(lme4)

## Loading required package: Matrix

18 subjects (3hrs sleep) 10 observations on each, p64 of bates the “new” random effects (vs nlme) see Bates book 2010 sleep deprivation 3hrs/night truckers

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

#dimension
dim(sleepstudy) 

## [1] 180   3

#see it
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

#show subject table
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 lattice xyplot

#lattice xyplot
lattice::xyplot(Reaction ~ Days | Subject, data=sleepstudy, 
       xlab="Days of sleep deprivation",
       ylab="Average reaction time (ms)",
       type=c("g","p","r"),
       layout = c(6,3)) #以6x3方式呈現

大部分睡眠天數越少，反應時間增加

3 linear model

sleeplmList <- lmList(Reaction ~ Days |Subject, data = sleepstudy)

#
sleeplmList    

## Call: lmList(formula = Reaction ~ Days | Subject, data = sleepstudy) 
## Coefficients:
##     (Intercept)      Days
## 308    244.1927 21.764702
## 309    205.0549  2.261785
## 310    203.4842  6.114899
## 330    289.6851  3.008073
## 331    285.7390  5.266019
## 332    264.2516  9.566768
## 333    275.0191  9.142045
## 334    240.1629 12.253141
## 335    263.0347 -2.881034
## 337    290.1041 19.025974
## 349    215.1118 13.493933
## 350    225.8346 19.504017
## 351    261.1470  6.433498
## 352    276.3721 13.566549
## 369    254.9681 11.348109
## 370    210.4491 18.056151
## 371    253.6360  9.188445
## 372    267.0448 11.298073
## 
## Degrees of freedom: 180 total; 144 residual
## Residual standard error: 25.59182

#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

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

## [1] 838.3423

var(coef(sleeplmList)[,2]) 

## [1] 43.01034

#quantile of  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

# stem leaf plot (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

#correlation coefficient
cor(coef(sleeplmList)[,1], coef(sleeplmList)[,2])  

## [1] -0.1375534

4 Comparing models

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

#lmer
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.4633  5.1793 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  Subject  (Intercept) 611.90   24.737       
##           Days         35.08    5.923   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.824  36.843
## Days          10.467      1.546   6.771
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.138

#
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

5 Reaction time and day plot

#轉換順序
sleepstudy$Subject = factor(sleepstudy$Subject, levels=c("309", "310", "335", "349", "351", "370", "371", "334", "330", "369", "332", "331", "350", "333", "372", "352", "308", "337"))
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) #一列排序 +

## `geom_smooth()` using formula 'y ~ x'

 labs(x = "Deprivation(in days)", y="Reaction time")+
 theme(axis.text.x = element_text(size = 6))

## NULL

ID 309的駕駛反應時間較少

ggplot(sleepstudy, aes(Days, Reaction)) +
 geom_point() +
 stat_smooth(method = "lm", aes(Days, Reaction), se = TRUE) +
 labs(x = "Deprivation (in Days)", 
      y = "Reaction time (in ms)") +
 theme_bw()

## `geom_smooth()` using formula 'y ~ x'