r - R : glm() vs rxGlm() 中的逻辑回归

Question

我在 R 中安装了很多 GLM。通常我使用 revoScaleR::rxGlm()因为我使用大型数据集并使用相当复杂的模型公式 - 和 glm()就是应付不过来。

在过去，这些都是基于泊松或 Gamma 错误结构和日志链接功能。这一切都很好。

今天我正在尝试建立一个逻辑回归模型，这是我以前在 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 响应，然后此扩展数据在不带权重参数的 glm 调用中运行。

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

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

A formula typically consists of a response, which in most RevoScaleR functions can be a single variable or multiple variables combined using cbind, the "~" operator, and one or more predictors,typically separated by the "+" operator. The rxSummary function typically requires a formula with no response.

是否 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 ???