Multivariable GLMMs

What is a multivariable model? Consider a classical experiment paradigm psychophysics. Usually, the experimenter manipulates a continuous stimulus (for example time, velocity, weight, etc.) and categorical or factorial predictors (gender, handiness, animacy, etc.). A multivariable model relates the response variable (which in psychophysics has usually a Binomial distribution, i.e. #success) to a linear combination of categorical and factorial predictors, by means of a specific link function (e.g. the probit or logit function).

A quick explanation of multivariable linear models can be found in Hidalgo and Goodman (2013). Refer to our article, Moscatelli et al. (2012) , for a more detailed explanation of multivariable models within GLMM framework.

Multivariable GLMM can be extremely useful to fit psychophysical data and test our experimental hypothesis. This quick tutorial is to introduce the R syntax to multivariable GLMM. First, load a dataset including factorial and continuous predictor. Otherwise you can simulate the dataset using twice MERsimulate:

#In two steps: Simulate a dataset with a factorial predictor ("condition")
datafr.1 <- MERsimulate(fixeff = c(-7.5, 0.0875), nsubject = 6,
 constant = T)                               
levels(datafr.1$Subject) = c("S1", "S2", "S3", "S4", "S5", "S6")
datafr.1$condition = rep("A", 54)

#note the difference in intercept (PSE) in the second dataset
datafr.2 <- MERsimulate(fixeff = c(-6.5, 0.0875),nsubject = 6,
 constant = T)                              
levels(datafr.2$Subject) = c("S1", "S2", "S3", "S4", "S5", "S6")
datafr.2$condition = rep("B", 54)

datafr = merge(datafr.1, datafr.2, all = T)
datafr$condition = as.factor(datafr$condition)

Now, use glmer{lme4.1} to fit the model. For clarity I first created a formula object that I called “formula.mod” and then passed to the glmer function. The model takes one of the two levels of the factorial predictor (In this case “A”) as baseline. This is the so-called “dummy coding”, see Moscatelli et al. (2012) for details. As we specified in our model formula, the model tests for a difference in intercept with respect to the baseline (the factorial predictor named “conditionB” in model summary) and for a difference in slope (“conditionB:X”).

#Here I assumed as random predictors the intercept, "X" and the "condition". 
#That is, I assumed that the effect of each predictor can vary randomly between subjects.
#This is not necessarily the case, you should try alternative models including for e.g. only
#a random intercept, or a random intercept and "X" etc.
formula.mod = cbind(Longer, Total - Longer) ~ X + condition + X * condition + (1 + condition + X| Subject)
mod1 <- glmer(formula = formula.mod, family = binomial(link = "probit"), data = datafr)

The model is coded in terms of intercept and slope, whereas in psychophysics we express our estimates in terms of JND and PSE. You can get the estimated PSE and JND for the “A” and “B” condition using MERpsychophycs. For a quick estimate, use the delta method as follows:

xplode.mod1 = xplode.mer(model = mod1, name.cont = "X", name.fact = "condition")
MERtreatment(xplode.mod1, datafr)

Delta method is based on asymptotic assumptions of a Gaussian distribution of the parameters, we cannot fully trust this assumption within GLMM framework. For a more correct estimate, use a bootstrap method as follows. First, create your function. Something like:

fun2mod1 = function(mer.obj){
  #allocate space: 4 parameters (jnd_A, jnd_B, pse_A, pse_B)
  jndpse = vector(mode = "numeric", length = 4)
  names(jndpse) = c("jnd_A","jnd_B", "pse_A", "pse_B")
  jndpse[1] = qnorm(0.75)/fixef(mer.obj)[2] #jnd_A
  jndpse[2] = qnorm(0.75)/(fixef(mer.obj)[2] + fixef(mer.obj)[4]) #jnd_B
  jndpse[3] = -fixef(mer.obj)[1]/fixef(mer.obj)[2] #pse_A
  jndpse[4] = -(fixef(mer.obj)[1] + fixef(mer.obj)[3])/(fixef(mer.obj)[2] + fixef(mer.obj)[4]) #pse_B

Then, run the bootstrap simulation (this takes time – up to 30 min, depending on the complexity of the model)

BootEstim = pseMer(mod1, B = 200, FUN = fun2mod1)

If you use the MERpsychophycs code or you would like to acknowledge this tutorial, please cite our article:

Moscatelli A, Mezzetti M, Lacquaniti F (2012). Modelling Psychophysical Data at the Population-Level: The Generalized Linear Mixed Model. Journal of Vision, 12(11):26, 1-17.


4 Comments to “Multivariable GLMMs”

  1. Hello! I am confused as to where you define “X”. working on plotting a model right now actually but struggling to follow your instructions because I’m unclear what X is. Obviously its some sort of variable in your df but I never saw it declared.


    • Hi Paula,

      sorry for late reply, busy time. X is the name of the predictor variable in the model formula. For example, if your formula is:

      y ~ speed,

      then replace “X” with “speed” in the examples. If you use gglplot2, you may create a new data frame for a more dense predictions, then you need to add a variable having the same name as the predictor in the formula in the new data frame.

      PS: with my colleague Priscilla, we are working on a tutorial to be submitted as a journal paper. If you want I can send you a preliminary draft privately.

  2. Hello!

    I’m currently working on a time perception project, and so far the MixedPsy package has been great. However, we are not sure about some of the details. If you could send me a preliminary draft of the tutorial, that would help a lot!

    Kind regards,

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 )

Google+ photo

You are commenting using your Google+ 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 )


Connecting to %s

%d bloggers like this: