Reordering forested model variables in sjPlot with plot_models (r)

Question

I have a group of models that I am plotting together in R using the plot_models function from the sjPlot package. The models are plotting confidence intervals for similar felm regressions that use different datasets. All the models share the same variables of interest, and their confidence intervals occupy the same ranges. I have successly plotted all of the models, however, I now cannot order them how would be ideal.

My code is the following:

placebo_models <- paste0("outcome", 1:13, " ~ Variable1 + Variable2 | fe_variable |0| fe_variable") |> lapply(\(x) felm(as.formula(x), data = df))


models_list <- list(placebo_models[[13]], placebo_models[[12]], placebo_models[[11]], placebo_models[[10]], placebo_models[[9]], placebo_models[[8]], placebo_models[[7]], placebo_models[[6]], placebo_models[[5]], placebo_models[[4]], placebo_models[[3]], placebo_models[[2]], placebo_models[[1]])
                  
testplot <- plot_models(models_list, colors = c("RED"), show.values = TRUE, p.threshold = 0.05, spacing = 0.8, dot.size = 1, digits = 5, title = "Placebo Models", axis.labels = c("Variable1", "Variable2")) + scale_color_discrete(labels = c("m1", "m2", "m3", "m4", "m5", "m6", "m7", "m8", "m9", "m10", "m11", "m12", "m13"))

My current plot output is the following:

Ideally, I would like to change this model to report in order (descending from the top) m1's Variable 1 confidence interval, followed by m1's Variable 2 confidence interval, before going to m2 (Variable1, then Variable2), m3 (Variable1, then Variable2), etc.

Is there a way to reorder the plotted lines such that this is possible? The variables are identical across all the models, so it should fit, as long as it can be restructured. Can the variables of different models be reordered and inter-mixed within the same plot like so?

I am unfortunately not able to share the original data, but this visualization should explain the problem. Many thanks!

Answer 1

The order in which elements appear in your legend is based on the order you feed them to plot_models. The order in which the y-axis appears is the order in which they are placed in your first model.

Check it out. This is the example published by sjPlot::plot_models(). First, as it's presented in help, then how you can change the order of the legend items or the plotted items. As for the y-axis, the ability to rearrange it may be impacted by your model. You can rearrange the order in which you call the independent variables in many regression models without changing the results. I run a few checks here to show that this lm is not impacted by variable order.

library(sjPlot)

data(efc)

# fit three models
fit1 <- lm(barthtot ~ c160age + c12hour + c161sex + c172code, 
           data = efc)
fit2 <- lm(neg_c_7 ~ c160age + c12hour + c161sex + c172code, 
           data = efc)
fit3 <- lm(tot_sc_e ~ c160age + c12hour + c161sex + c172code, 
           data = efc)

# plot multiple models
plot_models(fit1, fit2, fit3, grid = TRUE)

# plot multiple models with legend labels and
# point shapes instead of value labels
plot_models(
  fit1, fit2, fit3,
  axis.labels = c(
    "Carer's Age", "Hours of Care", "Carer's Sex", 
    "Educational Status"
  ),
  m.labels = c("Barthel Index", "Negative Impact", 
               "Services used"),
  show.values = FALSE, show.p = FALSE, p.shape = TRUE
)

If I wanted to rearrange the legend Dependent Variables, I need to reorder the plots in the call. So I'll reference the model order. The labels need to be in the same order as the plots.

plot_models(
      fit3, fit2, fit1,
      axis.labels = c(
        "Carer's Age", "Hours of Care", "Carer's Sex", 
        "Educational Status"
      ),
      m.labels = rev("Barthel Index", "Negative Impact", 
                   "Services used"),
      show.values = FALSE, show.p = FALSE, p.shape = TRUE
    )

Next, I've rearranged the call for fit1 and returned to the original plotting order, you'll see that this plot has a different order on the y-axis, but the legend matches the legend in the first plot.

fit1a <- lm(barthtot ~ c172code + c161sex + c12hour + c160age, 
           data = efc)
plot_models(
  fit1a, fit2, fit3,
  axis.labels = rev(c(
    "Carer's Age", "Hours of Care", "Carer's Sex", 
    "Educational Status")),
  m.labels = c("Barthel Index", "Negative Impact", 
               "Services used"),
  show.values = FALSE, show.p = FALSE, p.shape = TRUE
)

# validate models are identical with R2 and Fstat
summary(fit1)[[8]] # [1] 0.2695598 
summary(fit1a)[[8]] # [1] 0.2695598 

summary(fit1)[[10]][1] # value # 75.28364  
summary(fit1a)[[10]][1] # value # 75.28364