Подбор логистической модели с несколькими зависимыми переменными/LHS

Есть ли способ выполнить многомерную логистическую регрессию с помощью 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()

ищите полиномиальную аппроксимацию, например. ?nnet::multinom

Ben Bolker 19.04.2024 23:05

@BenBolker, похоже, OP не интересуется MNL, а просто имеет одни и те же предикторы, но разные переменные ответа. т.е. 2 разных glms, но с одинаковыми предикторами

Onyambu 20.04.2024 00:49

Уточню: MNL меня не интересует. Как правильно отмечают Оньямбу и @GRowInG, я заинтересован в создании предсказателейbinary_outcome_1~, предсказателейbinary_outcome_2~, предсказателейbinary_outcome_3~ и т. д.

allen.joseph 20.04.2024 09:12
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
2
3
65
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Посмотри на

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...)

Ben Bolker 20.04.2024 01:46

@BenBolker Первая часть моего ответа должна объяснить, почему результаты отличаются от того, что ожидает найти ОП. Я тоже не думаю, что речь идет о MNL, но, надеюсь, выясним, верно ли это предположение.

GRowInG 20.04.2024 08:01

Большое спасибо! Ваше решение сработало, хотя жаль, что 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()

allen.joseph 20.04.2024 09:22

Это тоже работает и помогает мне подумать о том, как это отобразить. 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()

allen.joseph 20.04.2024 09:26

Другие вопросы по теме