Class 11B

Repeated Measures

Study on sleep deprivation and reaction times.

This data comes from a sleep study in the lme4 package.

Let’s first visualize the data

library(lme4)
data(sleepstudy)
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

library(tidyverse)
ggplot(sleepstudy, aes(Days, Reaction))+
  geom_point()+
  geom_smooth(method="lm")+
  facet_wrap(~Subject)

It appears that both the slope and intercept vary across subjects

# model with random slope
mod<-lmer(Reaction~Days+(1|Subject), data = sleepstudy)
summary(mod)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1786.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2257 -0.5529  0.0109  0.5188  4.2506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 1378.2   37.12   
##  Residual              960.5   30.99   
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 251.4051     9.7467   25.79
## Days         10.4673     0.8042   13.02
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.371

# model with random slope and intercept
mod1<-lmer(Reaction~Days+(Days|Subject), data = sleepstudy)
summary(mod1)

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

Time Series and Correlation Structure

Data from calorie intake study and nutrition curriculum treatments over time.

calories<-read.csv("/Users/heatherkitada/Desktop/MATH138/calories.csv", 
                   header=TRUE)

Find the AR(1) parameter (using the nlme function)

library(nlme)

# correlation structure
# gls 
model.a = gls(Calories.per.day ~ Instruction + Month + Instruction*Month,
              data=calories)

ACF(model.a,
    form = ~ Month | Student)

##   lag       ACF
## 1   0 1.0000000
## 2   1 0.8989822
## 3   2 0.7462712
## 4   3 0.6249217
## 5   4 0.5365430
## 6   5 0.4564673

# lme
model.b = lme(Calories.per.day ~ Instruction + Month + Instruction*Month,
              random = ~1|Student,
              data=calories)

ACF(model.b)

##   lag        ACF
## 1   0  1.0000000
## 2   1  0.4286522
## 3   2 -0.1750233
## 4   3 -0.6283754
## 5   4 -0.8651088
## 6   5 -0.6957774

Fitting the model

modelGLS = gls(Calories.per.day ~ Instruction + Month + Instruction*Month,
            correlation = corAR1(form = ~ Month | Student,
                                 value = 0.8990),
            data=calories,
            method="REML")

modelLME = lme(Calories.per.day ~ Instruction + Month + Instruction*Month,
            random = ~1|Student,
            correlation = corAR1(form = ~ Month | Student,
                                 value = 0.4287),
            data=calories,
            method="REML")

Forcing R To Do Inference

# force R to do inference
library(car)

Anova(modelGLS)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Calories.per.day
##                   Df   Chisq Pr(>Chisq)    
## Instruction        2 10.4221   0.005456 ** 
## Month              1  0.0198   0.888045    
## Instruction:Month  2 38.6045  4.141e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(modelLME)

## Analysis of Deviance Table (Type II tests)
## 
## Response: Calories.per.day
##                     Chisq Df Pr(>Chisq)    
## Instruction       10.4221  2   0.005456 ** 
## Month              0.0198  1   0.888045    
## Instruction:Month 38.6045  2  4.141e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testing for Significance of Random Effects

Comparing a model with random effects to a model without!

# testing random effects
# use gls when we only have fixed effects
model.fixed = gls(Calories.per.day ~ Instruction + Month + Instruction*Month,
                  data=calories,
                  method="REML")

anova(modelGLS,
      model.fixed)

##             Model df      AIC      BIC    logLik   Test L.Ratio p-value
## modelGLS        1  8 714.9693 732.4865 -349.4847                       
## model.fixed     2  7 813.6213 828.9489 -399.8106 1 vs 2 100.652  <.0001