如何在 R 中对 2x2x2 设计进行纵向分析？

Question

我有一个数据集，其中包含以下处理变量（营养、肥料），它记录了水中藻类随时间的生长情况（t0、t1...t10）。在与标记为“氮”的肥料系列中，在t5天之后添加氮。在标记为“无”的系列中，不添加氮。

nutrition <- c("good","good","bad","bad","good","good","bad","bad","good","good","bad","bad","good","good","bad","bad")

fertlizer <- c("none", "nitrogen","none","nitrogen","none", "nitrogen","none","nitrogen","none", "nitrogen","none","nitrogen","none", "nitrogen","none","nitrogen")

t0 <- c(7,  6,  3, 20, 13,  4, 14,  9, 15,  5, 18, 19,  8,  1, 10, 16)
t1 <- c(12,  9,  3, 20,  4,  7,  6, 17, 19,  5, 18,  8, 15, 16, 10,  2)
t2 <- c(12,  9,  3, 20,  4,  7,  6, 17,7,  6,  3, 20, 13,  4, 14,  9)
t3 <- c(15,  5, 18, 19,  8,  1, 10, 16,4,  7,  6, 17,7,  6,  3, 20)
t4 <- c(6,7,12,4,7,18,9,10,2,10,11,14,15,1,14,16)
t5 <- c(4,  7,  6, 17,7,  6,  3, 20,15,  5, 18, 19,  8,  1, 10, 16)
t6 <- c(70,5,16,31,61,14,22,23,80,13,24,32,90,16,28,29)
t7 <- c(56,16,7,8,78,26,28,30,91,5,8,19,67,16,18,19)
t8 <- c(88,21,20,19,90,16,18,19,57,3, 20,  4,  7,67,13,12)
t9 <- c(62,12,15,27,71,20,  4,  7,72,6,  3, 20,73,14,  9, 15)
t10 <- c(40,13,7,19,50,3, 20, 7,66,14,  9, 15,80,16,18,19)


replicates <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16)

data <- data.frame(nutrition, fertlizer,replicates, t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)

data$nutrition <- as.factor(data$nutrition)
data$fertlizer <- as.factor(data$fertlizer)

我想比较各组之间的坡度，看看施肥干预后坡度是否发生变化。我想不使用肥料作为对照（即（好，无）或（坏，无）作为对照）

我将此数据转换为长格式，列标题为“营养”、“肥料”、“时间”、“重复”和“生长” 我创建了一个名为添加的新列，以区分 t5 之前和 t5 之后的时间段。 t5 之前的时间 -> 0，t5 之后的时间 ->1

   nutrition fertilizer time replicate growth addition
   good      none      t0   1         6        0
   good      none      t1   1         7        0
   ..
   ..
   good      none      t5   1         3        1

我运行以下纵向分析，每列具有以下结构：

nutrition: factor with 2 levels
fertlizer: factor with2 levels
time: factor with 10 levels
replicates: num 0,1,2,3...
growth: num 6, 7, 5 ...
addition: factor with 2 levels

lmer(增长 ~ 营养 + 肥料 + 时间 + 添加 + (1|replicates))

我收到一条错误消息，指出固定效应模型排名不足，因此删除了 x 个列。反正有这个问题吗？有没有改进模型编写方式的建议？

Answer 1

我对此不是 100% 确定，但我相信它应该能让你大致了解一下。

dtf <- structure(list(nutr = structure(c(2L, 2L, 1L, 1L, 2L, 2L, 
1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L), .Label = c("bad", "good"
), class = "factor"), fert = structure(c(2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L), .Label = c("N", 
"0"), class = "factor"), id = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 
10, 11, 12, 13, 14, 15, 16), t0 = c(7, 6, 3, 20, 13, 4, 14, 9, 
15, 5, 18, 19, 8, 1, 10, 16), t1 = c(12, 9, 3, 20, 4, 7, 6, 17, 
19, 5, 18, 8, 15, 16, 10, 2), t2 = c(12, 9, 3, 20, 4, 7, 6, 17, 
7, 6, 3, 20, 13, 4, 14, 9), t3 = c(15, 5, 18, 19, 8, 1, 10, 16, 
4, 7, 6, 17, 7, 6, 3, 20), t4 = c(6, 7, 12, 4, 7, 18, 9, 10, 
2, 10, 11, 14, 15, 1, 14, 16), t5 = c(4, 7, 6, 17, 7, 6, 3, 20, 
15, 5, 18, 19, 8, 1, 10, 16), t6 = c(70, 5, 16, 31, 61, 14, 22, 
23, 80, 13, 24, 32, 90, 16, 28, 29), t7 = c(56, 16, 7, 8, 78, 
26, 28, 30, 91, 5, 8, 19, 67, 16, 18, 19), t8 = c(88, 21, 20, 
19, 90, 16, 18, 19, 57, 3, 20, 4, 7, 67, 13, 12), t9 = c(62, 
12, 15, 27, 71, 20, 4, 7, 72, 6, 3, 20, 73, 14, 9, 15), t10 = c(40, 
13, 7, 19, 50, 3, 20, 7, 66, 14, 9, 15, 80, 16, 18, 19)), .Names = c("nutr", 
"fert", "id", "t0", "t1", "t2", "t3", "t4", "t5", "t6", "t7", 
"t8", "t9", "t10"), row.names = c(NA, -16L), class = "data.frame")

# Need to reshape into long format so that each column is a separate variable
library(reshape2)

dtf.long <- reshape2::melt(dtf, id.vars=1:3, variable.name="time")
dtf.long$time <- as.integer(sub("t", "", dtf.long$time))
dtf.long$fert2 <- dtf.long$time > 5 & dtf.long$fert == "N"

library(lattice)
xyplot(value ~ time | nutr * fert, data=dtf.long )

library(lme4)
m1.1 <- lmer(value ~ nutr * fert2 * time + (1 | id), dtf.long, REML=FALSE)
m1.2 <- lmer(value ~ nutr * fert2 * time + (1 + time | id), dtf.long, REML=FALSE)

# The random slope term doesn't appear to be adding anything of value
anova(m1.1, m1.2)

These slideslme4 中的纵向建模可能有用。

Answer 2

据我所知，该错误消息意味着您的预测器没有足够的信息来运行您想要的模型。也许您应该考虑使用带有块的 ANOVA 或 ANCOVA？考虑到数据的实验性质，走这条路是完全可以的。您还应该使用您的治疗组创建一个新列。比如有肥好营养，没肥好营养可以叫GWF和GWOF，营养不好的组合也一样。这将为您提供 4 次治疗，您的障碍可能是“时间”。

发布您的数据集样本将有助于解决您的问题