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

Question

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

我正在 lme4 中运行一个线性混合模型，分为 3 组，每组 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 会对成对测试的含义感到困惑。 (关于您是否真的想用连续预测器来分析这个模型，即作为响应曲面设计而不是双向方差分析，需要进行一段稍长的对话......)这是您的分析的稍微更紧凑的版本:

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 维护者提及这一点，看看是否可以检测并标记这个(确实非常常见)错误。