lmer 模型的 lsmeans 结果中的膨胀 DF答案

【问题标题】：Inflated DF in lsmeans results for an lmer modellmer 模型的 lsmeans 结果中的膨胀 DF
【发布时间】：2017-09-15 21:43:27
【问题描述】：

我使用 lme4 包中的 lmer 来运行线性混合效果模型。我有未处理 (5) 和处理过的地块 (10) 的 3 年温度数据。型号：

modela<-lmer(ave~yr*tr+(1|pl), REML=FALSE, data=mydata)

模型检查残差的正态性； qq范数图我的数据：

'data.frame':   6966 obs. of  7 variables:

$ yr : Factor w/ 3 levels "yr1","yr2","yr3": 1 1 1 1 1 1 1 1 1 1 ...

$ pl : Factor w/ 15 levels "C02","C03","C05",..: 1 1 1 1 1 1 1 1 1 1 ...

$ tr : Factor w/ 2 levels "Cont","OTC": 1 1 1 1 1 1 1 1 1 1 ...

$ ave: num  14.8 16.1 11.6 10.3 11.6 ...

交互很重要，所以我用了lsmeans：

lsmeans(modela, pairwise~yr*tr, adjust="tukey")

在对比中，我得到（仅摘录）

contrast                estimate        SE      df t.ratio p.value

yr1,Cont - yr2,Cont -0.727102895 0.2731808 6947.24  -2.662  0.0832

yr1,OTC - yr2,OTC   -0.990574030 0.2015650 6449.10  -4.914  <.0001

yr1,Cont - yr1,OTC  -0.005312771 0.3889335   31.89  -0.014  1.0000

yr2,Cont - yr2,OTC  -0.268783907 0.3929332   32.97  -0.684  0.9825

我的问题是关于某些对比的高 dfs，以及相关但无意义的低 p 值。

这可能是由于：

-在我的数据集中存在 NA（移除后有一些改进）

-样本量不等（例如，一种治疗方法中的 5 个，另一种治疗方法中的 10 个 - 然而，那些 (yr1,Cont - yr1, OTC) 似乎不是问题。

其他问题？

我搜索了 stakoverflow 问题，并进行了交叉验证。

感谢任何答案、想法、cmets。

【问题讨论】：

如果yr 是in-pl 因子，那么这些df 可能是正确的。请注意，对于lmer 型号，lsmeans 提供了两个 d.f. 的选择。方法——一种是基于 lmerTest 包中例程的 Satterthwaite 方法（默认），以及使用 pbkrtest 包中例程的 Kenward-Roger 方法。尝试将参数 lmer.df = "k" 添加到 lsmeans 调用中，看看后一种方法的结果是什么。如果它们具有可比性，我相信这里没有问题。不相关的注释——我想知道你为什么选择REML=FALSE 来拟合模型。
顺便说一下，P 值不会比显示的值大很多，除非 d.f.真的很小。
从这个问题的上下文来看，我猜如果你为每个图分别绘制残差与运行顺序（天？？），你会发现自己对模型的适用性不太满意.我进一步想象每个情节都存在大量的序列相关性。
谢谢@rvl。我选择 REML 是因为我使用最大似然选择了我最好的模型（与其他 R 用户讨论的结果，以及在我开始分析时学习 R）。是的，我都试过了。方法和结果非常相似。如何解释 dfs 的巨大差异？年内治疗30-40，治疗年内6000？这些是我的样本量：Yr1,Cont n=810, Yr1,OTC n=1619, Yr2,Cont n=809, Yr2,OTC n=1458, Yr3,Cont n=809, Yr3,OTC n=1457。 OTC/Cont 是治疗方法。
因为地块是处理的实验单位。它来自模型中的 (1|pl) 部分。

标签： r dataframe lme4 lsmeans

【解决方案1】：

在此示例中，处理是通过实验分配给地块的。分配给处理的少量地块严重限制了可用于统计比较处理的信息。（如果每个处理只有一个地块，甚至无法比较处理，因为您无法从地块的效果中区分出处理的效果。）您有 10 个地块分配给一个治疗和5到另一个。就治疗的主效应而言，您因此有 (10-1)+(5-1) = 13 d.f.对于治疗的主要效果，如果你这样做了

lsmeans(modela, pairwise ~ tr)

您将看到大约 13 d.f. （可能由于不平衡和缺失而减少）这些统计数据。当您比较年份和治疗的组合时，您会得到大约 3 倍的 d.f.因为有3年。然而，在其中一些比较中，被比较的每个组合中的年份都是相同的，并且在这些比较中，地块的变化大部分被抵消了（这是一个地块内的比较）；在这些情况下，d.f.基本上来自模型的残差，它有数千个 d.f.由于数据的不平衡，这些比较有点受到图间变化的污染，使得 d.f.略小于残差 d.f.

您似乎对诸如treat1、year1 与treat2、year3 之类的交叉比较并不特别感兴趣。我建议使用“按”变量来减少测试比较的数量，因为当你全部测试它们时，多重性校正是不必要的保守。它会是这样的：

modela.lsm = lsmeans(modela, ~ tr * yr)
pairs(modela.lsm, by = "yr")   # compare tr for each yr
pairs(modela.lsm, by = "tr")   # compare yr for each tr

这些调用会将 Tukey 校正单独应用于每个“按”组。如果您想对整个家庭进行多重校正，请执行以下操作：

rbind(pairs(modela.lsm, by = "yr"))
rbind(pairs(modela.lsm, by = "tr"))

默认情况下，使用多元 t 校正（Tukey 在这里不是正确的方法）。你甚至可以做

rbind(pairs(modela.lsm, by = "yr"), pairs(modela.lsm, by = "tr"))

将所有比较归为一个族并应用多元 t 调整。

【讨论】：

谢谢，这很有帮助。我认为使用 lsmeans (model, ~ tr|yr) 也可以做到这一点。你说得对，我对年份比较不感兴趣。