Obtaining Standardized coefficients from “rstanarm” package in R?

I was wondering if it might be possible (and perhaps recommended) to obtain standardized coefficients from stan_glm() in the rstanarm package?

Example:

library("rstanarm")
fit <- stan_glm(wt ~ vs*gear, data = mtcars)

What’s a standardized coefficient? I know about standardizing outcomes and predictors, but not coefficients.

There’s not a function in rstanarm for that but it’s easy to do with the output from rstanarm. For a coefficient b I think the “standardized coefficient” traditionally refers to

b * sd(x)/sd(y)

If that’s what you’re referring to then if you want the posterior distribution of the standardized coefficient then you just need to multiply each posterior draw of b by that ratio of standard deviations. For the example model in the original post that would be something like this:

fit <- stan_glm(wt ~ vs*gear, data = mtcars)
X <- model.matrix(fit)[, -1]                 # -1 to drop column of 1s for intercept
sd_X_over_sd_y <- apply(X, 2, sd) / sd(mtcars$wt)
b <- as.matrix(fit, pars = colnames(X))      # posterior distribution of regression coefficients
b_std <- sweep(b, 2, sd_X_over_sd_y, "*")    # multiply each row of b by sd_X_over_sd_y

Comparing the posterior means of b and b_std:

cbind(b = colMeans(b), b_std = colMeans(b_std))

                  b       b_std
vs      -0.69804364 -0.05837536
gear    -0.64534477 -0.07900151
vs:gear -0.04621978 -0.01514321

Hope that helps.

3 Likes

Standardized regression coefficients are reported by default in SPSS regression output, this explains their popularity. They are used (a little bit naively) to compare effects of different predictors on some outcome variable (especially in cases when values of variables used don’t have any natural interpretation).
Perhaps, even simpler way to obtain them is to standardize both the outcome and predictors before running analyses with rstanarm. The drawback of this approach, is that in case you would like to report both raw and standardized coefficients, you would need to fit the model twice.

BTW. There is a small error in code provided by jonah, which may confuse some readers. I believe one of the lines in the example should be:

sd_X_over_sd_y <- apply(X, 2, sd) / sd(mtcars$wt)

Thanks for catching that typo. You’re right it should use the wt variable not the mpg variable.