有没有方法使用glm()进行多元逻辑回归?我有几个二元结果,我知道你可以用线性回归(lm())和cbind()来实现这一点,但我似乎不知道如何用glm()和cbind()来实现这一点

library(tidyverse)

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

这将返回一个整洁的小方块,包含您的结果(每加仑、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()

我知道这不是最好使用的数据集(可能是由于完全分离),但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()

推荐答案

看一下

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,如here所示,这将产生以下方法:

#--------------
# 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

R相关问答推荐

在Julia中调用R函数

ggplot2中的X轴显示数值,单位为百,而不是十

然后根据不同的列值有条件地执行函数

如果第一个列表中的元素等于第二个列表的元素,则替换为第三个列表的元素

标识R中多个列中缺少的唯一值

未识别时区

如何调整曲线图中的y轴标签?

如何从R ggplot图片中获取SVG字符串?

如何在R forestplot中为多条垂直线分配唯一的 colored颜色 ?

从圆到R中的多边形的标绘雷达图

更新R中的数据表(使用data.table)

使用R将简单的JSON解析为嵌套框架

如何显示准确的p值而不是<;0.001*?

如何在Quarto中使用美人鱼图表中的标记来加粗文本

Broom.Mixed::Augment不适用于Sample::分析

GOGPLATE geom_boxploy色彩疯狂

隐藏基于 case 总数的值

网络抓取NBA.com

如何为包创建自定义roxygen2标签?

子样本间系数检验的比较