Есть ли способ выполнить многомерную логистическую регрессию с помощью glm()? У меня есть несколько двоичных результатов, и я знаю, что вы можете сделать это с помощью линейной регрессии (lm()) и cbind(), но я не могу понять, как это сделать с помощью glm() и cbind().
library(tidyverse)
lm(cbind(mpg, cyl, disp) ~ hp + drat + wt + qsec, data = mtcars) %>%
broom::tidy(conf.int = T)
Это вернет аккуратную таблицу с вашими результатами (миль на галлон, баллон, дисп). Как я могу распространить это на логистическую регрессию.
df <- mtcars %>%
mutate(across(c(vs, am), factor))
glm(cbind(am, vs) ~ hp + drat + wt + qsec, family = binomial(link = "logit"), data = df) %>%
broom::tidy()
Я знаю, что это не лучший набор данных для использования (вероятно, из-за полного разделения), но cbind просто возвращает неожиданный результат, чем если бы вы запускали glms отдельно.
glm(am ~ hp + drat + wt + qsec, family = binomial(link = "logit"), data = df) %>%
broom::tidy()
glm(vs ~ hp + drat + wt + qsec, family = binomial(link = "logit"), data = df) %>%
broom::tidy()
@BenBolker, похоже, OP не интересуется MNL, а просто имеет одни и те же предикторы, но разные переменные ответа. т.е. 2 разных glms, но с одинаковыми предикторами
Уточню: MNL меня не интересует. Как правильно отмечают Оньямбу и @GRowInG, я заинтересован в создании предсказателейbinary_outcome_1~, предсказателейbinary_outcome_2~, предсказателейbinary_outcome_3~ и т. д.
Посмотри на
glm(cbind(am, vs) ~ hp + drat + wt + qsec, family = binomial(link = "logit"), data = df, method = "model.frame") %>%
broom::tidy()
который возвращает
# A tibble: 5 × 13
column n mean sd median trimmed mad min max range skew kurtosis se
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 cbind(am, vs) 64 1.42 0.498 1 1.40 0 1 2 1 0.316 1.10 0.0622
2 hp 32 147. 68.6 123 141. 52 52 335 283 0.761 3.05 12.1
3 drat 32 3.60 0.535 3.70 3.58 0.475 2.76 4.93 2.17 0.279 2.44 0.0945
4 wt 32 3.22 0.978 3.32 3.15 0.517 1.51 5.42 3.91 0.444 3.17 0.173
5 qsec 32 17.8 1.79 17.7 17.8 0.955 14.5 22.9 8.4 0.387 3.55 0.316
В первой строке результатов в приведенном выше выводе вы можете видеть, что glm
, похоже, интерпретирует cbind(am, vs)
как объединенную зависимую переменную am и vs, имеющую несколько категорий и не являющуюся дихотомической (в этом примере 0 или 1). Как отметил Бен Болкер, для такого рода зависимой переменной потребуется полиномиальная логистическая регрессия.
Как вы упомянули, если вы одновременно запускаете отдельные модели glm
с одной зависимой переменной, вы получаете разные результаты и предупреждающие сообщения относительно полного разделения, а также подобранных вероятностей, имеющих значения (очень близкие к) 0 или 1.
Если вы хотите запустить две модели glm
за один раз и представить результаты в одном объекте, мы можем использовать purrr:map
, как показано здесь, что приводит к следующему подходу:
#--------------
# Load packages
#--------------
library(tidyverse)
#--------------
# Define data, including independent variables & dependent variables
#--------------
df <- mtcars
independent.variables.formula <- "~ hp + drat + wt + qsec"
dependent.variables <- c("am", "vs")
#--------------
# create output for both models in a data.frame showing the results
#--------------
df_res <- map(dependent.variables, function(DV) {
paste(DV, independent.variables.formula) %>%
as.formula %>%
glm(family=binomial(link = "logit"), data = df) %>%
broom::tidy(conf.int = T)
}) %>%
bind_rows() %>%
as.data.frame () %>%
mutate (DV = c(rep ("am", 5), rep ("vs", 5)),
across(c(2:4, 6:7), .fns = function(x) {format(round(x, 2), nsmall = 2)})) %>%
relocate (DV, .before = term)
df_res[ ,6] <- format.pval(df_res[ ,6], eps = .001, digits = 3) # format small p-values < 0.001 nicely
df_res
DV term estimate std.error statistic p.value conf.low conf.high
1 am (Intercept) 206.22 1166788.76 0.00 1.000 NA 238011.23
2 am hp 0.24 687.34 0.00 1.000 -139.85 NA
3 am drat 103.37 132242.73 0.00 0.999 -10359.33 12068.44
4 am wt -94.56 84508.03 0.00 0.999 -54340.19 14115.43
5 am qsec -20.03 34601.06 0.00 1.000 -3473.64 3149.50
6 vs (Intercept) -1928.91 1845493.29 0.00 0.999 -124797.75 120939.92
7 vs hp 1.59 2408.34 0.00 0.999 -489.26 NA
8 vs drat 92.35 139959.57 0.00 0.999 -16047.24 16231.94
9 vs wt -82.30 77899.01 0.00 0.999 -5580.25 4963.40
10 vs qsec 91.59 75674.95 0.00 0.999 -3734.28 4115.71
Я думаю, что ваш первый пример делает совсем не то, чего я ожидал от ОП, и, возможно, не то, что вы думаете. Почему вы включили method = "model.frame"
? (чтобы быть более ясным, первый вызов glm()
создает каркас модели, а не подгоняет GLM...)
@BenBolker Первая часть моего ответа должна объяснить, почему результаты отличаются от того, что ожидает найти ОП. Я тоже не думаю, что речь идет о MNL, но, надеюсь, выясним, верно ли это предположение.
Большое спасибо! Ваше решение сработало, хотя жаль, что cbind() выполняет всю тяжелую работу в lm(), но не в glm(). Надеюсь, все в порядке, я упростил ваш код! df_res <- карта(dependent.variables, function(DV) {paste(DV,dependent.variables.formula) %>% as.formula %>% glm(family=binomial(link = "logit"), data = df) %>% broom::tidy(conf.int = T) %>% mutate(dv = DV,signif = gtools::stars.pval(p.value),through(where(is.numeric), round, 5) )%>% relocate(dv, .before = term) }) %>%bind_rows()
Это тоже работает и помогает мне подумать о том, как это отобразить. df_res2 <- карта(зависимые.переменные, функция(DV) { glm(формула = as.formula(вставить(DV,dependent.variables.formula)), семья=биномиальная(ссылка = "logit"), data = df) % >% broom::tidy(conf.int = T) %>% mutate(dv = DV,signif = gtools::stars.pval(p.value),through(where(is.numeric), round, 5) ) %>% relocate(dv, .before = term) }) %>%bind_rows()
ищите полиномиальную аппроксимацию, например.
?nnet::multinom