拟合具有多个依赖变量/LHS的逻辑模型

Question

有没有办法使用glm()进行多元逻辑回归？我有几个二进制结果，我知道你可以用线性回归（lm（））和cbind（）来做到这一点，但我似乎不知道如何用glm（）和cbind（）来做到这一点

library(tidyverse)

lm(cbind(mpg, cyl, disp) ~ hp + drat + wt + qsec, data = mtcars) %>% 
  broom::tidy(conf.int = T)

这会返回一个包含结果的整洁小标题（mpg、cyl、disp）。我如何将其扩展到逻辑回归。

df <- mtcars %>% 
  mutate(across(c(vs, am), factor))

glm(cbind(am, vs) ~ hp + drat + wt + qsec, family = binomial(link = "logit"), data = df) %>% 
  broom::tidy()

我知道这不是最好的数据集（可能是由于完全分离），但与单独运行 glms 相比，cbind 只是返回意外的输出。

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

Answer 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）结果。正如 Ben Bolker 指出的那样，此类因变量需要进行多项回归。

正如您所提到的，如果您一次运行具有一个因变量的单独

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