Получение стандартизированных коэффициентов из пакета "rstanarm" в R?
Мне было интересно, возможно ли (и, возможно, рекомендуется) получить стандартизированные коэффициенты из stan_glm()
в rstanarm
пакет? (не нашел ничего конкретного в документации)
Могу ли я просто стандартизировать все переменные, как в обычной регрессии? (см. ниже)
Пример:
library("rstanarm")
fit <- stan_glm(wt ~ vs*gear, data = mtcars)
Стандартизация:
design <- wt ~ vs*gear
vars <- all.vars(design)
stand.vars <- lapply(mtcars[, vars], scale)
fit <- stan_glm(stand.vars, data = mtcars)
1 ответ
Я бы не сказал, что это положительно рекомендуется, но я бы порекомендовал вам не вычитать среднее значение выборки и делить на стандартное отклонение выборки результата, потому что неопределенность оценки в этих двух статистических показателях не будет распространена на апостериорное распределение.
Стандартизация предсказателей является более спорным. Вы можете сделать это, но это затрудняет выполнение апостериорного прогнозирования с новыми данными, потому что вы должны помнить, что нужно вычесть старые средства из новых данных и поделить на старые стандартные отклонения.
Наиболее эффективный в вычислительном отношении подход - это оставить переменные такими, какие они есть, но указать аргумент не по умолчанию QR = TRUE
особенно если вы все равно не собираетесь изменять стандартные (нормальные) априоры для коэффициентов. Затем вы можете стандартизировать апостериорные коэффициенты после факта, если интерес представляют стандартизованные коэффициенты. Для этого вы можете сделать
X <- model.matrix(fit)
sd_X <- apply(X, MARGIN = 2, FUN = sd)[-1]
sd_Y <- apply(posterior_predict(fit), MARGIN = 1, FUN = sd)
beta <- as.matrix(fit)[ , 2:ncol(X), drop = FALSE]
b <- sweep(sweep(beta, MARGIN = 2, STATS = sd_X, FUN = `*`),
MARGIN = 1, STATS = sd_Y, FUN = `/`)
summary(b)
Однако стандартизация коэффициентов регрессии просто создает иллюзию сопоставимости между переменными и ничего не говорит о том, насколько уместна разница в одно стандартное отклонение, особенно для фиктивных переменных. Если ваш вопрос действительно заключается в том, будет ли манипулирование этим предиктором или этим предиктором иметь большее значение для переменной результата, то просто смоделируйте такие манипуляции, как
PPD_0 <- posterior_predict(fit)
nd <- model.frame(fit)
nd[ , 2] <- nd[ , 2] + 1 # for example
PPD_1 <- posterior_predict(fit, newdata = nd)
summary(c(PPD_1 - PPD_0))
и повторите этот процесс для других интересующих манипуляций.