很抱歉,如果这个答案之前已经回答过,但是这里的大多数答案(例如这里或这里)并没有真正解决我的问题(或者我可能只是没有正确地看到它们是如何做的。我想使用(线性)混合效果对我的数据进行分析。我的设计非常简单:
- 我有两种细菌(M1 和 M2),它们要么与其他 8 种细菌(S1)配对,要么不配对(S2)。
- 在 S1 的情况下,测量一式四份,在 S2 的情况下,测量一式三份。
- 输出变量在多个时间点(重复测量)测量大约 72 小时。
- 然后将规定数量的生物质用于接种新瓶子,依此类推,共进行 4 个循环(SB1 到 4)。因此,连续周期的测量不是独立的(相关)。
- 每个周期的时间点数量并不总是相同(不平衡)
好吧,也许不是那么简单,让我们用一张图来澄清一下:
根据我的一位在混合模型方面经验丰富的同事的说法,这是一个重复测量的分裂情节设计。可悲的是,她是 SAS 用户,我很想使用 R(lme 或 lmer)进行分析。但是,我不是来翻译的(我想这可能是堆栈溢出的问题,尽管如果这里有人可以回答会很好),而是为了了解我的设计是否确实是不平衡重复测量的分裂分裂以及如何我可以构建一个足够的模型来分析它的嵌套、固定和随机效应。
这是我到目前为止的一些虚拟数据(我尝试使用拆分图命名法):
output <- rnorm(266)
mainplot <- c(rep("M1",133),rep("M2",133))
subplot <- c(rep("S1",76),rep("S2",57),
rep("S1",76),rep("S2",57))
subsubplot <-c(rep("SB1",24),rep("SB2",16),rep("SB3",20),rep("SB4",16),
rep("SB1",18),rep("SB2",12),rep("SB3",15),rep("SB4",12),
rep("SB1",24),rep("SB2",16),rep("SB3",20),rep("SB4",16),
rep("SB1",18),rep("SB2",12),rep("SB3",15),rep("SB4",12))
time<-rnorm(266)
mockset <- data.frame(output,mainplot,subplot,subsubplot,time)
mock.lmefit1 <- lme(output ~ 0 + mainplot * subplot * subsubplot,
data = mockset,
random = ~ time|mainplot/subplot/subsubplot)
模拟集真正代表了我的“不平衡重复测量”数据结构。在上面的模型中,我没有考虑子子图(即循环 SB1 到 4)之间的“自相关”,也没有明确纠正重复测量的“不平衡性”。我对我的拦截不感兴趣。
我的同事建议使用以下 SAS 代码:
proc mixed data=rdata;
class mainplot subplot subsubplot time;
model output= mainplot subplot
subsubplot subplot*subsubplot
time subplot*time
subplot*subsubplot*time/ddfm=KENWARDROGER;
random mainplot*subplot mainplot*subplot*subsubplot;
lsmeans subsubplot subplot subplot*subsubplot
time subplot*time subsubplot*time
subplot*subsubplot*time;
run;
如果我做对了*in SAS 类似于 R :,即它显示了术语之间的交互。对我来说,如何将此模型从 SAS 移植到 R 并不完全清楚。我也不太明白为什么这里固定了特定效果而其他效果是随机的。我无法选择加入 R 的另一件事是 SAS 的 ddfm(固定效应测试的分母自由度)选项。这有什么重要性?有谁知道 lme 如何处理这个问题?
mock.lmefitSAS <- lme(output ~ mainplot + subplot +
subsubplot + subplot:subsubplot +
time + subplot:time +
subsubplot:time +
subplot:subsubplot:time,
data = mockset,
random = ~ 1|mainplot:subplot + mainplot:subplot:subsubplot)
这段代码不起作用,因为我猜这不是在 lme 中给出交互项的正确方法(?错误是“组的无效公式”)但如前所述,我不是在这里进行 SAS-to-R 翻译...我只想知道考虑到我的设计以及我的 R 实现缺乏什么,SAS 模型是否以及为什么更正确。