5

我在 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包我偶尔把自己画成这样的角落，因为似乎没有多少其他人使用它）

4

1 回答 1

0

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

有两件事要尝试：

  • 扩展数据，去掉权重声明
  • 在 rxLogit 或 rxGlm 中使用 cbind(y,xy)~1 不带权重且不扩展数据

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

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

是否rxLogit允许cbind表示，如 glm() 确实（我举了一个例子glm1b），因为这将允许数据保持相同的大小......从rxLogit 页面，我猜不是 rxLogit，但 rxGLM 可能允许它，给定以下公式页面中的注意事项：

公式通常由响应组成，在大多数 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 ???
于 2020-04-24T18:36:19.160 回答