R 中的 Lsmeans 包 - lme 模型的自由度
Lsmeans package in R - degrees of freedom with lme models
我对 lsmeans 包在使用 nlme 包构建的线性混合模型中使用的自由度有疑问。
这里有一个例子来说明我基于 Oats 数据集的问题。我不想讨论这个模型是否与给定的数据集相关，我只是想重现我在另一个数据集上遇到的问题 ;-)。
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。
现在有谁知道为什么会发生这种情况，是否可以改变它？还是我遗漏了什么？
更新 - 初始错误消息
我最初注意到在一个特定情况下 lsmeans 的自由度降低，其中我的随机 运行 效果只有 2 个水平，并且当我对 Dunnett 的调整感兴趣时。由于我对对比比对 lsmeans 更感兴趣，既然我了解它的来源，我仍然可以使用它，但我把它放在那里以防万一有人有同样的错误并想知道为什么。
我在下面用 Oats 数据示例复制了它。我得到的错误发生在 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。
如果在model语句中省略ddfm选项，则默认为df的包含方法，在两个表中都列出了64个df。但是，我认为 SAS 在实施遏制方面是不正确的；请参阅我之前在 CrossValidated 中关于此主题的 post：https://stats.stackexchange.com/questions/140156/degrees-of-freedom-using-containment-method
我对 lsmeans 包在使用 nlme 包构建的线性混合模型中使用的自由度有疑问。
这里有一个例子来说明我基于 Oats 数据集的问题。我不想讨论这个模型是否与给定的数据集相关，我只是想重现我在另一个数据集上遇到的问题 ;-)。
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。
现在有谁知道为什么会发生这种情况，是否可以改变它？还是我遗漏了什么？
更新 - 初始错误消息
我最初注意到在一个特定情况下 lsmeans 的自由度降低，其中我的随机 运行 效果只有 2 个水平，并且当我对 Dunnett 的调整感兴趣时。由于我对对比比对 lsmeans 更感兴趣，既然我了解它的来源，我仍然可以使用它，但我把它放在那里以防万一有人有同样的错误并想知道为什么。
我在下面用 Oats 数据示例复制了它。我得到的错误发生在 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。
如果在model语句中省略ddfm选项，则默认为df的包含方法，在两个表中都列出了64个df。但是，我认为 SAS 在实施遏制方面是不正确的；请参阅我之前在 CrossValidated 中关于此主题的 post：https://stats.stackexchange.com/questions/140156/degrees-of-freedom-using-containment-method