机器算法验证 - anova() 命令对 lmer 模型对象有什么作用？ - 吾爱随笔录

anova() 命令对 lmer 模型对象有什么作用？

机器算法验证 r 方差分析混合模式 lme4-nlme

2022-01-17 19:45:25

lmer希望这是一个问题，这里有人可以回答我从适合的混合效应模型（来自lme4 R 包）分解平方和的性质。

首先，我应该说我知道使用这种方法存在争议，在实践中，我更有可能使用自举 LRT 来比较模型（正如 Faraway，2006 所建议的那样）。但是，我对如何复制结果感到困惑，所以为了我自己的理智，我想我会在这里问。

基本上，我开始掌握使用适合lme4包装的混合效果模型。我知道您可以使用该anova()命令来总结顺序测试模型中的固定效应。据我所知，这就是 Faraway (2006) 所指的“预期均方”方法。我想知道的是如何计算平方和？

我知道我可以从特定模型中获取估计值（使用coef()），假设它们是固定的，然后使用模型残差的平方和（有和没有感兴趣的因素）进行测试。这对于包含单个主体内因素的模型来说很好。但是，在实施裂区设计时，我得到的平方和值等于 R 使用aov()适当Error()名称产生的值。但是，这与模型对象上的命令产生的平方和不同，尽管 F 比率是相同的。anova()

当然，这是完全有道理的，因为混合模型中不需要Error()地层。但是，这必须意味着在混合模型中以某种方式惩罚平方和，以便提供适当的 F 比率。这是如何实现的？以及模型如何以某种方式纠正图间平方和但不纠正图内平方和。显然，这是通过为不同的效应指定不同的误差值来实现的经典裂区方差分析所必需的，那么混合效应模型如何允许这一点呢？

基本上，我希望能够自己复制anova()应用于 lmer 模型对象的命令的结果，以验证结果和我的理解，但是，目前我可以通过正常的主题内设计实现这一点，但不能用于拆分-情节设计，我似乎无法找出为什么会这样。

举个例子：

library(faraway)
library(lme4)
data(irrigation)

anova(lmer(yield ~ irrigation + variety + (1|field), data = irrigation))

Analysis of Variance Table
           Df Sum Sq Mean Sq F value
irrigation  3 1.6605  0.5535  0.3882
variety     1 2.2500  2.2500  1.5782

summary(aov(yield ~ irrigation + variety + Error(field/irrigation), data = irrigation))

Error: field
           Df Sum Sq Mean Sq F value Pr(>F)
irrigation  3  40.19   13.40   0.388  0.769
Residuals   4 138.03   34.51               

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
variety    1   2.25   2.250   1.578  0.249
Residuals  7   9.98   1.426

从上面可以看出，所有 F 比率都一致。多样性的平方和也一致。然而，灌溉的平方和不一致，但似乎 lmer 的输出是按比例缩放的。那么 anova() 命令实际上做了什么？

1个回答

使用来源，卢克。我们可以通过执行来查看 ANOVA 函数的内部getAnywhere(anova.Mermod)。该函数的第一部分是用于比较两个不同的模型。固定效应上的方差分析在else下半场大块出现：

 dc <- getME(object, "devcomp")
        X <- getME(object, "X")
        asgn <- attr(X, "assign")
        stopifnot(length(asgn) == (p <- dc$dims[["p"]]))
            ss <- as.vector(object@pp$RX() %*% object@beta)^2
        names(ss) <- colnames(X)
        terms <- terms(object)
        nmeffects <- attr(terms, "term.labels")[unique(asgn)]
        if ("(Intercept)" %in% names(ss)) 
            nmeffects <- c("(Intercept)", nmeffects)
        ss <- unlist(lapply(split(ss, asgn), sum))
        stopifnot(length(ss) == length(nmeffects))
        df <- vapply(split(asgn, asgn), length, 1L)
        ms <- ss/df
        f <- ms/(sigma(object)^2)
        table <- data.frame(df, ss, ms, f)
        dimnames(table) <- list(nmeffects, c("Df", "Sum Sq", 
            "Mean Sq", "F value"))
        if ("(Intercept)" %in% nmeffects) 
            table <- table[-match("(Intercept)", nmeffects), 
                ]
        attr(table, "heading") <- "Analysis of Variance Table"
        class(table) <- c("anova", "data.frame")
        table

object是 lmer 输出。我们在第 5 行开始计算平方和：ss <- as.vector ...代码将固定参数 (in beta) 乘以上三角矩阵；然后平方每一项。这是灌溉示例的上三角矩阵。每行对应于五个固定效应参数之一（截距，灌溉的 3 个自由度，多样性的 1 个自由度）。

zapsmall(irrigation.lmer@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]  [,5]
[1,] 0.813 0.203  0.203  0.203 0.407
[2,] 0.000 0.352 -0.117 -0.117 0.000
[3,] 0.000 0.000  0.332 -0.166 0.000
[4,] 0.000 0.000  0.000  0.287 0.000
[5,] 0.000 0.000  0.000  0.000 2.000

第一行为您提供截距的平方和，最后一行为您提供场内变化效果的 SS。第 2-4 行仅涉及灌溉水平的 3 个参数，因此预乘为您提供了 SS 的三部分用于灌溉。

这些片段本身并不有趣，因为它们来自 R 中的默认处理对比，但在一行中，ss <- unlist(lapply(split ....Bates 根据级别数和它们引用的因素收集了一些平方和。这里有很多簿记工作。我们也得到了自由度（灌溉为 3）。然后，他得到显示在打印输出上的均方anova。最后，他用组内残差除以所有均方sigma(object)^2。

发生什么了？的哲学lmer与所使用的矩方法无关aov。的想法lmer是最大化通过整合看不见的随机效应获得的边际可能性。在这种情况下，每个领域的随机生育水平。Pinheiro 和 Bates 的第 2 章描述了这个过程的丑陋。RX用来求平方和的矩阵就是它们的矩阵 $R_{00}$ 来自正文第 70 页的公式 2.17。这个矩阵是从一堆 QR 分解中获得的（除其他外）随机效应的设计矩阵和 $\sqrt{\sigma^2/\sigma_f^2}$ ，在哪里 $\sigma_f^2$ 是场效应的方差。这是您所询问的缺失因素，但它并没有以透明或简单的方式进入解决方案。

渐近地，固定效应的估计具有分布：

\hat{β} \sim N (β, σ^{2} [R_{00}^{- 1} R_{00}^{- T}])

$\hat{\beta} \sim \mathcal{N}(\beta, \sigma^2[R_{00}^{-1}R_{00}^{-T}])$

这意味着 $R_{00}\hat{\beta}$ 是独立分布的。如果 $\beta=0$ ，这些项的平方（除以 $\sigma^2$ ) 服从卡方分布。除以估计 $\sigma^2$ （另一个卡方除以 $\sigma^2$ ) 给出 F 统计量。我们不会除以组间分析的均方误差，因为这与这里发生的事情无关。我们需要比较通过 $R_{00}$ 通过将它们与估计值进行比较 $\sigma^2$ .

请注意，如果数据不平衡，您将不会获得相同的 F 统计量。如果您使用 ML 而不是 REML，您也不会获得相同的 F 统计量。

背后的想法aov是灌溉的预期均方是 $\sigma^2$ , $\sigma_f^2$ 和灌溉效果。场残差的期望均方是 $\sigma^2$ 和 $\sigma_f^2$ . 当灌溉效果为 0 时，这两个量都估计相同的东西，并且它们的比率服从 F 分布。

有趣的是，Bates 和 Pinheiro 建议使用 ANOVA 来拟合两个模型并进行似然比检验。后者倾向于反保守主义。

当然，如果数据不平衡，则不再清楚 ANOVA 正在测试什么假设。我从灌溉数据中删除了第一个观察结果并进行了改装。这里是 $R_{00}$ 再次矩阵：

zapsmall(fit2@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]   [,5]
[1,] 0.816 0.205  0.205  0.205  0.457
[2,] 0.000 0.354 -0.119 -0.119 -0.029
[3,] 0.000 0.000  0.334 -0.168 -0.040
[4,] 0.000 0.000  0.000  0.288 -0.071
[5,] 0.000 0.000  0.000  0.000  1.874

如您所见，灌溉参数的平方和现在也包含一些variety影响。

其它你可能感兴趣的问题

上一篇关于统计学硕士课程需要考虑的事项下一篇层次模型、神经网络、图形模型、贝叶斯网络之间有什么关系？