看一下
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