r - 事后测试线性混合模型 lsmeans 错误

Question

我有一个关于在线性混合模型上运行事后测试的问题：

我正在运行一个线性混合模型，分为lme43 组，每组 5 条蛇，每组以不同的通风率 ( Vent)，在不同时间点进行测量 ( Time)，将蛇指定为随机效应 ( ID)

下面的数据子集：

subset1 <- structure(list(ID = structure(c(5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 
6L, 6L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 
9L, 18L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 19L, 10L, 10L, 
10L, 10L, 10L, 20L, 20L, 20L, 20L, 20L, 4L, 4L, 4L, 4L, 4L, 11L, 
11L, 11L, 11L, 11L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 12L, 
13L, 14L, 15L, 16L, 17L, 17L, 17L, 17L, 17L), .Label = c("", 
"1_1_2", "10", "10_1_1", "13_1_4", "14_2_4", "15_3_4", "16_1_4", 
"17_2_4", "2_2_1", "5", "5_2_2", "5_2_3", "5_2_4", "5_2_5", "5_2_6", 
"7_1_2", "8", "9", "9_3_1"), class = "factor"), Vent = c(30L, 
30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 
30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 125L, 
125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 
125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 125L, 
125L, 125L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 
250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 250L, 
250L, 250L, 250L, 250L, 250L), Time = c(60L, 80L, 180L, 720L, 
1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 
60L, 80L, 180L, 720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 
80L, 180L, 720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 80L, 
180L, 720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 80L, 180L, 
720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 80L, 180L, 720L, 
1440L, 60L, 80L, 180L, 720L, 1440L, 60L, 80L, 180L, 720L, 1440L, 
60L, 80L, 180L, 720L, 1440L), corr.pO2 = c(224.1388673, 233.9456874, 
239.1553778, 107.2373336, 76.71835625, 164.6293748, 243.8501858, 
234.8205544, 71.74240501, 62.23789874, 69.69478654, 62.23789874, 
152.1142885, 79.61325688, 63.33285001, 240.8713061, 231.304842, 
222.7743953, 95.7912966, 64.41744793, 241.7255035, 238.2936023, 
138.1188987, 43.00663696, 50.64392111, 265.4973967, 274.0599252, 
285.0144919, 83.37647392, NA, 292.3660214, 281.6533627, 275.9747984, 
63.33285001, 56.59660394, 254.2521631, 222.3180596, 208.736288, 
88.83223104, 114.1782867, 208.255285, 232.1878564, 193.3861802, 
72.75355024, 60.01517133, 209.6956308, 245.9596884, 200.4342522, 
75.73874562, 67.61194011, 240.0146049, 261.1278627, 166.9318704, 
74.75152919, 73.75652657, 270.9724687, 251.7882317, 245.9596884, 
147.1396383, 50.64392111, 294.179467, 296.3431178, 284.6426934, 
73.75652657, 75.73874562, 233.0681297, 234.3834557, 143.3247511, 
73.75652657, 66.55672391, 245.9596884, 249.3041163, 223.6847954, 
92.35383362, 78.65544784)), .Names = c("ID", "Vent", "Time", 
"corr.pO2"), row.names = c(NA, 75L), class = "data.frame")

代码：

attach(subset1)

require(lme4)

with.vent = lmer(corr.pO2 ~ Vent * Time + (1|ID),REML = FALSE, data = subset1)

with.vent.no.int = lmer(corr.pO2 ~ Vent + Time + (1|ID),REML = FALSE, data = subset1)

anova(with.vent, with.vent.no.int)
#no significant interaction

通风效果测试：

without.vent = lmer(corr.pO2 ~ Time + (1|ID), REML = FALSE, data = subset1)

与通风口比较：

anova(with.vent.no.int, without.vent)
#no significant effect of ventilation treatment p=0.09199

测试时间的影响：

without.time = lmer(corr.pO2 ~ Vent + (1|ID), data = subset1)

anova(with.vent.no.int, without.time)
# highly significant effect of time on pO2 < 2.2e-16 *** 

所以尝试事后测试：

require(lsmeans)
lsmeans(with.vent.no.int, pairwise ~ Time, adjust = "tukey", data = subset1)

这是我得到错误的地方：

Error in solve.default(L %*% V0 %*% t(L), L) : 
  Lapack routine dgesv: system is exactly singular: U[1,1] = 0

我可以使用以下方法进行校正并运行成对测试：

pairwise.t.test(corr.pO2, Time, p.adj = "BH", paired = T)

但是知道这在其他变量有相互作用的情况下不起作用（就像我的其他数据一样），并且我希望在每个时间点和通风制度之间进行成对比较。这可能lsmeans()吗？

感谢您的意见，我知道似然比测试本身是有争议的。我已经考虑过混合效应方差分析，但是有一些缺失的数据点使得这不可能。数据之前被另一名学生分析为双向方差分析，没有重复测量，但我的感觉是这不合适，因为每条蛇都是在重复的时间点测量的

Answer 1

答案很简单：你需要确保你的Vent和Time预测的都是因数。否则lsmeans会对成对测试的含义感到困惑。（关于您是否真的想用连续预测器分析这个模型，即作为响应面设计而不是 2 路方差分析，还有一个稍微长一点的对话……）这是您分析的稍微更紧凑的版本：

subset1 <- transform(subset1,Vent=factor(Vent), Time=factor(Time))
require(lme4)
with.vent = lmer(corr.pO2 ~ Vent * Time + (1|ID),
       REML = FALSE, data = subset1)
drop1(with.vent,test="Chisq")  ## test interaction
with.vent.no.int = update(with.vent, . ~ . - Vent:Time)
drop1(with.vent.no.int,test="Chisq")  ## test main effects
require(lsmeans)
lsmeans(with.vent.no.int, pairwise ~ Time)

输出子集：

$contrasts
 contrast    estimate       SE    df t.ratio p.value
 60 - 80     -6.99222 12.76886 63.45  -0.548  0.9819
 60 - 180    14.74281 12.76886 63.45   1.155  0.7768
 60 - 720   147.27139 12.76886 63.45  11.534  <.0001
...

我同意错误信息是难以理解的。可能值得向lsmeans维护者提及这一点，看看是否有可能检测和标记这个（非常常见的）错误。

Answer 2

lsmeans 的下一次更新（可能在 2016 年 2 月 1 日左右）将捕获这种错误：

> lsmeans(with.vent.no.int, pairwise ~ Vent)

$lsmeans
     Vent   lsmean       SE    df lower.CL upper.CL
 135.1351 167.4871 6.859275 18.63 153.1111  181.863

Confidence level used: 0.95 

$contrasts


Warning message:
In contrast.ref.grid(result, method = contr, by, ...) :
  No contrasts were generated! Perhaps only one lsmean is involved.
  This can happen, for example, when your predictors are not factors.

该ref.grid功能可以方便地了解您所拥有的：

> ref.grid(with.vent.no.int)
'ref.grid' object with variables:
    Vent = 135.14
    Time = 483.24

Vent和都是Time协变量，因此默认使用它们的平均值。要改变这一点，你不一定要改变数据集；您可以将预测变量强制为模型中的因子：

> repaired = lmer(corr.pO2 ~ factor(Vent) + factor(Time) + (1|ID), 
                  REML = FALSE, data = subset1)
> ref.grid(repaired)
'ref.grid' object with variables:
    Vent =  30, 125, 250
    Time =   60,   80,  180,  720, 1440

> lsmeans(repaired, pairwise ~ Vent)
$lsmeans
 Vent   lsmean       SE    df lower.CL upper.CL
   30 146.0967 12.19373 18.16 120.4952 171.6981
  125 177.0917 12.29417 18.66 151.3274 202.8559
  250 173.2568 11.12879 26.72 150.4111 196.1024

Results are averaged over the levels of: Time 
Confidence level used: 0.95 

$contrasts
 contrast    estimate       SE    df t.ratio p.value
 30 - 125  -30.994975 17.31570 18.41  -1.790  0.2005
 30 - 250  -27.160077 16.50870 21.52  -1.645  0.2490
 125 - 250   3.834898 16.58302 21.81   0.231  0.9710

Results are averaged over the levels of: Time 
P value adjustment: tukey method for comparing a family of 3 estimates