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)
summary(datafr)

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)
summary(mod1)

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
  return(jndpse)
}

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)
BootEstim

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.

Advertisements

Leave a Reply

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

WordPress.com Logo

You are commenting using your WordPress.com 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

%d bloggers like this: