我有一个关于在线性混合模型上运行事后测试的问题:
我正在运行一个线性混合模型,分为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()
吗?
感谢您的意见,我知道似然比测试本身是有争议的。我已经考虑过混合效应方差分析,但是有一些缺失的数据点使得这不可能。数据之前被另一名学生分析为双向方差分析,没有重复测量,但我的感觉是这不合适,因为每条蛇都是在重复的时间点测量的