線形混合モデルでの事後テストの実行について質問があります。
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()ですか？
ご意見ありがとうございます。尤度比検定自体が物議をかもしていることは承知しています。混合効果の ANOVA を検討しましたが、これを不可能にするいくつかの欠落データ ポイントがあります。データは以前に別の学生によって二元配置分散分析として分析され、繰り返し測定はありませんでしたが、各ヘビが繰り返し測定されたため、これは不適切であると感じました