For an example analysis, we collated a .csv file of data from a repeated measures lexical decision experiment, which can be downloaded here. The code I have used in running the analysis described in this post can be found here.
This file holds data on:
1. participant attributes;
2. word attributes;
3. participant responses to words in the lexical decision task.
In this post, I want to take a quick look at how you would actually run a mixed-effects model, using the lmer() function furnished by the lme4 package, written by Doug Bates (Bates, 2010 – appears to a preprint of his presumed new book, see also Pinheiro & Bates, 2000). A useful how-to guide, one I followed when I first started doing this kind of modelling, is in Baayen (2008). An excellent, and very highly cited, introduction to mixed-effects modelling for psycholinguistics, is found in Baayen, Davidson & Bates (2008).
First, then, we enter the following lines of code to load the libraries we are likely to need, to set the working directory to the folder where we have copied the full database, and to read that database as a dataframe in the session’s workspace:
Load libraries, define functions, load data
We start by loading the libraries of functions we are likely to need, defining a set of functions not in libraries, setting the working directory, and loading the dataframe we will work with into the R session workspace.
# load libraries of packages #############################################################################
# library(languageR)
library(lme4)
library(ggplot2)
library(grid)
library(rms)
library(plyr)
library(reshape2)
library(psych)
# library(gridExtra)
# load data files #######################################################################################
# set the working directory
setwd("/Users/robdavies/Downloads")
# work out where your files are and substitute for the path for my files location, above
# # read in the item norms data
#
# item.norms <- read.csv(file="item norms 270312 050513.csv", head=TRUE,sep=",", na.strings = "-999")
# word.norms <- read.csv(file="word norms 140812.csv", head=TRUE,sep=",", na.strings = "-999")
item.coding <- read.csv(file="item coding info 270312.csv", head=TRUE,sep=",", na.strings = "-999")
# read in the merged new dataframe
item.words.norms <- read.csv(file="item.words.norms 050513.csv", head=TRUE,sep=",", na.strings = "NA")
# read in the subjects scores data file
subjects <- read.csv("ML scores 080612 220413.csv", header=T, na.strings = "-999")
# read in the behavioural data file
behaviour <- read.csv("ML data 080612.csv", header=T, na.strings = "-999")
# notice that we collated the subjects, items and behavioural data into a single dataframe which we then output
# using write.csv as a .csv file, we can read that collated data file in here, to save time
# read in the full database
RTs.norms.subjects <- read.csv("RTs.norms.subjects 061113.csv", header=T, na.strings = "-999")
# read in the database of correct responses to words only
nomissing.full <- read.csv("RTs.norms.subjects nomissing.full 061113.csv", header=T, na.strings = "-999")
# examine the data ######################################################################################
summary(nomissing.full)
# note the presence of centred variables
If you run all that code, you will see the dataframes and functions appear listed in the workspace window of R-studio.
Run the models
What I will do next is perform three sets of mixed-effects model analyses. I will talk about the first set in this post.
In this first set, I shall compare models varying in the predictor variables included in the model specification. We shall say that models with more predictors are more complex than models with fewer predictors (a subset of the predictors specified for the more complex models).
I shall then compare models of differing in complexity to see if increased complexity is justified by improved fit to data.
This is in response to: 1. reading Baayen (2008) and Pinheiro & Bates (2000) 2. recommendations by a reviewer; though I note that the reviewer recommendations are, in their application, qualified in my mind by some cautionary comments made by Pinheiro & Bates (2000), and confess that any problems of application or understanding are surely expressions of my misunderstanding.
What I am concerned about is the level of complexity that can be justified in the ultimate model by the contribution to improved fit by the model of the observed data. Of course, I have collected a set of predictors guided by theory, custom and practice in my field. But, as you have seen, I might choose to draw on a large set of potential predictors. I know that one cannot – to take one option used widely by psycholinguistics researchers – simply inspect a correlation table, comparing predictors with the outcome variable, and exclude those predictors that do not correlate with the outcome. Harrell (2001) comments that that approach has a biasing effect on the likelihood of detecting an effect that is really there (consuming, ‘ghost degrees of freedom’, or something like that). So what should we do or how should we satisfy ourselves that we should or should not include one or more predictor variables in our models?
You could start complex and step backwards to the simplest model that appears justified by model fit to observed data. I will return to that option another time.
In this post, I will start with the simplest model and build complexity up. Simpler models can be said to be nested in more complex models, e.g. a simple model including just main effects vs. a more complex model also including interactions. I will add complexity by adding sets of predictors, and I will compare each proximate pair of simple model and more complex model using the Likelihood Ratio Test (Pinheiro & Bates, 2000; see, especially, chapter 2, pp. 75-).
We can compare models with the same random effects (e.g. of subjects or items) but varying in complexity with respect to the fixed effects, i.e. replicable or manipulated variables like item length or participant age, and, separately, models with the same fixed effects but varying random effects.
Note that in comparing models varying in random effects, one estimates those effects using restricted maximum likelihood: you will see the argument REML = TRUE in the lmer() function call, that is what REML means. Note that in comparing models varying in fixed effects one must estimate those effects using maximum likelihood: you will see the argument REML = FALSE in the lmer() function call, that is what FALSE means. See comments in Baayen (2008) and Pinheiro & Bates (2000); I will come back to those important terms to focus on understanding, another time, focusing in this post on just writing the code for running the models.
We can compare the more and less complex models using the Likelihood Ration Test (LRT) which is executed with the anova([simple model], [more complex model]) function call. A 
Note that I log transform the RTs but one might also consider examining the reciprocal RT.
# we can move to some multilevel modelling ############################################################### # we might start by examining the effects of predictors in a series of additions of sets of predictors # we can hold the random effects, of subjects and items, constant, while varying the fixed effects # we can compare the relative fit of models to data using the anova() likelihood ratio test # in this series of models, we will use ML fitting full.lmer0 <- lmer(log(RT) ~ # Trial.order + # # cAge + cTOWREW + cTOWRENW + cART + # # cLength + cOLD + cBG_Mean + cLgSUBTLCD + cAOA + cIMG + (1|subjectID) + (1|item_name), data = nomissing.full, REML = FALSE) print(full.lmer0,corr=F) # full.lmer1 <- lmer(log(RT) ~ Trial.order + # cAge + cTOWREW + cTOWRENW + cART + # # cLength + cOLD + cBG_Mean + cLgSUBTLCD + cAOA + cIMG + (1|subjectID) + (1|item_name), data = nomissing.full, REML = FALSE) print(full.lmer1,corr=F) # full.lmer2 <- lmer(log(RT) ~ Trial.order + cAge + cTOWREW + cTOWRENW + cART + # cLength + cOLD + cBG_Mean + cLgSUBTLCD + cAOA + cIMG + (1|subjectID) + (1|item_name), data = nomissing.full, REML = FALSE) print(full.lmer2,corr=F) # full.lmer3 <- lmer(log(RT) ~ Trial.order + cAge + cTOWREW + cTOWRENW + cART + cLength + cOLD + cBG_Mean + cLgSUBTLCD + cAOA + cIMG + (1|subjectID) + (1|item_name), data = nomissing.full, REML = FALSE) print(full.lmer3,corr=F) # full.lmer4 <- lmer(log(RT) ~ Trial.order + (cAge + cTOWREW + cTOWRENW + cART)* (cLength + cOLD + cBG_Mean + cLgSUBTLCD + cAOA + cIMG) + (1|subjectID) + (1|item_name), data = nomissing.full, REML = FALSE) print(full.lmer4,corr=F) # compare models of varying complexity anova(full.lmer0, full.lmer1) anova(full.lmer1, full.lmer2) anova(full.lmer2, full.lmer3) anova(full.lmer3, full.lmer4)
This code will first print out the summaries of each model, along the lines you have already seen. It will then give you this:
# > # compare models of varying complexity # > # > anova(full.lmer0, full.lmer1) # Data: nomissing.full # Models: # full.lmer0: log(RT) ~ (1 | subjectID) + (1 | item_name) # full.lmer1: log(RT) ~ Trial.order + (1 | subjectID) + (1 | item_name) # Df AIC BIC logLik Chisq Chi Df Pr(>Chisq) # full.lmer0 4 -977.58 -951.31 492.79 # full.lmer1 5 -975.83 -943.00 492.92 0.2545 1 0.6139 # > anova(full.lmer1, full.lmer2) # Data: nomissing.full # Models: # full.lmer1: log(RT) ~ Trial.order + (1 | subjectID) + (1 | item_name) # full.lmer2: log(RT) ~ Trial.order + cAge + cTOWREW + cTOWRENW + cART + (1 | # full.lmer2: subjectID) + (1 | item_name) # Df AIC BIC logLik Chisq Chi Df Pr(>Chisq) # full.lmer1 5 -975.83 -943.00 492.92 # full.lmer2 9 -977.09 -917.98 497.54 9.2541 4 0.05505 . # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # > anova(full.lmer2, full.lmer3) # Data: nomissing.full # Models: # full.lmer2: log(RT) ~ Trial.order + cAge + cTOWREW + cTOWRENW + cART + (1 | # full.lmer2: subjectID) + (1 | item_name) # full.lmer3: log(RT) ~ Trial.order + cAge + cTOWREW + cTOWRENW + cART + cLength + # full.lmer3: cOLD + cBG_Mean + cLgSUBTLCD + cAOA + cIMG + (1 | subjectID) + # full.lmer3: (1 | item_name) # Df AIC BIC logLik Chisq Chi Df Pr(>Chisq) # full.lmer2 9 -977.09 -917.98 497.54 # full.lmer3 15 -1115.42 -1016.91 572.71 150.34 6 < 2.2e-16 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # > anova(full.lmer3, full.lmer4) # Data: nomissing.full # Models: # full.lmer3: log(RT) ~ Trial.order + cAge + cTOWREW + cTOWRENW + cART + cLength + # full.lmer3: cOLD + cBG_Mean + cLgSUBTLCD + cAOA + cIMG + (1 | subjectID) + # full.lmer3: (1 | item_name) # full.lmer4: log(RT) ~ Trial.order + (cAge + cTOWREW + cTOWRENW + cART) * # full.lmer4: (cLength + cOLD + cBG_Mean + cLgSUBTLCD + cAOA + cIMG) + # full.lmer4: (1 | subjectID) + (1 | item_name) # Df AIC BIC logLik Chisq Chi Df Pr(>Chisq) # full.lmer3 15 -1115.4 -1016.91 572.71 # full.lmer4 39 -1125.7 -869.56 601.84 58.262 24 0.0001119 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # it looks like the full model, including interactions but not the trialorder variable, is justified
Notice:
— It is not very pretty but it does tell you what predictors are included in each model you are comparing, for each pair of models being compared, the outcome of the ratio test, giving the numbers I usually report: the 
— I am looking at this output to see if, for each pair of models, the more complex model fits the data significantly better (if p =< 0.05).
–I add predictors in blocks and assess the improvement of fit associated with the whole block.
I have made a lot of promises to explain things more another time. So I will add one more. Pinheiro & Bates (2000; p. 88) argue that Likelihood Ratio Test comparisons of models varying in fixed effects tend to be anticonservative i.e. will see you observe significant differences in model fit more often than you should. I think they are talking, especially, about situations in which the number of model parameter differences (differences between the complex model and the nested simpler model) is large relative to the number of observations. This is not really a worry for this dataset, but I will come back to the substance of this view, and alternatives to the approach taken here.
What next?
What I would seek to do in the next series of analyses is determine if the effects that are significant appear to be significant due:
1. to the presence of other, nonsignificant effects — making the check by removing the nonsignificant effects then rerunning the model;
2. to the presence of outlier observations — making the check by removing observations associated with large model residuals then rerunning the model.
What have we learnt?
(The code I have used in running the analysis described in this post can be found here.)
Firstly, I am going to have to do a lot more explaining, for my own sake, and have promised to do so for a wide range of ideas, including those concerning the methods of model fitting, and, especially, the discussion around how one determines the appropriate set of fixed effects.
Secondly, that the model code is both simple and – on its face – reasonably intuitive.
Key vocabulary (mostly, a list of things to be explained)
fixed effects
random effects
random intercepts
random slopes
interactions
multiplicative products
centering
maximum likelihood
restricted maximum likelihood
Markov chain Monte Carlo
Likelihood Ratio Test
lmer()
anova()
Reading
Baayen, R. H. (2008). Analyzing linguistic data. Cambridge University Press.
Bates, D. M. (2005). Fitting linear mixed models in R. R News, 5, 27-30.
Belsley, D. A., Kuh, E., & Welsch, R. E. (1980). Regression diagnostics: Identifying influential data and sources of collinearity. New York: Wiley.
Cohen, J., Cohen, P., West, S. G., & Aiken, L. S. (2003). Applied multiple regression/correlation analysis for the behavioural sciences (3rd. edition). Mahwah, NJ: Lawrence Erlbaum Associates.
Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models. New York,NY: Cambridge University Press.
Pinheiro, J. C., & Bates, D. M. (2000). Mixed-effects models in S and S-plus (Statistics and Computing). New York: Springer.
ayeimanolr 