2

与 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
4

0 回答 0