r - LME模型在R-自由度中的Lsmeans软件包

Question

我对在使用nlme软件包构建的线性混合模型的情况下lsmeans软件包所使用的自由度有疑问。

这是一个示例，用于说明基于燕麦数据集的问题。我不是要讨论该模型与给定的数据集是否相关，而是要重现另一个数据集;-)的问题。

Oats.lme <- lme(yield ~ Variety, random = ~1 | Block, data = Oats)
anova(Oats.lme)

使用方差分析，我可以获得预期的64自由度。

numDF denDF  F-value p-value
(Intercept)     1    64 245.1409  <.0001
Variety         2    64   1.6654  0.1972

然后，我使用lsmeans函数：

lsmeans(Oats.lme, list(poly ~ Variety))

我得到

$`lsmeans of Variety`
Variety       lsmean       SE df lower.CL upper.CL
Golden Rain 104.5000 7.680866  5 84.75571 124.2443
Marvellous  109.7917 7.680866  5 90.04737 129.5360
Victory      97.6250 7.680866  5 77.88071 117.3693

Confidence level used: 0.95

$`polynomial contrasts of contrast`
contrast   estimate       SE df t.ratio p.value
linear     -6.87500  6.68529 64  -1.028  0.3076
quadratic -17.45833 11.57926 64  -1.508  0.1365

对于对比，我获得了相同的64 df，但是对于lsmeans本身，我只有5 df。我也使用SAS，并且对于相同类型的模型，对于lsmeans和对比度，我具有相同的df数（在当前示例中为64）。

我已经看到，使用lme4软件包时可能会更改自由度，但是我的代码嵌入在基于nlme的内部开发工具中，因此我基本上对nlme有所了解。

现在有人会为什么会发生这种情况吗？如果可以改变的话？还是我错过了什么？

更新-初始错误消息

最初，在一种特定情况下，当我对Dunnett的调整感兴趣时，我的随机游走效果只有2个级别，因此我注意到lsmeans的自由度降低了。因为我对对比的兴趣比对lsmeans的兴趣更大，所以现在我了解了它的来源，我仍然可以使用它，但是我把它放在那儿，以防万一有人有相同的错误并想知道为什么。

我在下面用燕麦数据示例复制了它。我获得的错误发生在lsmeans :::。qdunnx函数中，这是由于lsmeans的df为1。

Oats.lme <- lme(yield ~ Variety, random = ~1 | Block, data = subset(Oats,Block %in% c("I","II")))
lsm <- lsmeans(Oats.lme, trt.vs.ctrl ~ Variety)
summary(lsm,adjust = "dunnettx", infer = c(T, T), level = 0.95)

这是结果

$lsmeans
Variety      lsmean       SE df  lower.CL upper.CL
Golden Rain 123.250 15.88642  1 -78.60608 325.1061
Marvellous  125.500 15.88642  1 -76.35608 327.3561
Victory     115.125 15.88642  1 -86.73108 316.9811
Confidence level used: 0.95

$contrasts
contrast                 estimate      SE df t.ratio p.value
Marvellous - Golden Rain    2.250 12.8697 20   0.175  0.9695
Victory - Golden Rain      -8.125 12.8697 20  -0.631  0.7482
P value adjustment: dunnettx method for 2 tests

Error in if (abs(diff(r[1:2])) < 5e-04) return(r[1]) : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In qtukey(p, (1 + sqrt(1 + 8 * k))/2, df) : production de NaN

Answer 1

该模型表示响应变量有两个随机变化：一个是由于块引起的，另一个是由于品种引起的。每种变体的均值包括这两种变体的来源；但是这些均值的比较排除了块的差异，因为在同一块上比较了品种。

您只有六个块，因此有5个自由度用于估计块的变化，这说明了变异均值的自由度。比较具有更大的自由度，因为您不必考虑块的变化。

这里要考虑的另一件事是对nlme软件包的支持使用一种包含方法来获得自由度。从本质上讲，这涉及在最坏的情况下查看每种效果的自由度。如果改用lme4包和lmer函数拟合模型，则lsmeans将使用Satterthwaite或Kendall-Roger方法获得自由度，这些结果可能会更大。但是，方法的自由度仍将大大小于比较的自由度。

附录：SAS结果

这是一些具有相同数据和模型的SAS代码：

proc mixed data = Oats;
  class Variety Block;
  model yield = Variety / ddfm = satterth;
  random Block;
  lsmeans Variety / tdiff;

...以及lsmeans结果：

                           Least Squares Means

                                   Standard
Effect     Variety     Estimate       Error      DF    t Value    Pr > |t|
Variety    Golden_R      104.50      7.6809    8.87      13.61      <.0001
Variety    Marvello      109.79      7.6809    8.87      14.29      <.0001
Variety    Victory      97.6250      7.6809    8.87      12.71      <.0001

                      Differences of Least Squares Means
                                           Standard
Effect    Variety    _Variety   Estimate      Error     DF   t Value   Pr > |t|

Variety   Golden_R   Marvello    -5.2917     6.6853     64     -0.79     0.4316
Variety   Golden_R   Victory      6.8750     6.6853     64      1.03     0.3076
Variety   Marvello   Victory     12.1667     6.6853     64      1.82     0.0734

请注意，当使用Satterthwaite方法进行自由度计算时，SAS的比较结果显示为64 df，而平均值显示为8.87 df。

如果在 ddfm语句中省略了 model选项，则默认为df的包含方法，并且在两个表中都列出了64 df。但是，我认为SAS在实施收容措施方面是不正确的。请参阅我在CrossValidated中关于此主题的较早帖子： https://stats.stackexchange.com/questions/140156/degrees-of-freedom-using-containment-method