我正在研究 2x2 交叉研究的数据集。我有 10 个受试者,他们每个人都接受了 A 和 B 治疗,但顺序不同。(这是一项平衡的研究。)
我想看看 A 和 B 治疗如何改善血脂水平。我的想法是创建一个线性混合模型,将受试者作为随机效应;治疗、顺序和周期作为固定效应;最后,性别和年龄作为协变量。
我的数据:
#Reproducible data
id <- rep(1:10,3)
age <- rep(c("59","59","70","67","66","70","70","68","71","57"),3)
sex <- rep(c("F","M","F","M","F","F","F","M","F","M"),3)
sequence <- rep(c("1","2","1","2","1","2","1","2","2","2"),3)
period <- c(rep(0,10),rep(1,10),rep(2,10))
Treatment <- c(rep("C",10), rep(c("A","B"),4), "B","B",rep(c("B","A"),4), "A","A") #C is baseline
lipid <- c(18,6,30,12,14,19,10,22,22,27,13,28,14,23,12,27,9,10,13,22,13,22,29,12,16,24,15,13,17,11)
DF <- data.frame(id,age,sex,sequence,period,Treatment,lipid)
> head(DF)
id age sex sequence period Treatment lipid
1 1 59 F 1 0 C 18
2 2 59 M 2 0 C 6
3 3 70 F 1 0 C 30
4 4 67 M 2 0 C 12
5 5 66 F 1 0 C 14
6 6 70 F 2 0 C 19
我的线性混合模型:
library(lmerTest)
lm1 <- lmer(lipid~Treatment + sequence + period + sex + age + (1|id), data = DF, REML = F)
> summary(lm1)
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 1.344 1.159
Residual 34.986 5.915
Number of obs: 30, groups: id, 10
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 23.7890 18.2664 10.1410 1.302 0.222
TreatmentA -2.8500 2.8572 20.0000 -0.997 0.330
TreatmentB 2.2750 3.1018 20.0000 0.733 0.472
sequence2 4.7080 3.3324 10.0000 1.413 0.188
period1 -1.1250 2.6998 20.0000 -0.417 0.681
sexM -3.8351 3.7742 10.0000 -1.016 0.334
age -0.1078 0.2734 10.0000 -0.394 0.702
建立线性混合模型后,我想做事后测试来比较治疗 A 和 B。我尝试了 emmeans 和 multcomp,但它们给了我不同的结果。
艾默斯:
library(emmeans)
emm <- emmeans(lm1,"Treatment")
pairs(emm, adjust = "fdr")
> pairs(emm, adjust = "fdr")
contrast estimate SE df t.ratio p.value
C - A nonEst NA NA NA NA
C - B nonEst NA NA NA NA
A - B -5.12 2.93 23.5 -1.750 0.2794
多压缩:
library(multcomp)
summary(glht(lm1, linfct = mcp(Treatment = "Tukey")), test = adjusted("fdr"))
> summary(glht(lm1, linfct = mcp(Treatment = "Tukey")), test = adjusted("fdr"))
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
A - C == 0 -2.850 2.857 -0.997 0.463
B - C == 0 2.275 3.102 0.733 0.463
B - A == 0 5.125 2.700 1.898 0.173
(Adjusted p values reported -- fdr method)
我想问题是,
1) 基于研究设计,我的线性混合模型lm1 <- lmer(lipid~Treatment + sequence + period + sex + age + (1|id), data = DF, REML = F)
看起来还可以吗?还是我应该考虑其他交互项(例如治疗*序列)?
2) 为什么当multcomp给我值时, emmeans给我 CA 和 CB 的 NA ?由于结果不同,您会推荐哪一个对 lmer 模型进行事后测试?
任何想法表示赞赏,谢谢!