Package 'parameters'

Description

Bootstrap a statistical model n times to return a data frame of estimates.

Usage

bootstrap_model(model, iterations = 1000, ...)

## Default S3 method:
bootstrap_model(
  model,
  iterations = 1000,
  type = "ordinary",
  parallel = c("no", "multicore", "snow"),
  n_cpus = 1,
  verbose = FALSE,
  ...
)

## S3 method for class 'merMod'
bootstrap_model(
  model,
  iterations = 1000,
  type = "parametric",
  parallel = c("no", "multicore", "snow"),
  n_cpus = 1,
  cluster = NULL,
  verbose = FALSE,
  ...
)

Arguments

Details

By default, boot::boot() is used to generate bootstraps from the model data, which are then used to update() the model, i.e. refit the model with the bootstrapped samples. For merMod objects (lme4) or models from glmmTMB, the lme4::bootMer() function is used to obtain bootstrapped samples. bootstrap_parameters() summarizes the bootstrapped model estimates.

Value

A data frame of bootstrapped estimates.

Using with emmeans

The output can be passed directly to the various functions from the emmeans package, to obtain bootstrapped estimates, contrasts, simple slopes, etc. and their confidence intervals. These can then be passed to model_parameter() to obtain standard errors, p-values, etc. (see example).

Note that that p-values returned here are estimated under the assumption of translation equivariance: that shape of the sampling distribution is unaffected by the null being true or not. If this assumption does not hold, p-values can be biased, and it is suggested to use proper permutation tests to obtain non-parametric p-values.

See Also

bootstrap_parameters(), simulate_model(), simulate_parameters()

Examples

model <- lm(mpg ~ wt + factor(cyl), data = mtcars)
b <- bootstrap_model(model)
print(head(b))

est <- emmeans::emmeans(b, consec ~ cyl)
print(model_parameters(est))

Description

Compute bootstrapped parameters and their related indices such as Confidence Intervals (CI) and p-values.

Usage

bootstrap_parameters(model, ...)

## Default S3 method:
bootstrap_parameters(
  model,
  iterations = 1000,
  centrality = "median",
  ci = 0.95,
  ci_method = "quantile",
  test = "p-value",
  ...
)

Arguments

Details

This function first calls bootstrap_model() to generate bootstrapped coefficients. The resulting replicated for each coefficient are treated as "distribution", and is passed to bayestestR::describe_posterior() to calculate the related indices defined in the "test" argument.

Note that that p-values returned here are estimated under the assumption of translation equivariance: that shape of the sampling distribution is unaffected by the null being true or not. If this assumption does not hold, p-values can be biased, and it is suggested to use proper permutation tests to obtain non-parametric p-values.

Value

A data frame summarizing the bootstrapped parameters.

Using with emmeans

The output can be passed directly to the various functions from the emmeans package, to obtain bootstrapped estimates, contrasts, simple slopes, etc. and their confidence intervals. These can then be passed to model_parameter() to obtain standard errors, p-values, etc. (see example).

Note that that p-values returned here are estimated under the assumption of translation equivariance: that shape of the sampling distribution is unaffected by the null being true or not. If this assumption does not hold, p-values can be biased, and it is suggested to use proper permutation tests to obtain non-parametric p-values.

References

Davison, A. C., & Hinkley, D. V. (1997). Bootstrap methods and their application (Vol. 1). Cambridge university press.

See Also

bootstrap_model(), simulate_parameters(), simulate_model()

Examples

set.seed(2)
model <- lm(Sepal.Length ~ Species * Petal.Width, data = iris)
b <- bootstrap_parameters(model)
print(b)

est <- emmeans::emmeans(b, trt.vs.ctrl ~ Species)
print(model_parameters(est))

Description

Approximation of degrees of freedom based on a "between-within" heuristic.

Usage

ci_betwithin(model, ci = 0.95, ...)

dof_betwithin(model)

p_value_betwithin(model, dof = NULL, ...)

Arguments

Details

Small Sample Cluster corrected Degrees of Freedom

Inferential statistics (like p-values, confidence intervals and standard errors) may be biased in mixed models when the number of clusters is small (even if the sample size of level-1 units is high). In such cases it is recommended to approximate a more accurate number of degrees of freedom for such inferential statistics (see Li and Redden 2015). The Between-within denominator degrees of freedom approximation is recommended in particular for (generalized) linear mixed models with repeated measurements (longitudinal design). dof_betwithin() implements a heuristic based on the between-within approach. Note that this implementation does not return exactly the same results as shown in Li and Redden 2015, but similar.

Degrees of Freedom for Longitudinal Designs (Repeated Measures)

In particular for repeated measure designs (longitudinal data analysis), the between-within heuristic is likely to be more accurate than simply using the residual or infinite degrees of freedom, because dof_betwithin() returns different degrees of freedom for within-cluster and between-cluster effects.

Value

A data frame.

References

• Elff, M.; Heisig, J.P.; Schaeffer, M.; Shikano, S. (2019). Multilevel Analysis with Few Clusters: Improving Likelihood-based Methods to Provide Unbiased Estimates and Accurate Inference, British Journal of Political Science.

• Li, P., Redden, D. T. (2015). Comparing denominator degrees of freedom approximations for the generalized linear mixed model in analyzing binary outcome in small sample cluster-randomized trials. BMC Medical Research Methodology, 15(1), 38. doi:10.1186/s12874-015-0026-x

See Also

dof_betwithin() is a small helper-function to calculate approximated degrees of freedom of model parameters, based on the "between-within" heuristic.

Examples

if (require("lme4")) {
  data(sleepstudy)
  model <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
  dof_betwithin(model)
  p_value_betwithin(model)
}

Description

An approximate F-test based on the Kenward-Roger (1997) approach.

Usage

ci_kenward(model, ci = 0.95)

dof_kenward(model)

p_value_kenward(model, dof = NULL)

se_kenward(model)

Arguments

Details

Inferential statistics (like p-values, confidence intervals and standard errors) may be biased in mixed models when the number of clusters is small (even if the sample size of level-1 units is high). In such cases it is recommended to approximate a more accurate number of degrees of freedom for such inferential statistics. Unlike simpler approximation heuristics like the "m-l-1" rule (dof_ml1), the Kenward-Roger approximation is also applicable in more complex multilevel designs, e.g. with cross-classified clusters. However, the "m-l-1" heuristic also applies to generalized mixed models, while approaches like Kenward-Roger or Satterthwaite are limited to linear mixed models only.

Value

A data frame.

References

Kenward, M. G., & Roger, J. H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 983-997.

See Also

dof_kenward() and se_kenward() are small helper-functions to calculate approximated degrees of freedom and standard errors for model parameters, based on the Kenward-Roger (1997) approach.

dof_satterthwaite() and dof_ml1() approximate degrees of freedom based on Satterthwaite's method or the "m-l-1" rule.

Examples

if (require("lme4", quietly = TRUE)) {
  model <- lmer(Petal.Length ~ Sepal.Length + (1 | Species), data = iris)
  p_value_kenward(model)
}

Description

Approximation of degrees of freedom based on a "m-l-1" heuristic as suggested by Elff et al. (2019).

Usage

ci_ml1(model, ci = 0.95, ...)

dof_ml1(model)

p_value_ml1(model, dof = NULL, ...)

Arguments

Details

Small Sample Cluster corrected Degrees of Freedom

Inferential statistics (like p-values, confidence intervals and standard errors) may be biased in mixed models when the number of clusters is small (even if the sample size of level-1 units is high). In such cases it is recommended to approximate a more accurate number of degrees of freedom for such inferential statistics (see Li and Redden 2015). The m-l-1 heuristic is such an approach that uses a t-distribution with fewer degrees of freedom (dof_ml1()) to calculate p-values (p_value_ml1()) and confidence intervals (ci(method = "ml1")).

Degrees of Freedom for Longitudinal Designs (Repeated Measures)

In particular for repeated measure designs (longitudinal data analysis), the m-l-1 heuristic is likely to be more accurate than simply using the residual or infinite degrees of freedom, because dof_ml1() returns different degrees of freedom for within-cluster and between-cluster effects.

Limitations of the "m-l-1" Heuristic

Note that the "m-l-1" heuristic is not applicable (or at least less accurate) for complex multilevel designs, e.g. with cross-classified clusters. In such cases, more accurate approaches like the Kenward-Roger approximation (dof_kenward()) is recommended. However, the "m-l-1" heuristic also applies to generalized mixed models, while approaches like Kenward-Roger or Satterthwaite are limited to linear mixed models only.

Value

A data frame.

References

• Elff, M.; Heisig, J.P.; Schaeffer, M.; Shikano, S. (2019). Multilevel Analysis with Few Clusters: Improving Likelihood-based Methods to Provide Unbiased Estimates and Accurate Inference, British Journal of Political Science.

• Li, P., Redden, D. T. (2015). Comparing denominator degrees of freedom approximations for the generalized linear mixed model in analyzing binary outcome in small sample cluster-randomized trials. BMC Medical Research Methodology, 15(1), 38. doi:10.1186/s12874-015-0026-x

See Also

dof_ml1() is a small helper-function to calculate approximated degrees of freedom of model parameters, based on the "m-l-1" heuristic.

Examples

if (require("lme4")) {
  model <- lmer(Petal.Length ~ Sepal.Length + (1 | Species), data = iris)
  p_value_ml1(model)
}

Description

An approximate F-test based on the Satterthwaite (1946) approach.

Usage

ci_satterthwaite(model, ci = 0.95, ...)

dof_satterthwaite(model)

p_value_satterthwaite(model, dof = NULL, ...)

se_satterthwaite(model)

Arguments

Details

Inferential statistics (like p-values, confidence intervals and standard errors) may be biased in mixed models when the number of clusters is small (even if the sample size of level-1 units is high). In such cases it is recommended to approximate a more accurate number of degrees of freedom for such inferential statistics. Unlike simpler approximation heuristics like the "m-l-1" rule (dof_ml1), the Satterthwaite approximation is also applicable in more complex multilevel designs. However, the "m-l-1" heuristic also applies to generalized mixed models, while approaches like Kenward-Roger or Satterthwaite are limited to linear mixed models only.

Value

A data frame.

References

Satterthwaite FE (1946) An approximate distribution of estimates of variance components. Biometrics Bulletin 2 (6):110–4.

See Also

dof_satterthwaite() and se_satterthwaite() are small helper-functions to calculate approximated degrees of freedom and standard errors for model parameters, based on the Satterthwaite (1946) approach.

dof_kenward() and dof_ml1() approximate degrees of freedom based on Kenward-Roger's method or the "m-l-1" rule.

Examples

if (require("lme4", quietly = TRUE)) {
  model <- lmer(Petal.Length ~ Sepal.Length + (1 | Species), data = iris)
  p_value_satterthwaite(model)
}

Description

ci() attempts to return confidence intervals of model parameters.

Usage

## Default S3 method:
ci(x, ci = 0.95, dof = NULL, method = NULL, ...)

## S3 method for class 'glmmTMB'
ci(
  x,
  ci = 0.95,
  dof = NULL,
  method = "wald",
  component = "all",
  verbose = TRUE,
  ...
)

## S3 method for class 'merMod'
ci(x, ci = 0.95, dof = NULL, method = "wald", iterations = 500, ...)

Arguments

Value

A data frame containing the CI bounds.

Confidence intervals and approximation of degrees of freedom

There are different ways of approximating the degrees of freedom depending on different assumptions about the nature of the model and its sampling distribution. The ci_method argument modulates the method for computing degrees of freedom (df) that are used to calculate confidence intervals (CI) and the related p-values. Following options are allowed, depending on the model class:

Classical methods:

Classical inference is generally based on the Wald method. The Wald approach to inference computes a test statistic by dividing the parameter estimate by its standard error (Coefficient / SE), then comparing this statistic against a t- or normal distribution. This approach can be used to compute CIs and p-values.

"wald":

• Applies to non-Bayesian models. For linear models, CIs computed using the Wald method (SE and a t-distribution with residual df); p-values computed using the Wald method with a t-distribution with residual df. For other models, CIs computed using the Wald method (SE and a normal distribution); p-values computed using the Wald method with a normal distribution.

"normal"

• Applies to non-Bayesian models. Compute Wald CIs and p-values, but always use a normal distribution.

"residual"

• Applies to non-Bayesian models. Compute Wald CIs and p-values, but always use a t-distribution with residual df when possible. If the residual df for a model cannot be determined, a normal distribution is used instead.

Methods for mixed models:

Compared to fixed effects (or single-level) models, determining appropriate df for Wald-based inference in mixed models is more difficult. See the R GLMM FAQ for a discussion.

Several approximate methods for computing df are available, but you should also consider instead using profile likelihood ("profile") or bootstrap ("⁠boot"⁠) CIs and p-values instead.

"satterthwaite"

• Applies to linear mixed models. CIs computed using the Wald method (SE and a t-distribution with Satterthwaite df); p-values computed using the Wald method with a t-distribution with Satterthwaite df.

"kenward"

• Applies to linear mixed models. CIs computed using the Wald method (Kenward-Roger SE and a t-distribution with Kenward-Roger df); p-values computed using the Wald method with Kenward-Roger SE and t-distribution with Kenward-Roger df.

"ml1"

• Applies to linear mixed models. CIs computed using the Wald method (SE and a t-distribution with m-l-1 approximated df); p-values computed using the Wald method with a t-distribution with m-l-1 approximated df. See ci_ml1().

"betwithin"

• Applies to linear mixed models and generalized linear mixed models. CIs computed using the Wald method (SE and a t-distribution with between-within df); p-values computed using the Wald method with a t-distribution with between-within df. See ci_betwithin().

Likelihood-based methods:

Likelihood-based inference is based on comparing the likelihood for the maximum-likelihood estimate to the the likelihood for models with one or more parameter values changed (e.g., set to zero or a range of alternative values). Likelihood ratios for the maximum-likelihood and alternative models are compared to a $\chi$-squared distribution to compute CIs and p-values.

"profile"

• Applies to non-Bayesian models of class glm, polr, merMod or glmmTMB. CIs computed by profiling the likelihood curve for a parameter, using linear interpolation to find where likelihood ratio equals a critical value; p-values computed using the Wald method with a normal-distribution (note: this might change in a future update!)

"uniroot"

• Applies to non-Bayesian models of class glmmTMB. CIs computed by profiling the likelihood curve for a parameter, using root finding to find where likelihood ratio equals a critical value; p-values computed using the Wald method with a normal-distribution (note: this might change in a future update!)

Methods for bootstrapped or Bayesian models:

Bootstrap-based inference is based on resampling and refitting the model to the resampled datasets. The distribution of parameter estimates across resampled datasets is used to approximate the parameter's sampling distribution. Depending on the type of model, several different methods for bootstrapping and constructing CIs and p-values from the bootstrap distribution are available.

For Bayesian models, inference is based on drawing samples from the model posterior distribution.

"quantile" (or "eti")

• Applies to all models (including Bayesian models). For non-Bayesian models, only applies if bootstrap = TRUE. CIs computed as equal tailed intervals using the quantiles of the bootstrap or posterior samples; p-values are based on the probability of direction. See bayestestR::eti().

"hdi"

• Applies to all models (including Bayesian models). For non-Bayesian models, only applies if bootstrap = TRUE. CIs computed as highest density intervals for the bootstrap or posterior samples; p-values are based on the probability of direction. See bayestestR::hdi().

"bci" (or "bcai")

• Applies to all models (including Bayesian models). For non-Bayesian models, only applies if bootstrap = TRUE. CIs computed as bias corrected and accelerated intervals for the bootstrap or posterior samples; p-values are based on the probability of direction. See bayestestR::bci().

"si"

• Applies to Bayesian models with proper priors. CIs computed as support intervals comparing the posterior samples against the prior samples; p-values are based on the probability of direction. See bayestestR::si().

"boot"

• Applies to non-Bayesian models of class merMod. CIs computed using parametric bootstrapping (simulating data from the fitted model); p-values computed using the Wald method with a normal-distribution) (note: this might change in a future update!).

For all iteration-based methods other than "boot" ("hdi", "quantile", "ci", "eti", "si", "bci", "bcai"), p-values are based on the probability of direction (bayestestR::p_direction()), which is converted into a p-value using bayestestR::pd_to_p().

Examples

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

# regular confidence intervals
ci(model)

# using heteroscedasticity-robust standard errors
ci(model, vcov = "HC3")


library(parameters)
data(Salamanders, package = "glmmTMB")
model <- glmmTMB::glmmTMB(
  count ~ spp + mined + (1 | site),
  ziformula = ~mined,
  family = poisson(),
  data = Salamanders
)

ci(model)
ci(model, component = "zi")

Description

Compute hierarchical or kmeans cluster analysis and return the group assignment for each observation as vector.

Usage

cluster_analysis(
  x,
  n = NULL,
  method = "kmeans",
  include_factors = FALSE,
  standardize = TRUE,
  verbose = TRUE,
  distance_method = "euclidean",
  hclust_method = "complete",
  kmeans_method = "Hartigan-Wong",
  dbscan_eps = 15,
  iterations = 100,
  ...
)

Arguments

Details

The print() and plot() methods show the (standardized) mean value for each variable within each cluster. Thus, a higher absolute value indicates that a certain variable characteristic is more pronounced within that specific cluster (as compared to other cluster groups with lower absolute mean values).

Clusters classification can be obtained via print(x, newdata = NULL, ...).

Value

The group classification for each observation as vector. The returned vector includes missing values, so it has the same length as nrow(x).

Note

There is also a plot()-method implemented in the see-package.

References

• Maechler M, Rousseeuw P, Struyf A, Hubert M, Hornik K (2014) cluster: Cluster Analysis Basics and Extensions. R package.

See Also

• n_clusters() to determine the number of clusters to extract.

• cluster_discrimination() to determine the accuracy of cluster group classification via linear discriminant analysis (LDA).

• performance::check_clusterstructure() to check suitability of data for clustering.

• https://www.datanovia.com/en/lessons/

Examples

set.seed(33)
# K-Means ====================================================
rez <- cluster_analysis(iris[1:4], n = 3, method = "kmeans")
rez # Show results
predict(rez) # Get clusters
summary(rez) # Extract the centers values (can use 'plot()' on that)
if (requireNamespace("MASS", quietly = TRUE)) {
  cluster_discrimination(rez) # Perform LDA
}

# Hierarchical k-means (more robust k-means)
if (require("factoextra", quietly = TRUE)) {
  rez <- cluster_analysis(iris[1:4], n = 3, method = "hkmeans")
  rez # Show results
  predict(rez) # Get clusters
}

# Hierarchical Clustering (hclust) ===========================
rez <- cluster_analysis(iris[1:4], n = 3, method = "hclust")
rez # Show results
predict(rez) # Get clusters

# K-Medoids (pam) ============================================
if (require("cluster", quietly = TRUE)) {
  rez <- cluster_analysis(iris[1:4], n = 3, method = "pam")
  rez # Show results
  predict(rez) # Get clusters
}

# PAM with automated number of clusters
if (require("fpc", quietly = TRUE)) {
  rez <- cluster_analysis(iris[1:4], method = "pamk")
  rez # Show results
  predict(rez) # Get clusters
}

# DBSCAN ====================================================
if (require("dbscan", quietly = TRUE)) {
  # Note that you can assimilate more outliers (cluster 0) to neighbouring
  # clusters by setting borderPoints = TRUE.
  rez <- cluster_analysis(iris[1:4], method = "dbscan", dbscan_eps = 1.45)
  rez # Show results
  predict(rez) # Get clusters
}

# Mixture ====================================================
if (require("mclust", quietly = TRUE)) {
  library(mclust) # Needs the package to be loaded
  rez <- cluster_analysis(iris[1:4], method = "mixture")
  rez # Show results
  predict(rez) # Get clusters
}

Description

For each cluster, computes the mean (or other indices) of the variables. Can be used to retrieve the centers of clusters. Also returns the within Sum of Squares.

Usage

cluster_centers(data, clusters, fun = mean, ...)

Arguments

Value

A dataframe containing the cluster centers. Attributes include performance statistics and distance between each observation and its respective cluster centre.

Examples

k <- kmeans(iris[1:4], 3)
cluster_centers(iris[1:4], clusters = k$cluster)
cluster_centers(iris[1:4], clusters = k$cluster, fun = median)

Description

Computes linear discriminant analysis (LDA) on classified cluster groups, and determines the goodness of classification for each cluster group. See MASS::lda() for details.

Usage

cluster_discrimination(x, cluster_groups = NULL, ...)

Arguments

See Also

n_clusters() to determine the number of clusters to extract, cluster_analysis() to compute a cluster analysis and performance::check_clusterstructure() to check suitability of data for clustering.

Examples

# Retrieve group classification from hierarchical cluster analysis
clustering <- cluster_analysis(iris[, 1:4], n = 3)

# Goodness of group classification
cluster_discrimination(clustering)

Description

One of the core "issue" of statistical clustering is that, in many cases, different methods will give different results. The metaclustering approach proposed by easystats (that finds echoes in consensus clustering; see Monti et al., 2003) consists of treating the unique clustering solutions as a ensemble, from which we can derive a probability matrix. This matrix contains, for each pair of observations, the probability of being in the same cluster. For instance, if the 6th and the 9th row of a dataframe has been assigned to a similar cluster by 5 our of 10 clustering methods, then its probability of being grouped together is 0.5.

Usage

cluster_meta(list_of_clusters, rownames = NULL, ...)

Arguments

Details

Metaclustering is based on the hypothesis that, as each clustering algorithm embodies a different prism by which it sees the data, running an infinite amount of algorithms would result in the emergence of the "true" clusters. As the number of algorithms and parameters is finite, the probabilistic perspective is a useful proxy. This method is interesting where there is no obvious reasons to prefer one over another clustering method, as well as to investigate how robust some clusters are under different algorithms.

This metaclustering probability matrix can be transformed into a dissimilarity matrix (such as the one produced by dist()) and submitted for instance to hierarchical clustering (hclust()). See the example below.

Value

A matrix containing all the pairwise (between each observation) probabilities of being clustered together by the methods.

Examples

data <- iris[1:4]

rez1 <- cluster_analysis(data, n = 2, method = "kmeans")
rez2 <- cluster_analysis(data, n = 3, method = "kmeans")
rez3 <- cluster_analysis(data, n = 6, method = "kmeans")

list_of_clusters <- list(rez1, rez2, rez3)

m <- cluster_meta(list_of_clusters)

# Visualize matrix without reordering
heatmap(m, Rowv = NA, Colv = NA, scale = "none") # Without reordering
# Reordered heatmap
heatmap(m, scale = "none")

# Extract 3 clusters
predict(m, n = 3)

# Convert to dissimilarity
d <- as.dist(abs(m - 1))
model <- hclust(d)
plot(model, hang = -1)

Description

Compute performance indices for clustering solutions.

Usage

cluster_performance(model, ...)

## S3 method for class 'kmeans'
cluster_performance(model, ...)

## S3 method for class 'hclust'
cluster_performance(model, data, clusters, ...)

## S3 method for class 'dbscan'
cluster_performance(model, data, ...)

## S3 method for class 'parameters_clusters'
cluster_performance(model, ...)

Arguments

Examples

# kmeans
model <- kmeans(iris[1:4], 3)
cluster_performance(model)
# hclust
data <- iris[1:4]
model <- hclust(dist(data))
clusters <- cutree(model, 3)

rez <- cluster_performance(model, data, clusters)
rez

# DBSCAN
model <- dbscan::dbscan(iris[1:4], eps = 1.45, minPts = 10)

rez <- cluster_performance(model, iris[1:4])
rez

# Retrieve performance from parameters
params <- model_parameters(kmeans(iris[1:4], 3))
cluster_performance(params)

Description

Compute and extract model parameters of multiple regression models. See model_parameters() for further details.

Usage

compare_parameters(
  ...,
  ci = 0.95,
  effects = "fixed",
  component = "conditional",
  standardize = NULL,
  exponentiate = FALSE,
  ci_method = "wald",
  p_adjust = NULL,
  select = NULL,
  column_names = NULL,
  pretty_names = TRUE,
  coefficient_names = NULL,
  keep = NULL,
  drop = NULL,
  include_reference = FALSE,
  groups = NULL,
  verbose = TRUE
)

compare_models(
  ...,
  ci = 0.95,
  effects = "fixed",
  component = "conditional",
  standardize = NULL,
  exponentiate = FALSE,
  ci_method = "wald",
  p_adjust = NULL,
  select = NULL,
  column_names = NULL,
  pretty_names = TRUE,
  coefficient_names = NULL,
  keep = NULL,
  drop = NULL,
  include_reference = FALSE,
  groups = NULL,
  verbose = TRUE
)

Arguments

Details

This function is in an early stage and does not yet cope with more complex models, and probably does not yet properly render all model components. It should also be noted that when including models with interaction terms, not only do the values of the parameters change, but so does their meaning (from main effects, to simple slopes), thereby making such comparisons hard. Therefore, you should not use this function to compare models with interaction terms with models without interaction terms.

Value

A data frame of indices related to the model's parameters.

Examples

data(iris)
lm1 <- lm(Sepal.Length ~ Species, data = iris)
lm2 <- lm(Sepal.Length ~ Species + Petal.Length, data = iris)
compare_parameters(lm1, lm2)

# custom style
compare_parameters(lm1, lm2, select = "{estimate}{stars} ({se})")


# custom style, in HTML
result <- compare_parameters(lm1, lm2, select = "{estimate}<br>({se})|{p}")
print_html(result)


data(mtcars)
m1 <- lm(mpg ~ wt, data = mtcars)
m2 <- glm(vs ~ wt + cyl, data = mtcars, family = "binomial")
compare_parameters(m1, m2)

# exponentiate coefficients, but not for lm
compare_parameters(m1, m2, exponentiate = "nongaussian")

# change column names
compare_parameters("linear model" = m1, "logistic reg." = m2)
compare_parameters(m1, m2, column_names = c("linear model", "logistic reg."))

# or as list
compare_parameters(list(m1, m2))
compare_parameters(list("linear model" = m1, "logistic reg." = m2))

Description

Enables a conversion between Exploratory Factor Analysis (EFA) and Confirmatory Factor Analysis (CFA) lavaan-ready structure.

Usage

convert_efa_to_cfa(model, ...)

## S3 method for class 'fa'
convert_efa_to_cfa(
  model,
  threshold = "max",
  names = NULL,
  max_per_dimension = NULL,
  ...
)

efa_to_cfa(model, ...)

Arguments

Value

Converted index.

Examples

library(parameters)
data(attitude)
efa <- psych::fa(attitude, nfactors = 3)

model1 <- efa_to_cfa(efa)
model2 <- efa_to_cfa(efa, threshold = 0.3)
model3 <- efa_to_cfa(efa, max_per_dimension = 2)

suppressWarnings(anova(
  lavaan::cfa(model1, data = attitude),
  lavaan::cfa(model2, data = attitude),
  lavaan::cfa(model3, data = attitude)
))

Description

Estimate or extract degrees of freedom of models parameters.

Usage

degrees_of_freedom(model, method = "analytical", ...)

dof(model, method = "analytical", ...)

Arguments

Note

In many cases, degrees_of_freedom() returns the same as df.residuals(), or n-k (number of observations minus number of parameters). However, degrees_of_freedom() refers to the model's parameters degrees of freedom of the distribution for the related test statistic. Thus, for models with z-statistic, results from degrees_of_freedom() and df.residuals() differ. Furthermore, for other approximation methods like "kenward" or "satterthwaite", each model parameter can have a different degree of freedom.

Examples

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

model <- glm(vs ~ mpg * cyl, data = mtcars, family = "binomial")
dof(model)

model <- lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)
dof(model)

if (require("rstanarm", quietly = TRUE)) {
  model <- stan_glm(
    Sepal.Length ~ Petal.Length * Species,
    data = iris,
    chains = 2,
    refresh = 0
  )
  dof(model)
}

Description

Prints tables (i.e. data frame) in different output formats. print_md() is an alias for display(format = "markdown"), print_html() is an alias for display(format = "html"). print_table() is for specific use cases only, and currently only works for compare_parameters() objects.

Usage

## S3 method for class 'parameters_model'
display(
  object,
  format = "markdown",
  pretty_names = TRUE,
  split_components = TRUE,
  select = NULL,
  caption = NULL,
  subtitle = NULL,
  footer = NULL,
  align = NULL,
  digits = 2,
  ci_digits = digits,
  p_digits = 3,
  footer_digits = 3,
  ci_brackets = c("(", ")"),
  show_sigma = FALSE,
  show_formula = FALSE,
  zap_small = FALSE,
  font_size = "100%",
  line_padding = 4,
  column_labels = NULL,
  include_reference = FALSE,
  verbose = TRUE,
  ...
)

## S3 method for class 'parameters_sem'
display(
  object,
  format = "markdown",
  digits = 2,
  ci_digits = digits,
  p_digits = 3,
  ci_brackets = c("(", ")"),
  ...
)

## S3 method for class 'parameters_efa_summary'
display(object, format = "markdown", digits = 3, ...)

## S3 method for class 'parameters_efa'
display(
  object,
  format = "markdown",
  digits = 2,
  sort = FALSE,
  threshold = NULL,
  labels = NULL,
  ...
)

## S3 method for class 'equivalence_test_lm'
display(object, format = "markdown", digits = 2, ...)

print_table(x, digits = 2, p_digits = 3, theme = "default", ...)

Arguments

Details

display() is useful when the table-output from functions, which is usually printed as formatted text-table to console, should be formatted for pretty table-rendering in markdown documents, or if knitted from rmarkdown to PDF or Word files. See vignette for examples.

print_table() is a special function for compare_parameters() objects, which prints the output as a formatted HTML table. It is still somewhat experimental, thus, only a fixed layout-style is available at the moment (columns for estimates, confidence intervals and p-values). However, it is possible to include other model components, like zero-inflation, or random effects in the table. See 'Examples'. An alternative is to set engine = "tt" in print_html() to use the tinytable package for creating HTML tables.

Value

If format = "markdown", the return value will be a character vector in markdown-table format. If format = "html", an object of class gt_tbl. For print_table(), an object of class tinytable is returned.

See Also

print.parameters_model() and print.compare_parameters()

Examples

model <- lm(mpg ~ wt + cyl, data = mtcars)
mp <- model_parameters(model)
display(mp)


data(iris)
lm1 <- lm(Sepal.Length ~ Species, data = iris)
lm2 <- lm(Sepal.Length ~ Species + Petal.Length, data = iris)
lm3 <- lm(Sepal.Length ~ Species * Petal.Length, data = iris)
out <- compare_parameters(lm1, lm2, lm3)

print_html(
  out,
  select = "{coef}{stars}|({ci})",
  column_labels = c("Estimate", "95% CI")
)

# line break, unicode minus-sign
print_html(
  out,
  select = "{estimate}{stars}<br>({ci_low} \u2212 {ci_high})",
  column_labels = c("Est. (95% CI)")
)




data(iris)
data(Salamanders, package = "glmmTMB")
m1 <- lm(Sepal.Length ~ Species * Petal.Length, data = iris)
m2 <- lme4::lmer(
  Sepal.Length ~ Petal.Length + Petal.Width + (1 | Species),
  data = iris
)
m3 <- glmmTMB::glmmTMB(
  count ~ spp + mined + (1 | site),
  ziformula = ~mined,
  family = poisson(),
  data = Salamanders
)
out <- compare_parameters(m1, m2, m3, effects = "all", component = "all")
print_table(out)

Description

Computes Dominance Analysis Statistics and Designations

Usage

dominance_analysis(
  model,
  sets = NULL,
  all = NULL,
  conditional = TRUE,
  complete = TRUE,
  quote_args = NULL,
  contrasts = model$contrasts,
  ...
)

Arguments

Details

Computes two decompositions of the model's R2 and returns a matrix of designations from which predictor relative importance determinations can be obtained.

Note in the output that the "constant" subset is associated with a component of the model that does not directly contribute to the R2 such as an intercept. The "all" subset is apportioned a component of the fit statistic but is not considered a part of the dominance analysis and therefore does not receive a rank, conditional dominance statistics, or complete dominance designations.

The input model is parsed using insight::find_predictors(), does not yet support interactions, transformations, or offsets applied in the R formula, and will fail with an error if any such terms are detected.

The model submitted must accept an formula object as a formula argument. In addition, the model object must accept the data on which the model is estimated as a data argument. Formulas submitted using object references (i.e., lm(mtcars$mpg ~ mtcars$vs)) and functions that accept data as a non-data argument (e.g., survey::svyglm() uses design) will fail with an error.

Models that return TRUE for the insight::model_info() function's values "is_bayesian", "is_mixed", "is_gam", is_multivariate", "is_zero_inflated", or "is_hurdle" are not supported at current.

When performance::r2() returns multiple values, only the first is used by default.

Value

Object of class "parameters_da".

An object of class "parameters_da" is a list of data.frames composed of the following elements:

General

A data.frame which associates dominance statistics with model parameters. The variables in this data.frame include:

Parameter

Parameter names.

General_Dominance

Vector of general dominance statistics. The R2 ascribed to variables in the all argument are also reported here though they are not general dominance statistics.

Percent

Vector of general dominance statistics normalized to sum to 1.

Ranks

Vector of ranks applied to the general dominance statistics.

Subset

Names of the subset to which the parameter belongs in the dominance analysis. Each other data.frame returned will refer to these subset names.

Conditional

A data.frame of conditional dominance statistics. Each observation represents a subset and each variable represents an the average increment to R2 with a specific number of subsets in the model. NULL if conditional argument is FALSE.

Complete

A data.frame of complete dominance designations. The subsets in the observations are compared to the subsets referenced in each variable. Whether the subset in each variable dominates the subset in each observation is represented in the logical value. NULL if complete argument is FALSE.

Author(s)

Joseph Luchman

References

• Azen, R., & Budescu, D. V. (2003). The dominance analysis approach for comparing predictors in multiple regression. Psychological Methods, 8(2), 129-148. doi:10.1037/1082-989X.8.2.129

• Budescu, D. V. (1993). Dominance analysis: A new approach to the problem of relative importance of predictors in multiple regression. Psychological Bulletin, 114(3), 542-551. doi:10.1037/0033-2909.114.3.542

• Groemping, U. (2007). Estimators of relative importance in linear regression based on variance decomposition. The American Statistician, 61(2), 139-147. doi:10.1198/000313007X188252

See Also

domir::domin()

Examples

data(mtcars)

# Dominance Analysis with Logit Regression
model <- glm(vs ~ cyl + carb + mpg, data = mtcars, family = binomial())

performance::r2(model)
dominance_analysis(model)

# Dominance Analysis with Weighted Logit Regression
model_wt <- glm(vs ~ cyl + carb + mpg,
  data = mtcars,
  weights = wt, family = quasibinomial()
)

dominance_analysis(model_wt, quote_args = "weights")

Description

Compute the (conditional) equivalence test for frequentist models.

Usage

## S3 method for class 'lm'
equivalence_test(
  x,
  range = "default",
  ci = 0.95,
  rule = "classic",
  vcov = NULL,
  vcov_args = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'merMod'
equivalence_test(
  x,
  range = "default",
  ci = 0.95,
  rule = "classic",
  effects = c("fixed", "random"),
  vcov = NULL,
  vcov_args = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'ggeffects'
equivalence_test(
  x,
  range = "default",
  rule = "classic",
  test = "pairwise",
  verbose = TRUE,
  ...
)

Arguments

Details

In classical null hypothesis significance testing (NHST) within a frequentist framework, it is not possible to accept the null hypothesis, H0 - unlike in Bayesian statistics, where such probability statements are possible. "... one can only reject the null hypothesis if the test statistics falls into the critical region(s), or fail to reject this hypothesis. In the latter case, all we can say is that no significant effect was observed, but one cannot conclude that the null hypothesis is true." (Pernet 2017). One way to address this issues without Bayesian methods is Equivalence Testing, as implemented in equivalence_test(). While you either can reject the null hypothesis or claim an inconclusive result in NHST, the equivalence test - according to Pernet - adds a third category, "accept". Roughly speaking, the idea behind equivalence testing in a frequentist framework is to check whether an estimate and its uncertainty (i.e. confidence interval) falls within a region of "practical equivalence". Depending on the rule for this test (see below), statistical significance does not necessarily indicate whether the null hypothesis can be rejected or not, i.e. the classical interpretation of the p-value may differ from the results returned from the equivalence test.

Calculation of equivalence testing

• "bayes" - Bayesian rule (Kruschke 2018)

  This rule follows the "HDI+ROPE decision rule" (Kruschke, 2014, 2018) used for the Bayesian counterpart(). This means, if the confidence intervals are completely outside the ROPE, the "null hypothesis" for this parameter is "rejected". If the ROPE completely covers the CI, the null hypothesis is accepted. Else, it's undecided whether to accept or reject the null hypothesis. Desirable results are low proportions inside the ROPE (the closer to zero the better).

• "classic" - The TOST rule (Lakens 2017)

  This rule follows the "TOST rule", i.e. a two one-sided test procedure (Lakens 2017). Following this rule...

  • practical equivalence is assumed (i.e. H0 "accepted") when the narrow confidence intervals are completely inside the ROPE, no matter if the effect is statistically significant or not;

  • practical equivalence (i.e. H0) is rejected, when the coefficient is statistically significant, both when the narrow confidence intervals (i.e. 1-2*alpha) include or exclude the the ROPE boundaries, but the narrow confidence intervals are not fully covered by the ROPE;

  • else the decision whether to accept or reject practical equivalence is undecided (i.e. when effects are not statistically significant and the narrow confidence intervals overlaps the ROPE).

• "cet" - Conditional Equivalence Testing (Campbell/Gustafson 2018)

  The Conditional Equivalence Testing as described by Campbell and Gustafson 2018. According to this rule, practical equivalence is rejected when the coefficient is statistically significant. When the effect is not significant and the narrow confidence intervals are completely inside the ROPE, we accept (i.e. assume) practical equivalence, else it is undecided.

Levels of Confidence Intervals used for Equivalence Testing

For rule = "classic", "narrow" confidence intervals are used for equivalence testing. "Narrow" means, the the intervals is not 1 - alpha, but 1 - 2 * alpha. Thus, if ci = .95, alpha is assumed to be 0.05 and internally a ci-level of 0.90 is used. rule = "cet" uses both regular and narrow confidence intervals, while rule = "bayes" only uses the regular intervals.

p-Values

The equivalence p-value is the area of the (cumulative) confidence distribution that is outside of the region of equivalence. It can be interpreted as p-value for rejecting the alternative hypothesis and accepting the "null hypothesis" (i.e. assuming practical equivalence). That is, a high p-value means we reject the assumption of practical equivalence and accept the alternative hypothesis.

Second Generation p-Value (SGPV)

Second generation p-values (SGPV) were proposed as a statistic that represents the proportion of data-supported hypotheses that are also null hypotheses (Blume et al. 2018, Lakens and Delacre 2020). It represents the proportion of the full confidence interval range (assuming a normally or t-distributed, equal-tailed interval, based on the model) that is inside the ROPE. The SGPV ranges from zero to one. Higher values indicate that the effect is more likely to be practically equivalent ("not of interest").

Note that the assumed interval, which is used to calculate the SGPV, is an estimation of the full interval based on the chosen confidence level. For example, if the 95% confidence interval of a coefficient ranges from -1 to 1, the underlying full (normally or t-distributed) interval approximately ranges from -1.9 to 1.9, see also following code:

# simulate full normal distribution
out <- bayestestR::distribution_normal(10000, 0, 0.5)
# range of "full" distribution
range(out)
# range of 95% CI
round(quantile(out, probs = c(0.025, 0.975)), 2)

This ensures that the SGPV always refers to the general compatible parameter space of coefficients, independent from the confidence interval chosen for testing practical equivalence. Therefore, the SGPV of the full interval is similar to the ROPE coverage of Bayesian equivalence tests, see following code:

library(bayestestR)
library(brms)
m <- lm(mpg ~ gear + wt + cyl + hp, data = mtcars)
m2 <- brm(mpg ~ gear + wt + cyl + hp, data = mtcars)
# SGPV for frequentist models
equivalence_test(m)
# similar to ROPE coverage of Bayesian models
equivalence_test(m2)
# similar to ROPE coverage of simulated draws / bootstrap samples
equivalence_test(simulate_model(m))

ROPE range

Some attention is required for finding suitable values for the ROPE limits (argument range). See 'Details' in bayestestR::rope_range() for further information.

Value

A data frame.

Statistical inference - how to quantify evidence

There is no standardized approach to drawing conclusions based on the available data and statistical models. A frequently chosen but also much criticized approach is to evaluate results based on their statistical significance (Amrhein et al. 2017).

A more sophisticated way would be to test whether estimated effects exceed the "smallest effect size of interest", to avoid even the smallest effects being considered relevant simply because they are statistically significant, but clinically or practically irrelevant (Lakens et al. 2018, Lakens 2024).

A rather unconventional approach, which is nevertheless advocated by various authors, is to interpret results from classical regression models either in terms of probabilities, similar to the usual approach in Bayesian statistics (Schweder 2018; Schweder and Hjort 2003; Vos 2022) or in terms of relative measure of "evidence" or "compatibility" with the data (Greenland et al. 2022; Rafi and Greenland 2020), which nevertheless comes close to a probabilistic interpretation.

A more detailed discussion of this topic is found in the documentation of p_function().

The parameters package provides several options or functions to aid statistical inference. These are, for example:

• equivalence_test(), to compute the (conditional) equivalence test for frequentist models

• p_significance(), to compute the probability of practical significance, which can be conceptualized as a unidirectional equivalence test

• p_function(), or consonance function, to compute p-values and compatibility (confidence) intervals for statistical models

• the pd argument (setting pd = TRUE) in model_parameters() includes a column with the probability of direction, i.e. the probability that a parameter is strictly positive or negative. See bayestestR::p_direction() for details. If plotting is desired, the p_direction() function can be used, together with plot().

• the s_value argument (setting s_value = TRUE) in model_parameters() replaces the p-values with their related S-values (Rafi and Greenland 2020)

• finally, it is possible to generate distributions of model coefficients by generating bootstrap-samples (setting bootstrap = TRUE) or simulating draws from model coefficients using simulate_model(). These samples can then be treated as "posterior samples" and used in many functions from the bayestestR package.

Most of the above shown options or functions derive from methods originally implemented for Bayesian models (Makowski et al. 2019). However, assuming that model assumptions are met (which means, the model fits well to the data, the correct model is chosen that reflects the data generating process (distributional model family) etc.), it seems appropriate to interpret results from classical frequentist models in a "Bayesian way" (more details: documentation in p_function()).

Note

There is also a plot()-method implemented in the see-package.

References

• Amrhein, V., Korner-Nievergelt, F., and Roth, T. (2017). The earth is flat (p > 0.05): Significance thresholds and the crisis of unreplicable research. PeerJ, 5, e3544. doi:10.7717/peerj.3544

• Blume, J. D., D'Agostino McGowan, L., Dupont, W. D., & Greevy, R. A. (2018). Second-generation p-values: Improved rigor, reproducibility, & transparency in statistical analyses. PLOS ONE, 13(3), e0188299. https://doi.org/10.1371/journal.pone.0188299

• Campbell, H., & Gustafson, P. (2018). Conditional equivalence testing: An alternative remedy for publication bias. PLOS ONE, 13(4), e0195145. doi: 10.1371/journal.pone.0195145

• Greenland S, Rafi Z, Matthews R, Higgs M. To Aid Scientific Inference, Emphasize Unconditional Compatibility Descriptions of Statistics. (2022) https://arxiv.org/abs/1909.08583v7 (Accessed November 10, 2022)

• Kruschke, J. K. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press

• Kruschke, J. K. (2018). Rejecting or accepting parameter values in Bayesian estimation. Advances in Methods and Practices in Psychological Science, 1(2), 270-280. doi: 10.1177/2515245918771304

• Lakens, D. (2017). Equivalence Tests: A Practical Primer for t Tests, Correlations, and Meta-Analyses. Social Psychological and Personality Science, 8(4), 355–362. doi: 10.1177/1948550617697177

• Lakens, D. (2024). Improving Your Statistical Inferences (Version v1.5.1). Retrieved from https://lakens.github.io/statistical_inferences/. doi:10.5281/ZENODO.6409077

• Lakens, D., and Delacre, M. (2020). Equivalence Testing and the Second Generation P-Value. Meta-Psychology, 4. https://doi.org/10.15626/MP.2018.933

• Lakens, D., Scheel, A. M., and Isager, P. M. (2018). Equivalence Testing for Psychological Research: A Tutorial. Advances in Methods and Practices in Psychological Science, 1(2), 259–269. doi:10.1177/2515245918770963

• Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., and Lüdecke, D. (2019). Indices of Effect Existence and Significance in the Bayesian Framework. Frontiers in Psychology, 10, 2767. doi:10.3389/fpsyg.2019.02767

• Pernet, C. (2017). Null hypothesis significance testing: A guide to commonly misunderstood concepts and recommendations for good practice. F1000Research, 4, 621. doi: 10.12688/f1000research.6963.5

• Rafi Z, Greenland S. Semantic and cognitive tools to aid statistical science: replace confidence and significance by compatibility and surprise. BMC Medical Research Methodology (2020) 20:244.

• Schweder T. Confidence is epistemic probability for empirical science. Journal of Statistical Planning and Inference (2018) 195:116–125. doi:10.1016/j.jspi.2017.09.016

• Schweder T, Hjort NL. Frequentist analogues of priors and posteriors. In Stigum, B. (ed.), Econometrics and the Philosophy of Economics: Theory Data Confrontation in Economics, pp. 285-217. Princeton University Press, Princeton, NJ, 2003

• Vos P, Holbert D. Frequentist statistical inference without repeated sampling. Synthese 200, 89 (2022). doi:10.1007/s11229-022-03560-x

See Also

For more details, see bayestestR::equivalence_test(). Further readings can be found in the references. See also p_significance() for a unidirectional equivalence test.

Examples

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

# default rule
equivalence_test(model)

# using heteroscedasticity-robust standard errors
equivalence_test(model, vcov = "HC3")

# conditional equivalence test
equivalence_test(model, rule = "cet")

# plot method
if (require("see", quietly = TRUE)) {
  result <- equivalence_test(model)
  plot(result)
}

Description

The functions principal_components() and factor_analysis() can be used to perform a principal component analysis (PCA) or a factor analysis (FA). They return the loadings as a data frame, and various methods and functions are available to access / display other information (see the Details section).

Usage

factor_analysis(
  x,
  n = "auto",
  rotation = "none",
  sort = FALSE,
  threshold = NULL,
  standardize = TRUE,
  cor = NULL,
  ...
)

principal_components(
  x,
  n = "auto",
  rotation = "none",
  sparse = FALSE,
  sort = FALSE,
  threshold = NULL,
  standardize = TRUE,
  ...
)

rotated_data(pca_results, verbose = TRUE)

## S3 method for class 'parameters_efa'
predict(
  object,
  newdata = NULL,
  names = NULL,
  keep_na = TRUE,
  verbose = TRUE,
  ...
)

## S3 method for class 'parameters_efa'
print(x, digits = 2, sort = FALSE, threshold = NULL, labels = NULL, ...)

## S3 method for class 'parameters_efa'
sort(x, ...)

closest_component(pca_results)

Arguments

Details

Methods and Utilities

• n_components() and n_factors() automatically estimates the optimal number of dimensions to retain.

• performance::check_factorstructure() checks the suitability of the data for factor analysis using the sphericity (see performance::check_sphericity_bartlett()) and the KMO (see performance::check_kmo()) measure.

• performance::check_itemscale() computes various measures of internal consistencies applied to the (sub)scales (i.e., components) extracted from the PCA.

• Running summary() returns information related to each component/factor, such as the explained variance and the Eivenvalues.

• Running get_scores() computes scores for each subscale.

• Running closest_component() will return a numeric vector with the assigned component index for each column from the original data frame.

• Running rotated_data() will return the rotated data, including missing values, so it matches the original data frame.

• Running plot() visually displays the loadings (that requires the see-package to work).

Complexity

Complexity represents the number of latent components needed to account for the observed variables. Whereas a perfect simple structure solution has a complexity of 1 in that each item would only load on one factor, a solution with evenly distributed items has a complexity greater than 1 (Hofman, 1978; Pettersson and Turkheimer, 2010).

Uniqueness

Uniqueness represents the variance that is 'unique' to the variable and not shared with other variables. It is equal to ⁠1 – communality⁠ (variance that is shared with other variables). A uniqueness of 0.20 suggests that ⁠20%⁠ or that variable's variance is not shared with other variables in the overall factor model. The greater 'uniqueness' the lower the relevance of the variable in the factor model.

MSA

MSA represents the Kaiser-Meyer-Olkin Measure of Sampling Adequacy (Kaiser and Rice, 1974) for each item. It indicates whether there is enough data for each factor give reliable results for the PCA. The value should be > 0.6, and desirable values are > 0.8 (Tabachnick and Fidell, 2013).

PCA or FA?

There is a simplified rule of thumb that may help do decide whether to run a factor analysis or a principal component analysis:

• Run factor analysis if you assume or wish to test a theoretical model of latent factors causing observed variables.

• Run principal component analysis If you want to simply reduce your correlated observed variables to a smaller set of important independent composite variables.

(Source: CrossValidated)

Computing Item Scores

Use get_scores() to compute scores for the "subscales" represented by the extracted principal components. get_scores() takes the results from principal_components() and extracts the variables for each component found by the PCA. Then, for each of these "subscales", raw means are calculated (which equals adding up the single items and dividing by the number of items). This results in a sum score for each component from the PCA, which is on the same scale as the original, single items that were used to compute the PCA. One can also use predict() to back-predict scores for each component, to which one can provide newdata or a vector of names for the components.

Explained Variance and Eingenvalues

Use summary() to get the Eigenvalues and the explained variance for each extracted component. The eigenvectors and eigenvalues represent the "core" of a PCA: The eigenvectors (the principal components) determine the directions of the new feature space, and the eigenvalues determine their magnitude. In other words, the eigenvalues explain the variance of the data along the new feature axes.

Value

A data frame of loadings.

References

• Kaiser, H.F. and Rice. J. (1974). Little jiffy, mark iv. Educational and Psychological Measurement, 34(1):111–117

• Hofmann, R. (1978). Complexity and simplicity as objective indices descriptive of factor solutions. Multivariate Behavioral Research, 13:2, 247-250, doi:10.1207/s15327906mbr1302_9

• Pettersson, E., & Turkheimer, E. (2010). Item selection, evaluation, and simple structure in personality data. Journal of research in personality, 44(4), 407-420, doi:10.1016/j.jrp.2010.03.002

• Tabachnick, B. G., and Fidell, L. S. (2013). Using multivariate statistics (6th ed.). Boston: Pearson Education.

Examples

library(parameters)


# Principal Component Analysis (PCA) -------------------
principal_components(mtcars[, 1:7], n = "all", threshold = 0.2)

# Automated number of components
principal_components(mtcars[, 1:4], n = "auto")

# labels can be useful if variable names are not self-explanatory
print(
  principal_components(mtcars[, 1:4], n = "auto"),
  labels = c(
    "Miles/(US) gallon",
    "Number of cylinders",
    "Displacement (cu.in.)",
    "Gross horsepower"
  )
)

# Sparse PCA
principal_components(mtcars[, 1:7], n = 4, sparse = TRUE)
principal_components(mtcars[, 1:7], n = 4, sparse = "robust")

# Rotated PCA
principal_components(mtcars[, 1:7],
  n = 2, rotation = "oblimin",
  threshold = "max", sort = TRUE
)
principal_components(mtcars[, 1:7], n = 2, threshold = 2, sort = TRUE)

pca <- principal_components(mtcars[, 1:5], n = 2, rotation = "varimax")
pca # Print loadings
summary(pca) # Print information about the factors
predict(pca, names = c("Component1", "Component2")) # Back-predict scores

# which variables from the original data belong to which extracted component?
closest_component(pca)


# Factor Analysis (FA) ------------------------

factor_analysis(mtcars[, 1:7], n = "all", threshold = 0.2)
factor_analysis(mtcars[, 1:7], n = 2, rotation = "oblimin", threshold = "max", sort = TRUE)
factor_analysis(mtcars[, 1:7], n = 2, threshold = 2, sort = TRUE)

efa <- factor_analysis(mtcars[, 1:5], n = 2)
summary(efa)
predict(efa, verbose = FALSE)


# Automated number of components
factor_analysis(mtcars[, 1:4], n = "auto")

Description

A sample data set, used in tests and some examples.

Description

Format the name of the degrees-of-freedom adjustment methods.

Usage

format_df_adjust(
  method,
  approx_string = "-approximated",
  dof_string = " degrees of freedom"
)

Arguments

Value

A formatted string.

Examples

library(parameters)

format_df_adjust("kenward")
format_df_adjust("kenward", approx_string = "", dof_string = " DoF")

Description

Format order.

Usage

format_order(order, textual = TRUE, ...)

Arguments

Value

A formatted string.

Examples

format_order(2)
format_order(8)
format_order(25, textual = FALSE)

Description

Format the name of the p-value adjustment methods.

Usage

format_p_adjust(method)

Arguments

Value

A string with the full surname(s) of the author(s), including year of publication, for the adjustment-method.

Examples

library(parameters)

format_p_adjust("holm")
format_p_adjust("bonferroni")

Description

This functions formats the names of model parameters (coefficients) to make them more human-readable.

Usage

format_parameters(model, ...)

## Default S3 method:
format_parameters(model, brackets = c("[", "]"), ...)

Arguments

Value

A (names) character vector with formatted parameter names. The value names refer to the original names of the coefficients.

Interpretation of Interaction Terms

Note that the interpretation of interaction terms depends on many characteristics of the model. The number of parameters, and overall performance of the model, can differ or not between a * b, a : b, and a / b, suggesting that sometimes interaction terms give different parameterizations of the same model, but other times it gives completely different models (depending on a or b being factors of covariates, included as main effects or not, etc.). Their interpretation depends of the full context of the model, which should not be inferred from the parameters table alone - rather, we recommend to use packages that calculate estimated marginal means or marginal effects, such as modelbased, emmeans, ggeffects, or marginaleffects. To raise awareness for this issue, you may use print(...,show_formula=TRUE) to add the model-specification to the output of the print() method for model_parameters().

Examples

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

model <- lm(Sepal.Length ~ Petal.Length + (Species / Sepal.Width), data = iris)
format_parameters(model)

model <- lm(Sepal.Length ~ Species + poly(Sepal.Width, 2), data = iris)
format_parameters(model)

model <- lm(Sepal.Length ~ Species + poly(Sepal.Width, 2, raw = TRUE), data = iris)
format_parameters(model)

Description

A print()-method for objects from compare_parameters().

Usage

## S3 method for class 'compare_parameters'
format(
  x,
  split_components = TRUE,
  select = NULL,
  digits = 2,
  ci_digits = digits,
  p_digits = 3,
  ci_width = NULL,
  ci_brackets = NULL,
  zap_small = FALSE,
  format = NULL,
  groups = NULL,
  engine = NULL,
  ...
)

## S3 method for class 'compare_parameters'
print(
  x,
  split_components = TRUE,
  caption = NULL,
  subtitle = NULL,
  footer = NULL,
  digits = 2,
  ci_digits = digits,
  p_digits = 3,
  zap_small = FALSE,
  groups = NULL,
  column_width = NULL,
  ci_brackets = c("(", ")"),
  select = NULL,
  ...
)

## S3 method for class 'compare_parameters'
print_html(
  x,
  caption = NULL,
  subtitle = NULL,
  footer = NULL,
  digits = 2,
  ci_digits = digits,
  p_digits = 3,
  zap_small = FALSE,
  groups = NULL,
  select = NULL,
  ci_brackets = c("(", ")"),
  font_size = "100%",
  line_padding = 4,
  column_labels = NULL,
  engine = "gt",
  ...
)

## S3 method for class 'compare_parameters'
print_md(
  x,
  digits = 2,
  ci_digits = digits,
  p_digits = 3,
  caption = NULL,
  subtitle = NULL,
  footer = NULL,
  select = NULL,
  split_components = TRUE,
  ci_brackets = c("(", ")"),
  zap_small = FALSE,
  groups = NULL,
  engine = "tt",
  ...
)

Arguments

Value

Invisibly returns the original input object.

Global Options to Customize Messages and Tables when Printing

The verbose argument can be used to display or silence messages and warnings for the different functions in the parameters package. However, some messages providing additional information can be displayed or suppressed using options():

• parameters_summary: options(parameters_summary = TRUE) will override the summary argument in model_parameters() and always show the model summary for non-mixed models.

• parameters_mixed_summary: options(parameters_mixed_summary = TRUE) will override the summary argument in model_parameters() for mixed models, and will then always show the model summary.

• parameters_cimethod: options(parameters_cimethod = TRUE) will show the additional information about the approximation method used to calculate confidence intervals and p-values. Set to FALSE to hide this message when printing model_parameters() objects.

• parameters_exponentiate: options(parameters_exponentiate = TRUE) will show the additional information on how to interpret coefficients of models with log-transformed response variables or with log-/logit-links when the exponentiate argument in model_parameters() is not TRUE. Set this option to FALSE to hide this message when printing model_parameters() objects.

There are further options that can be used to modify the default behaviour for printed outputs:

• parameters_labels: options(parameters_labels = TRUE) will use variable and value labels for pretty names, if data is labelled. If no labels available, default pretty names are used.

• parameters_interaction: ⁠options(parameters_interaction = <character>)⁠ will replace the interaction mark (by default, *) with the related character.

• parameters_select: ⁠options(parameters_select = <value>)⁠ will set the default for the select argument. See argument's documentation for available options.

• easystats_html_engine: options(easystats_html_engine = "gt") will set the default HTML engine for tables to gt, i.e. the gt package is used to create HTML tables. If set to tt, the tinytable package is used.

• insight_use_symbols: options(insight_use_symbols = TRUE) will try to print unicode-chars for symbols as column names, wherever possible (e.g., ω instead of Omega).

Examples

data(iris)
lm1 <- lm(Sepal.Length ~ Species, data = iris)
lm2 <- lm(Sepal.Length ~ Species + Petal.Length, data = iris)

# custom style
result <- compare_parameters(lm1, lm2, select = "{estimate}{stars} ({se})")
print(result)

# custom style, in HTML
result <- compare_parameters(lm1, lm2, select = "{estimate}<br>({se})|{p}")
print_html(result)

Description

A print()-method for objects from model_parameters().

Usage

## S3 method for class 'parameters_model'
format(
  x,
  pretty_names = TRUE,
  split_components = TRUE,
  select = NULL,
  digits = 2,
  ci_digits = digits,
  p_digits = 3,
  ci_width = NULL,
  ci_brackets = NULL,
  zap_small = FALSE,
  format = NULL,
  groups = NULL,
  include_reference = FALSE,
  ...
)

## S3 method for class 'parameters_model'
print(
  x,
  pretty_names = TRUE,
  split_components = TRUE,
  select = NULL,
  caption = NULL,
  footer = NULL,
  digits = 2,
  ci_digits = digits,
  p_digits = 3,
  footer_digits = 3,
  show_sigma = FALSE,
  show_formula = FALSE,
  zap_small = FALSE,
  groups = NULL,
  column_width = NULL,
  ci_brackets = c("[", "]"),
  include_reference = FALSE,
  ...
)

## S3 method for class 'parameters_model'
summary(object, ...)

## S3 method for class 'parameters_model'
print_html(
  x,
  pretty_names = TRUE,
  split_components = TRUE,
  select = NULL,
  caption = NULL,
  subtitle = NULL,
  footer = NULL,
  align = NULL,
  digits = 2,
  ci_digits = digits,
  p_digits = 3,
  footer_digits = 3,
  ci_brackets = c("(", ")"),
  show_sigma = FALSE,
  show_formula = FALSE,
  zap_small = FALSE,
  groups = NULL,
  font_size = "100%",
  line_padding = 4,
  column_labels = NULL,
  include_reference = FALSE,
  verbose = TRUE,
  ...
)

## S3 method for class 'parameters_model'
print_md(
  x,
  pretty_names = TRUE,
  split_components = TRUE,
  select = NULL,
  caption = NULL,
  subtitle = NULL,
  footer = NULL,
  align = NULL,
  digits = 2,
  ci_digits = digits,
  p_digits = 3,
  footer_digits = 3,
  ci_brackets = c("(", ")"),
  show_sigma = FALSE,
  show_formula = FALSE,
  zap_small = FALSE,
  groups = NULL,
  include_reference = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Details

summary() is a convenient shortcut for print(object, select = "minimal", show_sigma = TRUE, show_formula = TRUE).

Value

Invisibly returns the original input object.

Global Options to Customize Messages and Tables when Printing

The verbose argument can be used to display or silence messages and warnings for the different functions in the parameters package. However, some messages providing additional information can be displayed or suppressed using options():

• parameters_summary: options(parameters_summary = TRUE) will override the summary argument in model_parameters() and always show the model summary for non-mixed models.

• parameters_mixed_summary: options(parameters_mixed_summary = TRUE) will override the summary argument in model_parameters() for mixed models, and will then always show the model summary.

• parameters_cimethod: options(parameters_cimethod = TRUE) will show the additional information about the approximation method used to calculate confidence intervals and p-values. Set to FALSE to hide this message when printing model_parameters() objects.

• parameters_exponentiate: options(parameters_exponentiate = TRUE) will show the additional information on how to interpret coefficients of models with log-transformed response variables or with log-/logit-links when the exponentiate argument in model_parameters() is not TRUE. Set this option to FALSE to hide this message when printing model_parameters() objects.

There are further options that can be used to modify the default behaviour for printed outputs:

• parameters_labels: options(parameters_labels = TRUE) will use variable and value labels for pretty names, if data is labelled. If no labels available, default pretty names are used.

• parameters_interaction: ⁠options(parameters_interaction = <character>)⁠ will replace the interaction mark (by default, *) with the related character.

• parameters_select: ⁠options(parameters_select = <value>)⁠ will set the default for the select argument. See argument's documentation for available options.

• easystats_html_engine: options(easystats_html_engine = "gt") will set the default HTML engine for tables to gt, i.e. the gt package is used to create HTML tables. If set to tt, the tinytable package is used.

• insight_use_symbols: options(insight_use_symbols = TRUE) will try to print unicode-chars for symbols as column names, wherever possible (e.g., ω instead of Omega).

Interpretation of Interaction Terms

Note that the interpretation of interaction terms depends on many characteristics of the model. The number of parameters, and overall performance of the model, can differ or not between a * b, a : b, and a / b, suggesting that sometimes interaction terms give different parameterizations of the same model, but other times it gives completely different models (depending on a or b being factors of covariates, included as main effects or not, etc.). Their interpretation depends of the full context of the model, which should not be inferred from the parameters table alone - rather, we recommend to use packages that calculate estimated marginal means or marginal effects, such as modelbased, emmeans, ggeffects, or marginaleffects. To raise awareness for this issue, you may use print(...,show_formula=TRUE) to add the model-specification to the output of the print() method for model_parameters().

Labeling the Degrees of Freedom

Throughout the parameters package, we decided to label the residual degrees of freedom df_error. The reason for this is that these degrees of freedom not always refer to the residuals. For certain models, they refer to the estimate error - in a linear model these are the same, but in - for instance - any mixed effects model, this isn't strictly true. Hence, we think that df_error is the most generic label for these degrees of freedom.

See Also

See also display().

Examples

library(parameters)
model <- glmmTMB::glmmTMB(
  count ~ spp + mined + (1 | site),
  ziformula = ~mined,
  family = poisson(),
  data = Salamanders
)
mp <- model_parameters(model)

print(mp, pretty_names = FALSE)

print(mp, split_components = FALSE)

print(mp, select = c("Parameter", "Coefficient", "SE"))

print(mp, select = "minimal")


# group parameters ------

data(iris)
model <- lm(
  Sepal.Width ~ Sepal.Length + Species + Petal.Length,
  data = iris
)
# don't select "Intercept" parameter
mp <- model_parameters(model, parameters = "^(?!\\(Intercept)")
groups <- list(
  "Focal Predictors" = c("Speciesversicolor", "Speciesvirginica"),
  "Controls" = c("Sepal.Length", "Petal.Length")
)
print(mp, groups = groups)

# or use row indices
print(mp, groups = list(
  "Focal Predictors" = c(1, 4),
  "Controls" = c(2, 3)
))

# only show coefficients, CI and p,
# put non-matched parameters to the end

data(mtcars)
mtcars$cyl <- as.factor(mtcars$cyl)
mtcars$gear <- as.factor(mtcars$gear)
model <- lm(mpg ~ hp + gear * vs + cyl + drat, data = mtcars)

# don't select "Intercept" parameter
mp <- model_parameters(model, parameters = "^(?!\\(Intercept)")
print(mp, groups = list(
  "Engine" = c("cyl6", "cyl8", "vs", "hp"),
  "Interactions" = c("gear4:vs", "gear5:vs")
))



# custom column layouts ------

data(iris)
lm1 <- lm(Sepal.Length ~ Species, data = iris)
lm2 <- lm(Sepal.Length ~ Species + Petal.Length, data = iris)

# custom style
result <- compare_parameters(lm1, lm2, select = "{estimate}{stars} ({se})")
print(result)


# custom style, in HTML
result <- compare_parameters(lm1, lm2, select = "{estimate}<br>({se})|{p}")
print_html(result)

Description

get_scores() takes n_items amount of items that load the most (either by loading cutoff or number) on a component, and then computes their average.

Usage

get_scores(x, n_items = NULL)

Arguments

Details

get_scores() takes the results from principal_components() and extracts the variables for each component found by the PCA. Then, for each of these "subscales", row means are calculated (which equals adding up the single items and dividing by the number of items). This results in a sum score for each component from the PCA, which is on the same scale as the original, single items that were used to compute the PCA.

Value

A data frame with subscales, which are average sum scores for all items from each component.

Examples

if (require("psych")) {
  pca <- principal_components(mtcars[, 1:7], n = 2, rotation = "varimax")

  # PCA extracted two components
  pca

  # assignment of items to each component
  closest_component(pca)

  # now we want to have sum scores for each component
  get_scores(pca)

  # compare to manually computed sum score for 2nd component, which
  # consists of items "hp" and "qsec"
  (mtcars$hp + mtcars$qsec) / 2
}

Description

Compute and extract model parameters. The available options and arguments depend on the modeling package and model class. Follow one of these links to read the model-specific documentation:

• Default method: lm, glm, stats, censReg, MASS, survey, ...

• Additive models: bamlss, gamlss, mgcv, scam, VGAM, Gam, gamm, ...

• ANOVA: afex, aov, anova, ...

• Bayesian: BayesFactor, blavaan, brms, MCMCglmm, posterior, rstanarm, bayesQR, bcplm, BGGM, blmrm, blrm, mcmc.list, MCMCglmm, ...

• Clustering: hclust, kmeans, mclust, pam, ...

• Correlations, t-tests, etc.: lmtest, htest, pairwise.htest, ...

• Meta-Analysis: metaBMA, metafor, metaplus, ...

• Mixed models: cplm, glmmTMB, lme4, lmerTest, nlme, ordinal, robustlmm, spaMM, mixed, MixMod, ...

• Multinomial, ordinal and cumulative link: brglm2, DirichletReg, nnet, ordinal, mlm, ...

• Multiple imputation: mice

• PCA, FA, CFA, SEM: FactoMineR, lavaan, psych, sem, ...

• Zero-inflated and hurdle: cplm, mhurdle, pscl, ...

• Other models: aod, bbmle, betareg, emmeans, epiR, ggeffects, glmx, ivfixed, ivprobit, JRM, lmodel2, logitsf, marginaleffects, margins, maxLik, mediation, mfx, multcomp, mvord, plm, PMCMRplus, quantreg, selection, systemfit, tidymodels, varEST, WRS2, bfsl, deltaMethod, fitdistr, mjoint, mle, model.avg, ...

Usage

model_parameters(model, ...)

parameters(model, ...)

Arguments

Value

A data frame of indices related to the model's parameters.

Standardization of model coefficients

Standardization is based on standardize_parameters(). In case of standardize = "refit", the data used to fit the model will be standardized and the model is completely refitted. In such cases, standard errors and confidence intervals refer to the standardized coefficient. The default, standardize = "refit", never standardizes categorical predictors (i.e. factors), which may be a different behaviour compared to other R packages or other software packages (like SPSS). To mimic behaviour of SPSS or packages such as lm.beta, use standardize = "basic".

Standardization Methods

• refit: This method is based on a complete model re-fit with a standardized version of the data. Hence, this method is equal to standardizing the variables before fitting the model. It is the "purest" and the most accurate (Neter et al., 1989), but it is also the most computationally costly and long (especially for heavy models such as Bayesian models). This method is particularly recommended for complex models that include interactions or transformations (e.g., polynomial or spline terms). The robust (default to FALSE) argument enables a robust standardization of data, i.e., based on the median and MAD instead of the mean and SD. See datawizard::standardize() for more details. Note that standardize_parameters(method = "refit") may not return the same results as fitting a model on data that has been standardized with standardize(); standardize_parameters() used the data used by the model fitting function, which might not be same data if there are missing values. see the remove_na argument in standardize().

• posthoc: Post-hoc standardization of the parameters, aiming at emulating the results obtained by "refit" without refitting the model. The coefficients are divided by the standard deviation (or MAD if robust) of the outcome (which becomes their expression 'unit'). Then, the coefficients related to numeric variables are additionally multiplied by the standard deviation (or MAD if robust) of the related terms, so that they correspond to changes of 1 SD of the predictor (e.g., "A change in 1 SD of x is related to a change of 0.24 of the SD of y). This does not apply to binary variables or factors, so the coefficients are still related to changes in levels. This method is not accurate and tend to give aberrant results when interactions are specified.

• basic: This method is similar to method = "posthoc", but treats all variables as continuous: it also scales the coefficient by the standard deviation of model's matrix' parameter of factors levels (transformed to integers) or binary predictors. Although being inappropriate for these cases, this method is the one implemented by default in other software packages, such as lm.beta::lm.beta().

• smart (Standardization of Model's parameters with Adjustment, Reconnaissance and Transformation - experimental): Similar to method = "posthoc" in that it does not involve model refitting. The difference is that the SD (or MAD if robust) of the response is computed on the relevant section of the data. For instance, if a factor with 3 levels A (the intercept), B and C is entered as a predictor, the effect corresponding to B vs. A will be scaled by the variance of the response at the intercept only. As a results, the coefficients for effects of factors are similar to a Glass' delta.

• pseudo (for 2-level (G)LMMs only): In this (post-hoc) method, the response and the predictor are standardized based on the level of prediction (levels are detected with performance::check_heterogeneity_bias()): Predictors are standardized based on their SD at level of prediction (see also datawizard::demean()); The outcome (in linear LMMs) is standardized based on a fitted random-intercept-model, where sqrt(random-intercept-variance) is used for level 2 predictors, and sqrt(residual-variance) is used for level 1 predictors (Hoffman 2015, page 342). A warning is given when a within-group variable is found to have access between-group variance.

See also package vignette.

Labeling the Degrees of Freedom

Throughout the parameters package, we decided to label the residual degrees of freedom df_error. The reason for this is that these degrees of freedom not always refer to the residuals. For certain models, they refer to the estimate error - in a linear model these are the same, but in - for instance - any mixed effects model, this isn't strictly true. Hence, we think that df_error is the most generic label for these degrees of freedom.

Confidence intervals and approximation of degrees of freedom

There are different ways of approximating the degrees of freedom depending on different assumptions about the nature of the model and its sampling distribution. The ci_method argument modulates the method for computing degrees of freedom (df) that are used to calculate confidence intervals (CI) and the related p-values. Following options are allowed, depending on the model class:

Classical methods:

Classical inference is generally based on the Wald method. The Wald approach to inference computes a test statistic by dividing the parameter estimate by its standard error (Coefficient / SE), then comparing this statistic against a t- or normal distribution. This approach can be used to compute CIs and p-values.

"wald":

• Applies to non-Bayesian models. For linear models, CIs computed using the Wald method (SE and a t-distribution with residual df); p-values computed using the Wald method with a t-distribution with residual df. For other models, CIs computed using the Wald method (SE and a normal distribution); p-values computed using the Wald method with a normal distribution.

"normal"

• Applies to non-Bayesian models. Compute Wald CIs and p-values, but always use a normal distribution.

"residual"

• Applies to non-Bayesian models. Compute Wald CIs and p-values, but always use a t-distribution with residual df when possible. If the residual df for a model cannot be determined, a normal distribution is used instead.

Methods for mixed models:

Compared to fixed effects (or single-level) models, determining appropriate df for Wald-based inference in mixed models is more difficult. See the R GLMM FAQ for a discussion.

Several approximate methods for computing df are available, but you should also consider instead using profile likelihood ("profile") or bootstrap ("⁠boot"⁠) CIs and p-values instead.

"satterthwaite"

• Applies to linear mixed models. CIs computed using the Wald method (SE and a t-distribution with Satterthwaite df); p-values computed using the Wald method with a t-distribution with Satterthwaite df.

"kenward"

• Applies to linear mixed models. CIs computed using the Wald method (Kenward-Roger SE and a t-distribution with Kenward-Roger df); p-values computed using the Wald method with Kenward-Roger SE and t-distribution with Kenward-Roger df.

"ml1"

• Applies to linear mixed models. CIs computed using the Wald method (SE and a t-distribution with m-l-1 approximated df); p-values computed using the Wald method with a t-distribution with m-l-1 approximated df. See ci_ml1().

"betwithin"

• Applies to linear mixed models and generalized linear mixed models. CIs computed using the Wald method (SE and a t-distribution with between-within df); p-values computed using the Wald method with a t-distribution with between-within df. See ci_betwithin().

Likelihood-based methods:

Likelihood-based inference is based on comparing the likelihood for the maximum-likelihood estimate to the the likelihood for models with one or more parameter values changed (e.g., set to zero or a range of alternative values). Likelihood ratios for the maximum-likelihood and alternative models are compared to a $\chi$-squared distribution to compute CIs and p-values.

"profile"

• Applies to non-Bayesian models of class glm, polr, merMod or glmmTMB. CIs computed by profiling the likelihood curve for a parameter, using linear interpolation to find where likelihood ratio equals a critical value; p-values computed using the Wald method with a normal-distribution (note: this might change in a future update!)

"uniroot"

• Applies to non-Bayesian models of class glmmTMB. CIs computed by profiling the likelihood curve for a parameter, using root finding to find where likelihood ratio equals a critical value; p-values computed using the Wald method with a normal-distribution (note: this might change in a future update!)

Methods for bootstrapped or Bayesian models:

Bootstrap-based inference is based on resampling and refitting the model to the resampled datasets. The distribution of parameter estimates across resampled datasets is used to approximate the parameter's sampling distribution. Depending on the type of model, several different methods for bootstrapping and constructing CIs and p-values from the bootstrap distribution are available.

For Bayesian models, inference is based on drawing samples from the model posterior distribution.

"quantile" (or "eti")

• Applies to all models (including Bayesian models). For non-Bayesian models, only applies if bootstrap = TRUE. CIs computed as equal tailed intervals using the quantiles of the bootstrap or posterior samples; p-values are based on the probability of direction. See bayestestR::eti().

"hdi"

• Applies to all models (including Bayesian models). For non-Bayesian models, only applies if bootstrap = TRUE. CIs computed as highest density intervals for the bootstrap or posterior samples; p-values are based on the probability of direction. See bayestestR::hdi().

"bci" (or "bcai")

• Applies to all models (including Bayesian models). For non-Bayesian models, only applies if bootstrap = TRUE. CIs computed as bias corrected and accelerated intervals for the bootstrap or posterior samples; p-values are based on the probability of direction. See bayestestR::bci().

"si"

• Applies to Bayesian models with proper priors. CIs computed as support intervals comparing the posterior samples against the prior samples; p-values are based on the probability of direction. See bayestestR::si().

"boot"

• Applies to non-Bayesian models of class merMod. CIs computed using parametric bootstrapping (simulating data from the fitted model); p-values computed using the Wald method with a normal-distribution) (note: this might change in a future update!).

For all iteration-based methods other than "boot" ("hdi", "quantile", "ci", "eti", "si", "bci", "bcai"), p-values are based on the probability of direction (bayestestR::p_direction()), which is converted into a p-value using bayestestR::pd_to_p().

Statistical inference - how to quantify evidence

There is no standardized approach to drawing conclusions based on the available data and statistical models. A frequently chosen but also much criticized approach is to evaluate results based on their statistical significance (Amrhein et al. 2017).

A more sophisticated way would be to test whether estimated effects exceed the "smallest effect size of interest", to avoid even the smallest effects being considered relevant simply because they are statistically significant, but clinically or practically irrelevant (Lakens et al. 2018, Lakens 2024).

A rather unconventional approach, which is nevertheless advocated by various authors, is to interpret results from classical regression models either in terms of probabilities, similar to the usual approach in Bayesian statistics (Schweder 2018; Schweder and Hjort 2003; Vos 2022) or in terms of relative measure of "evidence" or "compatibility" with the data (Greenland et al. 2022; Rafi and Greenland 2020), which nevertheless comes close to a probabilistic interpretation.

A more detailed discussion of this topic is found in the documentation of p_function().

The parameters package provides several options or functions to aid statistical inference. These are, for example:

• equivalence_test(), to compute the (conditional) equivalence test for frequentist models

• p_significance(), to compute the probability of practical significance, which can be conceptualized as a unidirectional equivalence test

• p_function(), or consonance function, to compute p-values and compatibility (confidence) intervals for statistical models

• the pd argument (setting pd = TRUE) in model_parameters() includes a column with the probability of direction, i.e. the probability that a parameter is strictly positive or negative. See bayestestR::p_direction() for details. If plotting is desired, the p_direction() function can be used, together with plot().

• the s_value argument (setting s_value = TRUE) in model_parameters() replaces the p-values with their related S-values (Rafi and Greenland 2020)

• finally, it is possible to generate distributions of model coefficients by generating bootstrap-samples (setting bootstrap = TRUE) or simulating draws from model coefficients using simulate_model(). These samples can then be treated as "posterior samples" and used in many functions from the bayestestR package.

Most of the above shown options or functions derive from methods originally implemented for Bayesian models (Makowski et al. 2019). However, assuming that model assumptions are met (which means, the model fits well to the data, the correct model is chosen that reflects the data generating process (distributional model family) etc.), it seems appropriate to interpret results from classical frequentist models in a "Bayesian way" (more details: documentation in p_function()).

Interpretation of Interaction Terms

Note that the interpretation of interaction terms depends on many characteristics of the model. The number of parameters, and overall performance of the model, can differ or not between a * b, a : b, and a / b, suggesting that sometimes interaction terms give different parameterizations of the same model, but other times it gives completely different models (depending on a or b being factors of covariates, included as main effects or not, etc.). Their interpretation depends of the full context of the model, which should not be inferred from the parameters table alone - rather, we recommend to use packages that calculate estimated marginal means or marginal effects, such as modelbased, emmeans, ggeffects, or marginaleffects. To raise awareness for this issue, you may use print(...,show_formula=TRUE) to add the model-specification to the output of the print() method for model_parameters().

Global Options to Customize Messages and Tables when Printing

The verbose argument can be used to display or silence messages and warnings for the different functions in the parameters package. However, some messages providing additional information can be displayed or suppressed using options():

• parameters_summary: options(parameters_summary = TRUE) will override the summary argument in model_parameters() and always show the model summary for non-mixed models.

• parameters_mixed_summary: options(parameters_mixed_summary = TRUE) will override the summary argument in model_parameters() for mixed models, and will then always show the model summary.

• parameters_cimethod: options(parameters_cimethod = TRUE) will show the additional information about the approximation method used to calculate confidence intervals and p-values. Set to FALSE to hide this message when printing model_parameters() objects.

• parameters_exponentiate: options(parameters_exponentiate = TRUE) will show the additional information on how to interpret coefficients of models with log-transformed response variables or with log-/logit-links when the exponentiate argument in model_parameters() is not TRUE. Set this option to FALSE to hide this message when printing model_parameters() objects.

There are further options that can be used to modify the default behaviour for printed outputs:

• parameters_labels: options(parameters_labels = TRUE) will use variable and value labels for pretty names, if data is labelled. If no labels available, default pretty names are used.

• parameters_interaction: ⁠options(parameters_interaction = <character>)⁠ will replace the interaction mark (by default, *) with the related character.

• parameters_select: ⁠options(parameters_select = <value>)⁠ will set the default for the select argument. See argument's documentation for available options.

• easystats_html_engine: options(easystats_html_engine = "gt") will set the default HTML engine for tables to gt, i.e. the gt package is used to create HTML tables. If set to tt, the tinytable package is used.

• insight_use_symbols: options(insight_use_symbols = TRUE) will try to print unicode-chars for symbols as column names, wherever possible (e.g., ω instead of Omega).

Note

The print() method has several arguments to tweak the output. There is also a plot()-method implemented in the see-package, and a dedicated method for use inside rmarkdown files, print_md().

For developers, if speed performance is an issue, you can use the (undocumented) pretty_names argument, e.g. model_parameters(..., pretty_names = FALSE). This will skip the formatting of the coefficient names and make model_parameters() faster.

References

• Amrhein, V., Korner-Nievergelt, F., and Roth, T. (2017). The earth is flat (p > 0.05): Significance thresholds and the crisis of unreplicable research. PeerJ, 5, e3544. doi:10.7717/peerj.3544

• Greenland S, Rafi Z, Matthews R, Higgs M. To Aid Scientific Inference, Emphasize Unconditional Compatibility Descriptions of Statistics. (2022) https://arxiv.org/abs/1909.08583v7 (Accessed November 10, 2022)

• Hoffman, L. (2015). Longitudinal analysis: Modeling within-person fluctuation and change. Routledge.

• Lakens, D. (2024). Improving Your Statistical Inferences (Version v1.5.1). Retrieved from https://lakens.github.io/statistical_inferences/. doi:10.5281/ZENODO.6409077

• Lakens, D., Scheel, A. M., and Isager, P. M. (2018). Equivalence Testing for Psychological Research: A Tutorial. Advances in Methods and Practices in Psychological Science, 1(2), 259–269. doi:10.1177/2515245918770963

• Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., and Lüdecke, D. (2019). Indices of Effect Existence and Significance in the Bayesian Framework. Frontiers in Psychology, 10, 2767. doi:10.3389/fpsyg.2019.02767

• Neter, J., Wasserman, W., and Kutner, M. H. (1989). Applied linear regression models.

• Rafi Z, Greenland S. Semantic and cognitive tools to aid statistical science: replace confidence and significance by compatibility and surprise. BMC Medical Research Methodology (2020) 20:244.

• Schweder T. Confidence is epistemic probability for empirical science. Journal of Statistical Planning and Inference (2018) 195:116–125. doi:10.1016/j.jspi.2017.09.016

• Schweder T, Hjort NL. Frequentist analogues of priors and posteriors. In Stigum, B. (ed.), Econometrics and the Philosophy of Economics: Theory Data Confrontation in Economics, pp. 285-217. Princeton University Press, Princeton, NJ, 2003

• Vos P, Holbert D. Frequentist statistical inference without repeated sampling. Synthese 200, 89 (2022). doi:10.1007/s11229-022-03560-x

See Also

insight::standardize_names() to rename columns into a consistent, standardized naming scheme.

Description

Parameters from ANOVAs

Usage

## S3 method for class 'aov'
model_parameters(
  model,
  type = NULL,
  df_error = NULL,
  ci = NULL,
  alternative = NULL,
  test = NULL,
  power = FALSE,
  es_type = NULL,
  keep = NULL,
  drop = NULL,
  table_wide = FALSE,
  verbose = TRUE,
  ...
)

## S3 method for class 'afex_aov'
model_parameters(
  model,
  es_type = NULL,
  df_error = NULL,
  type = NULL,
  keep = NULL,
  drop = NULL,
  verbose = TRUE,
  ...
)

Arguments

Details

• For an object of class htest, data is extracted via insight::get_data(), and passed to the relevant function according to:

  • A t-test depending on type: "cohens_d" (default), "hedges_g", or one of "p_superiority", "u1", "u2", "u3", "overlap".

    • For a Paired t-test: depending on type: "rm_rm", "rm_av", "rm_b", "rm_d", "rm_z".

  • A Chi-squared tests of independence or Fisher's Exact Test, depending on type: "cramers_v" (default), "tschuprows_t", "phi", "cohens_w", "pearsons_c", "cohens_h", "oddsratio", "riskratio", "arr", or "nnt".

  • A Chi-squared tests of goodness-of-fit, depending on type: "fei" (default) "cohens_w", "pearsons_c"

  • A One-way ANOVA test, depending on type: "eta" (default), "omega" or "epsilon" -squared, "f", or "f2".

  • A McNemar test returns Cohen's g.

  • A Wilcoxon test depending on type: returns "rank_biserial" correlation (default) or one of "p_superiority", "vda", "u2", "u3", "overlap".

  • A Kruskal-Wallis test depending on type: "epsilon" (default) or "eta".

  • A Friedman test returns Kendall's W. (Where applicable, ci and alternative are taken from the htest if not otherwise provided.)

• For an object of class BFBayesFactor, using bayestestR::describe_posterior(),

  • A t-test depending on type: "cohens_d" (default) or one of "p_superiority", "u1", "u2", "u3", "overlap".

  • A correlation test returns r.

  • A contingency table test, depending on type: "cramers_v" (default), "phi", "tschuprows_t", "cohens_w", "pearsons_c", "cohens_h", "oddsratio", or "riskratio", "arr", or "nnt".

  • A proportion test returns p.

• Objects of class anova, aov, aovlist or afex_aov, depending on type: "eta" (default), "omega" or "epsilon" -squared, "f", or "f2".

• Other objects are passed to parameters::standardize_parameters().

For statistical models it is recommended to directly use the listed functions, for the full range of options they provide.

Value

A data frame of indices related to the model's parameters.

Note

For ANOVA-tables from mixed models (i.e. anova(lmer())), only partial or adjusted effect sizes can be computed. Note that type 3 ANOVAs with interactions involved only give sensible and informative results when covariates are mean-centred and factors are coded with orthogonal contrasts (such as those produced by contr.sum, contr.poly, or contr.helmert, but not by the default contr.treatment).

Examples

df <- iris
df$Sepal.Big <- ifelse(df$Sepal.Width >= 3, "Yes", "No")

model <- aov(Sepal.Length ~ Sepal.Big, data = df)
model_parameters(model)

model_parameters(model, es_type = c("omega", "eta"), ci = 0.9)

model <- anova(lm(Sepal.Length ~ Sepal.Big, data = df))
model_parameters(model)
model_parameters(
  model,
  es_type = c("omega", "eta", "epsilon"),
  alternative = "greater"
)

model <- aov(Sepal.Length ~ Sepal.Big + Error(Species), data = df)
model_parameters(model)



df <- iris
df$Sepal.Big <- ifelse(df$Sepal.Width >= 3, "Yes", "No")
mm <- lme4::lmer(Sepal.Length ~ Sepal.Big + Petal.Width + (1 | Species), data = df)
model <- anova(mm)

# simple parameters table
model_parameters(model)

# parameters table including effect sizes
model_parameters(
  model,
  es_type = "eta",
  ci = 0.9,
  df_error = dof_satterthwaite(mm)[2:3]
)

Description

Format Bayesian Exploratory Factor Analysis objects from the BayesFM package.

Usage

## S3 method for class 'befa'
model_parameters(
  model,
  sort = FALSE,
  centrality = "median",
  dispersion = FALSE,
  ci = 0.95,
  ci_method = "eti",
  test = NULL,
  verbose = TRUE,
  ...
)

Arguments

Value

A data frame of loadings.

Examples

library(parameters)

if (require("BayesFM")) {
  efa <- BayesFM::befa(mtcars, iter = 1000)
  results <- model_parameters(efa, sort = TRUE, verbose = FALSE)
  results
  efa_to_cfa(results, verbose = FALSE)
}

Description

Parameters from BFBayesFactor objects from {BayesFactor} package.

Usage

## S3 method for class 'BFBayesFactor'
model_parameters(
  model,
  centrality = "median",
  dispersion = FALSE,
  ci = 0.95,
  ci_method = "eti",
  test = "pd",
  rope_range = "default",
  rope_ci = 0.95,
  priors = TRUE,
  es_type = NULL,
  include_proportions = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Details

The meaning of the extracted parameters:

• For BayesFactor::ttestBF(): Difference is the raw difference between the means.

• For BayesFactor::correlationBF(): rho is the linear correlation estimate (equivalent to Pearson's r).

• For BayesFactor::lmBF() / BayesFactor::generalTestBF() / BayesFactor::regressionBF() / BayesFactor::anovaBF(): in addition to parameters of the fixed and random effects, there are: mu is the (mean-centered) intercept; sig2 is the model's sigma; g / ⁠g_*⁠ are the g parameters; See the Bayes Factors for ANOVAs paper (doi:10.1016/j.jmp.2012.08.001).

Value

A data frame of indices related to the model's parameters.

Examples

# Bayesian t-test
model <- BayesFactor::ttestBF(x = rnorm(100, 1, 1))
model_parameters(model)
model_parameters(model, es_type = "cohens_d", ci = 0.9)

# Bayesian contingency table analysis
data(raceDolls)
bf <- BayesFactor::contingencyTableBF(
  raceDolls,
  sampleType = "indepMulti",
  fixedMargin = "cols"
)
model_parameters(bf,
  centrality = "mean",
  dispersion = TRUE,
  verbose = FALSE,
  es_type = "cramers_v"
)

Description

Extract and compute indices and measures to describe parameters of generalized additive models (GAM(M)s).

Usage

## S3 method for class 'cgam'
model_parameters(
  model,
  ci = 0.95,
  ci_method = "residual",
  bootstrap = FALSE,
  iterations = 1000,
  standardize = NULL,
  exponentiate = FALSE,
  p_adjust = NULL,
  keep = NULL,
  drop = NULL,
  verbose = TRUE,
  ...
)

## S3 method for class 'gamm'
model_parameters(
  model,
  ci = 0.95,
  bootstrap = FALSE,
  iterations = 1000,
  verbose = TRUE,
  ...
)

## S3 method for class 'Gam'
model_parameters(
  model,
  es_type = NULL,
  df_error = NULL,
  type = NULL,
  table_wide = FALSE,
  verbose = TRUE,
  ...
)

## S3 method for class 'scam'
model_parameters(
  model,
  ci = 0.95,
  ci_method = "residual",
  bootstrap = FALSE,
  iterations = 1000,
  standardize = NULL,
  exponentiate = FALSE,
  p_adjust = NULL,
  keep = NULL,
  drop = NULL,
  verbose = TRUE,
  ...
)

Arguments