In easystats/parameters: Processing of Model Parameters

library(knitr)
knitr::opts_chunk$set(
  comment = ">",
  warning = FALSE,
  message = FALSE
)
options(digits = 2)
options(knitr.kable.NA = "")

pkgs <- c("effectsize", "parameters", "correlation")
if (!all(sapply(pkgs, requireNamespace, quietly = TRUE))) {
  knitr::opts_chunk$set(eval = FALSE)
}

set.seed(333)

Introduction

Standardizing parameters (i.e., coefficients) can allow for their comparison within and between models, variables and studies. Moreover, as it returns coefficients expressed in terms of change of variance (for instance, coefficients expressed in terms of SD of the response variable), it can allow for the usage of effect size interpretation guidelines, such as Cohen's (1988) famous rules of thumb.

However, standardizing a model's parameters should not be automatically and mindlessly done: for some research fields, particular variables or types of studies (e.g., replications), it sometimes makes more sense to keep, use and interpret the original parameters, especially if they are well known or easily understood.

Critically, parameters standardization is not a trivial process. Different techniques exist, that can lead to drastically different results. Thus, it is critical that the standardization method is explicitly documented and detailed.

Standardizing Parameters of Simple Models

Standardized Associations

library(parameters)
library(effectsize)

m <- lm(rating ~ complaints, data = attitude)

standardize_parameters(m)

Standardizing the coefficient of this simple linear regression gives a value of

round(standardize_parameters(m)[2, 2], 2)

But did you know that for a simple regression this is actually the same as a correlation? Thus, you can eventually apply some (in)famous interpretation guidelines (e.g., Cohen's rules of thumb).

library(correlation)
correlation(attitude, select = c("rating", "complaints"))

Standardized Differences

How does it work in the case of differences, when factors are entered and differences between a given level and a reference level? You might have heard that it is similar to a Cohen's *d*. Well, let's see.

mtcars <- datasets::mtcars

# Select portion of data containing the two levels of interest
mtcars$am <- factor(mtcars$am, labels = c("Manual", "Automatic"))

m <- lm(mpg ~ am, data = mtcars)
standardize_parameters(m)

This linear model suggests that the standardized difference between Manual (the reference level - the model's intercept) and Automatic is of 1.20 standard deviation of mpg (because the response variable was standardized, right?). Let's compute the Cohen's *d* between these two levels:

library(effectsize)
cohens_d(mpg ~ am, data = mtcars)

It is larger! Why? How? Both differences should be expressed in units of SD! But which SDs? Different SDs!

When looking at the difference between groups as a slope, the standardized parameter is the difference between the means in $SD_{mpg}$. That is, the slope between Manual and Automatic is a change of 1.20 $SD_{mpg}$s.

However, when looking a the difference as a distance between two populations, Cohen's d is the distance between the means in units of pooled SDs. That is, the distance between Manual and Automatic is of 1.48 SDs of each of the groups (here assumed to be equal).

In this simple model, the pooled SD is the residual SD, so we can also estimate Cohen's d as:

coef(m)[2] / sigma(m)

And we can also get an approximation of Cohen's d by converting the $t$-statistic from the regression model via t_to_d():

model_parameters(m)

t_to_d(4.11, df_error = 30)

It is also interesting to note that using the smart method (explained in detail below) when standardizing parameters will give you indices equivalent to Glass' *delta*, which is a standardized difference expressed in terms of SD of the reference group.

m <- lm(mpg ~ am, data = mtcars)

standardize_parameters(m, method = "smart")

glass_delta(mpg ~ am, data = mtcars)

... So note that some standardized differences are different than others! :)

Standardizing Parameters of Linear Models

As mentioned above, standardization of parameters can also be used to compare among parameters within the same model. Essentially, what prevents us from normally being able to compare among different parameters is that their underlying variables are on different scales.[^But also as noted above, this is not always an issue. For example, when the variables scale is important for the interpretation of results, standardization might in fact hinder interpretation!]

For example, in the following example, we use a liner regression model to predict a worker's salary (in Shmekels) from their age (years), seniority (years), overtime (xtra_hours) and how many compliments they give their boss (n_comps).

Let us explore the different parameter standardization methods provided by parameters.

Standardized Slopes are Not (Always) Correlations

We saw that in simple linear models, the standardized slope is equal to the correlation between the outcome and predictor - does this hold for multiple regression as well? As in each effect in a regression model is "adjusted" for the other ones, we might expect coefficients to be somewhat alike to partial correlations. Let's first start by computing the partial correlation between numeric predictors and the outcome.

data("hardlyworking", package = "effectsize")
head(hardlyworking)

correlation(
  hardlyworking,
  select = "salary",
  select2 = c("xtra_hours", "n_comps", "age", "seniority"),
  partial = TRUE # get partial correlations
)

Let's compare these to the standardized slopes:

mod <- lm(salary ~ xtra_hours + n_comps + age + seniority,
  data = hardlyworking
)

standardize_parameters(mod)

They are quite different! It seems then that standardized slopes in multiple linear regressions are not the same a correlations or partial correlations :(

However, not all hope is lost yet - we can still try and recover the partial correlations from our model, in another way: by converting the t-statistics (and their degrees of freedom, df) into a partial correlation coefficient r.

params <- model_parameters(mod)

t_to_r(params$t[-1], df_error = params$df_error[-1])

Wow, the retrieved correlations coefficients from the regression model are exactly the same as the partial correlations we estimated above! So these "r" effect sizes can also be used.

Methods of Standardizing Parameters

Let's convert age into a 3-level factor:

hardlyworking$age_g <- cut(hardlyworking$age,
  breaks = c(25, 30, 35, 45)
)

mod <- lm(salary ~ xtra_hours + n_comps + age_g + seniority,
  data = hardlyworking
)

model_parameters(mod)

It seems like the best or most important predictor is n_comps as it has the coefficient. However, it is hard to compare among predictors, as they are on different scales. To address this issue, we must have all the predictors on the same scale - usually in the arbitrary unit of standard deviations.

"refit": Re-fitting the model with standardized data

This method is based on a complete model re-fit with a standardized version of data. Hence, this method is equal to standardizing the variables before fitting the model. It is the "purest" and the most accurate [@neter1989applied], but it is also the most computationally costly and long (especially for heavy models such as Bayesian models, or complex mixed models). This method is particularly recommended for models that include interactions or transformations (e.g., exponentiation, log, polynomial or spline terms).

standardize_parameters(mod, method = "refit")

standardize_parameters also has a robust argument (default to FALSE), which enables a robust standardization of the data, i.e., based on the median and MAD instead of the mean and SD:

standardize_parameters(mod, method = "refit", robust = TRUE)

Note that since age_g is a factor, it is not numerically standardized, and so it standardized parameter is still not directly comparable to those of numeric variables. To address this, we can set two_sd = TRUE, thereby scaling parameters on 2 SDs (or MADs) of the predictors [@gelman2008scaling].

standardize_parameters(mod, method = "refit", two_sd = TRUE)

parameters also comes with a helper function that returns the re-fit model, without summarizing it, which can then be used as the original model would:

mod_z <- standardize(mod, two_sd = FALSE, robust = FALSE)
mod_z

model_parameters(mod_z)

"posthoc": Refit without refitting

Post-hoc standardization of the parameters aims 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.

But why is interaction for posthoc-standardization problematic? A regression model estimates coefficient between two variables when the other predictors are at 0 (are fixed at 0, that people interpret as "adjusted for"). When a standardized data is passed (in the refit method), the effects and interactions are estimated at the means of the other predictors (because 0 is the mean for a standardized variable). Whereas in posthoc-standardization, this coefficient corresponds to something different (because the 0 has different meanings in standardized and non-standardized data).

standardize_parameters(mod, method = "posthoc")

"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 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.

standardize_parameters(mod, method = "smart")

"basic": Raw scaling of the model frame

This method is similar to method = "posthoc", but treats all variables as continuous: it scales the coefficient by the standard deviation of model's matrix' parameter of factors levels (transformed to integers) or binary predictors. Although it can be argued that this might be inappropriate for these cases, this method allows for easier importance judgment across all predictor type (numeric, factor, interactions...). It is also the type of standardization implemented by default in other software packages (also lm.beta::lm.beta()), and, such as can be used for reproducibility and replication purposes.

standardize_parameters(mod, method = "basic")

Standardizing Parameters In Mixed Models

Linear mixed models (LMM/HLM/MLM) offer an additional conundrum to standardization - how does one even calculate the SDs of the various predictors? Or of the response - is it the deviations within each group? Or perhaps between them?

The solution: standardize according to level of the predictor [@hoffman2015longitudinal, page 342]! Level 1 parameters are standardized according to variance within groups, while level 2 parameters are standardized according to variance between groups. The resulting standardized coefficient are also called pseudo-standardized coefficients.[^Note that like method "basic", these are based on the model matrix.]

m <- lme4::lmer(Reaction ~ Days + (Days | Subject), data = lme4::sleepstudy)

standardize_parameters(m, method = "pseudo", ci_method = "satterthwaite")

# compare to:
standardize_parameters(m, method = "basic", ci_method = "satterthwaite")

Standardizing Parameters In Generalized Linear Models

Unlike linear (/mixed) models, in generalized linear (/mixed) models (GLMs) there is less of a need for standardization. Why? Because in many GLMs the estimated coefficients are themselves measures of effect size, such as odds-ratios (OR) in logistic regression, or incidence rate ratios (IRR) in Poisson regressions. This is because in such model the outcome is not on an arbitrary scale - that is, the meaning of rates and probabilities are changed by arbitrary linear transformations.

But still, some standardization is sometimes needed, for the predictors. Luckily, standardize_parameters() (and standardize()) are smart enough to know when GLMs are passed so as to only standardize according to the predictors:

mod_b <- glm(am ~ mpg + factor(cyl),
  data = mtcars,
  family = binomial()
)

standardize_parameters(mod_b, method = "refit", two_sd = TRUE)
# standardize_parameters(mod_b, method = "posthoc", two_sd = TRUE)
# standardize_parameters(mod_b, method = "basic")

These can then be converted to OR (with exp()) and discussed as the "change in Odds as a function of a change in one SD of x".

std <- standardize_parameters(mod_b, method = "refit", two_sd = TRUE)
exp(std$Std_Coefficient)

Or we can directly ask for the coefficients to be exponentiated:

standardize_parameters(mod_b, method = "refit", two_sd = TRUE, exponentiate = TRUE)

References

easystats/parameters documentation built on April 27, 2024, 7:28 p.m.