【问题标题】:How to predict terms of merMod objects (lme4)?如何预测 merMod 对象(lme4)的术语?
【发布时间】:2015-07-14 22:42:22
【问题描述】:

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

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

【问题讨论】:

    标签: r predict lme4


    【解决方案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
    

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

    【讨论】:

      猜你喜欢
      • 2015-02-08
      • 1970-01-01
      • 1970-01-01
      • 2021-01-13
      • 2019-10-17
      • 2013-05-16
      • 1970-01-01
      • 2014-02-17
      • 1970-01-01
      相关资源
      最近更新 更多