Fitting GLMMs with lme4

This page it is NOT meant to be a tutorial for Mixed Models – if you need such a tutorial, you may read our article (Moscatelli, Mezzetti and Lacquaniti, 2012), or one of the other books and tutorials suggested in our Link Page. A comprehensive and detailed article on the package lme4, written by the package’s authors and maintainers, is Bates et al (2014) , now available on arxiv.

The main purpose of this page is to provide a quick introduction to the syntax of the lme4 package, in order to fit a GLMM to the data. Load the MERpsychophysics package before running the code. If you installed the lme4 package, it will load automatically when loading MERpsychophysics. Examples based on lme4.1 and MERpsychophysics.1.

Here I assume that you have already imported your data set in the R workspace. Otherwise, you can simulate a data set using the  function MERsimulate (all MER- function are not from lme4 package, but from our MERpsychophysics package). For example, the line below generates a data frame named “datafr”:

datafr = MERsimulate(nsubjects = 10)

It is possible to change the number of subjects, trials, stimulus values etc. (see the Reference Manual for details). Check the structure of the data frame with the R function str(). The variable datafr$X contains the continuous predictor (i.e., the values of the experimental stimulus), while the variables datafr$Longer and datafr$Total contain the response of the subject. We assume a Binomial response, defined by the number of Yes (or Longer, in our case) responses and the total number of trials. Now we are ready to fit the model.

#fit the GLMM(probit link function):: one random effect
mod1 = glmer(formula = cbind(Longer, Total - Longer) ~ X + (1 | Subject),
family = binomial(link = "probit"), data = datafr)

The fitted model is named as mod1. We can check the model using the function summary(). Let’s have a closer look at the syntax. We use the glmer function, from package lme4, in order to fit the model. Type:


in the R workspace to open the help page of the function.

The formula argument contain the model formula. The tilde “~” divides the formula in two sides, the binomial response variable on the left side, and the predictors on the right side. In Mixed Model framework, we define two types of predictor, the fixed effects and the random effects predictor. In our example, the values of the experimental stimulus X is the fixed-effect predictor, while the term (1 | Subject) is the random-effect predictor. The random-effect predictor is defined between parentheses, in this case we wrote 1 to include only a random intercept in the model. The “|” character divides the random effect from the grouping factor. To clarify the issue: The model assumes that in each different subjects, the response has a different distribution. In this first model, we assume that the response distribution of different subjects changes ONLY for the intercept value – that is, all the underling psychometric functions of our 10 subjects are shifted one to the other, but they have the same slope.

If this is not the case, and we assume that different subjects have different slopes? Then, we can fit a second model including 3 random effects: The random intercept, the random slope, and the correlation of the two. The syntax is here:

#fit the GLMM(probit link function)::two random effects
mod2 = glmer(formula = cbind(Longer, Total - Longer) ~ X + (1 + X| Subject),
family = binomial(link = "probit"), data = datafr)

The family argument specifies that our (conditional) response follows a binomial distribution, and that we are using a Probit link function. This is the function you have to use, if you want your sigmoid function to be a cumulative Gaussian distribution. The argument data specifies the data frame where we stored our variables.

Now we compare mod1 and mod2, in order to decide which of the two gives the best fitting of the data. We can compare the AIC of the two models – the model with the lowest AIC provides the best fitting to the data. Otherwise, we can use the likelihood ratio test (R function anova) to compare the two models.

#compare the two models: AIC, BIC, LR test
anova(mod1, mod2, test = "Chisq")

Once we have chose our model, we can use the functions of the MERpsychophysics in order to estimate the PSE and JND.

#Estimate the PSE with bootstrap method or delta method
#Increase the size to make the estimate more reliable
#(although, this will increase the duration of the bootstrap)
boot.mod2 = pseMer(mod2, B = 200)

#Using delta method:
xplode.mod2 = xplode.mer(model = mod2, name.cont = "X")

This is for a simple GLMM having a single fixed-effect, continuous predictor. In most of psychophysical experiment, we manipulate both continuous (e.g., time, length, force, etc.) and factorial predictors (e.g. gender, handiness, etc.). These type of date can be modelled using a multivariable GLMM. Next page introduces the basic R syntax to fit multivariable GLMM.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s