Logisti c регрессия в R: glm () против rxGlm () - PullRequest
Новая
       40

Logisti c регрессия в R: glm () против rxGlm ()

4 голосов
/ 15 апреля 2020

Я подхожу ко многим GLM в R. Обычно я использовал revoScaleR::rxGlm() для этого, потому что я работаю с большими наборами данных и использую довольно сложные формулы модели - а glm() просто не справится.

В прошлом все они основывались на пуассоновских или гамма-структурах ошибок и функциях логарифмической связи. Все это хорошо работает.

Сегодня я пытаюсь построить регрессионную модель логистики c, чего раньше не делал в 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, я иногда рисую себя в таком углу, потому что не многие другие используют его)

1 Ответ

0 голосов
/ 24 апреля 2020

Я тут летаю вслепую, потому что сам не могу проверить это в RevoScaleR - но не могли бы вы попробовать запустить приведенный ниже код и оставить комментарий о результатах? Затем я могу отредактировать / удалить эту запись соответственно

Две вещи, которые нужно попробовать:

  • Развернуть данные, избавиться от утверждения о весах
  • использовать cbind (y, xy) ~ 1 в rxLo git или rxGlm без весов и без расширения данных

Если зависимая переменная должна быть двоичной, то данные необходимо развернуть так, чтобы каждая строка соответствовала каждому ответу 1 или 0, а затем эти расширенные данные выполняются в вызове glm без аргумента весов.

Я попытался продемонстрировать это на вашем примере, применив метки к df_reprex и затем сделав соответствующий df_reprex_expanded - я знаю, что это неудачно, потому что вы сказали, что данные, с которыми вы работали, уже были большими.

Разрешает ли rxLogit представление cbind, как это делает glm () (я поставил пример как glm1b), потому что это позволило бы данным оставаться одинакового размера ... от rxLo git page , я предполагаю, что не для rxLo git, но 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 ???
...