Random effects plots

DATE: 2020-10-25
AUTHOR: John L. Godlee


I was forwarded a question about how to customise sjPlot::plot_model()[1]. The colleague wanted to plot slopes of fixed effects at different levels of a random effect, in a random slope linear mixed effects model.

1: https://cran.r-project.org/web/packages/sjPlot/index.html (https://cran.r-project.org)

Load some packages:

library(ggplot2)
library(lmerTest)
library(sjPlot)
library(tibble)
library(purrr)

The example model was:

model <- lmerTest::lmer(mpg ~ cyl + disp + hp + drat + wt + qsec + carb +
  (1|am) + (cyl + disp + hp + drat + wt + qsec + carb|gear),
  data = mtcars, REML = FALSE)

This model doesn't even converge properly, but resembles the structure of their actual model.

The random effects are plotted like this:

sjPlot::plot_model(model, type="re", vline.color="#A9A9A9", dot.size=1.5,
  show.values=F, value.offset=.2)[[1]]

This produces a pretty ugly plot, which compresses the x axis and makes it difficult to compare slopes.

Random effects plotted using plot_model() (https://johngodlee.xyz)

Specifically, they wanted to plot the random effects for gear only, rearrange the figure to have 3 rows and 3 columns, and order the facets with intercept at the start.

The code below basically simplifies the source code of sjPlot:::plot_type_ranef() for our specific situation.

First, extract random effects for gear:

rand_ef <- ranef(model)[[1]]

Then, extract the standard errors:

vars.m <- attr(rand_ef, "postVar")
K <- dim(vars.m)[1]
J <- dim(vars.m)[3]
names.full <- dimnames(rand_ef)
rand_se <- array(NA, c(J, K))
for (j in 1:J) {
  rand_se[j, ] <- sqrt(diag(as.matrix(vars.m[, , j])))
  }
dimnames(rand_se) <- list(names.full[[1]], names.full[[2]])

Define some parameters for confidence intervals:

ci.lvl <- 0.95
ci <- 1 - ((1 - ci.lvl) / 2)

Make rownames a column for both effects and SEs:

rand_ef <- rownames_to_column(rand_ef)
rand_se <- rownames_to_column(as.data.frame(rand_se))

Define group names:

grp.names <- colnames(rand_ef)
alabels <- rand_ef[["rowname"]]

For each effect, calculate values per level of gear:

mydf <- map_df(2:ncol(rand_ef), function(i) {
  out <- data_frame(estimate = rand_ef[[i]])

  # Calculate confidence intervals
  out$conf.low = rand_ef[[i]] - (stats::qnorm(ci) * rand_se[[i]])
  out$conf.high = rand_ef[[i]] + (stats::qnorm(ci) * rand_se[[i]])

  # set column names (variable / coefficient name)
  # as group indicator, and save axis labels and title in variable
  out$facet <- grp.names[i]
  out$term <- factor(alabels)
  # create default grouping, depending on the effect:
  # split positive and negative associations with outcome
  # into different groups
  out$group <- dplyr::if_else(out$estimate > 0, "pos", "neg")

  return(out)
})

Make facet a factor to control plotting order:

mydf$facet <- factor(mydf$facet,
  levels = c("(Intercept)", "wt", "cyl", "carb", "disp", "drat",
    "hp", "qsec"))

Create plot:

ggplot() +
  geom_point(data = mydf,
    aes(x = estimate, y = term, colour = group)) +
  geom_errorbarh(data = mydf,
    aes(xmin = conf.low, xmax = conf.high, y = term, colour = group),
    height = 0) +
  scale_colour_manual(values = c("red", "blue")) +
  facet_wrap(~facet)

New plot (https://johngodlee.xyz)