MRRI MLR Demo Examples

Dan Mirman

These are examples of two different kinds of multilevel regression (MLR) analyses that I presented as part of a very basic introduction to MLR at MRRI in January, 2015. The main point of these examples was just to give people a sense of what running a MLR in R involves, not as full educational examples.

Example 1: Longitudinal data

Read in data based on a study of effects of transitional probability (TP) on novel word learning (Mirman et al., 2008):

##     Subject         TP          Block         Accuracy        Correct    
##  Min.   :244.0   High:280   Min.   : 1.0   Min.   :0.000   Min.   :0.00  
##  1st Qu.:319.5   Low :280   1st Qu.: 3.0   1st Qu.:0.667   1st Qu.:4.00  
##  Median :346.5              Median : 5.5   Median :0.833   Median :5.00  
##  Mean   :344.4              Mean   : 5.5   Mean   :0.805   Mean   :4.83  
##  3rd Qu.:370.5              3rd Qu.: 8.0   3rd Qu.:1.000   3rd Qu.:6.00  
##  Max.   :395.0              Max.   :10.0   Max.   :1.000   Max.   :6.00  
##    Incorrect   
##  Min.   :0.00  
##  1st Qu.:0.00  
##  Median :1.00  
##  Mean   :1.17  
##  3rd Qu.:2.00  
##  Max.   :6.00
WordLearnEx <- read.delim("http://www.danmirman.org/gca/WordLearnEx.txt")
summary(WordLearnEx)

Here's a graph of those data, made using the ggplot2 package:

library(ggplot2)
ggplot(WordLearnEx, aes(Block, Accuracy, color=TP)) + 
  stat_summary(fun.data=mean_se, geom="pointrange")

plot of chunk unnamed-chunk-3

We're going to analyze these data using a second-order (quadratic) orthogonal polynomial, so let's build that and add to the data frame

#make orthogonal polynomial using the Block values
t <- poly(1:max(WordLearnEx$Block), 2)
t
##                 1           2
##  [1,] -0.49543369  0.52223297
##  [2,] -0.38533732  0.17407766
##  [3,] -0.27524094 -0.08703883
##  [4,] -0.16514456 -0.26111648
##  [5,] -0.05504819 -0.34815531
##  [6,]  0.05504819 -0.34815531
##  [7,]  0.16514456 -0.26111648
##  [8,]  0.27524094 -0.08703883
##  [9,]  0.38533732  0.17407766
## [10,]  0.49543369  0.52223297
## attr(,"degree")
## [1] 1 2
## attr(,"coefs")
## attr(,"coefs")$alpha
## [1] 5.5 5.5
## 
## attr(,"coefs")$norm2
## [1]   1.0  10.0  82.5 528.0
## 
## attr(,"class")
## [1] "poly"   "matrix"
#add it into the data frame
WordLearnEx[, paste("ot", 1:2, sep="")] <- t[WordLearnEx$Block, 1:2]
#re-check data
summary(WordLearnEx)
##     Subject         TP          Block         Accuracy        Correct    
##  Min.   :244.0   High:280   Min.   : 1.0   Min.   :0.000   Min.   :0.00  
##  1st Qu.:319.5   Low :280   1st Qu.: 3.0   1st Qu.:0.667   1st Qu.:4.00  
##  Median :346.5              Median : 5.5   Median :0.833   Median :5.00  
##  Mean   :344.4              Mean   : 5.5   Mean   :0.805   Mean   :4.83  
##  3rd Qu.:370.5              3rd Qu.: 8.0   3rd Qu.:1.000   3rd Qu.:6.00  
##  Max.   :395.0              Max.   :10.0   Max.   :1.000   Max.   :6.00  
##    Incorrect         ot1               ot2          
##  Min.   :0.00   Min.   :-0.4954   Min.   :-0.34816  
##  1st Qu.:0.00   1st Qu.:-0.2752   1st Qu.:-0.26112  
##  Median :1.00   Median : 0.0000   Median :-0.08704  
##  Mean   :1.17   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.:2.00   3rd Qu.: 0.2752   3rd Qu.: 0.17408  
##  Max.   :6.00   Max.   : 0.4954   Max.   : 0.52223

We'll start by fitting a base model that has just the time terms (i.e., a model of the overall learning curve), then add the TP condition terms.

m.base <- lmer(Accuracy ~ ot1 + ot2 + 
                  (ot1 + ot2 | Subject), data=WordLearnEx, REML=FALSE)
m.0 <- lmer(Accuracy ~ ot1 + ot2 + TP +
                  (ot1 + ot2 | Subject), data=WordLearnEx, REML=FALSE)
m.1 <- lmer(Accuracy ~ ot1 + ot2 + TP + ot1:TP +
               (ot1 + ot2 | Subject), data=WordLearnEx, REML=FALSE)
m.2 <- lmer(Accuracy ~ ot1 + ot2 + 
               TP + ot1:TP + ot2:TP +
               (ot1 + ot2 | Subject), data=WordLearnEx, REML=FALSE)
anova(m.base, m.0, m.1, m.2)
## Data: WordLearnEx
## Models:
## m.base: Accuracy ~ ot1 + ot2 + (ot1 + ot2 | Subject)
## m.0: Accuracy ~ ot1 + ot2 + TP + (ot1 + ot2 | Subject)
## m.1: Accuracy ~ ot1 + ot2 + TP + ot1:TP + (ot1 + ot2 | Subject)
## m.2: Accuracy ~ ot1 + ot2 + TP + ot1:TP + ot2:TP + (ot1 + ot2 | Subject)
##        Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m.base 10 -330.77 -287.50 175.39  -350.77                           
## m.0    11 -330.32 -282.72 176.16  -352.32 1.5506      1    0.21304  
## m.1    12 -328.68 -276.75 176.34  -352.68 0.3576      1    0.54984  
## m.2    13 -332.63 -276.37 179.32  -358.63 5.9506      1    0.01471 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Looks like adding the TP effect on the quadratic term significantly improved model fit. In other words, there was a significant effect of TP on the (quadratic) curvature of the learning curves. We can get more detailed output for the full model, including parameter estimates, etc.

summary(m.2)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: 
## Accuracy ~ ot1 + ot2 + TP + ot1:TP + ot2:TP + (ot1 + ot2 | Subject)
##    Data: WordLearnEx
## 
##      AIC      BIC   logLik deviance df.resid 
##   -332.6   -276.4    179.3   -358.6      547 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6182 -0.5358  0.1263  0.5665  2.6157 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  Subject  (Intercept) 0.01076  0.10374             
##           ot1         0.01542  0.12419  -0.33      
##           ot2         0.00628  0.07925  -0.28 -0.82
##  Residual             0.02456  0.15672             
## Number of obs: 560, groups:  Subject, 56
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.831486   0.021727   38.27
## ot1          0.287391   0.037788    7.61
## ot2         -0.167304   0.033188   -5.04
## TPLow       -0.052961   0.030727   -1.72
## ot1:TPLow   -0.001075   0.053441   -0.02
## ot2:TPLow    0.116455   0.046935    2.48
## 
## Correlation of Fixed Effects:
##           (Intr) ot1    ot2    TPLow  o1:TPL
## ot1       -0.183                            
## ot2       -0.114 -0.229                     
## TPLow     -0.707  0.129  0.081              
## ot1:TPLow  0.129 -0.707  0.162 -0.183       
## ot2:TPLow  0.081  0.162 -0.707 -0.114 -0.229

and we can plot the model fit

ggplot(WordLearnEx, aes(Block, Accuracy, color=TP)) + 
  stat_summary(fun.data=mean_se, geom="pointrange") +
  stat_summary(aes(y=fitted(m.2)), fun.y=mean, geom="line")

plot of chunk unnamed-chunk-7

Example 2: Crossed random effects of subjects and items

These data are from a pilot study of phonological neighborhood density effects in younger and older adults.

##     Subject         Group           Item      Density           RT       
##  S103   : 115   Older  :1494   bag    :  30   High:1397   Min.   :  440  
##  S122   : 109   Younger:1370   blonde :  30   Low :1467   1st Qu.: 1161  
##  S105   : 107                  chain  :  30               Median : 1536  
##  S102   : 105                  chair  :  30               Mean   : 1730  
##  S108   : 102                  child  :  30               3rd Qu.: 2040  
##  S110   : 102                  couch  :  30               Max.   :12753  
##  (Other):2224                  (Other):2684
density <- read.delim("http://www.danmirman.org/gca/Density.txt")
summary(density)
ggplot(density, aes(Group, RT, fill=Density)) + geom_violin() + 
  scale_y_log10(breaks=c(400, 800, 1600, 3200, 6400, 12800))

plot of chunk unnamed-chunk-10

These are trial-level data (error trials have already been removed). The old-fashioned thing do would be to run separate by-subject and by-item ANOVAs. To do that, we need to compute by-subject and by-item averages (using the plyr package)

library(plyr)
#by-subject means
density.subj <- ddply(density, .(Subject, Group, Density), summarize, RT = mean(RT))
summary(density.subj)
##     Subject       Group    Density         RT      
##  S101   : 2   Older  :32   High:30   Min.   :1027  
##  S102   : 2   Younger:28   Low :30   1st Qu.:1383  
##  S103   : 2                          Median :1603  
##  S104   : 2                          Mean   :1743  
##  S105   : 2                          3rd Qu.:1960  
##  S106   : 2                          Max.   :3601  
##  (Other):48
#by-item means
density.item <- ddply(density, .(Item, Density, Group), summarize, RT = mean(RT))
summary(density.item)
##       Item     Density        Group           RT        
##  bag    :  2   High:120   Older  :120   Min.   : 942.2  
##  ball   :  2   Low :120   Younger:120   1st Qu.:1366.0  
##  band   :  2                            Median :1686.4  
##  bean   :  2                            Mean   :1767.4  
##  bear   :  2                            3rd Qu.:2040.0  
##  bell   :  2                            Max.   :3866.8  
##  (Other):228

Now we can run the ANOVAs (repeated measures and mixed within-between ANOVAs are easier to run using the ez package)

library(ez)
## Warning: package 'ez' was built under R version 3.1.2
#by-subjects mixed ANOVA with Density as within-subject variable and Group as between-subjects variable
ezANOVA(density.subj, dv=RT, wid=Subject, within=Density, between=Group, type=3)
## Warning: Data is unbalanced (unequal N per group). Make sure you specified
## a well-considered value for the type argument to ezANOVA().
## $ANOVA
##          Effect DFn DFd         F            p p<.05         ges
## 2         Group   1  28 17.470709 0.0002590132     * 0.374070220
## 3       Density   1  28 12.455240 0.0014610742     * 0.018425767
## 4 Group:Density   1  28  5.256639 0.0295826013     * 0.007860159
#by-items mixed ANOVA with Group as within-item variable and Density as between-items variable
ezANOVA(density.item, dv=RT, wid=Item, within=Group, between=Density, type=3)
## $ANOVA
##          Effect DFn DFd          F            p p<.05         ges
## 2       Density   1 118   2.437651 1.211315e-01       0.014696412
## 3         Group   1 118 223.461729 5.359396e-29     * 0.344870084
## 4 Density:Group   1 118   1.964657 1.636411e-01       0.004606873

Here's how we could analyze these data using MLR with crossed random effects of subjects and items

m <- lmer(RT ~ Density*Group + (Density | Subject) + (Group | Item), data=density, 
          contrasts = list(Density = "contr.sum", Group = "contr.sum"))

There are no p-values due to some ambiguity about how the degrees of freedom should be computed, but for reasonably large data sets, we can just use a normal approximation

coefs <- data.frame(coef(summary(m)))
coefs$p.normal <- 2*(1-pnorm(abs(coefs$t.value)))
coefs
##                   Estimate Std..Error   t.value     p.normal
## (Intercept)     1753.62267   79.49054 22.060772 0.000000e+00
## Density1          58.35157   29.62191  1.969879 4.885224e-02
## Group1           325.28831   75.89407  4.286083 1.818507e-05
## Density1:Group1   36.15172   17.85479  2.024763 4.289168e-02

In the MLR model, each of the effects is statistically significant. We can also use the model to get more reasonable estimates for making a graph (using the effects package)

library(effects)
## Loading required package: lattice
## Loading required package: grid
## Loading required package: colorspace
ef <- as.data.frame(effect("Density:Group", m))
ef
##   Density   Group      fit       se    lower    upper
## 1    High   Older 2173.414 123.9369 1930.400 2416.429
## 2     Low   Older 1984.408 106.8746 1774.849 2193.967
## 3    High Younger 1450.534 124.1431 1207.115 1693.953
## 4     Low Younger 1406.135 104.4298 1201.369 1610.900
ggplot(ef, aes(Group, fit, fill=Density)) + 
  geom_bar(position="dodge", stat="identity") +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), 
                width=0.4, position=position_dodge(width=0.9)) +
  ylab("RT (ms)")

plot of chunk unnamed-chunk-15