Смещение логистической регрессии редких событий: как смоделировать недооцененные p с минимальным примером?

У CrossValidated есть несколько вопросов о том, когда и как применять коррекцию смещения редкого события, разработанную King and Zeng (2001) . Я ищу что-то другое: минимальную демонстрацию, основанную на симуляции, которая существует.

В частности, король и дзенг

«... в данных по редким событиям смещения вероятностей могут быть существенно значимыми с размерами выборки в тысячах и имеют предсказуемое направление: оценочные вероятности событий слишком малы».

Вот моя попытка симулировать такой уклон в R:

# FUNCTIONS
do.one.sim = function(p){
    N = length(p)
    # Draw fake data based on probabilities p  
    y = rbinom(N, 1, p)  
    # Extract the fitted probability.
    #    If p is constant, glm does y ~ 1, the intercept-only model.
    #    If p is not constant, assume its smallest value is p[1]:
    glm(y ~ p, family = 'binomial')$fitted[1]
}
mean.of.K.estimates = function(p, K){
    mean(replicate(K, do.one.sim(p) ))
}

# MONTE CARLO
N = 100
p = rep(0.01, N)
reps = 100
# The following line may take about 30 seconds
sim = replicate(reps, mean.of.K.estimates(p, K=100))
# Z-score:
abs(p[1]-mean(sim))/(sd(sim)/sqrt(reps))
# Distribution of average probability estimates:
hist(sim)

Когда я запускаю это, я, как правило, получаю очень маленькие z-оценки, и гистограмма оценок очень близка к центру по истине p = 0,01.

Чего мне не хватает? Неужели моя симуляция недостаточно велика, чтобы показать истинный (и, очевидно, очень маленький) уклон? Требует ли смещение некоторого ковариата (больше, чем перехват), который будет включен?

Обновление 1: Кинг и Цзэн включают грубое приближение для смещения в уравнении 12 своей статьи. Отмечая в знаменателе, я резко сократил быть и повторно запускал моделирование, но до сих пор нет смещения расчетных вероятностей событий не очевидно. (Я использовал это только для вдохновения. Обратите внимание, что мой вопрос выше касается оценочных вероятностей событий, а не .)$\beta_0$NN5$\hat \beta_0$

Обновление 2: Следуя предложению в комментариях, я включил в регрессию независимую переменную, что привело к эквивалентным результатам:

p.small = 0.01
p.large = 0.2
p = c(rep(p.small, round(N/2) ), rep(p.large, N- round(N/2) ) )
sim = replicate(reps, mean.of.K.estimates(p, K=100))

Объяснение: Я использовал pсебя в качестве независимой переменной, где pесть вектор с повторениями небольшого значения (0,01) и большего значения (0,2). В конце simсохраняются только оценочные вероятности, соответствующие и нет признаков смещения.$p = 0.01$

Обновление 3 (5 мая 2016 г.): это заметно не меняет результаты, но моя новая функция внутреннего моделирования

do.one.sim = function(p){
    N = length(p)
    # Draw fake data based on probabilities p  
    y = rbinom(N, 1, p)
    if(sum(y) == 0){ # then the glm MLE = minus infinity to get p = 0
        return(0)
    }else{
        # Extract the fitted probability.
        #    If p is constant, glm does y ~ 1, the intercept only model.
        #    If p is not constant, assume its smallest value is p[1]:
        return(glm(y ~ p, family = 'binomial')$fitted[1])
    }
}

Объяснение: MLE, когда y тождественно равен нулю, не существует ( спасибо за комментарии здесь за напоминание ). R не может выдать предупреждение, потому что его " положительный допуск сходимости " фактически удовлетворен. Более свободно говоря, MLE существует и является минус бесконечность, что соответствует ; отсюда и моя функция обновления. Единственная другая связная вещь, которую я могу придумать, - это отбросить те прогоны симуляции, где y тождественно равен нулю, но это явно приведет к результатам, еще более противоречащим первоначальному утверждению, что «оценочные вероятности события слишком малы».$p=0$

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

Прежде всего, несколько общих комментариев:

• В статье, которую вы цитируете, речь идет о предвзятости. То, что мне не было ясно раньше (также в отношении комментариев, которые были сделаны выше), это что-то особенное в случаях, когда у вас 10/10000, а не 10/30 наблюдений. Однако после некоторых симуляций я бы согласился, что есть.

• Проблема, которую я имел в виду (я часто сталкивался с этим, и недавно в статье «Методы в области экологии и эволюции» я писал об этом, но я не смог найти ссылку на нее) заключается в том, что вы можете получить вырожденные случаи с GLM в небольших данных. ситуации, когда MLE находится на расстоянии FAAAR от истины или даже на бесконечности - / + (полагаю, из-за нелинейной связи). Мне не ясно, как следует относиться к этим случаям при оценке смещения, но из моих моделей я бы сказал, что они кажутся ключевыми для смещения редких событий. Моя интуиция заключается в том, чтобы удалить их, но тогда не совсем ясно, как далеко они должны быть удалены. Может быть, что-то иметь в виду для исправления предвзятости.

• Кроме того, эти вырожденные случаи кажутся склонными вызывать численные проблемы (поэтому я увеличил максимальное значение в функции glm, но можно подумать и об увеличении эпсилона, чтобы удостовериться, что кто-то действительно сообщает истинное MLE).

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

set.seed(123)
replicates = 1000
N= 40
slope = 2 # slope (linear scale)
intercept = - 1 # intercept (linear scale)

bias <- matrix(NA, nrow = replicates, ncol = 3)
incidencePredBias <- rep(NA, replicates)

for (i in 1:replicates){
  pred = runif(N,min=-1,max=1) 
  linearResponse = intercept + slope*pred
  data = rbinom(N, 1, plogis(linearResponse))  
  fit <- glm(data ~ pred, family = 'binomial', control = list(maxit = 300))
  bias[i,1:2] = fit$coefficients - c(intercept, slope)
  bias[i,3] = mean(predict(fit,type = "response")) - mean(plogis(linearResponse))
}

par(mfrow = c(1,3))
text = c("Bias intercept", "Bias slope", "Bias prediction")

for (i in 1:3){
  hist(bias[,i], breaks = 100, main = text[i])
  abline(v=mean(bias[,i]), col = "red", lwd = 3)  
}

apply(bias, 2, mean)
apply(bias, 2, sd) / sqrt(replicates)

Результирующее смещение и стандартные ошибки для перехвата, наклона и прогнозирования

-0.120429315  0.296453122 -0.001619793
 0.016105833  0.032835468  0.002040664

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

Если я устанавливаю параметры в ситуации редкого события

N= 4000
slope = 2 # slope (linear scale)
intercept = - 10 # intercept (linear scale)

Я получаю больший уклон для перехвата, но все еще ничего в прогнозе

   -1.716144e+01  4.271145e-01 -3.793141e-06
    5.039331e-01  4.806615e-01  4.356062e-06

На гистограмме оценочных значений мы видим явление вырожденных оценок параметров (если их так называть)

Давайте удалим все строки, для которых оценки перехвата <20

apply(bias[bias[,1] > -20,], 2, mean)
apply(bias[bias[,1] > -20,], 2, sd) / sqrt(length(bias[,1] > -10))

Смещение уменьшается, и на рисунках все становится немного яснее - оценки параметров явно не распределяются нормально. Интересно, что это означает для достоверности CI, о которых сообщают.

-0.6694874106  1.9740437782  0.0002079945
1.329322e-01 1.619451e-01 3.242677e-06

Я бы пришел к выводу, что смещение редкого события на перехвате обусловлено самими редкими событиями, а именно теми редкими, чрезвычайно малыми оценками. Не уверен, хотим мы их удалить или нет, не уверен, какой будет отсечка.

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

Смещение редких событий происходит только при наличии регрессоров. Это не произойдет в модели только для перехвата, подобной моделируемой здесь. См. Этот пост для деталей: http://statisticalhorizons.com/linear-vs-logistic#comment-276108

Рисунок 7 в документе, кажется, наиболее непосредственно касается вопроса о предвзятости в прогнозах. Я не совсем понимаю эту цифру (в частности, интерпретация «оценочные вероятности событий слишком малы» кажется чрезмерным упрощением), но мне удалось воспроизвести нечто подобное, основываясь на кратком описании их моделирования в разделе 6.1:

n_grid = 40
x_grid = seq(0, 7, length.out = n_grid)
beta0 = -6
beta1 = 1

inverse_logit = function(x) 1/(1 + exp(-x))

do.one.sim = function(){
    N = 5000
    x = rnorm(N)
    p = inverse_logit(beta0 + beta1*x)
    # Draw fake data based on probabilities p
    y = rbinom(N, 1, p)
    if(sum(y) == 0){ # then the glm MLE = minus infinity to get p = 0
        return(rep(0, n_grid))
    }else{
        # Extract the error
        mod = glm(y ~ x, family = 'binomial')
        truth = inverse_logit(beta0 + beta1*x_grid)
        pred = predict(mod, newdata = data.frame(x = x_grid),
            type = 'response')
        return(pred - truth)
    }
}
mean.of.K.estimates = function(K){
    rowMeans(replicate(K, do.one.sim()))
}

set.seed(1)
bias = replicate(10, mean.of.K.estimates(100))
maxes = as.numeric(apply(bias, 1, max))
mins = as.numeric(apply(bias, 1, min))

par(mfrow = c(3, 1), mar = c(4,4,2,2))
plot(x_grid, rowMeans(bias), type = 'l',
    ylim = c(min(bias), max(bias)),
    xlab = 'x', ylab = 'bias')
lines(x_grid, maxes, lty = 2)
lines(x_grid, mins, lty = 2)
plot(x_grid, dnorm(x_grid), type = 'l',
    xlab = 'x', ylab = 'standard normal density')
plot(x_grid, inverse_logit(beta0 + beta1*x_grid),
    xlab = 'x', ylab = 'true simulation P(Y = 1)',
    type = 'l')

Первый график - моя репликация их фигуры 7 с добавлением пунктирных кривых, представляющих весь спектр результатов за 10 испытаний.

Согласно статье, xздесь есть предикторная переменная в регрессии, взятая из стандартной нормали. Таким образом, как показано на втором графике, относительная частота наблюдений для x > 3(где наиболее заметное отклонение происходит на первом графике) уменьшается все меньше.

Третий график показывает «истинные» вероятности моделирования в процессе генерации как функцию от x. Похоже, что наибольший уклон происходит там, где xон редок или отсутствует.

Взятые вместе, они предполагают, что ФП полностью неверно истолковал центральное утверждение статьи, путая «редкое событие» (т.е. x > 3 ) с «событием, для которого P(Y = 1)очень мало». Предположительно, статья касается первого, а не второго.