Setup and Model Fitting

library(parameters)
library(effectsize)
library(insight)
library(see)
library(glmmTMB)
library(lme4)
library(lavaan)
library(metafor)
library(ggplot2)
library(splines)
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)
#> # Fixed Effects
#> 
#> Parameter   | Log-Mean |   SE |         95% CI |     z |      p
#> ---------------------------------------------------------------
#> (Intercept) |    -0.36 | 0.28 | [-0.90,  0.18] | -1.30 | 0.194 
#> spp [PR]    |    -1.27 | 0.24 | [-1.74, -0.80] | -5.27 | < .001
#> spp [DM]    |     0.27 | 0.14 | [ 0.00,  0.54] |  1.95 | 0.051 
#> spp [EC-A]  |    -0.57 | 0.21 | [-0.97, -0.16] | -2.75 | 0.006 
#> spp [EC-L]  |     0.67 | 0.13 | [ 0.41,  0.92] |  5.20 | < .001
#> spp [DES-L] |     0.63 | 0.13 | [ 0.38,  0.87] |  4.96 | < .001
#> spp [DF]    |     0.12 | 0.15 | [-0.17,  0.40] |  0.78 | 0.435 
#> mined [no]  |     1.27 | 0.27 | [ 0.74,  1.80] |  4.72 | < .001
#> 
#> # Zero-Inflated
#> 
#> Parameter   | Log-Odds |   SE |         95% CI |     z |      p
#> ---------------------------------------------------------------
#> (Intercept) |     0.79 | 0.27 | [ 0.26,  1.32] |  2.90 | 0.004 
#> mined [no]  |    -1.84 | 0.31 | [-2.46, -1.23] | -5.87 | < .001


model2 <- lm(Sepal.Length ~ Species * bs(Petal.Width, degree = 2), data = iris)
model_parameters(model2)
#> Parameter                                       | Coefficient |    SE |          95% CI | t(141) |      p
#> ---------------------------------------------------------------------------------------------------------
#> (Intercept)                                     |        4.79 |  0.17 | [  4.45,  5.13] |  27.66 | < .001
#> Species [versicolor]                            |       -3.73 |  2.14 | [ -7.96,  0.50] |  -1.74 | 0.083 
#> Species [virginica]                             |       -2.67 |  2.88 | [ -8.36,  3.03] |  -0.93 | 0.356 
#> Petal.Width [1st degree]                        |        2.53 |  2.36 | [ -2.13,  7.20] |   1.07 | 0.285 
#> Petal.Width [2nd degree]                        |      -11.18 | 21.14 | [-52.98, 30.62] |  -0.53 | 0.598 
#> Species [versicolor] * Petal.Width [1st degree] |        5.48 |  4.84 | [ -4.09, 15.05] |   1.13 | 0.260 
#> Species [virginica] * Petal.Width [1st degree]  |        2.37 |  4.35 | [ -6.22, 10.96] |   0.54 | 0.587 
#> Species [versicolor] * Petal.Width [2nd degree] |       14.84 | 21.16 | [-26.99, 56.68] |   0.70 | 0.484 
#> Species [virginica] * Petal.Width [2nd degree]  |       15.81 | 21.32 | [-26.35, 57.96] |   0.74 | 0.460


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)

plot(result)

When size_text is given, coefficients and confidence intervals are added to the plot.

plot(result, size_text = 4)

This also works for exponentiated coefficients.

result <- model_parameters(model1, exponentiate = TRUE)
plot(result, 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)

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 <- parameters::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")

plot(result)

Including group levels of random effects

result <- model_parameters(brms_model,
  effects = "all",
  component = "all", group_level = TRUE
)
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
#> Meta-analysis using 'metafor'
#> 
#> Parameter | Coefficient |   SE |         95% CI |     z |      p | Weight
#> -------------------------------------------------------------------------
#> Study 1   |       -0.39 | 0.32 | [-1.01,  0.23] | -1.24 | 0.215  |   9.95
#> Study 2   |        0.68 | 0.32 | [ 0.05,  1.30] |  2.13 | 0.033  |   9.95
#> Study 3   |        0.28 | 0.13 | [ 0.03,  0.54] |  2.17 | 0.030  |  59.17
#> Study 4   |       -1.40 | 0.36 | [-2.10, -0.69] | -3.88 | < .001 |   7.72
#> Overall   |       -0.18 | 0.44 | [-1.05,  0.68] | -0.42 | 0.676  |
plot(result)

If size_text is not NULL, estimates and confidence intervals are included in the plot.

plot(result, 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)
model <- rma(yi, vi, mods = ~disease, data = dat, digits = 3, slab = author)
result <- model_parameters(model)

result
#> # Depression
#> 
#> Parameter            | Coefficient |   SE |         95% CI |     z |      p | Weight
#> ------------------------------------------------------------------------------------
#> Aronson              |       -0.89 | 0.57 | [-2.01,  0.23] | -1.56 | 0.119  |   3.07
#> Ferguson & Simes     |       -1.59 | 0.44 | [-2.45, -0.72] | -3.59 | < .001 |   5.14
#> Rosenthal et al      |       -1.35 | 0.64 | [-2.61, -0.08] | -2.09 | 0.036  |   2.41
#> Frimodt-Moller et al |       -0.22 | 0.23 | [-0.66,  0.23] | -0.96 | 0.336  |  19.53
#> Coetzee & Berjak     |       -0.47 | 0.24 | [-0.94,  0.00] | -1.98 | 0.048  |  17.72
#> Overall              |       -0.12 | 0.59 | [-1.28,  1.04] | -0.20 | 0.841  |       
#> 
#> # CVD
#> 
#> Parameter          | Coefficient |   SE |         95% CI |      z |      p | Weight
#> -----------------------------------------------------------------------------------
#> Hart & Sutherland  |       -1.44 | 0.14 | [-1.72, -1.16] | -10.19 | < .001 |  49.97
#> Stein & Aronson    |       -0.79 | 0.08 | [-0.95, -0.62] |  -9.46 | < .001 | 144.81
#> Vandiviere et al   |       -1.62 | 0.47 | [-2.55, -0.70] |  -3.43 | < .001 |   4.48
#> TPT Madras         |        0.01 | 0.06 | [-0.11,  0.14] |   0.19 | 0.849  | 252.42
#> Comstock et al     |       -0.34 | 0.11 | [-0.56, -0.12] |  -3.05 | 0.002  |  80.57
#> Comstock & Webster |        0.45 | 0.73 | [-0.98,  1.88] |   0.61 | 0.541  |   1.88
#> Overall            |        0.03 | 0.56 | [-1.08,  1.13] |   0.05 | 0.963  |       
#> 
#> # Cancer
#> 
#> Parameter         | Coefficient |   SE |         95% CI |     z |      p | Weight
#> ---------------------------------------------------------------------------------
#> Rosenthal et al.1 |       -1.37 | 0.27 | [-1.90, -0.84] | -5.07 | < .001 |  13.69
#> Comstock et al.1  |       -0.02 | 0.27 | [-0.54,  0.51] | -0.06 | 0.948  |  14.00
#> Overall           |       -0.69 | 0.49 | [-1.65,  0.26] | -1.42 | 0.155  |
plot(result)

Bayesian Meta-Analysis using brms

# We download the model to save computation time. 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))
#
# 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
#> # Studies
#> 
#> Parameter              | Median |       89% CI | Weight |     pd | % in ROPE |  Rhat |      ESS
#> -----------------------------------------------------------------------------------------------
#> Call et al             |   0.63 | [0.33, 0.95] |   3.83 | 99.95% |     0.26% | 1.000 | 10283.00
#> Cavanagh et al         |   0.43 | [0.17, 0.68] |   5.09 | 99.45% |     2.35% | 1.000 |  7644.00
#> DanitzOrsillo          |   1.05 | [0.59, 1.51] |   2.89 |   100% |        0% | 1.000 |  4733.00
#> de Vibe et al          |   0.25 | [0.07, 0.44] |   8.49 | 98.39% |     9.78% | 1.000 |  4947.00
#> Frazier et al          |   0.45 | [0.25, 0.66] |   6.91 | 99.91% |     0.46% | 1.000 |  7127.00
#> Frogeli et al          |   0.60 | [0.35, 0.86] |   5.10 | 99.99% |     0.14% | 1.000 |  9028.00
#> Gallego et al          |   0.65 | [0.38, 0.96] |   4.45 |   100% |     0.04% | 1.000 |  9883.00
#> Hazlett Stevens   Oren |   0.54 | [0.27, 0.81] |   4.75 | 99.85% |     0.57% | 1.000 |  9039.00
#> Hintz et al            |   0.36 | [0.13, 0.59] |   5.95 | 99.36% |     3.51% | 1.000 |  7263.00
#> Kang et al             |   0.84 | [0.47, 1.27] |   2.97 |   100% |     0.01% | 1.000 |  7505.00
#> Kuhlmann et al         |   0.27 | [0.01, 0.55] |   5.14 | 93.21% |    14.96% | 1.000 |  6314.00
#> Lever Taylor et al     |   0.46 | [0.19, 0.76] |   4.33 | 99.34% |     2.15% | 1.000 |  9747.00
#> Phang et al            |   0.55 | [0.26, 0.85] |   4.09 | 99.83% |     0.61% | 1.000 | 10419.00
#> Rasanen et al          |   0.49 | [0.18, 0.80] |   3.88 | 99.26% |     2.20% | 1.000 | 10918.00
#> Ratanasiripong         |   0.54 | [0.19, 0.92] |   2.85 | 99.01% |     2.51% | 1.000 | 11984.00
#> Shapiro et al          |   0.96 | [0.57, 1.42] |   3.17 |   100% |        0% | 1.000 |  6546.00
#> SongLindquist          |   0.59 | [0.29, 0.86] |   4.41 | 99.95% |     0.39% | 1.000 |  9709.00
#> Warnecke et al         |   0.58 | [0.28, 0.89] |   4.02 | 99.81% |     0.78% | 1.000 |  9922.00
#> Overall                |   0.57 | [0.42, 0.71] |        |   100% |        0% | 1.000 |  3921.00
#> 
#> # Tau
#> 
#> Parameter | Median |       89% CI |   pd | % in ROPE
#> ----------------------------------------------------
#> tau       |   0.29 | [0.15, 0.46] | 100% |     1.36%
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
#> # TOST-test for Practical Equivalence
#> 
#>   ROPE: [-1.99 1.99]
#> 
#>         Parameter        H0 inside ROPE        90% CI
#>       (Intercept)  Rejected      0.00 % [59.33 68.41]
#>              time Undecided     83.52 % [-0.76  2.53]
#>               age  Accepted    100.00 % [-0.26  0.32]
#>   education [mid]  Rejected      0.00 % [ 5.13 12.39]
#>  education [high]  Rejected      0.00 % [10.14 18.57]

plot(result)

result <- equivalence_test(model4, rule = "cet")
result
#> # Conditional Equivalence Testing
#> 
#>   ROPE: [-1.99 1.99]
#> 
#>         Parameter        H0 inside ROPE        90% CI
#>       (Intercept)  Rejected      0.00 % [59.33 68.41]
#>              time Undecided     83.52 % [-0.76  2.53]
#>               age  Accepted    100.00 % [-0.26  0.32]
#>   education [mid]  Rejected      0.00 % [ 5.13 12.39]
#>  education [high]  Rejected      0.00 % [10.14 18.57]

plot(result)

For random effects

result <- equivalence_test(model3, effects = "random")
result
#> # TOST-test for Practical Equivalence
#> 
#>   ROPE: [-5.63 5.63]
#> 
#> Group: grp
#> 
#>  Parameter       H0 inside ROPE       90% CI
#>          1 Accepted    100.00 % [-1.91 4.28]
#>          2 Accepted    100.00 % [-3.27 2.84]
#>          3 Accepted    100.00 % [-3.33 2.83]
#>          4 Accepted    100.00 % [-3.56 2.60]
#>          5 Accepted    100.00 % [-3.33 2.85]
#> 
#> Group: Subject
#> 
#>  Parameter        H0 inside ROPE          90% CI
#>        308  Rejected      0.00 % [ 25.22  56.50]
#>        309  Rejected      0.00 % [-93.51 -62.24]
#>        310  Rejected      0.00 % [-78.91 -47.65]
#>        330 Undecided     36.01 % [-11.32  19.96]
#>        331 Undecided     35.46 % [ -5.45  25.80]
#>        332 Undecided     36.02 % [ -7.17  24.10]
#>        333  Rejected     15.44 % [  0.81  32.07]
#>        334 Undecided     36.00 % [-18.65  12.65]
#>        335  Rejected      0.00 % [-61.11 -29.85]
#>        337  Rejected      0.00 % [ 56.54  87.79]
#>        349  Rejected      1.22 % [-36.57  -5.25]
#>        350 Undecided     22.52 % [ -1.41  29.88]
#>        351 Undecided     35.97 % [-23.74   7.59]
#>        352  Rejected      0.00 % [ 20.76  52.12]
#>        369 Undecided     36.05 % [ -8.66  22.60]
#>        370 Undecided     35.99 % [-21.96   9.34]
#>        371 Undecided     36.05 % [-18.96  12.29]
#>        372  Rejected      9.80 % [  2.57  33.84]

plot(result)

From simulated model parameters

simulated_parameters <- simulate_model(model1, component = "conditional")
result <- equivalence_test(simulated_parameters)
result
#> # Test for Practical Equivalence
#> 
#>   ROPE: [-0.10 0.10]
#> 
#> Parameter   |        H0 | inside ROPE |       95% HDI
#> -----------------------------------------------------
#> (Intercept) | Undecided |     11.99 % | [-0.86  0.24]
#> sppPR       |  Rejected |      0.00 % | [-1.74 -0.84]
#> sppDM       | Undecided |     10.30 % | [-0.00  0.52]
#> sppEC-A     |  Rejected |      0.00 % | [-0.96 -0.19]
#> sppEC-L     |  Rejected |      0.00 % | [ 0.41  0.90]
#> sppDES-L    |  Rejected |      0.00 % | [ 0.40  0.89]
#> sppDF       | Undecided |     38.80 % | [-0.15  0.38]
#> minedno     |  Rejected |      0.00 % | [ 0.72  1.81]

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)

data(iris)
result <- cluster_analysis(iris[, 1:4], n_clusters = 3)
result
#> # Cluster Analysis (mean z-score by cluster)
#> 
#>          Term Group 1 Group 2 Group 3
#>  Sepal.Length   -1.01    0.09    1.24
#>   Sepal.Width    0.85   -0.70    0.07
#>  Petal.Length   -1.30    0.38    1.14
#>   Petal.Width   -1.25    0.31    1.19
#> 
#> # Accuracy of Cluster Group Classification
#> 
#>  Group Accuracy
#>      1  100.00%
#>      2   95.31%
#>      3   97.22%
#> 
#> Overall accuracy of classification: 97.33%

plot(result)


result <- cluster_analysis(iris[, 1:4], n_clusters = 4)
plot(result, n_columns = 2)

Number of Components/Factors to Retain

(related function documentation)

data(mtcars)
result <- n_factors(mtcars, type = "PCA")
result
#> # Method Agreement Procedure:
#> 
#> The choice of 3 dimensions is supported by 5 (29.41%) methods out of 17 (Bartlett, CNG, SE Scree, R2, Velicer's MAP).

plot(result)

plot(result, type = "line")

Number of Clusters to Retain

(related function documentation)

data(iris)
result <- n_clusters(standardize(iris[, 1:4]))
result
#> # Method Agreement Procedure:
#> 
#> The choice of 2 clusters is supported by 12 (42.86%) methods out of 28 (CH, Cindex, DB, Silhouette, Beale, Ratkowsky, PtBiserial, McClain, Dunn, SDindex, Mixture, Tibs2001SEmax).

plot(result)

plot(result, type = "line")

Description of Variable Distributions

(related function documentation)

Histogram for Numbers with Fractional Part

data(iris)
result <- describe_distribution(iris$Sepal.Length)
result
#> Mean |   SD |  IQR |        Range | Skewness | Kurtosis |   n | n_Missing
#> -------------------------------------------------------------------------
#> 5.84 | 0.83 | 1.30 | [4.30, 7.90] |     0.31 |    -0.55 | 150 |         0

plot(result)

Add Range of Dispersion (SD or MAD)

plot(result, dispersion = TRUE)

Thin Bars for Integer Values

set.seed(333)
x <- sample(1:100, 1000, replace = TRUE)
result <- describe_distribution(x)
result
#>  Mean |    SD |   IQR |          Range | Skewness | Kurtosis |    n | n_Missing
#> -------------------------------------------------------------------------------
#> 50.18 | 28.66 | 48.75 | [1.00, 100.00] |     0.02 |    -1.16 | 1000 |         0

plot(result)

Use a Normal Curve instead of Ribbon

plot(result, dispersion = TRUE, dispersion_style = "curve")

Highlighting Categories

set.seed(123)
result <- describe_distribution(sample(LETTERS[1:10], 1000, TRUE))

# highlight one category
plot(result, highlight = "D")


# highlight multiple categories
plot(result, highlight = c("D", "H"), size_bar = .4)


# own color scales - pass a named vector to 'scale_fill_manual()'
# the name of the non-highlighted color is "no_highlight".
plot(result, highlight = c("D", "H", "A"), size_bar = .4) +
  scale_fill_manual(values = c(D = "red", H = "green", A = "gold", no_highlight = "steelblue"))