3
$\begingroup$

I am running a probit regression with a random effect:

m1<-glmer(Binary~Explan+(1|Random),family=binomial(link="probit"))

where Explan is a three-level categorical variable.

I want to calculate the mean predicted probabilities for each level of Explan. I tried doing so using this code:

newdata=data.frame(Explan="First")
predict(m1,newdata,type="response")

where First is a level of the categorical Explan variable.

However I get the following error message:

Error: (p <- ncol(X)) == ncol(Y) is not TRUE

Were this a logit model, I would simply strip the model of the intercept and then back-transform the model summary coefficients to get the predicted values that I'm after, but I am unsure of how I would go about this with a mixed-effects probit model.

Any help in extracting the predicted probabilities would be greatly appreciated.

Improve this question
$\endgroup$
2

1 Answer 1

Reset to default
2
$\begingroup$

This question and excellent exchange was the impetus for creating the predictInterval function in the merTools package. predictInterval is designed to use the arm::sim functions to generate distributions of parameters in the model and then to use those distributions to generate simulated values of the response given the newdata provided by the user. It's simple to use -- all you would need to do is:

library(merTools)
preds <- predictInterval(m1, newdata = newdata, n.sims = 999)

You can specify a whole host of other values to predictInterval including setting the interval for the prediction intervals, choosing whether to report the mean or median of the distribution, and choosing whether or not to include the residual variance from the model.

It's not a full prediction interval because the variability of the theta parameters in the lmer object are not included, but all of the other variation is captured through this method, giving a pretty decent approximation.

merTools also allows you to easily construct simulated scenarios for analysis. For example, if you wanted to compute the effect of changing a level of a factor variable for the average observation you could do this:

subExample <- list(Explan= "First")
example5 <- draw(m1, type = 'average', varList = subExample) 
predictInterval(m1, newdata = example5, n.sims = 1000)

Check out the package vignettes here and here to learn more.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you very much for this, however unfortunately it would appear that the predictInterval function doesn't work for probit models yet, output: Prediction for this NLMMs or GLMMs that are not mixed logistic regressions is not yet implemented. $\endgroup$
    – milvus_milvus
    Commented Aug 24, 2015 at 20:30
  • $\begingroup$ Ah, yeah. I think I can easily work around that within the day. We have had some debate about less standard link functions, but we should definitely be able to support probit. I thought we had changed that but haven't yet. I'll put a patch up on the GitHub repo in the next hour or so. Then you can use devtools::install_github("jknowles/merTools") to get the development version and try it out. $\endgroup$
    – jknowles
    Commented Aug 24, 2015 at 21:38
  • $\begingroup$ Thanks for such a quick patch. I have tried it and it seems to work - in as much as I have a new error: Error in [.data.frame(newdata, , j) : undefined columns selected. I suspect this is human error on my part, but I would be grateful if you had any advice. I have checked for spelling/spacing errors etc., as they seem to occasionally throw up this error. $\endgroup$
    – milvus_milvus
    Commented Aug 25, 2015 at 0:26
  • $\begingroup$ I think the key is that newdata has to have the same number of columns as the original dataset. We provide some functionality to make it easier to do this in the package. Can you confirm though that if you try: predictInterval(m1, newdata = m1@frame[1, ]) it works? $\endgroup$
    – jknowles
    Commented Aug 25, 2015 at 13:47

Not the answer you're looking for? Browse other questions tagged or ask your own question.