R-用lme4分析重复测量不平衡设计？

Question

在我的实验中，我修剪了植物并测量了它们在季节结束时的反应，例如产生的叶片质量。我操纵了剪裁强度和剪裁时间，并跨越了这两种治疗方法。我还包括了一个控制剪裁处理，产生了 5 种不同的剪裁处理组合。每次处理 12 株植物，我在两年内跟踪了总共 60 株植物。也就是说，我在第 1 年收集了这 60 株植物的测量值，并在第 2 年再次收集了相同的植物。

最简单的方法是分别分析 5 种不同的治疗方法。然而，我想获得时间和强度的影响以及它们的相互作用，但是由于控制处理没有与时间或强度完全交叉，这使得我的实验设计不平衡并且在统计上很棘手。为了使这更复杂一点，我也想在我的模型中加入年份的影响。

理想情况下，我可以使用 lme4 执行此操作，之后使用 lsmeans 包进行多重比较变得轻而易举。

当我尝试运行我的模型时

     m1<-lmer(log(plant.leaf.g+1)~timing*intensity*year+(1|id), data=cmv) #not significant

我收到警告“固定效应模型矩阵秩不足，因此删除 8 列/系数”。

有谁知道我可以让这种不平衡的混合模型与 lme4 一起使用吗？

这是我的数据的一个子集，其中“从不”在时间和强度下“零”任意替换了“控制”治疗：

id  year    timing  intensity   treatment   plant.leaf.g
91  2015    early   low early-low   315.944
92  2015    never   zero    control 99.28
93  2015    late    high    late-high   663.936
94  2015    early   low early-low   25.488
95  2015    early   high    early-high  453.57
96  2015    late    low late-low    90.804
97  2015    never   zero    control 1312.098
98  2015    late    high    late-high   959.82
99  2015    late    low late-low    28.014
100 2015    late    high    late-high   178.56
91  2014    early   low early-low   289.14
92  2014    never   zero    control 61.774
93  2014    late    high    late-high   639.936
94  2014    early   low early-low   138.39
95  2014    early   high    early-high  168.216
96  2014    late    low late-low    51.008
97  2014    never   zero    control 966.112
98  2014    late    high    late-high   279.048
99  2014    late    low late-low    23.936
100 2014    late    high    late-high   169.344

cmv<-structure(list(id = c(91L, 92L, 93L, 94L, 95L, 96L, 97L, 98L, 
99L, 100L, 101L, 102L, 103L, 105L, 106L, 107L, 108L, 109L, 110L, 
91L, 92L, 93L, 94L, 95L, 96L, 97L, 98L, 99L, 100L, 101L, 102L, 
103L, 104L, 105L, 106L, 107L, 108L, 109L, 110L), year = c(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), timing = structure(c(1L, 3L, 2L, 1L, 1L, 2L, 3L, 
2L, 2L, 2L, 2L, 1L, 1L, 2L, 3L, 1L, 1L, 3L, 2L, 1L, 3L, 2L, 1L, 
1L, 2L, 3L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 3L, 1L, 1L, 3L, 2L
), .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, 2L, 3L, 1L, 2L, 1L, 2L, 3L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 
2L, 3L, 2L, 1L, 3L, 1L), .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, 3L, 1L, 4L, 3L, 2L, 
5L, 1L, 4L, 5L, 4L, 5L, 2L, 2L, 5L, 5L, 1L, 3L, 2L, 1L, 4L
), .Label = c("control", "early-high", "early-low", "late-high", 
"late-low"), class = "factor"), plant.stem.g = c(315.944, 
99.28, 663.936, 25.488, 453.57, 90.804, 1312.098, 959.82, 
28.014, 178.56, 158.12, 387.528, 288.75, 327.348, 770.44, 
835.05, 457.188, 942.002, 229.194, 289.14, 61.774, 639.936, 
138.39, 168.216, 51.008, 966.112, 279.048, 23.936, 169.344, 
154.14, 703.04, 836.4, 511.92, 463.524, 245.226, 267.41, 
439.392, 714.85, 68.012)), .Names = c("id", "year", "timing", 
"intensity", "treatment", "plant.stem.g"), class = "data.frame", row.names =     c(NA, 
-39L))

注意：我已经让 m1=aov(plant.leaf.g~intensity*timing*year+Error(id), data=cmv) 运行，但我读到我应该使用 car 包中的 Anova type="3" 函数来获取我的 p 值，但我无法做到这与 Error(id) 术语。我也无法与TukeyHSD 函数或multcomp 包进行多重比较。

Answer 1

本质上并没有错

 m1<-lmer(log(plant.leaf.g+1)~timing*intensity*year+(1|id), 
          data=cmv)

（除了其中带有零的对数转换数据很棘手；您确定加 1 是正确的吗？只有叶质量是无单位的才有意义。您可以考虑添加 min(plant.leaf.g[plant.leaf.g>0])/2 代替...）

出现警告（不是错误）是因为您的数据集中没有时间、强度和年份的所有组合，但您要求 R 估计每个组合的参数。一些合理的选择是：

• 忽略警告（在比较每个因素的整体影响时，您可能会得到合理的答案）
• 降低模型的复杂性，特别是通过消除三向交互（即使用(timing+intensity+year)^2）（我假设这会起作用，但如果有组合，您可能需要进一步简化模型数据中缺失的时间和强度）
• 从 3 向交互作用构建单向 ANOVA，例如cmv$int <- with(cmv,interaction(timing,intensity,year,drop=TRUE))（但你将无法分离主要效果和交互）