我对在使用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
最佳答案
该模型表示响应变量有两个随机变化:一个是由于块引起的,另一个是由于品种引起的。每种变体的均值包括这两种变体的来源;但是这些均值的比较排除了块的差异,因为在同一块上比较了品种。
您只有六个块,因此有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关于r - LME模型在R-自由度中的Lsmeans软件包,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/42815322/