R 中的逻辑回归：glm() 与 rxGlm()

Question

我在 R 中拟合了很多 GLM。通常我为此使用 revoScaleR::rxGlm()，因为我处理大型数据集并使用非常复杂的模型公式 - 而 glm() 就是无法应付。

过去，这些都是基于泊松或伽马错误结构和日志链接函数。这一切都很好。

今天我正在尝试构建一个逻辑回归模型，这是我以前在 R 中没有做过的，我偶然发现了一个问题。我正在使用 revoScaleR::rxLogit() 尽管 revoScaleR::rxGlm() 产生相同的输出 - 并且有同样的问题。

考虑这个代表：

df_reprex <- data.frame(x = c(1, 1, 2, 2), # number of trials
                        y = c(0, 1, 0, 1)) # number of successes

df_reprex$p <- df_reprex$y / df_reprex$x # success rate

# overall average success rate is 2/6 = 0.333, so I hope the model outputs will give this number

glm_1 <- glm(p ~ 1,
             family = binomial,
             data = df_reprex,
             weights = x)

exp(glm_1$coefficients[1]) / (1 + exp(glm_1$coefficients[1])) # overall fitted average 0.333 - correct

glm_2 <- rxLogit(p ~ 1,
                 data = df_reprex,
                 pweights = "x")

exp(glm_2$coefficients[1]) / (1 + exp(glm_2$coefficients[1])) # overall fitted average 0.167 - incorrect

第一次调用glm() 会产生正确的答案。对rxLogit() 的第二次调用没有。阅读rxLogit() 的文档：https://docs.microsoft.com/en-us/machine-learning-server/r-reference/revoscaler/rxlogit 它声明“因变量必须是二进制的”。

所以看起来rxLogit() 需要我使用y 作为因变量而不是p。但是，如果我运行

glm_2 <- rxLogit(y ~ 1,
                 data = df_reprex,
                 pweights = "x")

我得到一个总体平均水平

exp(glm_2$coefficients[1]) / (1 + exp(glm_2$coefficients[1]))

取而代之的是 0.5，这也不是正确答案。

有谁知道我该如何解决这个问题？我是否需要在模型公式中使用offset() 术语，或者更改权重，或者...

（通过使用revoScaleR 包，我偶尔会将自己画到这样的角落，因为似乎没有多少其他人使用它）

Answer 1

我在这里瞎了眼，因为我自己无法在 RevoScaleR 中验证这些 - 但是您会尝试运行下面的代码并就结果发表评论吗？然后我可以相应地编辑/删除这篇文章

尝试两件事：

• 扩展数据，去掉权重声明
• 在 rxLogit 或 rxGlm 中使用 cbind(y,x-y)~1，无需权重且无需扩展数据

---

如果因变量需要是二进制的，则必须扩展数据，以便每一行对应于每个 1 或 0 响应，然后在不带 weights 参数的 glm 调用中运行此扩展数据。

我尝试通过将标签应用于df_reprex 然后制作相应的df_reprex_expanded 来通过您的示例来证明这一点——我知道这很不幸，因为您说您正在使用的数据已经很大。

rxLogit 是否允许 cbind 表示，就像 glm() 一样（我举了一个例子为glm1b），因为这将允许数据保持相同的大小......来自rxLogit page，我猜不适用于 rxLogit，但 rxGLM 可能允许它，因为formula page 中有以下注释：

一个公式通常由一个响应组成，在大多数 RevoScaleR 函数可以是单个变量或多个变量组合 使用 cbind、“~”运算符和一个或多个预测变量，通常 由“+”运算符分隔。 rxSummary 函数通常 需要一个没有响应的公式。

下面示例中的glm_2b 或glm_2c 是否有效？

df_reprex <- data.frame(x = c(1, 1, 2, 2), # number of trials
                        y = c(0, 1, 0, 1), # number of successes
                        trial=c("first", "second", "third", "fourth")) # trial label

df_reprex$p <- df_reprex$y / df_reprex$x # success rate

# overall average success rate is 2/6 = 0.333, so I hope the model outputs will give this number

glm_1 <- glm(p ~ 1,
             family = binomial,
             data = df_reprex,
             weights = x)

exp(glm_1$coefficients[1]) / (1 + exp(glm_1$coefficients[1])) # overall fitted average 0.333 - correct


df_reprex_expanded <- data.frame(y=c(0,1,0,0,1,0),
                                trial=c("first","second","third", "third", "fourth", "fourth"))

## binary dependent variable
## expanded data
## no weights
glm_1a <- glm(y ~ 1,
              family = binomial,
              data = df_reprex_expanded)


exp(glm_1a$coefficients[1]) / (1 + exp(glm_1a$coefficients[1])) # overall fitted average 0.333 - correct

## cbind(success, failures) dependent variable
## compressed data
## no weights
glm_1b <- glm(cbind(y,x-y)~1,
              family=binomial,
              data=df_reprex)

exp(glm_1b$coefficients[1]) / (1 + exp(glm_1b$coefficients[1])) # overall fitted average 0.333 - correct


glm_2 <- rxLogit(p ~ 1,
                 data = df_reprex,
                 pweights = "x")

exp(glm_2$coefficients[1]) / (1 + exp(glm_2$coefficients[1])) # overall fitted average 0.167 - incorrect

glm_2a <- rxLogit(y ~ 1,
                 data = df_reprex_expanded)

exp(glm_2a$coefficients[1]) / (1 + exp(glm_2a$coefficients[1])) # overall fitted average ???

# try cbind() in rxLogit.  If no, then try rxGlm below
glm_2b <- rxLogit(cbind(y,x-y)~1,
              data=df_reprex)

exp(glm_2b$coefficients[1]) / (1 + exp(glm_2b$coefficients[1])) # overall fitted average ???

# cbind() + rxGlm + family=binomial FTW(?)
glm_2c <- rxGlm(cbind(y,x-y)~1,
              family=binomial,
              data=df_reprex)

exp(glm_2c$coefficients[1]) / (1 + exp(glm_2c$coefficients[1])) # overall fitted average ???