I am fitting a model with somewhat collinear predictors, and then predicting out of it. This works fine when using the original model object, but for performance reasons, I want to save a "chopped"/reduced-size version of the model, and then load that reduced model again later to make predictions.
Despite the collinearity, I am able to predict()
using the full model, but when I attempt the same operation with the reduced-size model, I get
Error in X %*% fixef(object) : non-conformable arguments
What can I add to the glmer_chop()
function that will let me predict out of the reduced_model
with no issues, just as I can with original_model
?
Reprex:
library(lme4)
# This is intended to reduce the size of the (g)lmer object
glmer_chop <- function(object) {
newobj <- object
newobj@frame <- model.frame(object)[0, 1:ncol(model.frame(object))]
attr(newobj@frame, "terms") <- attr(object@frame, "terms")
newobj@pp <- with(object@pp,
new("merPredD",
Lambdat=Lambdat,
Lind=Lind,
theta=theta,
delu=delu,
u=u,u0=u0,
n=nrow(X),
X=matrix(1,nrow=nrow(X)),
Zt=Zt))
newobj@resp <- new("glmResp",family=binomial(),y=numeric(0))
return(newobj)
}
original_model <- glmer(vs ~
disp +
hp +
I(cyl > 6) + # These two vars are correlated
factor(cyl) + # These two vars are correlated
(1 | am),
family = "binomial",
data = mtcars)
summary(original_model)
reduced_model <- glmer_chop(original_model)
mtcars$pred_original <- mtcars$pred_reduced <- NA
# Predicting from full model works
mtcars$pred_original <- predict(original_model, newdata = mtcars, type = "response")
# Predicting from reduced model errors with:
# Error in X %*% fixef(object) : non-conformable arguments
mtcars$pred_reduced <- predict(reduced_model, newdata = mtcars, type = "response")
It turns out only a small change is needed: in the @pp
slot, keep all of the original X
, not just a reduced matrix version of it.
glmer_chop <- function(object) {
newobj <- object
newobj@frame <- model.frame(object)[0, 1:ncol(model.frame(object))]
attr(newobj@frame, "terms") <- attr(object@frame, "terms")
newobj@pp <- with(object@pp,
new("merPredD",
Lambdat=Lambdat,
Lind=Lind,
theta=theta,
delu=delu,
u=u,u0=u0,
n=nrow(X),
X=X, # THIS IS THE CHANGE
Zt=Zt))
newobj@resp <- new("glmResp",family=binomial(),y=numeric(0))
return(newobj)
}
Here's why: predict()
calls fixef()
, which in turn calls getME()
, which retrieves the entire X
. fixef()
uses that object to assign dimnames
to the fixed effects; without them, the matrix multiplication fails. (The reason it fails is not clear to me, but that's how it is.) Compare the output of fixef()
for the original and reduced models:
> fixef(original_model)
(Intercept) disp hp I(cyl > 6)TRUE factor(cyl)6
335.504295 -1.398626 1.096582 -23793.196306 -76.846436
> fixef(reduced_model)
[1] 335.504295 -1.398626 1.096582 -23793.196306 -76.846436
With the edit above, the reduced model looks like the original model.
> fixef(revised_reduced_model)
(Intercept) disp hp I(cyl > 6)TRUE factor(cyl)6
335.504295 -1.398626 1.096582 -23793.196306 -76.846436