LT17
LT17

Reputation: 353

Add legend to sjPlot model outputs modified with ggplot_build

I have two models, and I want to plot their estimates (using sjPlot), combining them in one plot. I create my models:

data(iris)

mod1 <- glmmTMB::glmmTMB(Sepal.Length ~ Petal.Length + Petal.Width + (1|Species),
                         data = iris,
                         family = "gaussian")

mod2 <- glmmTMB::glmmTMB(Sepal.Width ~ Petal.Length + Petal.Width + (1|Species),
                         data = iris,
                         family = "gaussian")

And I create my plot objects, with some small personalization:

plot_mod1 <- plot_model(model = mod1,
                        type = "est", # compute the estimates for linear models (standardized beta coefficients)
                        ci.lvl = 0.95, # adds 95% CI
                        ci.style = "whisker", # adds the horizontal lines
                        transform = NULL, # avoid automatic exponential tranformation of coefficients
                        axis.labels = "", # avoid automatic labelling
                        show.legend = T,
                        colors = "blue",
                        vline.color = "grey50") 

plot_mod2 <- plot_model(model = mod2,
                        type = "est", # compute the estimates for linear models (standardized beta coefficients)
                        ci.lvl = 0.95, # adds 95% CI
                        ci.style = "whisker", # adds the horizontal lines
                        transform = NULL, # avoid automatic exponential tranformation of coefficients
                        axis.labels = "", # avoid automatic labelling
                        show.legend = T,
                        colors = "red",
                        vline.color = "grey50") 

p1 <- plot_mod1 + 
  # add large x axis limits - no idea why it needs ylim and not xlim
  scale_y_continuous(breaks = seq(-0.5, 1.5, 0.5), limits = c(-0.5, 1.5)) + 
  # change y axis labels - no idea why it needs the scale_x and not scale_y
  scale_x_discrete(labels = list(
    Petal.Length = "Petal length",
    Petal.Width = "Petal width")) + 
  # add horizontal line to separate variables
  geom_vline(xintercept = 1.5, linetype = "dashed", color = "black") + 
  ylab(c("Estimate")) + 
  labs(title = "My plot")

p2 <- plot_mod2 + 
  # add large x axis limits - no idea why it needs ylim and not xlim
  scale_y_continuous(breaks = seq(-0.5, 1.5, 0.5), limits = c(-0.5, 1.5)) + 
  # change y axis labels - no idea why it needs the scale_x and not scale_y
  scale_x_discrete(labels = list(
    Petal.Length = "Petal length",
    Petal.Width = "Petal width")) + 
  ylab(c("Estimate")) + 
  labs(title = "My plot 2")

Even if I don't understand why the axes are reversed, i.e., why I need to specify discrete names on the x-axis if they are plotted on the y-axis, the plots look fine:

cowplot::plot_grid(p1,p2)

enter image description here

Now I need some deeper manipulation (change in transparency, plotting the estimates from the two models in the same plot...), so I extract the plot information with ggplot_build and play around with that:

qq1 <- ggplot_build(p1)
qq2 <- ggplot_build(p2)

# change plot 1 point size and transparency
qq1$data[[2]]$size <- 5  
qq1$data[[2]]$alpha <- 0.4
qq1$data[[3]]$alpha <- 0.4

# move mod1 estimates up 
qq1$data[[2]]$x <- (qq1$data[[2]]$x) + 0.25
qq1$data[[3]]$x <- (qq1$data[[3]]$x) + 0.25

# move mod2 estimates down
qq2$data[[1]]$x <- (qq2$data[[1]]$x) - 0.25
qq2$data[[2]]$x <- (qq2$data[[2]]$x) - 0.25

# merge everything
qq1$data[[2]] <- bind_rows(qq1$data[[2]], 
                            qq2$data[[1]])

qq1$data[[3]] <- bind_rows(qq1$data[[3]], 
                            qq2$data[[2]])

plot(ggplot_gtable(qq1))

And the output is what I want:

enter image description here

Finally, I would like to have a legend, with blue for mod1 and red for mod2. However, I have no clue where to map my legend. I tried sjPlot::plot_model(show.legend = T) but nothing appears (as they work for marginal effects plots and not for coefficients). I tried adding scale_color_manual(values = c("blue", "red") to p1 <- plot_mod1 + ..., but it changes the colours of the points without creating a legend. I am getting a bit lost as I don't know how to add a legend and where is it being mapped. Any clarification would be appreciated. I am using R version 4.2.2 and sjPlot_2.8.14.

Upvotes: 0

Views: 252

Answers (1)

stefan
stefan

Reputation: 125238

An easier approach would be to extract the model data using sjPlot::get_model_data. Then use ggplot2 to create your plot from scratch:

library(ggplot2)

mod_data <- lapply(list(mod1, mod2), \(mod) sjPlot::get_model_data(
  model = mod,
  type = "est",
  ci.lvl = 0.95,
  ci.style = "whisker",
  transform = NULL
))
names(mod_data) <- c("mod1", "mod2")
mod_data <- dplyr::bind_rows(mod_data, .id = "model")
mod_data$model <- forcats::fct_rev(mod_data$model)

ggplot(mod_data, aes(x = estimate, y = term, color = model)) +
  geom_vline(xintercept = 0, linewidth = .25) +
  geom_pointrange(
    aes(xmax = conf.high, xmin = conf.low, size = model, alpha = model),
    position = position_dodge(width = .75)
  ) +
  scale_y_discrete(labels = c(
    Petal.Length = "Petal length",
    Petal.Width = "Petal width"
  )) +
  scale_color_manual(values = c(mod1 = "blue", mod2 = "red"), breaks = rev) +
  scale_size_manual(values = c(mod1 = 1.25, mod2 = .25), breaks = rev) +
  scale_alpha_manual(values = c(mod1 = .4, mod2 = 1), breaks = rev) +
  geom_hline(yintercept = 1.5, linetype = "dashed", color = "black") +
  labs(title = "My plot", x = "Estimate", y = NULL) +
  scale_x_continuous(
    breaks = seq(-0.5, 1.5, 0.5), limits = c(-0.5, 1.5)
  )

enter image description here

Upvotes: 0

Related Questions