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.
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")
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")
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))
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)")
