Удаление полностью разделенных наблюдений из glm ()

Question

Я немного анализирую данные, используя данные HMDA из пакета AER; однако переменные, которые я использовал для подгонки к модели, по-видимому, содержат некоторые наблюдения, которые идеально определяют результаты, проблема, известная как «разделение». Поэтому я попытался исправить это, используя решение, рекомендованное этим потоком , но когда я попытался выполнить первый набор исходного кода из glm.fit(), R вернул сообщение об ошибке:

Error in family$family : object of type 'closure' is not subsettable

, поэтому я не мог продолжить, чтобы удалить эти полностью определенные наблюдения из моих данных с помощью этого кода. Мне интересно, может ли кто-нибудь помочь мне исправить это?

Мой текущий код указан ниже для справки.

# load the AER package and HMDA data
library(AER)
data(HMDA)

# fit a 2-degree olynomial probit model 
probit.fit <- glm(deny ~ poly(hirat, 2), family = binomial, data = HMDA)

# using the revised source code from that stackexchage thread to find out observations that received a warning message
library(tidyverse)
library(dplyr)
library(broom)

eps <- 10 * .Machine$double.eps
if (family$family == "binomial") {
  if (any(mu > 1 - eps) || any(mu < eps)) 
    warning("glm.fit: fitted probabilities numerically 0 or 1 occurred", 
            call. = FALSE)
}

# this return the following error message
# Error in family$family : object of type 'closure' is not subsettable


probit.resids <- augment(probit.fit) %>%
  mutate(p = 1 / (1 + exp(-.fitted)),
         warning = p > 1-eps)

arrange(probit.resids, desc(.fitted)) %>%  
  select(2:5, p, warning) %>% 
  slice(1:10)


HMDA.nwarning <- filter(HMDA, !probit.resids$warning)

# using HMDA.nwarning should solve the problem...
probit.fit <- glm(deny ~ poly(hirat, 2), family = binomial, data = HMDA.nwarning)

Answer 1

В этом фрагменте кода

if (family$family == "binomial") {
  if (any(mu > 1 - eps) || any(mu < eps)) 
    warning("glm.fit: fitted probabilities numerically 0 or 1 occurred", 
            call. = FALSE)
}

есть функция binomial(), вызываемая при запуске glm с family == "binomial". Если вы посмотрите под glm (просто введите glm):

if (is.character(family)) 
        family <- get(family, mode = "function", envir = parent.frame())
    if (is.function(family)) 
        family <- family()
    if (is.null(family$family)) {
        print(family)
        stop("'family' not recognized")
    }

И функция glm проверяет binomial () $ family во время подгонки, и если любое из прогнозируемых значений отличается от 1 или 0 к EPS это предупреждение выдается.

Вам не нужно запускать эту часть, и да, вам нужно установить eps <- 10 * .Machine$double.eps. Итак, давайте запустим приведенный ниже код, и если вы запустите пробит, вам нужно указать link="probit" в биномиальном виде, в противном случае по умолчанию будет lo git:

library(AER)
library(tidyverse)
library(dplyr)
library(broom)

data(HMDA)

probit.fit <- glm(deny ~ poly(hirat, 2), family = binomial(link="probit"), data = HMDA)

eps <- 10 * .Machine$double.eps

probit.resids <- augment(probit.fit) %>%
  mutate(p = 1 / (1 + exp(-.fitted)),
         warning = p > 1-eps)

В столбце предупреждение указывает, если наблюдения повышаются предупреждение, в этом наборе данных есть одно:

table(probit.resids$warning)

FALSE  TRUE 
 2379     1

Мы можем использовать следующий шаг, чтобы отфильтровать его

HMDA.nwarning <- filter(HMDA, !probit.resids$warning)
dim(HMDA.nwarning)
[1] 2379   14

и повторно запустить регрессию:

probit.fit <- glm(deny ~ poly(hirat, 2), family = binomial(link="probit"), data = HMDA.nwarning)
coefficients(probit.fit)
(Intercept) poly(hirat, 2)1 poly(hirat, 2)2 
      -1.191292        8.708494        6.884404