[R-sig-ME] Error with predict.merMod and random slopes

[Ben Bolker]

On 13-11-12 07:50 AM, Nick Donnelly wrote:
> Hello All,
> 
> This is my first time posting, and I'm not a very experienced user of R and
> lme4, so I apologise in advance if I've missed something important or this
> issue has been addressed previously.
> 
> I'm trying to use a binomial mixed model to predict whether subjects make a
> certain type of response in trials in a behavioural task. I've been fitting
> a binomial mixed model using lme4, and have been using the ID of my subjects
> as a random intercept, and using the predict.merMod method to predict the
> outcome of new trial data, which has been working. 
> 
> However, having recently read Barr et al 2013
> (http://idiom.ucsd.edu/~rlevy/papers/barr-etal-2013-jml.pdf) I've been
> trying to make my models maximal by adding random slopes for my
> within-subjects factors. When I add random slopes to my models and use
> predict.merMod I now get an error:
> 
> "Error: sum(nb) == q is not TRUE"
> 
> Here is the code, and a link to the data I am working with, which is
> hopefully fully self-contained and working:
> 
> library(lme4)
> 
> latencyData <-
> read.csv("https://dl.dropboxusercontent.com/s/uhn86f5nae662i6
> /latencyDataTest.dat",header = TRUE)
> 
> #latencyData$prem is the outcome of each behavioural trial (premature or
> #not premature) - the dependent variable
> #latencyData$group is a between subjects factor 
> #latencyData$latency is a reaction time that is measured as part of each
> #trial
> #latencyData$PO is the outcome of the previous trial
> #latencyData$subjectID is a code for the subject which made each trial
> 
> #take the full dataset apart from one trial, which I will use for
> #prediction later
> x <- 1
> ix <- 1:dim(latencyData)[1] != x
> 
> dataTrain <- latencyData[ix,]
> dataTest  <- latencyData[x ,]
> 
> 
> #fit a logistic regression model with the subject ID as random intercept
> m1 <- glmer(prem ~ group*latency*PO + (1|subjectID),data = dataTrain, 
> family = binomial)
> 
> #fit a logistic regression model that attempts to be "maximal" after Barr
> #et al 2013, including random slopes for the subject's latency, the
> #previous trial outcome, and their interaction(the highest order
> within-#subjects factors)
> 
> m2 <- glmer(prem ~ group*latency*PO + (latency:PO|subjectID),
> data = dataTrain, family = binomial)
> 
> #Try a model with just 1 random slope
> m3 <- glmer(prem ~ group*latency*PO + (PO|subjectID),
> data = dataTrain,family = binomial)
> 
> 
> #try and predict the left out trial with the random intercept only model-
> #gives a result
> predict(m1,dataTest,REform = NULL,type = "response")
> 
> #try and predict the left out trial with the random intercepts and slopes 
> #model
> predict(m2,dataTest,REform = NULL,type = "response")
> 
> #this gives an error:"Error: sum(nb) == q is not TRUE"
> 
> #try and predict the left out trial with the random intercepts and single
> #random slope model
> predict(m3,dataTest,REform = NULL,type = "response")
> 
> #this gives the same error:"Error: sum(nb) == q is not TRUE"
> 
> #If I remove the RE from the prediction, I can get a result rather than an
> #error
> predict(m3,dataTest,REform = NA,type = "response")
> 
> However, I would like to be able to make predictions using the maximal model
> and the random effects structure. Is there anything anyone can suggest I can
> do to get this working, or is it not yet an available feature in lme4 (I
> know the predict.merMod method is new)?


  Thanks for the very detailed report.  This *should* work, the fact
that it doesn't constitutes a bug.

  I'm guessing that you might be able to work around the problems by
defining the interactions up front, e.g.

transform(dataTrain,latencyPO=interaction(latency,PO))

and then using (latencyPO|subjectID) as the random effect.

  I will take a closer look (maybe with a smaller example so that it
runs faster ...)

  Ben Bolker