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)
I have a dataframe as above and produce a glmer model (with the variable "animal" devation-coded):
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)
I wish to visualize this model using ggcoef_model. However, when I use
ggcoef_model(testM, include = !contains(c("owner","Residual.sd")))
It produces the following visualization
In the bottom of the graph, under the "animal * treat" header we see "cat * A2" and "cow * A2", but no "dog * A2" - why not? How can I display the "dog * A2" level as well?
With deviation coding (also known as sum contrasts), the sum of the contrasts for each factor sums to zero. Each level is compared to the overall mean rather than a specific reference level. This means that each level's coefficient represents the deviation from the mean effect of the factor levels.
For your model, the animal
factor is deviation-coded, and thus the baseline for the interaction terms will involve the reference level of animal
(in your case, likely dog
if dog
is the first level).
You can verify the reference level using levels(df$animal)
.
The estimate for dog:A2
is implicitly represented in the intercept and the overall mean, hence it is not shown explicitly. However we can make use of the emmeans
package which can be used to compute and visualise the estimated marginal means for each combination of factors, which includes the interaction effects. You could try something like this:
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)
These can then be visualised:
plot(emmeans_results, comparisons = TRUE)
which now includes the "missing" Dog:A2
interaction
Edit to address the 2nd question in the first comment, adding a combined plot of all main effects and interactions:
library(ggplot2)
emmeans_main_effects_animal <- emmeans(testM, ~ animal)
emmeans_main_effects_treat <- emmeans(testM, ~ treat)
emmeans_interactions <- emmeans(testM, ~ animal * treat)
emm_main_animal_df <- as.data.frame(emmeans_main_effects_animal)
emm_main_treat_df <- as.data.frame(emmeans_main_effects_treat)
emm_interactions_df <- as.data.frame(emmeans_interactions)
ggplot() +
# Main effects for animal
geom_point(data = emm_main_animal_df, aes(x = animal, y = emmean), color = "blue") +
geom_errorbar(data = emm_main_animal_df, aes(x = animal, ymin = asymp.LCL, ymax = asymp.UCL), width = 0.2, color = "blue") +
# Main effects for treat
geom_point(data = emm_main_treat_df, aes(x = treat, y = emmean), color = "green") +
geom_errorbar(data = emm_main_treat_df, aes(x = treat, ymin = asymp.LCL, ymax = asymp.UCL), width = 0.2, color = "green") +
# Interaction effects
geom_point(data = emm_interactions_df, aes(x = interaction(animal, treat), y = emmean), color = "red") +
geom_errorbar(data = emm_interactions_df, aes(x = interaction(animal, treat), ymin = asymp.LCL, ymax = asymp.UCL), width = 0.2, color = "red") +
labs(title = "Main Effects and Interaction Effects", x = "Levels", y = "Estimated Marginal Means") +
theme_minimal()