如何预测 merMod 对象（lme4）的术语？

Posted 2023-02-16

【中文标题】如何预测 merMod 对象（lme4）的术语？

【英文标题】：How to predict terms of merMod objects (lme4)?

【发布时间】：2015-07-14 22:42:22

【问题描述】：

对于简单的glm 对象，我可以使用predict(fit, type = "terms") 检索具有每个项的拟合值的矩阵。

lmer 的等价物是什么？ glmer 适合的型号？据我所知，predict.merMod 函数不支持type = terms。

【参考方案1】：

lmer 的等价物是什么？ glmer合身模特？

我认为没有。不过，您可以按如下方式轻松制作一个

#####
# fit model with one terms which is a matrix
library(lme4)
fit <- lmer(Reaction ~ cbind(Days, (Days > 3) * Days) + (Days | Subject), 
            sleepstudy)

#####
# very similar code to `predict.lm`
pred_terms_merMod <- function(fit, newdata)
  tt <- terms(fit)
  beta <- fixef(fit)

  mm <- model.matrix(tt, newdata)
  aa <- attr(mm, "assign")
  ll <- attr(tt, "term.labels")
  hasintercept <- attr(tt, "intercept") > 0L
  if (hasintercept) 
    ll <- c("(Intercept)", ll)
  aaa <- factor(aa, labels = ll)
  asgn <- split(order(aa), aaa)
  if (hasintercept) 
    asgn$"(Intercept)" <- NULL
    avx <- colMeans(mm)
    termsconst <- sum(avx * beta)
  
  nterms <- length(asgn)
  if (nterms > 0) 
    predictor <- matrix(ncol = nterms, nrow = NROW(mm))
    dimnames(predictor) <- list(rownames(mm), names(asgn))
    if (hasintercept) 
      mm <- sweep(mm, 2L, avx, check.margin = FALSE)
    for (i in seq.int(1L, nterms, length.out = nterms)) 
      idx <- asgn[[i]]
      predictor[, i] <- mm[, idx, drop = FALSE] %*% beta[idx]
    
   else 
    predictor <- ip <- matrix(0, n, 0L)
  
  attr(predictor, "constant") <- if (hasintercept) termsconst else 0
  predictor


# use the function
newdata <- data.frame(Days = c(1, 5), Reaction = c(0, 0))
(out <- pred_terms_merMod(fit, newdata))
#R>   cbind(Days, (Days > 3) * Days)
#R> 1                        -21.173
#R> 2                         21.173
#R> attr(,"constant")
#R> [1] 283.24

#####
# confirm results
beta. <-  fixef(fit)
beta.[1] + beta.[2]
#R> (Intercept) 
#R>      262.07 
out[1] + attr(out, "constant")
#R> [1] 262.07

beta.[1] + (beta.[2] + beta.[3]) * 5
#R> (Intercept) 
#R>      304.41 
out[2] + attr(out, "constant")
#R> [1] 304.41

据我所知，将上述内容扩展到包括标准错误应该很简单。