在我的实验中,我剪下植物并测量它们的反应,例如在季节结束时产生的叶片质量。我操纵了剪裁强度和剪裁时间,并跨越了这两种治疗方法。我还包括了一个控制剪裁处理,产生了 5 种不同的剪裁处理组合。每次处理 12 株植物,我在两年内跟踪了总共 60 株植物。也就是说,我在第 1 年收集了这 60 株植物的测量值,并在第 2 年再次收集了相同的植物。
这是我的设计,其中时间下的“从不”和强度下的“零”任意替换了“控制”处理:
Year Timing intensity treatments
2015 early high early-high
2015 early low early-low
2015 late high late-high
2015 late low late-low
2015 never zero control
2014 early high early-high
2014 early low early-low
2014 late high late-high
2014 late low late-low
2014 never zero control
我按照 Ben Bolker 的一项建议忽略了运行 lme4 的警告,然后对模型进行了 F 检验(R- 分析重复测量不平衡设计与 lme4?):
m1<-lmer(log(plant.leaf.g)~timing*intensity*year+(1|id), data=cmv)
Anova(m1, type="III", test="F")
方差分析输出给了我时间和强度之间的显着相互作用(p = 0.006),然后我使用以下方法进行了多重比较测试:
cmv$SHD<-interaction(cmv$timing, cmv$intensity)
m2<-lmer(log(plant.leaf.g)~-1+SHD+(1|id),data=cmv, na.action=na.exclude)
summary(glht(m2, linfct=mcp(SHD="Tukey")))
这是我输出的剪辑,其中 p=0.08 的唯一重要对:
Estimate Std. Error z value Pr(>|z|)
late.2014 - early.2014 == 0 -0.6584 0.3448 -1.910 0.3844
never.2014 - early.2014 == 0 0.1450 0.4102 0.354 0.9992
early.2015 - early.2014 == 0 -0.4906 0.2786 -1.761 0.4788
late.2015 - early.2014 == 0 -0.1687 0.3494 -0.483 0.9965
never.2015 - early.2014 == 0 0.4201 0.4079 1.030 0.9032
never.2014 - late.2014 == 0 0.8034 0.4119 1.951 0.3597
early.2015 - late.2014 == 0 0.1678 0.3419 0.491 0.9963
late.2015 - late.2014 == 0 0.4897 0.2724 1.797 0.4553
never.2015 - late.2014 == 0 1.0785 0.4119 2.618 0.0885 .
early.2015 - never.2014 == 0 -0.6356 0.4074 -1.560 0.6133
为什么 Anova 认为时间*强度非常显着,但在我的多重比较测试中没有显着性?还有另一种方法我应该进行多重比较吗?
在其他多重比较输出中,我得到的 p 值高达 1.00000,这正常吗?
data<-structure(list(id = c(91L, 92L, 93L, 94L, 95L, 96L, 97L, 98L,
99L, 100L, 101L, 102L, 103L, 105L, 106L, 107L, 108L, 109L, 110L,
111L, 112L, 113L, 114L, 115L, 116L, 117L, 119L, 120L, 121L, 122L,
123L, 124L, 125L, 126L, 127L, 128L, 129L, 130L, 131L, 132L, 133L,
134L, 135L, 136L, 137L, 138L, 139L, 140L, 141L, 142L, 143L, 144L,
146L, 147L, 148L, 149L, 150L, 91L, 92L, 93L, 94L, 95L, 96L, 97L,
98L, 99L, 100L, 101L, 102L, 103L, 105L, 106L, 107L, 108L, 109L,
110L, 111L, 112L, 113L, 114L, 115L, 116L, 117L, 119L, 120L, 121L,
122L, 123L, 124L, 125L, 126L, 127L, 128L, 129L, 130L, 131L, 132L,
133L, 134L, 135L, 136L, 137L, 138L, 139L, 140L, 141L, 142L, 143L,
144L, 146L, 147L, 148L, 149L, 150L), quad = c(2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), year = c(2015L,
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L,
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L,
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L,
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L,
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L,
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L,
2015L, 2015L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L,
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L,
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L,
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L,
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L,
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L,
2014L, 2014L, 2014L, 2014L, 2014L), timing = structure(c(1L,
3L, 2L, 1L, 1L, 2L, 3L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 3L, 1L, 1L,
3L, 2L, 3L, 1L, 3L, 2L, 3L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L,
1L, 3L, 3L, 2L, 2L, 1L, 2L, 3L, 2L, 1L, 1L, 2L, 2L, 3L, 1L, 2L,
2L, 2L, 1L, 1L, 2L, 1L, 1L, 3L, 1L, 3L, 2L, 1L, 1L, 2L, 3L, 2L,
2L, 2L, 2L, 1L, 1L, 2L, 3L, 1L, 1L, 3L, 2L, 3L, 1L, 3L, 2L, 3L,
1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 2L,
3L, 2L, 1L, 1L, 2L, 2L, 3L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 3L,
1L), .Label = c("early", "late", "never"), class = "factor"),
intensity = structure(c(2L, 3L, 1L, 2L, 1L, 2L, 3L, 1L, 2L,
1L, 2L, 1L, 1L, 2L, 3L, 2L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L,
1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 3L, 3L, 2L, 1L, 1L,
1L, 3L, 1L, 1L, 2L, 2L, 1L, 3L, 2L, 1L, 2L, 1L, 1L, 2L, 2L,
2L, 1L, 3L, 2L, 3L, 1L, 2L, 1L, 2L, 3L, 1L, 2L, 1L, 2L, 1L,
1L, 2L, 3L, 2L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 2L, 1L,
1L, 2L, 2L, 2L, 2L, 1L, 2L, 3L, 3L, 2L, 1L, 1L, 1L, 3L, 1L,
1L, 2L, 2L, 1L, 3L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 3L, 2L
), .Label = c("high", "low", "zero"), class = "factor"),
treatment = structure(c(3L, 1L, 4L, 3L, 2L, 5L, 1L, 4L, 5L,
4L, 5L, 2L, 2L, 5L, 1L, 3L, 2L, 1L, 4L, 1L, 2L, 1L, 4L, 1L,
2L, 3L, 2L, 4L, 3L, 5L, 5L, 3L, 2L, 3L, 1L, 1L, 5L, 4L, 2L,
4L, 1L, 4L, 2L, 3L, 5L, 4L, 1L, 3L, 4L, 5L, 4L, 2L, 3L, 5L,
3L, 2L, 1L, 3L, 1L, 4L, 3L, 2L, 5L, 1L, 4L, 5L, 4L, 5L, 2L,
2L, 5L, 1L, 3L, 2L, 1L, 4L, 1L, 2L, 1L, 4L, 1L, 2L, 3L, 2L,
4L, 3L, 5L, 5L, 3L, 2L, 3L, 1L, 1L, 5L, 4L, 2L, 4L, 1L, 4L,
2L, 3L, 5L, 4L, 1L, 3L, 4L, 5L, 4L, 4L, 3L, 5L, 2L, 1L, 3L
), .Label = c("control", "early-high", "early-low", "late-high",
"late-low"), class = "factor"), plant.leaf.g = c(846.216,
382.704, 2393.088, 61.832, 1315.86, 275.816, 3705.862, 3500.52,
67.482, 432, 487.492, 1228.618, 776.16, 1575, 735.9, 2417.75,
1342.92, 2359.046, 686.726, 1385.856, 343.684, 2277.312,
465.528, 2314.584, 508.4, 1243.644, 1064.448, 1020.646, NA,
494.832, 1318.248, 1516.4, 1271.218, 512.512, 157.878, 3753.992,
586.032, 1042.176, 889.632, 651.052, 498.042, 625.872, 16.28,
497.51, 593.75, 706.84, 2238.742, 232.584, 671.532, 90.72,
1412.442, 902.728, 3077.184, 619.106, 0.576, 400.452, 684.522,
849.852, 152.76, 1280.448, 274.47, 387.614, 98.496, 2304.504,
644.952, 35.392, 250.56, 267.33, 2212.08, 2392.596, 751.944,
629.418, 731.544, 1013.196, 1516.4, 130.536, 2910.6, 554.4,
2163.35, 223.86, 2369.376, 551.976, 985.6, 1482.24, 815.386,
1664.132, 596.376, 1581.432, 217.128, 1041.656, 951.168,
256.172, 1587.148, 359.448, 546.48, 1226.544, 371.64, 293.504,
177.726, 343.26, 691.24, 207.604, 588.924, 1405.258, 136.17,
451.432, 576.18, 424.804, 884.534, 2466.45, 1524.432, 973.208,
369.474, 410.048)), .Names = c("id", "quad", "year", "timing",
"intensity", "treatment", "plant.leaf.g"), class = "data.frame", row.names = c(NA,
-114L))
PS。我一辈子都无法使用 lsmeans 来处理这种不平衡的设计。输出中报告了很多 NA。