Different methods for model parameters standardization

[DominiqueMakowski]

Currently, model standardization does the full standardization, meaning that it actually directly standardizes the outcome and predictors (omitting the factors and binary variables) and then refit the model.

While this is IMHO the most meaningful method to obtain standardized coefs (that can be interpreted, for numeric predictors, as how an increase of 1 SD in the predictor changes the SD of the outcome), it is also slow (necessitating to refit the model), which can be a problem, for instance in Bayesian analyses.

Although a posteriori standardization of parameters does not always make sense (especially for factors and interactions effects), it would be good to allow for different (and faster) standardization methods, such as the simple scaling of parameters based on the outcome's SD.

I am not sure what design would be the best, and several different methods exist.

Some references:

• Bring, J. (1994). How to standardize regression coefficients. The American Statistician, 48(3), 209-213.
• Menard, S. (2004). Six approaches to calculating standardized logistic regression coefficients. The American Statistician, 58(3), 218-223.
• Gelman, A. (2008). Scaling regression inputs by dividing by two standard deviations. Statistics in medicine, 27(15), 2865-2873.
• Schielzeth, H. (2010). Simple means to improve the interpretability of regression coefficients. Methods in Ecology and Evolution, 1(2), 103-113.
• Menard, S. (2011). Standards for standardized logistic regression coefficients. Social Forces, 89(4), 1409-1428.

[strengejacke]

I hardly dare to say it, but there is of course also a function in sjstats:
https://github.com/strengejacke/sjstats/blob/master/R/std_b.R

The way how parameters are scaled is the "default" that is used in many statistical packages, so results are not unexpected. The function also includes the standardization from Gelman (by 2 SD), and for mixed models and stanreg objects.

[strengejacke]

But as you can see, I'm removing a lot of functions in favor of the easystats packages, where they're much better findable.

[DominiqueMakowski]

This time I knew it was in sjstats ☺️ but nevertheless I think this is an important feature for parameters.

Here's a brief overview (to be improved):

Refit method

• Standardize data and refit
• Bring et al. (1994): This standardization is sometimes recommended to improve computational accuracy (see Neter, Wasserman, and Kutner, 1989).

Full

• First, calculate the regression coefficients by using unstandardized variables and then
  multiply them by the ratio between the standard deviation of the respective independent variable and the standard deviation of the dependent variable (Bring et al., 1994)
• Menard et al. 2011: For logistic models, b*M = b(sX)R/slogit(Ŷ) where b*M is the standardized logistic regression coefficient, b is the unstandardized logistic regression coefficient, sX is he standard deviation of the independent variable X, slogit(Ŷ) is the standard deviation of logit(Ŷ), that is, the standard deviation of the predicted values of logit(Y), and R is the correlation between the observed values of Y (either 0 or 1) and the predicted values of Y (predicted probabilities for each case).

Partial

• Menard et al. 2011: For logistic models, Goodman (1972) describes a “standardized” coefficient that is really the normal distribution version of the Wald test statistic for logistic regression, b = b/SE(b), where b is the unstandardized (logit or logistic regression) coefficient and SE(b) is the standard error of that coefficient (Menard 2010)
• Menard et al. 2011: A second approach, suggested for logistic regression models by Menard (1995) and Agresti (1996), is to ignore the variance in logit(Y) and simply standardize the coefficients based on the variance in the predictors, using b = (b)(sX).

2 SD

• Gelman, A. (2008): For these reasons, we recommend the general practice of scaling numeric inputs by dividing by two standard deviations, which allows the coefficients to be interpreted in the same way as with binary inputs. (We leave binary inputs unscaled because their coefficients can already be interpreted directly.)
• Implemented here (and sjstats)

[strengejacke]

We should definitely include such a function, and also offer different ways of standardizing, including model refit.

[pdwaggoner]

For whatever value it may have, I think standardization is a great idea. Most packages do these things (e.g., car, MASS, scales and even the scale function in base R). How are you thinking to go about incorporating these tests. Just trying to get a sense of how you guys develop and expand packages.

[strengejacke]

You mean the technical approach? First, writing a generic

standardize_coefficients <- function(x, ...) {
  UseMethod("standardize_coefficients")
}

and then the related methods for objects, like standardize_coefficients.lm etc. I would start with what we have already done, so copying the code from sjstats. Then add step by step methods for other models / objects.

[DominiqueMakowski]

I gave it a try (on the dev branch).

I also added a small testing suite code to help to compare the new solutions and methods to refit (that can be considered as the baseline):

library(parameters)
library(insight)
library(dplyr)
library(sjstats)
library(MuMIn)
library(lm.beta)

comparison <- function(model){
  out <- standardize_parameters(model, method="refit", robust=FALSE)
  robust_true <- standardize_parameters(model, method="refit", robust=FALSE)$estimate

  out$full <- out$estimate - standardize_parameters(model, method="full", robust=FALSE)$estimate
  out$full_robust <- robust_true - standardize_parameters(model, method="full", robust=FALSE)$estimate


  # out$MuMin <- out$estimate - MuMIn::std.coef(model, partial.sd=FALSE)[, 1]
  out$sjstats <- out$estimate - c(NA, sjstats::std_beta(model)[, 2])
  out$lm.beta <- out$estimate - lm.beta::lm.beta(model)$standardized.coefficients
  out <- dplyr::select(out, -estimate)
  return(out)
}

data <- dplyr::mutate(iris, Group_Sepal.Width = as.factor(ifelse(Sepal.Width > 3, "High", "Low")))

Simple Models with Numerics

model <- lm(Sepal.Length ~ Petal.Width * Sepal.Width, data=data)
knitr::kable(comparison(model), digits=2)

parameter	full	full_robust	sjstats	lm.beta
(Intercept)	-0.01	-0.01	NA	-0.01
Petal.Width	-0.28	-0.28	-0.28	-0.28
Sepal.Width	-0.06	-0.06	-0.06	-0.06
Petal.Width:Sepal.Width	0.01	0.01	0.24	0.24

Simple Models with Factors

model <- lm(Sepal.Length ~ Species * Group_Sepal.Width, data=data)
knitr::kable(comparison(model), digits=2)

parameter	full	full_robust	sjstats	lm.beta
(Intercept)	-0.93	-0.93	NA	-0.93
Speciesversicolor	0.00	0.00	0.90	0.90
Speciesvirginica	0.00	0.00	1.10	1.10
Group_Sepal.WidthLow	0.00	0.00	-0.27	-0.27
Speciesversicolor:Group_Sepal.WidthLow	0.00	0.00	-0.14	-0.14
Speciesvirginica:Group_Sepal.WidthLow	0.00	0.00	0.08	0.08

Model with Factors and Numerics

model <- lm(Sepal.Length ~ Species * Petal.Width, data=data)
knitr::kable(comparison(model), digits=2)

parameter	full	full_robust	sjstats	lm.beta
(Intercept)	0.06	0.06	NA	0.06
Speciesversicolor	0.72	0.72	0.25	0.25
Speciesvirginica	-0.40	-0.40	-0.09	-0.09
Petal.Width	0.00	0.00	0.00	0.00
Speciesversicolor:Petal.Width	0.00	0.00	0.07	0.07
Speciesvirginica:Petal.Width	0.00	0.00	0.07	0.07

Model with Polynomials

model <- lm(Sepal.Length ~ Species * poly(Petal.Width, 2), data=data)
knitr::kable(comparison(model), digits=2)

parameter	full	full_robust	sjstats	lm.beta
(Intercept)	-4.55	-4.55	NA	-4.55
Speciesversicolor	0.00	0.00	1.59	1.59
Speciesvirginica	0.00	0.00	1.87	1.87
poly(Petal.Width, 2)1	0.00	0.00	-39.95	-39.95
poly(Petal.Width, 2)2	0.00	0.00	-17.87	-17.87
Speciesversicolor:poly(Petal.Width, 2)1	0.00	0.00	66.45	66.45
Speciesvirginica:poly(Petal.Width, 2)1	0.00	0.00	64.90	64.90
Speciesversicolor:poly(Petal.Width, 2)2	0.00	0.00	4.46	4.46
Speciesvirginica:poly(Petal.Width, 2)2	0.00	0.00	12.58	12.58

[DominiqueMakowski]

Currently, all "external" packages (sjstats MuMin and lm.beta) give the same results. Some of these results differ a lot from the currently implemented methods. Importantly, the a posteriori standardisation is similar to a refit standardization only in a very few cases (for simple models with numeric predictors).

[DominiqueMakowski]

Here's a clearer version with only the difference from refit (the baseline):

library(parameters)
library(insight)
library(dplyr)
library(sjstats)
library(MuMIn)
library(lm.beta)

data <- dplyr::mutate(iris, Group_Sepal.Width = ifelse(Sepal.Width > 3, "High", "Low"))
model <- lm(Sepal.Length ~ Species * Petal.Width * Sepal.Width + poly(Petal.Length, 2), data=data)

out <- standardize_parameters(model, method="refit", robust=FALSE)
robust_true <- standardize_parameters(model, method="refit", robust=FALSE)$estimate

out$full <- out$estimate - standardize_parameters(model, method="full", robust=FALSE)$estimate
out$full_robust <- robust_true - standardize_parameters(model, method="full", robust=FALSE)$estimate


# out$MuMin <- out$estimate - MuMIn::std.coef(model, partial.sd=FALSE)[, 1]
out$sjstats <- out$estimate - c(NA, sjstats::std_beta(model)[, 2])
out$lm.beta <- out$estimate - lm.beta::lm.beta(model)$standardized.coefficients
out <- dplyr::select(out, -estimate)
knitr::kable(out, digits=2)

parameter	full	full_robust	sjstats	lm.beta
(Intercept)	1.01	1.01	NA	1.01
Speciesversicolor	0.18	0.18	-0.37	-0.37
Speciesvirginica	3.49	3.49	0.94	0.94
Petal.Width	-1.73	-1.73	-1.73	-1.73
Sepal.Width	-0.39	-0.39	-0.39	-0.39
poly(Petal.Length, 2)1	0.00	0.00	16.26	16.26
poly(Petal.Length, 2)2	0.00	0.00	1.88	1.88
Speciesversicolor:Petal.Width	-0.08	-0.08	-0.23	-0.23
Speciesvirginica:Petal.Width	-0.59	-0.59	-0.57	-0.57
Speciesversicolor:Sepal.Width	-0.02	-0.02	-0.53	-0.53
Speciesvirginica:Sepal.Width	-0.13	-0.13	-1.57	-1.57
Petal.Width:Sepal.Width	0.08	0.08	1.49	1.49
Speciesversicolor:Petal.Width:Sepal.Width	0.00	0.00	0.05	0.05
Speciesvirginica:Petal.Width:Sepal.Width	0.03	0.03	0.66	0.66

^{Created on 2019-04-20 by the reprex package (v0.2.1)}

[DominiqueMakowski]

Sorry for this compulsive editing, but it seems like the discrepancies are highly volatile, especially when factors are involed...

[strengejacke]

what is full and full_robust?

[DominiqueMakowski]

full is the "classic" solution (using SD of response and of predictor), but with particular (and exploratory) treatment of interactions and factors. robust scales by mad instead of SD

[DominiqueMakowski]

I think we have a clearer view now of the general design.