Summary of Model Parameters

Source: vignettes/model_parameters.Rmd

The model_parameters() function (also accessible via the shortcut parameters()) allows you to extract the parameters and their characteristics from various models in a consistent way. It can be considered as a lightweight alternative to broom::tidy(), with some notable differences:

  • The names of the returned data frame are specific to their content. For instance, the column containing the statistic is named following the statistic name, i.e., t, z, etc., instead of a generic name such as statistic (however, you can get standardized (generic) column names using standardize_names()).
  • It is able to compute or extract indices not available by default, such as p-values, CIs, etc.
  • It includes feature engineering capabilities, including parameters bootstrapping.

Correlations and t-tests

Frequentist

cor.test(iris$Sepal.Length, iris$Sepal.Width) %>%
  parameters()
#> Parameter1        |       Parameter2 |     r |     t |  df |     p |        95% CI |  Method
#> --------------------------------------------------------------------------------------------
#> iris$Sepal.Length | iris$Sepal.Width | -0.12 | -1.44 | 148 | 0.152 | [-0.27, 0.04] | Pearson
t.test(mpg ~ vs, data = mtcars) %>%
  parameters()
#> Parameter | Group | Mean_Group1 | Mean_Group2 | Difference |     t |    df |      p |          95% CI |                  Method
#> -------------------------------------------------------------------------------------------------------------------------------
#> mpg       |    vs |       16.62 |       24.56 |       7.94 | -4.67 | 22.72 | < .001 | [-11.46, -4.42] | Welch Two Sample t-test

Bayesian

library(BayesFactor)

BayesFactor::correlationBF(iris$Sepal.Length, iris$Sepal.Width) %>%
  parameters()
#> Parameter | Median |        89% CI |     pd | % in ROPE |              Prior | Effects |   Component |   BF
#> -----------------------------------------------------------------------------------------------------------
#> rho       |  -0.11 | [-0.23, 0.02] | 92.90% |    43.13% | Cauchy (0 +- 0.33) |   fixed | conditional | 0.51
BayesFactor::ttestBF(formula = mpg ~ vs, data = mtcars) %>%
  parameters()
#> Parameter  | Median |          89% CI |     pd | % in ROPE |              Prior | Effects |   Component |     BF
#> ----------------------------------------------------------------------------------------------------------------
#> Difference |  -7.30 | [-10.15, -4.58] | 99.98% |        0% | Cauchy (0 +- 0.71) |   fixed | conditional | 529.27

ANOVAs

Indices of effect size for ANOVAs, such as partial and non-partial versions of eta_squared(), epsilon_sqared() or omega_squared(), were moved to the effectsize-package. However, parameters uses these function to compute such indices for parameters summaries.

Simple

aov(Sepal.Length ~ Species, data = iris) %>%
  parameters(omega_squared = "partial", eta_squared = "partial", epsilon_squared = "partial")
#> Parameter | Sum_Squares |  df | Mean_Square |      F |      p | Omega_Sq (partial) | Eta_Sq (partial) | Epsilon_Sq (partial)
#> ----------------------------------------------------------------------------------------------------------------------------
#> Species   |       63.21 |   2 |       31.61 | 119.26 | < .001 |               0.61 |             0.62 |                 0.61
#> Residuals |       38.96 | 147 |        0.27 |        |        |                    |                  |

Repeated measures

parameters() (resp. its alias model_parameters()) also works on repeated measures ANOVAs, whether computed from aov() or from a mixed model.

aov(mpg ~ am + Error(gear), data = mtcars) %>%
  parameters()
#> Group  | Parameter | Sum_Squares | df | Mean_Square |    F |     p
#> ------------------------------------------------------------------
#> gear   |        am |      259.75 |  1 |      259.75 |      |      
#> Within |        am |      145.45 |  1 |      145.45 | 5.85 | 0.022
#> Within | Residuals |      720.85 | 29 |       24.86 |      |

Regressions (GLMs, Mixed Models, GAMs, …)

parameters() (resp. its alias model_parameters()) was mainly built with regression models in mind. It works for many types of models and packages, including mixed models and Bayesian models.

GLMs

glm(vs ~ poly(mpg, 2) + cyl, data = mtcars) %>%
  parameters()
#> Parameter        | Coefficient |   SE |         95% CI |     t | df |      p
#> ----------------------------------------------------------------------------
#> (Intercept)      |        2.04 | 0.39 | [ 1.27,  2.80] |  5.22 | 28 | < .001
#> mpg [1st degree] |       -0.33 | 0.61 | [-1.53,  0.87] | -0.53 | 28 | 0.599 
#> mpg [2nd degree] |        0.10 | 0.32 | [-0.54,  0.74] |  0.31 | 28 | 0.762 
#> cyl              |       -0.26 | 0.06 | [-0.38, -0.14] | -4.14 | 28 | < .001

Mixed Models

library(lme4)

lmer(Sepal.Width ~ Petal.Length + (1|Species), data = iris) %>%
  parameters()
#> Parameter    | Coefficient |   SE |       95% CI |    t |  df |      p
#> ----------------------------------------------------------------------
#> (Intercept)  |        2.00 | 0.56 | [0.90, 3.10] | 3.56 | 146 | < .001
#> Petal.Length |        0.28 | 0.06 | [0.17, 0.40] | 4.75 | 146 | < .001

Mixed Model with Zero-Inflation Model

library(GLMMadaptive)
library(glmmTMB)
data("Salamanders")
model <- mixed_model(
  count ~ spp + mined,
  random = ~1 | site,
  zi_fixed = ~spp + mined,
  family = zi.negative.binomial(),
  data = Salamanders
)
parameters(model)
#> # Fixed Effects component
#> 
#> Parameter   | Coefficient |   SE |        95% CI |     z |      p
#> -----------------------------------------------------------------
#> (Intercept) |       -0.63 | 0.40 | [-1.42, 0.16] | -1.56 | 0.118 
#> spp [PR]    |       -0.99 | 0.70 | [-2.35, 0.38] | -1.41 | 0.157 
#> spp [DM]    |        0.17 | 0.24 | [-0.29, 0.63] |  0.72 | 0.469 
#> spp [EC-A]  |       -0.39 | 0.35 | [-1.07, 0.29] | -1.13 | 0.258 
#> spp [EC-L]  |        0.49 | 0.24 | [ 0.02, 0.96] |  2.03 | 0.043 
#> spp [DES-L] |        0.59 | 0.23 | [ 0.14, 1.04] |  2.57 | 0.010 
#> spp [DF]    |       -0.11 | 0.24 | [-0.59, 0.37] | -0.46 | 0.642 
#> mined [no]  |        1.45 | 0.37 | [ 0.73, 2.17] |  3.95 | < .001
#> 
#> # Zero-Inflated component
#> 
#> Parameter   | Coefficient |   SE |         95% CI |     z |      p
#> ------------------------------------------------------------------
#> (Intercept) |        0.90 | 0.64 | [-0.35,  2.15] |  1.41 | 0.159 
#> spp [PR]    |        1.12 | 1.50 | [-1.82,  4.06] |  0.74 | 0.456 
#> spp [DM]    |       -0.95 | 0.82 | [-2.56,  0.65] | -1.17 | 0.244 
#> spp [EC-A]  |        1.04 | 0.72 | [-0.38,  2.46] |  1.44 | 0.150 
#> spp [EC-L]  |       -0.58 | 0.74 | [-2.03,  0.88] | -0.77 | 0.439 
#> spp [DES-L] |       -0.91 | 0.78 | [-2.43,  0.61] | -1.18 | 0.239 
#> spp [DF]    |       -2.63 | 2.37 | [-7.27,  2.02] | -1.11 | 0.268 
#> mined [no]  |       -2.56 | 0.63 | [-3.80, -1.32] | -4.06 | < .001

Mixed Models with Dispersion Model

library(glmmTMB)
sim1 <- function(nfac = 40, nt = 100, facsd = 0.1, tsd = 0.15, mu = 0, residsd = 1) {
  dat <- expand.grid(fac = factor(letters[1:nfac]), t = 1:nt)
  n <- nrow(dat)
  dat$REfac <- rnorm(nfac, sd = facsd)[dat$fac]
  dat$REt <- rnorm(nt, sd = tsd)[dat$t]
  dat$x <- rnorm(n, mean = mu, sd = residsd) + dat$REfac + dat$REt
  dat
}
set.seed(101)
d1 <- sim1(mu = 100, residsd = 10)
d2 <- sim1(mu = 200, residsd = 5)
d1$sd <- "ten"
d2$sd <- "five"
dat <- rbind(d1, d2)
model <- glmmTMB(x ~ sd + (1 | t), dispformula =  ~ sd, data = dat)

parameters(model)
#> # Fixed Effects component
#> 
#> Parameter   | Coefficient |   SE |            95% CI |       z |   df |      p
#> ------------------------------------------------------------------------------
#> (Intercept) |      200.03 | 0.10 | [ 199.84, 200.22] | 2056.35 | 5195 | < .001
#> sd [ten]    |      -99.71 | 0.22 | [-100.14, -99.29] | -458.39 | 5195 | < .001
#> 
#> # Dispersion component
#> 
#> Parameter   | Coefficient |   SE |       95% CI |      z |   df |      p
#> ------------------------------------------------------------------------
#> (Intercept) |        3.20 | 0.03 | [3.15, 3.26] | 115.48 | 5195 | < .001
#> sd [ten]    |        1.39 | 0.04 | [1.31, 1.46] |  35.35 | 5195 | < .001

Bayesian Models

model_parameters() works fine with Bayesian models from the rstanarm package…

library(rstanarm)

stan_glm(mpg ~ wt * cyl, data = mtcars) %>%
  parameters()
#> # Fixed effects 
#> 
#> Parameter   | Median |          89% CI |     pd | % in ROPE |  Rhat | ESS |               Prior
#> -----------------------------------------------------------------------------------------------
#> (Intercept) |  52.34 | [ 45.00, 62.69] |   100% |        0% | 1.026 | 184 | Normal (0 +- 60.27)
#> wt          |  -7.95 | [-11.64, -4.86] |   100% |     0.20% | 1.028 | 189 | Normal (0 +- 15.40)
#> cyl         |  -3.54 | [ -5.10, -2.24] | 99.80% |     0.80% | 1.027 | 185 |  Normal (0 +- 8.44)
#> wt:cyl      |   0.71 | [  0.27,  1.22] | 97.80% |    32.80% | 1.034 | 168 |  Normal (0 +- 1.36)

… as well as for (more complex) models from the brms package. For more complex models, other model components can be printed using the arguments effects and component arguments.

library(brms)
data(fish)
set.seed(123)
model <- brm(bf(
   count ~ persons + child + camper + (1 | persons),
   zi ~ child + camper + (1 | persons)
 ),
 data = fish,
 family = zero_inflated_poisson()
)
parameters(model, component = "conditional")
#> Parameter   | Median |         89% CI |     pd | % in ROPE | ESS |  Rhat
#> ------------------------------------------------------------------------
#> b_Intercept |  -0.87 | [-1.49, -0.08] | 96.80% |     4.80% |  78 | 1.000
#> b_persons   |   0.84 | [ 0.60,  1.06] |   100% |        0% |  75 | 0.997
#> b_child     |  -1.16 | [-1.32, -1.00] |   100% |        0% | 107 | 1.027
#> b_camper1   |   0.74 | [ 0.52,  0.91] |   100% |        0% | 224 | 0.993

parameters(model, effects = "all", component = "all")
#> # Fixed Effects (Count Model) 
#> 
#> Parameter   | Median |         89% CI |     pd | % in ROPE | ESS |  Rhat
#> ------------------------------------------------------------------------
#> (Intercept) |  -0.87 | [-1.49, -0.08] | 96.80% |     4.80% |  78 | 1.000
#> persons     |   0.84 | [ 0.60,  1.06] |   100% |        0% |  75 | 0.997
#> child       |  -1.16 | [-1.32, -1.00] |   100% |        0% | 107 | 1.027
#> camper1     |   0.74 | [ 0.52,  0.91] |   100% |        0% | 224 | 0.993
#> 
#> # Fixed Effects (Zero-Inflated Model) 
#> 
#> Parameter   | Median |         89% CI |     pd | % in ROPE | ESS |  Rhat
#> ------------------------------------------------------------------------
#> (Intercept) |  -0.76 | [-1.66,  0.51] | 87.20% |    10.40% |  98 | 0.992
#> child       |   1.87 | [ 1.37,  2.43] |   100% |        0% | 262 | 0.999
#> camper1     |  -0.83 | [-1.44, -0.22] | 99.20% |     0.80% | 168 | 0.997
#> 
#> # Random Effects (Count Model) 
#> 
#> Parameter | Median |        89% CI |     pd | % in ROPE | ESS |  Rhat
#> ---------------------------------------------------------------------
#> persons.1 |  -0.01 | [-0.40, 0.35] | 59.20% |    57.60% |  80 | 1.012
#> persons.2 |   0.03 | [-0.15, 0.33] | 61.60% |    60.80% |  88 | 0.994
#> persons.3 |  -0.02 | [-0.38, 0.11] | 63.20% |    64.80% |  66 | 1.008
#> persons.4 |   0.00 | [-0.51, 0.29] | 51.20% |    62.40% |  76 | 0.992
#> 
#> # Random Effects (Zero-Inflated Model) 
#> 
#> Parameter | Median |         89% CI |     pd | % in ROPE | ESS |  Rhat
#> ----------------------------------------------------------------------
#> persons.1 |   1.38 | [ 0.58,  2.66] | 97.60% |     1.60% | 108 | 0.992
#> persons.2 |   0.27 | [-0.62,  1.40] | 68.80% |    13.60% | 100 | 1.002
#> persons.3 |  -0.11 | [-1.36,  0.86] | 60.80% |    16.80% |  96 | 0.993
#> persons.4 |  -1.19 | [-2.62, -0.31] | 95.20% |     0.80% | 115 | 0.992

Structural Models (PCA, EFA, CFA, SEM…)

The parameters package extends the support to structural models.

Principal Component Analysis (PCA) and Exploratory Factor Analysis (EFA)

library(psych)

psych::pca(mtcars, nfactors = 3) %>%
  parameters()
#> # Rotated loadings from Principal Component Analysis (varimax-rotation)
#> 
#> Variable |   RC2 |   RC3 |   RC1 | Complexity | Uniqueness
#> ----------------------------------------------------------
#> mpg      |  0.66 | -0.41 | -0.54 |       2.63 |       0.10
#> cyl      | -0.62 |  0.67 |  0.34 |       2.49 |       0.05
#> disp     | -0.72 |  0.52 |  0.35 |       2.33 |       0.10
#> hp       | -0.30 |  0.64 |  0.63 |       2.40 |       0.10
#> drat     |  0.85 | -0.26 | -0.05 |       1.19 |       0.21
#> wt       | -0.78 |  0.21 |  0.51 |       1.90 |       0.08
#> qsec     | -0.18 | -0.91 | -0.28 |       1.28 |       0.06
#> vs       |  0.28 | -0.86 | -0.23 |       1.36 |       0.12
#> am       |  0.92 |  0.14 | -0.11 |       1.08 |       0.12
#> gear     |  0.91 | -0.02 |  0.26 |       1.16 |       0.10
#> carb     |  0.11 |  0.44 |  0.85 |       1.53 |       0.07
#> 
#> The 3 principal components (varimax rotation) accounted for 89.87% of the total variance of the original data (RC2 = 41.43%, RC3 = 29.06%, RC1 = 19.39%).
library(FactoMineR)

FactoMineR::FAMD(iris, ncp = 3) %>%
  parameters()
#> # Loadings from Factor Analysis (no rotation)
#> 
#> Variable     | Dim.1 | Dim.2 | Dim.3 | Complexity
#> -------------------------------------------------
#> Sepal.Length |  0.75 |  0.07 |  0.10 |       1.05
#> Sepal.Width  |  0.23 |  0.51 |  0.23 |       1.86
#> Petal.Length |  0.98 |  0.00 |  0.00 |       1.00
#> Petal.Width  |  0.94 |  0.01 |  0.00 |       1.00
#> Species      |  0.96 |  0.75 |  0.26 |       2.05
#> 
#> The 3 latent factors accounted for 96.73% of the total variance of the original data (Dim.1 = 64.50%, Dim.2 = 22.37%, Dim.3 = 9.86%).

Confirmatory Factor Analysis (CFA) and Structural Equation Models (SEM)

Frequentist

library(lavaan)

model <- lavaan::cfa(' visual  =~ x1 + x2 + x3
                       textual =~ x4 + x5 + x6
                       speed   =~ x7 + x8 + x9 ',
                       data = HolzingerSwineford1939)

model_parameters(model)
#> # Loading type
#> 
#> Link          | Coefficient |   SE |       95% CI |      p
#> ----------------------------------------------------------
#> visual =~ x1  |        1.00 | 0.00 | [1.00, 1.00] | < .001
#> visual =~ x2  |        0.55 | 0.10 | [0.36, 0.75] | < .001
#> visual =~ x3  |        0.73 | 0.11 | [0.52, 0.94] | < .001
#> textual =~ x4 |        1.00 | 0.00 | [1.00, 1.00] | < .001
#> textual =~ x5 |        1.11 | 0.07 | [0.98, 1.24] | < .001
#> textual =~ x6 |        0.93 | 0.06 | [0.82, 1.03] | < .001
#> speed =~ x7   |        1.00 | 0.00 | [1.00, 1.00] | < .001
#> speed =~ x8   |        1.18 | 0.16 | [0.86, 1.50] | < .001
#> speed =~ x9   |        1.08 | 0.15 | [0.79, 1.38] | < .001
#> 
#> # Correlation type
#> 
#> Link              | Coefficient |   SE |       95% CI |      p
#> --------------------------------------------------------------
#> visual ~~ textual |        0.41 | 0.07 | [0.26, 0.55] | < .001
#> visual ~~ speed   |        0.26 | 0.06 | [0.15, 0.37] | < .001
#> textual ~~ speed  |        0.17 | 0.05 | [0.08, 0.27] | < .001

Bayesian

blavaan to be done.

Meta-Analysis

parameters() also works for rma-objects from the metafor package.

library(metafor)

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

rma(yi = effectsize, sei = standarderror, method = "REML", data = mydat) %>%
  model_parameters()
#> 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  |

Plotting Model Parameters

There is a plot()-method implemented in the see-package. Several examples are shown in this vignette.