Как получить стандартизированные коэффициенты регрессии для загруженных моделей или моделей с устойчивой SE в R

Question

В R, когда мы подбираем регрессионную модель, стандартизированные коэффициенты, к сожалению, недоступны. К счастью, есть пакет lm.beta, который их извлекает. Например:

fit.model1.lm <- lm(Sepal.Length ~ Sepal.Width + Petal.Width, iris)
lm.beta(fit.model1.lm)

Однако, если мы используем начальную загрузку (например, из пакета simpleboot):

fit.model1.boot <- lm.boot(fit.model1.lm)

Или если мы используем устойчивые стандартные ошибки (например, из sandwich package):

coeftest(fit.model1.lm, vcov = vcovHC(fit.model1.lm))

Невозможно получить стандартизированные коэффициенты через lm.beta, поскольку они работают только с объектами класса lm.

Как их получить в таком случае?

Answer 1

Чтобы получить стандартизированные коэффициенты, мы также можем scale наши переменные до оценки:

library("lm.beta")
fit.model1.lm <- lm(Sepal.Length ~ Sepal.Width + Petal.Width, iris)
lm.beta(fit.model1.lm)
# Call:
#   lm(formula = Sepal.Length ~ Sepal.Width + Petal.Width, data = iris)
# 
# Standardized Coefficients::
# (Intercept) Sepal.Width Petal.Width 
#   0.0000000   0.2100575   0.8948486 
(fit.model1.lm.st <- lm(Sepal.Length ~ Sepal.Width + Petal.Width, 
                        cbind(lapply(iris[-5], scale), iris[5])))
# Call:
#   lm(formula = Sepal.Length ~ Sepal.Width + Petal.Width, data = cbind(lapply(iris[-5], 
#                                                                              scale), iris[5]))
# 
# Coefficients:
# (Intercept)  Sepal.Width  Petal.Width  
#  -1.422e-15    2.101e-01    8.948e-01  

Так что использование scale d lm fit (здесь fit.model1.lm.st) в lm.boot даст Вы нужных стандартизированных bootstrap коэффициентов.

library("simpleboot")
set.seed(42)  # for sake of reproducibility
(fit.model1.boot <- lm.boot(fit.model1.lm.st, R=200))
# BOOTSTRAP OF LINEAR MODEL  (method = rows)
# 
# Original Model Fit
# ------------------
# Call:
# lm(formula = Sepal.Length ~ Sepal.Width + Petal.Width, data = cbind(lapply(iris[-5], 
#                                                                              scale), iris[5]))
# 
# Coefficients:
# (Intercept)  Sepal.Width  Petal.Width  
#  -1.422e-15    2.101e-01    8.948e-01  

library(lmtest)
library(sandwich)
coeftest(fit.model1.lm.st, vcov=vcovHC(fit.model1.lm.st))
# t test of coefficients:
#   
#                Estimate  Std. Error t value  Pr(>|t|)    
# (Intercept) -1.4222e-15  4.4920e-02  0.0000         1    
# Sepal.Width  2.1006e-01  4.7055e-02  4.4641 1.591e-05 ***
# Petal.Width  8.9485e-01  5.2321e-02 17.1031 < 2.2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1