Search code examples
rggplot2statistics

ggcoef_model not displaying all interaction comparisons


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 Example

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?


Solution

  • 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)
    

    enter image description here

    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()
    

    enter image description here