In easystats/see: Model Visualisation Toolbox for 'easystats' and 'ggplot2'

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  warning = FALSE,
  message = FALSE
)

pkgs <- c("ggplot2", "mgcv", "modelbased", "lme4", "rstanarm")
successfully_loaded <- vapply(pkgs, requireNamespace, FUN.VALUE = logical(1L), quietly = TRUE)
can_evaluate <- all(successfully_loaded)

if (can_evaluate) {
  knitr::opts_chunk$set(eval = TRUE)
  vapply(pkgs, require, FUN.VALUE = logical(1L), quietly = TRUE, character.only = TRUE)
} else {
  knitr::opts_chunk$set(eval = FALSE)
}

Introduction

modelbased is a package in easystats ecosystem to help with model-based estimations, to easily compute of marginal means, contrast analysis and model predictions.

For more, see: https://easystats.github.io/modelbased/

This vignette can be referred to by citing the package:

citation("see")

Let's first load all the needed libraries and set a common ggplot theme for all plots:

library(modelbased)
library(rstanarm)
library(ggplot2)
library(see)
library(lme4)
library(mgcv)

theme_set(theme_modern())

Pairwise Contrasts

model <- stan_glm(Sepal.Width ~ Species, data = iris, refresh = 0)

contrasts <- estimate_contrasts(model)
means <- estimate_means(model)

plot(contrasts, means)

Estimate model-based predictions for the response

Interactions, with continuous interaction terms

model <- lm(mpg ~ wt * gear, data = mtcars)

result <- estimate_expectation(model, data = "grid")
plot(result)

Interactions, with continuous interaction terms

mtcars$gear <- as.factor(mtcars$gear)
model <- lm(mpg ~ wt * gear, data = mtcars)

result <- estimate_expectation(model, data = "grid")
plot(result)

# full range
result <- estimate_relation(model, by = c("wt", "gear"), preserve_range = FALSE)
plot(result)

Interactions between two continuous variables

model <- lm(mpg ~ hp * wt, data = mtcars)

slopes <- estimate_slopes(model, trend = "hp", by = "wt")

plot(slopes)

Group-level scores of mixed models

model <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)

result <- estimate_grouplevel(model)
plot(result)

model <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)

result <- estimate_grouplevel(model)
plot(result)

Estimate slopes

model <- lm(Sepal.Width ~ Species * Petal.Length, data = iris)

result <- estimate_slopes(model, trend = "Petal.Length", by = "Species")
plot(result)

model <- lm(Petal.Length ~ Species * poly(Sepal.Width, 3), data = iris)

result <- estimate_slopes(model, by = c("Sepal.Width", "Species"))
plot(result)

model <- lm(Petal.Length ~ Species * poly(Sepal.Width, 3), data = iris)

result <- estimate_slopes(model, by = c("Sepal.Width", "Species"))
suppressWarnings(print(plot(result)))

Estimate derivatives

Linear-model

model_lm <- lm(mpg ~ wt, data = mtcars)

plot(estimate_relation(model_lm))

Non-linear model

# Fit a non-linear General Additive Model (GAM)
model <- mgcv::gam(Sepal.Width ~ s(Petal.Length), data = iris)

# 1. Compute derivatives
deriv <- estimate_slopes(model,
  trend = "Petal.Length",
  by = "Petal.Length",
  length = 100
)

# 2. Visualize predictions and derivative
plots(
  plot(estimate_relation(model)),
  plot(deriv),
  n_rows = 2
)