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
)

# The following are needed by `{parameters}` for cluster analysis:
# `{factoextra}`, `{mclust}`, and `{NbClust}`
#
# `{httr}` is needed for `insight::download_model()`
can_evaluate <- FALSE
pkgs <- c(
  "brms", "curl", "factoextra", "ggplot2", "glmmTMB", "httr",
  "lavaan", "NbClust", "mclust", "metafor", "lme4", "splines"
)
successfully_loaded <- vapply(pkgs, requireNamespace, FUN.VALUE = logical(1L), quietly = TRUE)
all_deps_available <- all(successfully_loaded)

# even if all dependencies are available, evaluate only if internet access is available
if (all_deps_available) {
  can_evaluate <- curl::has_internet()
}

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

This vignette can be referred to by citing the package:

citation("see")

Introduction

parameters' primary goal is to provide utilities for processing the parameters of various statistical models (see here for a list of supported models). Beyond computing p-values, CIs, Bayesian indices and other measures for a wide variety of models, this package implements features like bootstrapping of parameters and models, feature reduction (feature extraction and variable selection), or tools for data reduction like functions to perform cluster, factor or principal component analysis.

Another important goal of the parameters package is to facilitate and streamline the process of reporting results of statistical models, which includes the easy and intuitive calculation of standardized estimates or robust standard errors and p-values. parameters therefor offers a simple and unified syntax to process a large variety of (model) objects from many different packages.

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

Setup and Model Fitting

library(parameters)
library(effectsize)
library(insight)
library(see)
library(glmmTMB)
library(lme4)
library(lavaan)
library(metafor)
library(ggplot2)
data("Salamanders")
data("iris")
data("sleepstudy")
data("qol_cancer")

set.seed(12345)
sleepstudy$grp <- sample(1:5, size = 180, replace = TRUE)

theme_set(theme_modern())

# fit three example model
model1 <- glmmTMB(
  count ~ spp + mined + (1 | site),
  ziformula = ~mined,
  family = poisson(),
  data = Salamanders
)
model_parameters(model1, effects = "fixed")


model2 <- lm(Sepal.Length ~ Species * splines::bs(Petal.Width, degree = 2), data = iris)
model_parameters(model2, effects = "fixed")


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

model4 <- lm(QoL ~ time + age + education, data = qol_cancer)

Model Parameters

(related function documentation)

The plot()-method for model_parameters() creates a so called “forest plot”. In case of models with multiple components, parameters are separated into facets by model component.

result <- model_parameters(model1, effects = "fixed")

plot(result)

When show_labels is TRUE, coefficients and confidence intervals are added to the plot. Adjust the text size with size_text.

plot(result, show_labels = TRUE, size_text = 4)

This also works for exponentiated coefficients.

result <- model_parameters(model1, exponentiate = TRUE, effects = "fixed")
plot(result, show_labels = TRUE, size_text = 4)

It is also possible to plot only the count-model component. This is done in model_parameters() via the component argument. In easystats-functions, the count-component has the more generic name "conditional".

result <- model_parameters(model1, exponentiate = TRUE, component = "conditional")
plot(result)

As compared to the classical summary()-output, model_parameters(), and hence the plot()-method, tries to create human readable, prettier parameters names.

result <- model_parameters(model2)

plot(result)

Including group levels of random effects

result <- model_parameters(model3, group_level = TRUE)
plot(result)

Changing parameter names in the plot

By default, model_parameters() returns a data frame, where the parameter names are found in the column Parameter. These names are used by default in the generated plot:

plot(model_parameters(model4))

However, there are several ways to change the the names of the parameters. Since plot() returns ggplot objects, these can be easily modified, e.g. by adding a scale_y_discrete() layer. Note the correct order of the labels!

library(ggplot2)
plot(model_parameters(model4)) +
  scale_y_discrete(
    labels = c("Eucation (high-level)", "Education (mid-level)", "Age", "Time")
  )

Another way would be changing the values of the Parameter column before calling plot():

mp <- model_parameters(model4)
mp$Parameter <- c(
  "(Intercept)", "Time", "Age", "Education (mid-level)",
  "Eucation (high-level)"
)
plot(mp)

Simulated Model Parameters

simulate_parameters() computes simulated draws of parameters and their related indices such as Confidence Intervals (CI) and p-values. Simulating parameter draws can be seen as a (computationally faster) alternative to bootstrapping.

As simulate_parameters() is based on simulate_model() and thus simulates many draws for each parameter, plot() will produce similar plots as the density estimation plots from Bayesian models.

result <<- simulate_parameters(model1)

plot(result)

plot(result, stack = FALSE)

To avoid vertical overlapping, use normalize_height.

plot(result, stack = FALSE, normalize_height = TRUE)

plot(result, n_columns = 2)

plot(result, n_columns = 2, stack = FALSE)

Model Parameters of SEM models

structure <- " visual  =~ x1 + x2 + x3
               textual =~ x4 + x5 + x6
               speed   =~ x7 + x8 + x9 "

model <- lavaan::cfa(structure, data = HolzingerSwineford1939)
result <- model_parameters(model)
plot(result)

Model Parameters of Bayesian models

model_parameters() for Bayesian models will produce "forest plots" (instead of density estimations).

# We download the model to save computation time. Here is the code
# to refit the exact model used below...

# zinb <- read.csv("http://stats.idre.ucla.edu/stat/data/fish.csv")
# set.seed(123)
# model <- brm(bf(
#     count ~ persons + child + camper + (1 | persons),
#     zi ~ child + camper + (1 | persons)
#   ),
#   data = zinb,
#   family = zero_inflated_poisson()
# )
brms_model <- insight::download_model("brms_zi_2")
result <- model_parameters(
  brms_model,
  effects = "all",
  component = "all",
  verbose = FALSE
)

plot(result)

Including group levels of random effects

result <- model_parameters(brms_model,
  effects = "all",
  component = "all",
  group_level = TRUE,
  verbose = FALSE
)
plot(result)

One-column Layout

plot(result, n_column = 1)

Including Intercepts and Variance Estimates for Random Intercepts

plot(result, show_intercept = TRUE)

Model Parameters of Meta-Analysis models

mydat <<- data.frame(
  effectsize = c(-0.393, 0.675, 0.282, -1.398),
  standarderror = c(0.317, 0.317, 0.13, 0.36)
)

ma <- rma(yi = effectsize, sei = standarderror, method = "REML", data = mydat)
result <- model_parameters(ma)

result
plot(result)

If show_labels is TRUE, estimates and confidence intervals are included in the plot. Adjust the text size with size_text.

plot(result, show_labels = TRUE, size_text = 4)

Funnel plots

If type = "funnel", a funnel plot is created.

plot(result, type = "funnel")

Model Parameters of Meta-Analysis Models with Subgroups

set.seed(123)
data(dat.bcg)
dat <- escalc(
  measure = "RR",
  ai = tpos,
  bi = tneg,
  ci = cpos,
  di = cneg,
  data = dat.bcg
)
dat$author <- make.unique(dat$author)
dat$disease <- sample(c("Cancer", "CVD", "Depression"), size = nrow(dat), replace = TRUE)
mydat <<- dat
model <- rma(yi, vi, mods = ~disease, data = mydat, digits = 3, slab = author)
result <- model_parameters(model)

result
plot(result)

Bayesian Meta-Analysis using brms

We download the model to save computation time, but here is the code to refit the exact model used below:

# Data from
# https://github.com/MathiasHarrer/Doing-Meta-Analysis-in-R/blob/master/_data/Meta_Analysis_Data.RData
set.seed(123)
priors <- c(
  prior(normal(0, 1), class = Intercept),
  prior(cauchy(0, 0.5), class = sd)
)

model <- brm(
  TE | se(seTE) ~ 1 + (1 | Author),
  data = Meta_Analysis_Data,
  prior = priors,
  iter = 4000
)

library(brms)
model <- insight::download_model("brms_meta_1")
result <- model_parameters(model)

result
plot(result)

Comparison of Models

(related function documentation)

data(iris)
# shorter variable name
iris$Length <- iris$Petal.Length
lm1 <- lm(Sepal.Length ~ Species, data = iris)
lm2 <- lm(Sepal.Length ~ Species + Length, data = iris)
lm3 <- lm(Sepal.Length ~ Species * Length, data = iris)

result <- compare_parameters(lm1, lm2, lm3)
plot(result)

plot(result, size_text = 3.8) +
  labs(y = NULL) +
  theme(legend.position = "bottom")

Equivalence Testing

(related function documentation)

For fixed effects

# default rules, like in bayestestR::equivalence_test()
result <- equivalence_test(model4)
result

plot(result)

result <- equivalence_test(model4, rule = "cet")
result

plot(result)

For random effects

result <- equivalence_test(model3, effects = "random")
result

plot(result)

p-value function and consonance/compatibility plot

(related function documentation)

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

result <- p_function(model)
plot(result)

result <- p_function(model, ci_levels = c(0.5, 0.7, emph = 0.9))
plot(result, size_line = c(0.6, 1.2))

Probability of Direction

(related function documentation)

data(qol_cancer)
model <- lm(QoL ~ time + age + education, data = qol_cancer)
result <- p_direction(model)

plot(result)

Practical Significance

(related function documentation)

data(qol_cancer)
model <- lm(QoL ~ time + age + education, data = qol_cancer)
result <- p_significance(model)

plot(result)

Principal Component Analysis

(related function documentation)

data(mtcars)
result <- principal_components(mtcars[, 1:7], n = "all", threshold = 0.2)
result

plot(result)

result <- principal_components(
  mtcars[, 1:7],
  n = 3,
  rotation = "varimax",
  threshold = "max",
  sort = TRUE
)

result

plot(result, type = "line", text_color = "white") +
  theme_abyss()

Cluster Analysis

(related function documentation)

Clustering solution

data(iris)
result <- cluster_analysis(iris[, 1:4], n = 3)
result
plot(result)

result <- cluster_analysis(iris[, 1:4], n = 4)
plot(result)

Cluster centers

s <- summary(result)
s

plot(s)

Number of Components/Factors to Retain

(related function documentation)

data(mtcars)
result <- n_factors(mtcars, type = "PCA")
result

plot(result)
plot(result, type = "line")

Number of Clusters to Retain

(related function documentation)

data(iris)
result <- n_clusters(standardize(iris[, 1:4]))
result

plot(result)
plot(result, type = "line")

easystats/see documentation built on March 1, 2025, 3:54 p.m.