library(tidyverse)
library(broom)
library(parameters)
library(tinytable)
model1 <- lm(body_mass ~ flipper_len * species, data = penguins){parameters} vs. {broom}
Here’s a basic model:
General printing
We can show it as a standard tibble with tidy()—this is all fine and normal
tidy(model1)
## # A tibble: 6 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -2536. 879. -2.88 4.19e- 3
## 2 flipper_len 32.8 4.63 7.10 7.69e-12
## 3 speciesChinstrap -501. 1523. -0.329 7.42e- 1
## 4 speciesGentoo -4251. 1427. -2.98 3.11e- 3
## 5 flipper_len:speciesChinstrap 1.74 7.86 0.222 8.25e- 1
## 6 flipper_len:speciesGentoo 21.8 6.94 3.14 1.84e- 3And we can see confidence intervals:
tidy(model1, conf.int = TRUE)
## # A tibble: 6 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -2536. 879. -2.88 4.19e- 3 -4266. -806.
## 2 flipper_len 32.8 4.63 7.10 7.69e-12 23.7 41.9
## 3 speciesChinstrap -501. 1523. -0.329 7.42e- 1 -3498. 2495.
## 4 speciesGentoo -4251. 1427. -2.98 3.11e- 3 -7059. -1444.
## 5 flipper_len:speciesChinstrap 1.74 7.86 0.222 8.25e- 1 -13.7 17.2
## 6 flipper_len:speciesGentoo 21.8 6.94 3.14 1.84e- 3 8.14 35.4We can use model_parameters() to get the same details, but with confidence intervals by default and nicer p-value formatting in a Markdown-esque table:
model_parameters(model1)
## Parameter | Coefficient | SE | 95% CI | t(336) | p
## --------------------------------------------------------------------------------------------------
## (Intercept) | -2535.84 | 879.47 | [-4265.79, -805.88] | -2.88 | 0.004
## flipper len | 32.83 | 4.63 | [ 23.73, 41.93] | 7.10 | < .001
## species [Chinstrap] | -501.36 | 1523.46 | [-3498.08, 2495.36] | -0.33 | 0.742
## species [Gentoo] | -4251.44 | 1427.33 | [-7059.08, -1443.81] | -2.98 | 0.003
## flipper len × species [Chinstrap] | 1.74 | 7.86 | [ -13.71, 17.19] | 0.22 | 0.825
## flipper len × species [Gentoo] | 21.79 | 6.94 | [ 8.14, 35.44] | 3.14 | 0.002
##
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.Manipulating the results
You can do stuff with the {broom} results, like this:
tidy(model1, conf.int = TRUE) |>
ggplot(aes(x = estimate, y = term)) +
geom_pointrange(aes(xmin = conf.low, xmax = conf.high))
Even though it’s printed by default in a non-tibble way, the {parameters} version is still a normal data frame behind the scenes, so you still plot it just fine:
model_parameters(model1) |>
ggplot(aes(x = Coefficient, y = Parameter)) +
geom_pointrange(aes(xmin = CI_low, xmax = CI_high))
Printing with Quarto
{broom} tables are tibbles and can be used with table-making functions like {tinytable}:
tidy(model1, conf.int = TRUE) |>
mutate(p.value = scales::label_pvalue()(p.value)) |>
tt(digits = 3) |>
style_tt(i = 2, background = "orange")| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | -2535.84 | 879.47 | -2.883 | 0.004 | -4265.79 | -805.9 |
| flipper_len | 32.83 | 4.63 | 7.095 | <0.001 | 23.73 | 41.9 |
| speciesChinstrap | -501.36 | 1523.46 | -0.329 | 0.742 | -3498.08 | 2495.4 |
| speciesGentoo | -4251.44 | 1427.33 | -2.979 | 0.003 | -7059.08 | -1443.8 |
| flipper_len:speciesChinstrap | 1.74 | 7.86 | 0.222 | 0.825 | -13.71 | 17.2 |
| flipper_len:speciesGentoo | 21.79 | 6.94 | 3.139 | 0.002 | 8.14 | 35.4 |
{parameters} tables can do the same thing:
model_parameters(model1) |>
select(-CI, -df_error) |>
mutate(p = scales::label_pvalue()(p)) |>
tt(digits = 3) |>
style_tt(i = 2, background = "orange")| Parameter | Coefficient | SE | CI_low | CI_high | t | p |
|---|---|---|---|---|---|---|
| (Intercept) | -2535.84 | 879.47 | -4265.79 | -805.9 | -2.883 | 0.004 |
| flipper_len | 32.83 | 4.63 | 23.73 | 41.9 | 7.095 | <0.001 |
| speciesChinstrap | -501.36 | 1523.46 | -3498.08 | 2495.4 | -0.329 | 0.742 |
| speciesGentoo | -4251.44 | 1427.33 | -7059.08 | -1443.8 | -2.979 | 0.003 |
| flipper_len:speciesChinstrap | 1.74 | 7.86 | -13.71 | 17.2 | 0.222 | 0.825 |
| flipper_len:speciesGentoo | 21.79 | 6.94 | 8.14 | 35.4 | 3.139 | 0.002 |
Or—even better—you can use display() to show it as a {tinytable} object:
model_parameters(model1) |>
display(format = "tt") |>
style_tt(i = 2, background = "orange")| Parameter | Coefficient | SE | 95% CI | t(336) | p |
|---|---|---|---|---|---|
| (Intercept) | -2535.84 | 879.47 | (-4265.79, -805.88) | -2.88 | 0.004 |
| flipper len | 32.83 | 4.63 | (23.73, 41.93) | 7.10 | < .001 |
| species (Chinstrap) | -501.36 | 1523.46 | (-3498.08, 2495.36) | -0.33 | 0.742 |
| species (Gentoo) | -4251.44 | 1427.33 | (-7059.08, -1443.81) | -2.98 | 0.003 |
| flipper len × species (Chinstrap) | 1.74 | 7.86 | (-13.71, 17.19) | 0.22 | 0.825 |
| flipper len × species (Gentoo) | 21.79 | 6.94 | (8.14, 35.44) | 3.14 | 0.002 |
Multilevel models
Multilevel models need {broom.mixed} to work; {parameters} supports them directly.
library(lme4)
model2 <- lmer(
body_mass ~ flipper_len * sex + (1 + flipper_len | species),
data = penguins
)library(broom.mixed)
tidy(model2, conf.int = TRUE)
## # A tibble: 8 × 8
## effect group term estimate std.error statistic conf.low conf.high
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) 43.6 694. 0.0628 -1316. 1404.
## 2 fixed <NA> flipper_len 19.0 4.25 4.47 10.7 27.3
## 3 fixed <NA> sexmale 83.3 485. 0.172 -867. 1033.
## 4 fixed <NA> flipper_len:sexmale 2.19 2.42 0.904 -2.56 6.93
## 5 ran_pars species sd__(Intercept) 524. NA NA NA NA
## 6 ran_pars species cor__(Intercept).flipper_len -1.000 NA NA NA NA
## 7 ran_pars species sd__flipper_len 4.90 NA NA NA NA
## 8 ran_pars Residual sd__Observation 294. NA NA NA NAmodel_parameters(model2)
## # Fixed Effects
##
## Parameter | Coefficient | SE | 95% CI | t(325) | p
## ---------------------------------------------------------------------------------------
## (Intercept) | 43.56 | 693.87 | [-1321.49, 1408.60] | 0.06 | 0.950
## flipper len | 19.00 | 4.25 | [ 10.64, 27.37] | 4.47 | < .001
## sex [male] | 83.26 | 484.74 | [ -870.36, 1036.89] | 0.17 | 0.864
## flipper len × sex [male] | 2.19 | 2.42 | [ -2.57, 6.94] | 0.90 | 0.367
##
## # Random Effects
##
## Parameter | Coefficient | SE | 95% CI
## ---------------------------------------------------------------------------------
## SD (Intercept: species) | 523.52 | 1082.80 | [ 9.09, 30164.15]
## SD (flipper_len: species) | 4.90 | 4.63 | [ 0.77, 31.29]
## Cor (Intercept~flipper_len: species) | -1.00 | 2.74 | [ -1.00, 1.00]
## SD (Residual) | 293.68 | 11.52 | [271.96, 317.14]
##
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation. Uncertainty intervals for random effect variances computed using a Wald z-distribution approximation.Standardized coefficients
{parameters} can automatically standardize model results (with lots of different methods—see here for more details)
model_parameters(model1, standardize = "refit")
## Parameter | Coefficient | SE | 95% CI | t(336) | p
## -----------------------------------------------------------------------------------------
## (Intercept) | -0.18 | 0.07 | [-0.32, -0.03] | -2.39 | 0.017
## flipper len | 0.58 | 0.08 | [ 0.42, 0.74] | 7.10 | < .001
## species [Chinstrap] | -0.19 | 0.10 | [-0.39, 0.01] | -1.87 | 0.062
## species [Gentoo] | 0.16 | 0.13 | [-0.11, 0.42] | 1.17 | 0.242
## flipper len × species [Chinstrap] | 0.03 | 0.14 | [-0.24, 0.30] | 0.22 | 0.825
## flipper len × species [Gentoo] | 0.38 | 0.12 | [ 0.14, 0.62] | 3.14 | 0.002
##
## Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.model_parameters(model1, standardize = "refit") |>
filter(Parameter != "(Intercept)") |>
ggplot(aes(x = Coefficient, y = Parameter)) +
geom_pointrange(aes(xmin = CI_low, xmax = CI_high))