lme4 deviant/tratment 对比编码与 R 中的交互 - 水平缺失
lme4 deviant/tratment contrast coding with interactions in R - levels are missing
我有一个带有 2 向交互项的混合效应模型(使用 lme4),每个项都有多个级别(每个 4 个),我想根据它们的总均值来研究它们的效果。我在这里展示汽车数据集中的这个例子,并省略了错误项,因为这个例子不需要它:
## shorten data frame for simplicity
df=Cars93[c(1:15),]
df=Cars93[is.element(Cars93$Make,c('Acura Integra', 'Audi 90','BMW 535i','Subaru Legacy')),]
df$Make=drop.levels(df$Make)
df$Model=drop.levels(df$Model)
## define contrasts (every factor has 4 levels)
contrasts(df$Make) = contr.treatment(4)
contrasts(df$Model) = contr.treatment(4)
## model
m1 <- lm(Price ~ Model*Make,data=df)
summary(m1)
如您所见,交互项中省略了第一级。我想在输出中包含所有 4 个级别,参考总均值(通常称为偏差编码)。这些是我查看的来源:https://marissabarlaz.github.io/portfolio/contrastcoding/#coding-schemes and 。不过,最后一个参考文献并未报告相互作用。
简单的回答是你想要的东西是不可能直接实现的。您必须使用稍微不同的方法。
在具有交互作用的模型中,您想使用平均值为零而不是特定水平的对比。否则,低阶效应(即主效应)不是主效应而是简单效应(当其他因子水平处于其参考水平时评估)。这在我关于混合模型的章节中有更详细的解释:
http://singmann.org/download/publications/singmann_kellen-introduction-mixed-models.pdf
为了得到你想要的结果,你必须以合理的方式拟合模型,然后将其传递给 emmeans
以与截距(即未加权的总均值)进行比较。这也适用于如下所示的交互(因为您的代码不起作用,我使用 warpbreaks
)。
afex::set_sum_contrasts() ## uses contr.sum globally
library("emmeans")
## model
m1 <- lm(breaks ~ wool * tension,data=warpbreaks)
car::Anova(m1, type = 3)
coef(m1)[1]
# (Intercept)
# 28.14815
## both CIs include grand mean:
emmeans(m1, "wool")
# wool emmean SE df lower.CL upper.CL
# A 31.0 2.11 48 26.8 35.3
# B 25.3 2.11 48 21.0 29.5
#
# Results are averaged over the levels of: tension
# Confidence level used: 0.95
## same using test
emmeans(m1, "wool", null = coef(m1)[1], infer = TRUE)
# wool emmean SE df lower.CL upper.CL null t.ratio p.value
# A 31.0 2.11 48 26.8 35.3 28.1 1.372 0.1764
# B 25.3 2.11 48 21.0 29.5 28.1 -1.372 0.1764
#
# Results are averaged over the levels of: tension
# Confidence level used: 0.95
emmeans(m1, "tension", null = coef(m1)[1], infer = TRUE)
# tension emmean SE df lower.CL upper.CL null t.ratio p.value
# L 36.4 2.58 48 31.2 41.6 28.1 3.196 0.0025
# M 26.4 2.58 48 21.2 31.6 28.1 -0.682 0.4984
# H 21.7 2.58 48 16.5 26.9 28.1 -2.514 0.0154
#
# Results are averaged over the levels of: wool
# Confidence level used: 0.95
emmeans(m1, c("tension", "wool"), null = coef(m1)[1], infer = TRUE)
# tension wool emmean SE df lower.CL upper.CL null t.ratio p.value
# L A 44.6 3.65 48 37.2 51.9 28.1 4.499 <.0001
# M A 24.0 3.65 48 16.7 31.3 28.1 -1.137 0.2610
# H A 24.6 3.65 48 17.2 31.9 28.1 -0.985 0.3295
# L B 28.2 3.65 48 20.9 35.6 28.1 0.020 0.9839
# M B 28.8 3.65 48 21.4 36.1 28.1 0.173 0.8636
# H B 18.8 3.65 48 11.4 26.1 28.1 -2.570 0.0133
#
# Confidence level used: 0.95
请注意,对于 coef()
,您可能希望将 fixef()
用于 lme4
模型。
我有一个带有 2 向交互项的混合效应模型(使用 lme4),每个项都有多个级别(每个 4 个),我想根据它们的总均值来研究它们的效果。我在这里展示汽车数据集中的这个例子,并省略了错误项,因为这个例子不需要它:
## shorten data frame for simplicity
df=Cars93[c(1:15),]
df=Cars93[is.element(Cars93$Make,c('Acura Integra', 'Audi 90','BMW 535i','Subaru Legacy')),]
df$Make=drop.levels(df$Make)
df$Model=drop.levels(df$Model)
## define contrasts (every factor has 4 levels)
contrasts(df$Make) = contr.treatment(4)
contrasts(df$Model) = contr.treatment(4)
## model
m1 <- lm(Price ~ Model*Make,data=df)
summary(m1)
如您所见,交互项中省略了第一级。我想在输出中包含所有 4 个级别,参考总均值(通常称为偏差编码)。这些是我查看的来源:https://marissabarlaz.github.io/portfolio/contrastcoding/#coding-schemes and
简单的回答是你想要的东西是不可能直接实现的。您必须使用稍微不同的方法。
在具有交互作用的模型中,您想使用平均值为零而不是特定水平的对比。否则,低阶效应(即主效应)不是主效应而是简单效应(当其他因子水平处于其参考水平时评估)。这在我关于混合模型的章节中有更详细的解释: http://singmann.org/download/publications/singmann_kellen-introduction-mixed-models.pdf
为了得到你想要的结果,你必须以合理的方式拟合模型,然后将其传递给 emmeans
以与截距(即未加权的总均值)进行比较。这也适用于如下所示的交互(因为您的代码不起作用,我使用 warpbreaks
)。
afex::set_sum_contrasts() ## uses contr.sum globally
library("emmeans")
## model
m1 <- lm(breaks ~ wool * tension,data=warpbreaks)
car::Anova(m1, type = 3)
coef(m1)[1]
# (Intercept)
# 28.14815
## both CIs include grand mean:
emmeans(m1, "wool")
# wool emmean SE df lower.CL upper.CL
# A 31.0 2.11 48 26.8 35.3
# B 25.3 2.11 48 21.0 29.5
#
# Results are averaged over the levels of: tension
# Confidence level used: 0.95
## same using test
emmeans(m1, "wool", null = coef(m1)[1], infer = TRUE)
# wool emmean SE df lower.CL upper.CL null t.ratio p.value
# A 31.0 2.11 48 26.8 35.3 28.1 1.372 0.1764
# B 25.3 2.11 48 21.0 29.5 28.1 -1.372 0.1764
#
# Results are averaged over the levels of: tension
# Confidence level used: 0.95
emmeans(m1, "tension", null = coef(m1)[1], infer = TRUE)
# tension emmean SE df lower.CL upper.CL null t.ratio p.value
# L 36.4 2.58 48 31.2 41.6 28.1 3.196 0.0025
# M 26.4 2.58 48 21.2 31.6 28.1 -0.682 0.4984
# H 21.7 2.58 48 16.5 26.9 28.1 -2.514 0.0154
#
# Results are averaged over the levels of: wool
# Confidence level used: 0.95
emmeans(m1, c("tension", "wool"), null = coef(m1)[1], infer = TRUE)
# tension wool emmean SE df lower.CL upper.CL null t.ratio p.value
# L A 44.6 3.65 48 37.2 51.9 28.1 4.499 <.0001
# M A 24.0 3.65 48 16.7 31.3 28.1 -1.137 0.2610
# H A 24.6 3.65 48 17.2 31.9 28.1 -0.985 0.3295
# L B 28.2 3.65 48 20.9 35.6 28.1 0.020 0.9839
# M B 28.8 3.65 48 21.4 36.1 28.1 0.173 0.8636
# H B 18.8 3.65 48 11.4 26.1 28.1 -2.570 0.0133
#
# Confidence level used: 0.95
请注意,对于 coef()
,您可能希望将 fixef()
用于 lme4
模型。