我有一个包含以下处理变量(营养、肥料)的数据集,它记录了水中藻类随时间(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)
我想比较各组之间的坡度,看看在施肥干预后坡度是否发生变化。我不想使用肥料作为对照(即(好,无)或(坏,无)作为对照)
我将这些数据转换为长格式,列标题为“营养”、“肥料”、“时间”、“复制”和“生长”
我创建了一个名为addition 的新列来区分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|重复))
我收到一条错误消息,说固定效应模型排名不足,因此删除了 x 列数。反正有这个问题吗?是否有改进模型编写方式的建议?
最佳答案
我不是 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 slides 在
lme4
中纵向建模可能有用。关于r - 如何对 R 中的 2x2x2 设计进行纵向分析?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/47739844/