hunger <- sample(x = 0:1, size = 50, replace = TRUE)
treat <- sample(x = c("A1","A2"), prob = c(.75, .25), size = 50, replace = TRUE)
owner <- sample(x = c("Alice","Ben"), prob = c(.5, .5), size = 50, replace = TRUE)
animal <- sample(x = c("dog","cat","cow"), prob = c(.33, .33, .33), size = 50, replace = TRUE)
df <- as.data.frame(cbind(animal,treat,owner,hunger))
df$hunger <- as.numeric(df$hunger)
df$animal <- as.factor(df$animal)
我有一个如上所述的数据框并生成一个 glmer 模型(带有变量“动物”偏差编码):
contrasts(df$animal) = contr.sum(3) #setup contrast for deviation coding (1-2, 1-3. 2-3, comparisons)
testM <- glmer(hunger ~ animal + treat + animal*treat +
(1 + animal + treat + animal*treat||owner), data = df, family="binomial", glmerControl(optimizer="bobyqa", optCtrl = list(maxfun=1e5)))
summary(testM)
我希望使用 ggcoef_model 可视化该模型。然而,当我使用
ggcoef_model(testM, include = !contains(c("owner","Residual.sd")))
在图表底部,在“动物 * 治疗”标题下,我们看到“猫 * A2”和“牛 * A2”,但没有“狗 * A2” - 为什么不呢?如何才能同时显示“dog * A2”级别?
如果没有最小的工作示例,很难确定,但从描述来看,它听起来很像它与模型中分类变量及其交互的编码方式以及在此类模型中如何解释截距有关。这可能是可以利用估计边际均值的情况。
当您对因子
animal
使用偏差编码(总和对比度)时,模型会将其中一个级别设置为参考(基线),其他级别的系数将解释为与该参考级别的偏差。因此,总会有“缺失”的估计,但它们实际上并没有缺失,它们只是构成截距估计的一部分。
对于您的模型,动物因素是偏差编码的,因此交互项的基线将涉及
animal
的参考水平(在您的情况下,如果 dog
是第一个水平,则可能是 dog
)。
您可以使用
levels(df$animal)
验证参考水平。
emmeans
包可用于计算和可视化每个因素组合的估计边际均值,其中包括交互效应。你可以尝试这样的事情:
library(emmeans)
library(lme4)
contrasts(df$animal) = contr.sum(3)
testM <- glmer(hunger ~ animal + treat + animal*treat +
(1 + animal + treat + animal*treat || owner), data = df, family = "binomial", glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e5)))
emmeans_results <- emmeans(testM, ~ animal * treat)
summary(emmeans_results)
pairwise_comparisons <- contrast(emmeans_results, method = "pairwise")
summary(pairwise_comparisons)
然后可以将这些可视化:
plot(emmeans_results, comparisons = TRUE)