stats::anova 和 car::Anova 在评估 lme4 的线性混合效应模型时的区别
Difference between stats::anova and car::Anova in evaluating a linear mixed effect model from lme4
我正在了解有关 lme4 软件包的更多信息,并且对两者都表示赞赏 Bodo Winter's tutorial and this guide on Tufts。但是,这两个指南在建议确定固定效应显着性的方法时有所不同。
Winters 建议使用 R 的 anova
函数来比较一个具有所讨论固定效应的模型和一个不具有固定效应的模型。
相比之下,塔夫茨首先建议使用car
包的Anova
功能(他们还建议使用anova
方法)。
但是,从下面的播放示例中可以看出,这两种方法 return 不同的卡方和 p 值。
library(lme4)
# meaningless models
lmer_wt_null = lmer(mpg ~ (1 + wt | cyl), data = mtcars, REML = FALSE)
lmer_wt_full = lmer(mpg ~ wt + (1 + wt | cyl), data = mtcars, REML = FALSE)
# stats::anova output (Winters)
anova(lmer_wt_null, lmer_wt_full)
# Data: mtcars
# Models:
# lmer_wt_null: mpg ~ (1 + wt | cyl)
# lmer_wt_full: mpg ~ wt + (1 + wt | cyl)
# Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
# lmer_wt_null 5 167.29 174.62 -78.647 157.29
# lmer_wt_full 6 163.14 171.93 -75.568 151.14 6.1563 1 0.01309 *
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
library(car)
# car::anova output (Tufts)
Anova(lmer_wt_full)
# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: mpg
# Chisq Df Pr(>Chisq)
# wt 19.213 1 1.169e-05 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
这两种方法有哪些不同之处,这些 p 值之间的差异有何意义?
我几乎可以肯定我遗漏了一些基本的东西。谢谢。
我打算投票将其迁移到 CrossValidated,但是:这是一个很好的例子,说明了 Wald 和似然比检验之间的区别 p-values。
anova()
的结果基于似然比检验,等同于(您可以检查)计算 chi-squared 分布的上尾区:
pchisq(deviance(lmer_wt_null)-deviance(lmer_wt_full), df=1, lower.tail=FALSE)
car::Anova()
的结果基于 Wald 检验,该检验做出了更强的假设(他们假设 log-likelihood 曲面是二次曲面)。这里的测试是基于正态分布1,2:
的上尾的two-tailed测试
(cc <- coef(summary(lmer_wt_full)))
2*pnorm(abs(cc["wt","t value"]),lower.tail=FALSE)
我们可以通过计算和绘制 可能性曲线 来更深入地了解这一点;与二次曲线的偏差显示 Wald 检验失败的地方。
pp <- profile(lmer_wt_full)
dd <- as.data.frame(pp)
est <- cc["wt","Estimate"]
se <- cc["wt","Std. Error"]
library(ggplot2)
ggplot(subset(dd,.par=="wt" & .zeta>-2.6 & .zeta<2.6),aes(x=.focal,y=.zeta))+
geom_point()+geom_line()+
geom_abline(intercept=-est*se,slope=se,colour="red")+
geom_hline(yintercept=c(-1,1)*1.96)
ggsave("lmerprof.png")
黑线表示可能性曲线。 y-axis 显示 2* 偏差差异的带符号平方根——这基本上是一个 Normal-deviate 尺度。在此尺度上,二次 log-likelihood 曲面的 Wald 假设对应于线性轮廓。红线表示 Wald 近似。
我们还可以将基于似然概况的置信区间与基于 Wald 近似的置信区间进行比较(图中 +/- 1.96 截止点之间的区域):
可能性概况:
confint(pp)["wt",]
## 2.5 % 97.5 %
##-7.042211 -1.561525
沃尔德:
confint(lmer_wt_full,method="Wald")["wt",]
## 2.5 % 97.5 %
##-6.018525 -2.299228
1在许多情况下,这是基于 t-distribution 的,但这让我们陷入了更难的问题,即如何估计自由
2如果我们愿意,我们可以找到等效的 chi-squared 测试,但这通常是使用正态统计
我正在了解有关 lme4 软件包的更多信息,并且对两者都表示赞赏 Bodo Winter's tutorial and this guide on Tufts。但是,这两个指南在建议确定固定效应显着性的方法时有所不同。
Winters 建议使用 R 的 anova
函数来比较一个具有所讨论固定效应的模型和一个不具有固定效应的模型。
相比之下,塔夫茨首先建议使用car
包的Anova
功能(他们还建议使用anova
方法)。
但是,从下面的播放示例中可以看出,这两种方法 return 不同的卡方和 p 值。
library(lme4)
# meaningless models
lmer_wt_null = lmer(mpg ~ (1 + wt | cyl), data = mtcars, REML = FALSE)
lmer_wt_full = lmer(mpg ~ wt + (1 + wt | cyl), data = mtcars, REML = FALSE)
# stats::anova output (Winters)
anova(lmer_wt_null, lmer_wt_full)
# Data: mtcars
# Models:
# lmer_wt_null: mpg ~ (1 + wt | cyl)
# lmer_wt_full: mpg ~ wt + (1 + wt | cyl)
# Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
# lmer_wt_null 5 167.29 174.62 -78.647 157.29
# lmer_wt_full 6 163.14 171.93 -75.568 151.14 6.1563 1 0.01309 *
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
library(car)
# car::anova output (Tufts)
Anova(lmer_wt_full)
# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: mpg
# Chisq Df Pr(>Chisq)
# wt 19.213 1 1.169e-05 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
这两种方法有哪些不同之处,这些 p 值之间的差异有何意义?
我几乎可以肯定我遗漏了一些基本的东西。谢谢。
我打算投票将其迁移到 CrossValidated,但是:这是一个很好的例子,说明了 Wald 和似然比检验之间的区别 p-values。
anova()
的结果基于似然比检验,等同于(您可以检查)计算 chi-squared 分布的上尾区:
pchisq(deviance(lmer_wt_null)-deviance(lmer_wt_full), df=1, lower.tail=FALSE)
car::Anova()
的结果基于 Wald 检验,该检验做出了更强的假设(他们假设 log-likelihood 曲面是二次曲面)。这里的测试是基于正态分布1,2:
(cc <- coef(summary(lmer_wt_full)))
2*pnorm(abs(cc["wt","t value"]),lower.tail=FALSE)
我们可以通过计算和绘制 可能性曲线 来更深入地了解这一点;与二次曲线的偏差显示 Wald 检验失败的地方。
pp <- profile(lmer_wt_full)
dd <- as.data.frame(pp)
est <- cc["wt","Estimate"]
se <- cc["wt","Std. Error"]
library(ggplot2)
ggplot(subset(dd,.par=="wt" & .zeta>-2.6 & .zeta<2.6),aes(x=.focal,y=.zeta))+
geom_point()+geom_line()+
geom_abline(intercept=-est*se,slope=se,colour="red")+
geom_hline(yintercept=c(-1,1)*1.96)
ggsave("lmerprof.png")
黑线表示可能性曲线。 y-axis 显示 2* 偏差差异的带符号平方根——这基本上是一个 Normal-deviate 尺度。在此尺度上,二次 log-likelihood 曲面的 Wald 假设对应于线性轮廓。红线表示 Wald 近似。
我们还可以将基于似然概况的置信区间与基于 Wald 近似的置信区间进行比较(图中 +/- 1.96 截止点之间的区域):
可能性概况:
confint(pp)["wt",]
## 2.5 % 97.5 %
##-7.042211 -1.561525
沃尔德:
confint(lmer_wt_full,method="Wald")["wt",]
## 2.5 % 97.5 %
##-6.018525 -2.299228
1在许多情况下,这是基于 t-distribution 的,但这让我们陷入了更难的问题,即如何估计自由
2如果我们愿意,我们可以找到等效的 chi-squared 测试,但这通常是使用正态统计