与 lme4:: lmer 对象相比,当我在 lmerTest:: lmer 上使用 anova 时,我得到的平方和和平均平方和高 10 倍。请参阅下面的 R 日志文件。请注意警告,当我附加 lmerTest 包时, stats::sigma 函数会覆盖 lme4::sigma 函数,我怀疑正是这导致了差异。此外,anova 报告现在说它是 III 型 anova,而不是预期的 I 型。这是 lmerTest 包中的错误,还是使用 Kenward-Roger 近似来改变自由度SumSQ 和 MSS 的计算以及我不明白的 anova 类型的规范?
我会附加测试文件,但它是机密的临床试验信息。如有必要,我可以看看我是否可以拼凑出一个测试用例。
提前感谢大家对此提供的任何建议。
> library(lme4)
Loading required package: Matrix
Attaching package: ‘lme4’
The following object is masked from ‘package:stats’:
sigma
> test100 <- lmer(log(value) ~ prepost * lowhi + (1|CID/LotNo/rep),
REML = F, data = GSIRlong, subset = !is.na(value))
> library(lmerTest)
Attaching package: ‘lmerTest’
The following object is masked from ‘package:lme4’:
lmer
The following object is masked from ‘package:stats’:
step
Warning message:
replacing previous import ‘lme4::sigma’ by ‘stats::sigma’ when loading
‘lmerTest’
> test200 <- lmer(log(value) ~ prepost * lowhi + (1|CID/LotNo/rep),
REML = F, data = GSIRlong, subset = !is.na(value))
> anova(test100)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
prepost 1 3.956 3.956 18.4825
lowhi 1 130.647 130.647 610.3836
prepost:lowhi 1 0.038 0.038 0.1758
> anova(test200, ddf = 'Ken')
Analysis of Variance Table of type III with Kenward-Roger
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
prepost 37.15 37.15 1 308.04 18.68 2.094e-05 ***
lowhi 1207.97 1207.97 1 376.43 607.33 < 2.2e-16 ***
prepost:lowhi 0.35 0.35 1 376.43 0.17 0.676
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
更新:谢谢,本。我在 lmerTest 上做了一点代码考古,看看我是否可以对上述异常现象做出解释。首先,事实证明 lmerTest::lmer 只是将模型提交给 lme4::lmer,然后将结果重新标记为“mermodLmerTest”对象。这样做的唯一效果是调用 lmerTest 包中的 summary() 和 anova() 版本,而不是来自 base 和 stats 的通常默认值。(这些 lmerTest 函数已编译,我还没有进一步查看 C++ 代码。) lmerTest::summary 只是在 base::summary 结果中添加了三列,给出了 df、t 值和 Pr。请注意,默认情况下,lmerTest::anova 计算类型 III 的 anova,而不是 stats::anova 中的类型 I。(解释我上面的第二个问题。)如果一个人的模型包括交互,这不是一个很好的选择。
然而,使用 summary 和 anova 的 nlmeTest 版本还有其他惊喜,如下面的 R 控制台文件所示。我使用了 lmerTest 包含的 sleepstudy 数据,因此该代码应该是可复制的。
一个。请注意,“sleepstudy”有 180 条记录(有 3 个变量)
湾。fm1 和 fm1a 的摘要是相同的,除了添加了固定效应列。但请注意,在 lmerTest::summary 中,截距和天数的 ddf 分别为 1371 和 1281;奇怪的是,“sleepstudy”中只有 180 条记录。
C。就像我上面的原始模型一样,nlm4 和 nlmrTest 版本的 anova 给出了非常不同的 Sum Sq 和 Mean Sq 值。(分别为 30031 和 446.65)。
d:有趣的是,使用 Satterthwaite 和 Kenward-Rogers 估计的 DenDF 的 nlmrTest 版本的 anova 差异很大(分别为 5794080 和 28)。KR 值似乎更合理。
鉴于上述问题,我现在不愿意依赖 lmerTest 来提供可靠的 p 值。根据您(Doug Bates 的)博客条目(https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html),我现在正在使用(并推荐)Dan 发布的方法米尔曼(http://mindingthebrain.blogspot.ch/2014/02/three-ways-to-get-parameter-specific-p.html) 在下面的最后一段代码中估计一个朴素的 t 检验 p 值(假设基本上是无限自由度)和一个 Kenward-Rogers 对 df 的估计(使用 R 包 'pbkrtest' - lmerTest 使用的相同包)。我找不到 R 代码来计算 Satterthwaite 估计值。朴素的 t 检验 p 值显然是反保守的,但 KR 估计值被认为是相当不错的。如果两者给出了相似的 p 估计值,那么我认为人们可以对 [naive t-test, KR 估计值] 范围内的 p 值感到满意。
> library(lme4); library(lmerTest); library(pbkrtest);
dim(sleepstudy)
[1] 180 3
>
> fm1 <- lme4::lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
> fm1a <- lmerTest::lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
>
> summary(fm1)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (Days | Subject)
Data: sleepstudy
REML criterion at convergence: 1743.6
Scaled residuals:
Min 1Q Median 3Q Max
-3.9536 -0.4634 0.0231 0.4634 5.1793
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 612.09 24.740
Days 35.07 5.922 0.07
Residual 654.94 25.592
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 251.405 6.825 36.84
Days 10.467 1.546 6.77
Correlation of Fixed Effects:
(Intr)
Days -0.138
> summary(fm1a)
Linear mixed model fit by REML t-tests use Satterthwaite approximations to
degrees of freedom [lmerMod]
Formula: Reaction ~ Days + (Days | Subject)
Data: sleepstudy
REML criterion at convergence: 1743.6
Scaled residuals:
Min 1Q Median 3Q Max
-3.9536 -0.4634 0.0231 0.4634 5.1793
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 612.09 24.740
Days 35.07 5.922 0.07
Residual 654.94 25.592
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 251.405 6.825 1371.100 302.06 <2e-16 ***
Days 10.467 1.546 1281.700 55.52 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
Days -0.138
Warning message:
In deviance.merMod(object, ...) :
deviance() is deprecated for REML fits; use REMLcrit for the REML
criterion or deviance(.,REML=FALSE) for deviance calculated at the REML fit
>
> anova(fm1)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
Days 1 30031 30031 45.853
> anova(fm1a, ddf = 'Sat', type = 1)
Analysis of Variance Table of type I with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
Days 446.65 446.65 1 5794080 45.853 1.275e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning message:
In deviance.merMod(object, ...) :
deviance() is deprecated for REML fits; use REMLcrit for the REML criterion or deviance(.,REML=FALSE) for deviance calculated at the REML fit
> anova(fm1a, ddf = 'Ken', type = 1)
Analysis of Variance Table of type I with Kenward-Roger
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
Days 446.65 446.65 1 27.997 45.853 2.359e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning message:
In deviance.merMod(object, ...) :
deviance() is deprecated for REML fits; use REMLcrit for the REML criterion or deviance(.,REML=FALSE) for deviance calculated at the REML fit
>
> # t.test
> coefs <- data.frame(coef(summary(fm1)))
> coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
> coefs
Estimate Std..Error t.value p.z
(Intercept) 251.40510 6.824556 36.838311 0.000000e+00
Days 10.46729 1.545789 6.771485 1.274669e-11
>
> # Kenward-Rogers
> df.KR <- get_ddf_Lb(fm1, fixef(fm1))
> df.KR
[1] 25.89366
> coefs$p.KR <- 2 * (1 - pt(abs(coefs$t.value), df.KR))
> coefs
Estimate Std..Error t.value p.z p.KR
(Intercept) 251.40510 6.824556 36.838311 0.000000e+00 0.0000e+00
Days 10.46729 1.545789 6.771485 1.274669e-11 3.5447e-07