从 R 中的 glm 和 lmer 提取的残差方差答案

【问题标题】：Residual variance extracted from glm and lmer in R从 R 中的 glm 和 lmer 提取的残差方差
【发布时间】：2012-09-27 07:59:59
【问题描述】：

我正在尝试将我所阅读的有关多级建模的内容与我对 R 中 glm 的了解合并。我现在使用来自 here 的高度增长数据。

我已经完成了如下所示的一些编码：

library(lme4)
library(ggplot2)

setwd("~/Documents/r_code/multilevel_modelling/")

rm(list=ls())

oxford.df <- read.fwf("oxboys/OXBOYS.DAT",widths=c(2,7,6,1))
names(oxford.df) <- c("stu_code","age_central","height","occasion_id")
oxford.df <- oxford.df[!is.na(oxford.df[,"age_central"]),]
oxford.df[,"stu_code"] <- factor(as.character(oxford.df[,"stu_code"]))
oxford.df[,"dummy"] <- 1

chart <- ggplot(data=oxford.df,aes(x=occasion_id,y=height))
chart <- chart + geom_point(aes(colour=stu_code))

# see if lm and glm give the same estimate
glm.01 <- lm(height~age_central+occasion_id,data=oxford.df)
glm.02 <- glm(height~age_central+occasion_id,data=oxford.df,family="gaussian")
summary(glm.02)
vcov(glm.02)
var(glm.02$residual)
(logLik(glm.01)*-2)-(logLik(glm.02)*-2)
1-pchisq(-2.273737e-13,1)
# lm and glm give the same estimation
# so glm.02 will be used from now on

# see if lmer without level2 variable give same result as glm.02
mlm.03 <- lmer(height~age_central+occasion_id+(1|dummy),data=oxford.df,REML=FALSE)
(logLik(glm.02)*-2)-(logLik(mlm.03)*-2)
# 1-pchisq(-3.408097e-07,1)
# glm.02 and mlm.03 give the same estimation, only if REML=FALSE

mlm.03 给我以下输出：

> mlm.03
Linear mixed model fit by maximum likelihood 
Formula: height ~ age_central + occasion_id + (1 | dummy) 
   Data: oxford.df 
  AIC  BIC logLik deviance REMLdev
 1650 1667 -819.9     1640    1633
Random effects:
 Groups   Name        Variance Std.Dev.
 dummy    (Intercept)  0.000   0.0000  
 Residual             64.712   8.0444  
Number of obs: 234, groups: dummy, 1

Fixed effects:
            Estimate Std. Error t value
(Intercept)  142.994     21.132   6.767
age_central    1.340     17.183   0.078
occasion_id    1.299      4.303   0.302

Correlation of Fixed Effects:
            (Intr) ag_cnt
age_central  0.999       
occasion_id -1.000 -0.999

您可以看到 random effect 部分中的残差存在差异，这是我从 Jos WR Twisk 的 Applied Multilevel Analysis - A Practical Guide 中读到的，这表示模型中“无法解释的差异”的数量。

我想知道是否可以从 glm.02 得出相同的剩余方差，所以我尝试了以下方法：

> var(resid(glm.01))
[1] 64.98952
> sd(resid(glm.01))
[1] 8.061608

结果与mlm.03 输出略有不同。这是否指的是 mlm.03 中所述的相同“残差”？

【问题讨论】：

这个问题在这里是题外话（应该在 stackoverflow 上）所以我投票关闭，但你可以访问与 attr(VarCorr(mlm.03),"sc")^2 匹配的 lmer 模型的剩余方差，因为你的模型被命名为 @ 987654334@。可以使用VarCorr() 函数访问其他方差分量。

标签： r variance

【解决方案1】：

您的glm.02 和glm.01 使用最小二乘法估计了一个简单的线性回归模型。另一方面，mlm.03 是通过最大似然估计的线性混合模型。我不知道你的数据集，但看起来你使用 dummy 变量在 2 级创建了一个零方差的集群结构。

所以您的问题基本上有两个答案，但在您的情况下，只有第二个答案很重要。 glm.02 和 mlm.03 模型确实不包含相同的残差方差估计，因为...

模型通常不同（混合效应与经典回归）。但是，在您的情况下，dummy 变量似乎抑制了混合模型中的附加方差分量。所以对我来说，模型似乎是平等的。
用于估计残差方差的方法不同。 glm 使用 LS，lmer 在您的代码中使用 ML。残差方差的 ML 估计略有偏差（导致方差估计更小）。这可以通过使用 REML 而不是 ML 来估计方差分量来解决。

但是，对于似然比检验，使用经典 ML（而不是 REML）仍然是必要且正确的。使用 REML 比较两种可能性是不正确的。

干杯！

【讨论】：