R:对 lmer 的事后测试。emmeans 和 multcomp 软件包

机器算法验证 r 交叉研究
2022-03-21 00:34:15

我正在研究 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 模型进行事后测试?

任何想法表示赞赏,谢谢!

3个回答

在建模中,您必须小心不要以不同的方式包含完全相同的情况。例如,您已经发现所有period = 0案例的设计都Treatment C无法获得有用的结果。在输出中,当 3 个级别意味着应该有 2 个系数时summary(lm1),这导致仅报告 1个系数(这就是我发现问题的方式),并且正确拒绝提供涉及何时无法摆脱的对比periodemmeansTreatment Cperiod = 0

同样,您不能在模型中同时包含sequenceperiod,因为它们代表相同的事物。价值sequence = 1意味着得到Treatment Athen ; 意思是相反的。因此,对于那些拥有者的所有测量都将在之后进行,对于那些拥有者的测量将在之后进行,反之亦然Treatment Bsequence = 2period = 1sequence = 1Treatment Asequence = 2Treatment Bperiod = 2

如果您只有 10 个个体,每个个体进行 3 次测量(基线、第一次治疗后、第二次治疗后),如果您尝试评估超过 1 个或 2 个预测变量,您将面临过度拟合的风险。对于具有连续结果的生物医学情况,通常的经验法则是每个评估的预测变量大约有 15 个观察值。

但是,如果这是一项随机试验,那么您可以利用随机化本身,理想情况下,这应该消除对结果的任何贡献,除了设计中包含的那些:在您的情况下,两种治疗方法及其应用顺序。最好证明随机化在平衡您所知道的事物(基线血脂值、性别、年龄)方面做了合理的工作,但您可能不必在模型中将性别或年龄作为协变量包括在内。协变量的校正在大型研究中可能会有所帮助,但是您没有足够的案例来充分地做到这一点,而只有 10 名参与者。

我喜欢@Isabella Ghement 在评论中的建议,将基线血脂值合并为协变量,而不是将它们建模为结果。似乎没有实际Treatment C的,或者如果在该治疗之前似乎没有血脂值的记录,所以真的没有什么可以建模的。遵循该建议可能会消除包含随机效应截距项的需要,这将对可能不合适的截距施加高斯结构。

然后,您的模型归结为对Treatment Avs的评估Treatment B以及给予治疗的 2 个顺序。这似乎需要一个交互项 ( lipid ~ baselineLipid + Treatment + sequence + Treatment:sequence),其中Treatment0 表示TreatmentA, 1 表示TreatmentB你的研究可能仍然不够大,无法处理这么多的预测变量,但我认为这个模型有更好的机会公平地代表你的结果。

我将尝试回答(2)。

首先,给出的代码不会产生显示的结果。为了获得这些结果,我需要在拟合模型之前对数据进行以下操作:

DF$Treatment = relevel(factor(DF$Treatment), ref = "C")
DF$age = as.numeric(DF$age)
DF$period = factor(DF$period)

在以后的帖子中,请实际测试提供的代码,以确保您准确地表示发生了什么。

现在,关于这个问题。与multcomp等许多(大多数)包不同,emmeans 包测试估计可估计性与秩不足模型引起的歧义有关。模型排名不足的事实是问题中显示的重要遗漏。注意警告信息:

> lm1 <- lmer(lipid~Treatment + sequence + period + sex + age + (1|id), data = DF, REML = F)
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

当存在秩不足时,这意味着某些模型预测不是唯一的;它们根据模型的参数化方式而有所不同。在这种情况下,模型lm1使用默认的对比编码 ,contr.treatment其中每个因素使用 0--1 指标,省略了参考水平的指标。通过这种对比编码,我们得到了 OP 中显示的结果:(省略了一些输出)

> library(emmeans)
> pairs(emmeans(lm1, "Treatment"))

 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.2078 

> library(multcomp)
> summary(glht(lm1, linfct = mcp(Treatment = "Tukey")), test = adjusted("fdr"))

           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)

contr.sum现在,让我们使用编码拟合相同的模型:

> lm2 = update(lm1, contrasts = list(Treatment = "contr.sum", 
+   sequence = "contr.sum", period = "contr.sum", sex = "contr.sum"))
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

> fixef(lm2)
(Intercept)  Treatment1  Treatment2   sequence1     period1        sex1         age 
 23.6587644   1.6916667  -3.4083333  -2.3539871  -1.1250000   1.9175647  -0.1077586 

请注意lm1,由于对比度编码不同,大多数回归系数与 中的不同。现在用这个模型做同样的分析:

> pairs(emmeans(lm2, "Treatment"))

 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.2078 

> summary(glht(lm2, linfct = mcp(Treatment = "Tukey")), test = adjusted("fdr"))

            Estimate Std. Error z value Pr(>|z|)
A - C == 0   -5.100      5.065  -1.007    0.471
B - C == 0    0.025      4.613   0.005    0.996
B - A == 0    5.125      2.700   1.898    0.173
(Adjusted p values reported -- fdr method)

两个模型的emmeans结果相同。但是,multcomp结果是不同的,尽管B - A对比度相同。这种对比是独一无二我希望这能解释为什么emmeans没有显示两个比较,以及为什么multcomp真的应该测试可估计性。

我无法回答为什么multcomp会为一个可估计的比较产生不同的标准误差。

按照目前删除的建议sequence,我建议还包括period嵌套在 ID 中并将其从固定效果中删除,即lmer(lipid~Treatment + sex + age + (1|id/period), data = DF, REML = F)

我发现本指南很有帮助: https ://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#model-specification