仅在使用“切碎的”(g)lmer 模型对象时才存在不合格参数

Question

我正在拟合一个具有一定共线预测变量的模型，然后根据它进行预测。使用原始模型对象时效果很好，但出于性能原因，我想保存模型的“切碎”/缩小尺寸版本，然后稍后再次加载该缩小模型以进行预测。

尽管存在共线性，我仍然能够

predict()

使用完整模型，但是当我尝试使用缩小尺寸的模型进行相同操作时，我得到

Error in X %*% fixef(object) : non-conformable arguments

我可以在

glmer_chop()

函数中添加什么，让我能够毫无问题地从

reduced_model

中进行预测，就像使用

original_model

一样？

代表：

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

Answer 1

事实证明只需要一个小小的改变：在

@pp

插槽中，保留所有原始

X

，而不仅仅是它的简化矩阵版本。

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

原因如下：

predict()

调用

fixef()

，又调用

getME()

，检索整个

X

。

fixef()

使用该对象将

dimnames

分配给固定效果；没有它们，矩阵乘法就会失败。 （我不清楚失败的原因，但事实就是如此。）比较原始模型和简化模型的

fixef()

的输出：

> 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

通过上述编辑，缩小后的模型看起来就像原始模型一样。

> fixef(revised_reduced_model)
   (Intercept)           disp             hp I(cyl > 6)TRUE   factor(cyl)6 
    335.504295      -1.398626       1.096582  -23793.196306     -76.846436