Summary of Model Parameters

• Correlations and t-tests
  • Frequentist
  • Bayesian
• ANOVAs
  • Simple
  • Repeated measures
• Regressions (GLMs, Mixed Models, GAMs, …)
  • GLMs
  • Mixed Models
  • Mixed Model with Zero-Inflation Model
  • Mixed Models with Dispersion Model
  • Bayesian Models
• Structural Models (PCA, EFA, CFA, SEM…)
  • Principal Component Analysis (PCA) and Exploratory Factor Analysis (EFA)
  • Confirmatory Factor Analysis (CFA) and Structural Equation Models (SEM)
• Meta-Analysis
• Plotting Model Parameters

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(148) |     p |        95% CI |  Method
#> ---------------------------------------------------------------------------------------
#> iris$Sepal.Length | iris$Sepal.Width | -0.12 |  -1.44 | 0.152 | [-0.27, 0.04] | Pearson

t.test(mpg ~ vs, data = mtcars) %>% 
  parameters()

#> Parameter | Group | Mean_Group1 | Mean_Group2 | Difference | t(22.72) |      p |          95% CI |                  Method
#> --------------------------------------------------------------------------------------------------------------------------
#> mpg       |    vs |       16.62 |       24.56 |       7.94 |    -4.67 | < .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 |    BF
#> ------------------------------------------------------------------------------------
#> rho       |  -0.11 | [-0.23, 0.02] | 92.25% |    44.23% | Cauchy (0 +- 0.33) | 0.509

BayesFactor::ttestBF(formula = mpg ~ vs, data = mtcars) %>% 
  parameters()

#> Parameter  | Median |        89% CI |     pd | % in ROPE |              Prior |     BF
#> --------------------------------------------------------------------------------------
#> Difference |   7.31 | [4.50, 10.12] | 99.98% |        0% | Cauchy (0 +- 0.71) | 529.27

ANOVAs

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

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 | Omega2 (partial) | Eta2 (partial) | Epsilon2 (partial)
#> ----------------------------------------------------------------------------------------------------------------------
#> Species   |       63.21 |   2 |       31.61 | 119.26 | < .001 |             0.61 |           0.62 |               0.61
#> Residuals |       38.96 | 147 |        0.27 |        |        |                  |                |

aov(Sepal.Length ~ Species * Sepal.Width, data = iris) %>%
  parameters(omega_squared = "partial", eta_squared = "partial", ci = .8)

#> Parameter           | Sum_Squares |  df | Mean_Square |      F |      p | Omega2 (partial) | Omega2 80% CI | Eta2 (partial) |  Eta2 80% CI
#> ------------------------------------------------------------------------------------------------------------------------------------------
#> Species             |       63.21 |   2 |       31.61 | 163.44 | < .001 |             0.68 |  [0.63, 0.72] |           0.69 | [0.64, 0.73]
#> Sepal.Width         |       10.95 |   1 |       10.95 |  56.64 | < .001 |             0.27 |  [0.19, 0.34] |           0.28 | [0.21, 0.36]
#> Species:Sepal.Width |        0.16 |   2 |        0.08 |   0.41 | 0.667  |        -7.98e-03 |  [0.00, 0.00] |       5.61e-03 | [0.00, 0.02]
#> Residuals           |       27.85 | 144 |        0.19 |        |        |                  |               |                |

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, family = binomial()) %>% 
  parameters()

#> Parameter        | Log-Odds |   SE |          95% CI |     z |     p
#> --------------------------------------------------------------------
#> (Intercept)      |    13.51 | 7.20 | [  2.56, 33.42] |  1.88 | 0.060
#> mpg [1st degree] |    -6.64 | 8.99 | [-27.81, 11.13] | -0.74 | 0.461
#> mpg [2nd degree] |     1.16 | 3.59 | [ -7.91,  8.56] |  0.32 | 0.746
#> cyl              |    -2.28 | 1.18 | [ -5.58, -0.51] | -1.92 | 0.055

glm(vs ~ poly(mpg, 2) + cyl, data = mtcars, family = binomial()) %>% 
  parameters(exponentiate = TRUE, df_method = "wald")

#> Parameter        | Odds Ratio |       SE |           95% CI |     z |     p
#> ---------------------------------------------------------------------------
#> (Intercept)      |   7.38e+05 | 5.31e+06 | [0.55, 9.87e+11] |  1.88 | 0.060
#> mpg [1st degree] |   1.31e-03 |     0.01 | [0.00, 59497.56] | -0.74 | 0.461
#> mpg [2nd degree] |       3.20 |    11.49 | [0.00,  3637.30] |  0.32 | 0.746
#> cyl              |       0.10 |     0.12 | [0.01,     1.05] | -1.92 | 0.055

Mixed Models

library(lme4)

lmer(Sepal.Width ~ Petal.Length + (1|Species), data = iris) %>% 
  parameters()

#> Parameter    | Coefficient |   SE |       95% CI | t(146) |      p
#> ------------------------------------------------------------------
#> (Intercept)  |        2.00 | 0.56 | [0.90, 3.10] |   3.56 | < .001
#> Petal.Length |        0.28 | 0.06 | [0.17, 0.40] |   4.75 | < .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
#> 
#> Parameter   | Log-Mean |   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
#> 
#> Parameter   | Log-Odds |   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
#> 
#> Parameter   | Coefficient |   SE |            95% CI |       z |      p
#> -----------------------------------------------------------------------
#> (Intercept) |      200.03 | 0.10 | [ 199.84, 200.22] | 2056.35 | < .001
#> sd [ten]    |      -99.71 | 0.22 | [-100.14, -99.29] | -458.39 | < .001
#> 
#> # Dispersion
#> 
#> Parameter   | Coefficient |   SE |       95% CI |      z |      p
#> -----------------------------------------------------------------
#> (Intercept) |        3.20 | 0.03 | [3.15, 3.26] | 115.48 | < .001
#> sd [ten]    |        1.39 | 0.04 | [1.31, 1.46] |  35.35 | < .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 |              CI |     pd | % in ROPE |  Rhat | ESS |                   Prior
#> ---------------------------------------------------------------------------------------------------
#> (Intercept) |  52.89 | [ 44.39, 62.73] |   100% |        0% | 1.043 |  36 | Normal (20.09 +- 15.07)
#> wt          |  -8.17 | [-11.75, -4.86] |   100% |        0% | 1.045 |  38 |  Normal (0.00 +- 15.40)
#> cyl         |  -3.61 | [ -5.10, -2.02] |   100% |     0.20% | 1.036 |  39 |   Normal (0.00 +- 8.44)
#> wt:cyl      |   0.73 | [  0.31,  1.25] | 99.20% |    31.80% | 1.047 |  34 |   Normal (0.00 +- 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")
#> # Fixed effects 
#> 
#> Parameter   | Median |         89% CI |     pd | % in ROPE |  Rhat |  ESS
#> -------------------------------------------------------------------------
#> (Intercept) |  -0.86 | [-1.40, -0.12] | 97.88% |     0.68% | 1.014 |  178
#> persons     |   0.84 | [ 0.64,  1.05] |   100% |        0% | 1.014 |  171
#> child       |  -1.16 | [-1.30, -1.00] |   100% |        0% | 1.001 | 1421
#> camper1     |   0.73 | [ 0.59,  0.88] |   100% |        0% | 1.000 | 3552
#> 
#> Using highest density intervals as credible intervals.

parameters(model, effects = "all", component = "all")
#> # Fixed effects (conditional) 
#> 
#> Parameter   | Median |         89% CI |     pd | % in ROPE |  Rhat |  ESS
#> -------------------------------------------------------------------------
#> (Intercept) |  -0.86 | [-1.40, -0.12] | 97.88% |     0.68% | 1.014 |  178
#> persons     |   0.84 | [ 0.64,  1.05] |   100% |        0% | 1.014 |  171
#> child       |  -1.16 | [-1.30, -1.00] |   100% |        0% | 1.001 | 1421
#> camper1     |   0.73 | [ 0.59,  0.88] |   100% |        0% | 1.000 | 3552
#> 
#> # Fixed effects (zero-inflated) 
#> 
#> Parameter   | Median |         89% CI |     pd | % in ROPE |  Rhat |  ESS
#> -------------------------------------------------------------------------
#> (Intercept) |  -0.65 | [-1.81,  0.37] | 84.72% |     7.25% | 1.007 |  988
#> child       |   1.86 | [ 1.33,  2.39] |   100% |        0% | 1.003 |  999
#> camper1     |  -0.81 | [-1.39, -0.23] | 98.75% |     1.65% | 1.001 | 2103
#> 
#> # Random effects (conditional) SD/Cor: persons
#> 
#> Parameter   | Median |       89% CI |   pd | % in ROPE |  Rhat | ESS
#> --------------------------------------------------------------------
#> (Intercept) |   0.12 | [0.00, 0.44] | 100% |    43.25% | 1.018 | 293
#> 
#> # Random effects (zero-inflated) SD/Cor: persons
#> 
#> Parameter    | Median |       89% CI |   pd | % in ROPE |  Rhat | ESS
#> ---------------------------------------------------------------------
#> zi_Intercept |   1.29 | [0.52, 2.56] | 100% |        0% | 1.009 | 484
#> 
#> Using highest density intervals as credible intervals.

To include information about the random effect parameters (group levels), set group_level = TRUE:

parameters(model, effects = "all", component = "conditional", group_level = TRUE)
#> # Fixed effects 
#> 
#> Parameter   | Median |         89% CI |     pd | % in ROPE |  Rhat |  ESS
#> -------------------------------------------------------------------------
#> (Intercept) |  -0.86 | [-1.40, -0.12] | 97.88% |     0.68% | 1.014 |  178
#> persons     |   0.84 | [ 0.64,  1.05] |   100% |        0% | 1.014 |  171
#> child       |  -1.16 | [-1.30, -1.00] |   100% |        0% | 1.001 | 1421
#> camper1     |   0.73 | [ 0.59,  0.88] |   100% |        0% | 1.000 | 3552
#> 
#> # Random effects Intercept: persons
#> 
#> Parameter |    Median |        89% CI |     pd | % in ROPE |  Rhat | ESS
#> ------------------------------------------------------------------------
#> persons.1 | -4.77e-03 | [-0.35, 0.36] | 54.00% |    61.05% | 1.019 | 182
#> persons.2 |      0.02 | [-0.15, 0.34] | 65.18% |    62.30% | 1.011 | 230
#> persons.3 |     -0.01 | [-0.23, 0.15] | 58.48% |    68.80% | 1.006 | 341
#> persons.4 |  2.72e-03 | [-0.22, 0.34] | 52.78% |    62.28% | 1.009 | 250
#> 
#> # Random effects SD/Cor: persons
#> 
#> Parameter   | Median |       89% CI |   pd | % in ROPE |  Rhat | ESS
#> --------------------------------------------------------------------
#> (Intercept) |   0.12 | [0.00, 0.44] | 100% |    43.25% | 1.018 | 293
#> 
#> Using highest density intervals as credible intervals.

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 | 1.32e-03 | 1.99e-03 |       1.00
#> Petal.Width  |  0.94 |     0.01 | 2.82e-05 |       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
#> 
#> 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
#> 
#> 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.