#r #lme4 #coefficients #rstanarm #glmmtmb
#r #lme4 #коэффициенты #rstanarm #glmmtmb
Вопрос:
Мне интересно, есть ли простой способ изменить, какие значения находятся в перехвате, возможно, математически, без повторного запуска больших моделей. В качестве примера:
mtcars$cyl<-as.factor(mtcars$cyl)
summary(
lm(mpg~cyl hp,data=mtcars)
)
Вывод:
Call:
lm(formula = mpg ~ cyl hp, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-4.818 -1.959 0.080 1.627 6.812
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.65012 1.58779 18.044 < 2e-16 ***
cyl6 -5.96766 1.63928 -3.640 0.00109 **
cyl8 -8.52085 2.32607 -3.663 0.00103 **
hp -0.02404 0.01541 -1.560 0.12995
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.146 on 28 degrees of freedom
Multiple R-squared: 0.7539, Adjusted R-squared: 0.7275
F-statistic: 28.59 on 3 and 28 DF, p-value: 1.14e-08
Теперь я могу изменить опорный уровень на 6 cyl и увидеть, как 8 cyl теперь сравнивается с 6 cyl, а не с 4 cyl:
mtcars$cyl<-relevel(mtcars$cyl,"6")
summary(
lm(mpg~cyl hp,data=mtcars)
)
Вывод:
Call:
lm(formula = mpg ~ cyl hp, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-4.818 -1.959 0.080 1.627 6.812
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.68246 2.22805 10.18 6.48e-11 ***
cyl4 5.96766 1.63928 3.64 0.00109 **
cyl8 -2.55320 1.97867 -1.29 0.20748
hp -0.02404 0.01541 -1.56 0.12995
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.146 on 28 degrees of freedom
Multiple R-squared: 0.7539, Adjusted R-squared: 0.7275
F-statistic: 28.59 on 3 and 28 DF, p-value: 1.14e-08
Что мне интересно, есть ли способ получить эти значения без повторного запуска модели? Вы можете видеть, что сравнение от 4 cyl до 6 cyl одинаково в каждой модели ( -5.96 и 5.96 ), но как мне получить оценку для «другого» коэффициента в любой модели (например -2.55 , из первой модели). Конечно, в этом случае для запуска другой модели требуется доля секунды. Но с очень большими моделями было бы удобно иметь возможность изменять опорный уровень без повторного запуска. Существуют ли относительно простые способы преобразовать все оценки и стандартные ошибки, основанные на другом опорном уровне, или это слишком сложно сделать?
Любые решения для lme4 , glmmTMB , или rstanarm моделей будут оценены.
Комментарии:
1. проверьте
emmeansпакет.2. Существует также
multcomp::glht, который может позволить вам создавать запросы вручную.3. потенциально полезно: stats.stackexchange.com/questions/168650 /…
Ответ №1:
Вот функция, которая выдаст вам коэффициенты для каждой перестановки заданной факторной переменной без необходимости повторного запуска модели или указания контрастов:
rearrange_model_factors <- function(model, var)
{
var <- deparse(substitute(var))
coefs <- coef(model)
level_coefs <- grep(paste0("^", var), names(coefs))
coefs[level_coefs] <- coefs[1] coefs[level_coefs]
used_levels <- gsub(var, "", names(coefs[level_coefs]))
all_levels <- levels(model$model[[var]])
names(coefs)[1] <- paste0(var, setdiff(all_levels, used_levels))
level_coefs <- grep(paste0("^", var), names(coefs))
levs <- coefs[level_coefs]
perms <- gtools::permutations(length(levs), length(levs))
perms <- lapply(seq(nrow(perms)), function(i) levs[perms[i,]])
lapply(perms, function(x) {
x[-1] <- x[-1] - x[1]
coefs[level_coefs] <- x
names(coefs)[level_coefs] <- names(x)
names(coefs)[1] <- "(Intercept)"
coefs
})
}
Предположим, у вас была такая модель:
iris_mod <- lm(Sepal.Width ~ Species Sepal.Length, data = iris)
Чтобы увидеть, как изменились бы ваши коэффициенты, если Species бы они были в другом порядке, вы бы просто сделали:
rearrange_model_factors(iris_mod, Species)
#> [[1]]
#> (Intercept) Speciesversicolor Speciesvirginica Sepal.Length
#> 1.6765001 -0.9833885 -1.0075104 0.3498801
#>
#> [[2]]
#> (Intercept) Speciesvirginica Speciesversicolor Sepal.Length
#> 1.6765001 -1.0075104 -0.9833885 0.3498801
#>
#> [[3]]
#> (Intercept) Speciessetosa Speciesvirginica Sepal.Length
#> 0.69311160 0.98338851 -0.02412184 0.34988012
#>
#> [[4]]
#> (Intercept) Speciesvirginica Speciessetosa Sepal.Length
#> 0.69311160 -0.02412184 0.98338851 0.34988012
#>
#> [[5]]
#> (Intercept) Speciessetosa Speciesversicolor Sepal.Length
#> 0.66898976 1.00751035 0.02412184 0.34988012
#>
#> [[6]]
#> (Intercept) Speciesversicolor Speciessetosa Sepal.Length
#> 0.66898976 0.02412184 1.00751035 0.34988012
Или на вашем собственном примере:
mtcars$cyl <- as.factor(mtcars$cyl)
rearrange_model_factors(lm(mpg ~ cyl hp, data = mtcars), cyl)
#> [[1]]
#> (Intercept) cyl6 cyl8 hp
#> 28.65011816 -5.96765508 -8.52085075 -0.02403883
#>
#> [[2]]
#> (Intercept) cyl8 cyl6 hp
#> 28.65011816 -8.52085075 -5.96765508 -0.02403883
#>
#> [[3]]
#> (Intercept) cyl4 cyl8 hp
#> 22.68246309 5.96765508 -2.55319567 -0.02403883
#>
#> [[4]]
#> (Intercept) cyl8 cyl4 hp
#> 22.68246309 -2.55319567 5.96765508 -0.02403883
#>
#> [[5]]
#> (Intercept) cyl4 cyl6 hp
#> 20.12926741 8.52085075 2.55319567 -0.02403883
#>
#> [[6]]
#> (Intercept) cyl6 cyl4 hp
#> 20.12926741 2.55319567 8.52085075 -0.02403883
Нам нужно немного изложения, чтобы понять, почему это работает.
Хотя приведенная выше функция запускает модель только один раз, давайте начнем с создания списка, содержащего 3 версии mtcars , в которых уровни базовых коэффициентов cyl различны.
df_list <- list(mtcars_4 = within(mtcars, cyl <- factor(cyl, c(4, 6, 8))),
mtcars_6 = within(mtcars, cyl <- factor(cyl, c(6, 4, 8))),
mtcars_8 = within(mtcars, cyl <- factor(cyl, c(8, 4, 6))))
Теперь мы можем извлечь коэффициенты вашей модели для всех трех версий сразу с помощью lapply . Для ясности мы удалим hp коэффициент, который все равно остается неизменным во всех трех версиях:
coefs <- lapply(df_list, function(x) coef(lm(mpg ~ cyl hp, data = x))[-4])
coefs
#> $mtcars_4
#> (Intercept) cyl6 cyl8
#> 28.650118 -5.967655 -8.520851
#>
#> $mtcars_6
#> (Intercept) cyl4 cyl8
#> 22.682463 5.967655 -2.553196
#>
#> $mtcars_8
#> (Intercept) cyl4 cyl6
#> 20.129267 8.520851 2.553196
Теперь мы напоминаем себе, что коэффициент для каждого уровня фактора задается относительно базового уровня. Это означает, что для коэффициентов без перехвата мы можем просто добавить значение перехвата к их коэффициентам, чтобы получить их абсолютное значение. Это означает, что эти числа представляют ожидаемое значение mpg , когда hp равно 0 для всех трех уровней cyl
coefs <- lapply(coefs, function(x) c(x[1], x[-1] x[1]))
coefs
#> $mtcars_4
#> (Intercept) cyl6 cyl8
#> 28.65012 22.68246 20.12927
#>
#> $mtcars_6
#> (Intercept) cyl4 cyl8
#> 22.68246 28.65012 20.12927
#>
#> $mtcars_8
#> (Intercept) cyl4 cyl6
#> 20.12927 28.65012 22.68246
Поскольку теперь у нас есть все три значения как абсолютные, давайте переименуем «Intercept» в соответствующий уровень фактора:
coefs <- mapply(function(x, y) { names(x)[1] <- y; x},
x = coefs, y = c("cyl4", "cyl6", "cyl8"), SIMPLIFY = FALSE)
coefs
#> $mtcars_4
#> cyl4 cyl6 cyl8
#> 28.65012 22.68246 20.12927
#>
#> $mtcars_6
#> cyl6 cyl4 cyl8
#> 22.68246 28.65012 20.12927
#>
#> $mtcars_8
#> cyl8 cyl4 cyl6
#> 20.12927 28.65012 22.68246
Наконец, давайте изменим порядок, чтобы мы могли сравнить абсолютные значения всех трех уровней факторов:
coefs <- lapply(coefs, function(x) x[order(names(x))])
coefs
#> $mtcars_4
#> cyl4 cyl6 cyl8
#> 28.65012 22.68246 20.12927
#>
#> $mtcars_6
#> cyl4 cyl6 cyl8
#> 28.65012 22.68246 20.12927
#>
#> $mtcars_8
#> cyl4 cyl6 cyl8
#> 28.65012 22.68246 20.12927
Мы видим, что все они одинаковы. Вот почему порядок факторов является произвольным в lm : изменение порядка уровней факторов в конечном итоге дает те же числовые прогнозы, даже если сводка выглядит иначе.
TL; DR
Итак, ответ на ваш вопрос о том, где вы получаете -2.55, если у вас есть только первая модель, — это найти разницу между коэффициентами без перехвата. В этом случае
(-8.520851) -(-5.967655)
#> [1] -2.553196
В качестве альтернативы, добавьте перехват к коэффициентам без перехвата, и вы сможете увидеть, каким было бы перехват, если бы какой-либо из уровней был базовым, и вы можете получить коэффициент для любого уровня относительно любого другого путем простого вычитания. Вот как rearrange_model_factors работает моя функция.
Создано 2020-10-05 пакетом reprex (версия 0.3.0)
Комментарии:
1. Я не думаю, что это делает то, что хочет OP: они хотят иметь возможность извлекать коэффициенты / ошибки std и т.д. для модели с повторными параметрами без переоборудования модели . Я думаю, что предварительная установка всех возможных уровней не удовлетворяет этому запросу.
2. @BenBolker это не то, что я предлагаю. Я просто показываю эквивалентность различных упорядочений факторов. ОП спрашивает: » Но как мне получить оценку для «другого» коэффициента в любой модели (например, -2,55 из первой модели)». — и я показываю, как
3. @BenBolker я имею в виду, как получить коэффициенты по-другому упорядоченного фактора без повторного запуска модели. Извините, мое изложение не было более ясным.
4. вы показываете, как это сделать вручную. Я думаю, что OP ищет оборудование…
5. @Dylan_Gomes нет, SE и CI сложнее. Однако теперь вы знаете, что вам нужно только добавить перехват к коэффициентам фактора, чтобы получить абсолютные значения на каждом уровне, и SE останется прежним — вы можете рассчитать CI для каждого из них, используя этот факт.