Sesgo de regresión logística de eventos raros: ¿cómo simular las p subestimadas con un ejemplo mínimo?

CrossValidated tiene varias preguntas sobre cuándo y cómo aplicar la corrección de sesgo de eventos raros por King y Zeng (2001) . Estoy buscando algo diferente: una demostración mínima basada en simulación de que existe el sesgo.

En particular, el estado de King y Zeng

"... en datos de eventos raros, los sesgos en las probabilidades pueden ser significativamente significativos con tamaños de muestra en miles y están en una dirección predecible: las probabilidades de eventos estimadas son demasiado pequeñas".

Aquí está mi intento de simular tal sesgo en 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)

Cuando ejecuto esto, tiendo a obtener puntuaciones z muy pequeñas, y el histograma de estimaciones está muy cerca de centrarse sobre la verdad p = 0.01.

¿Qué me estoy perdiendo? ¿Es que mi simulación no es lo suficientemente grande como para mostrar el sesgo verdadero (y evidentemente muy pequeño)? ¿El sesgo requiere algún tipo de covariable (más que la intercepción) para ser incluido?

Actualización 1: King y Zeng incluyen una aproximación aproximada del sesgo de en la ecuación 12 de su artículo. Observando el en el denominador, reduje drásticamente a ser y volví a ejecutar la simulación, pero aún no hay sesgo en las probabilidades estimadas de eventos. (Utilicé esto solo como inspiración. Tenga en cuenta que mi pregunta anterior es sobre las probabilidades estimadas de eventos, no ).$\beta_0$NN5$\hat \beta_0$

Actualización 2: Siguiendo una sugerencia en los comentarios, incluí una variable independiente en la regresión, lo que condujo a resultados equivalentes:

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))

Explicación: me utilicé pcomo la variable independiente, donde pes un vector con repeticiones de un valor pequeño (0.01) y un valor mayor (0.2). Al final, simalmacena solo las probabilidades estimadas correspondientes a y no hay signos de sesgo.$p = 0.01$

Actualización 3 (5 de mayo de 2016): Esto no cambia notablemente los resultados, pero mi nueva función de simulación interna es

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])
    }
}

Explicación: El MLE cuando y es idénticamente cero no existe ( gracias a los comentarios aquí para el recordatorio ). R no puede lanzar una advertencia porque su " tolerancia de convergencia positiva " en realidad se satisface. Hablando más liberalmente, el MLE existe y es menos infinito, que corresponde a ; de ahí la actualización de mi función. La única otra cosa coherente que se me ocurre hacer es descartar esas ejecuciones de la simulación donde y es idénticamente cero, pero eso claramente conduciría a resultados aún más contrarios a la afirmación inicial de que "las probabilidades estimadas de eventos son demasiado pequeñas".$p=0$

Esta es una pregunta interesante: he hecho algunas simulaciones que publico a continuación con la esperanza de que esto estimule una mayor discusión.

En primer lugar, algunos comentarios generales:

• El artículo que cita es sobre el sesgo de eventos raros. Lo que no estaba claro para mí antes (también con respecto a los comentarios que se hicieron anteriormente) es si hay algo especial sobre los casos en los que tiene 10/10000 en lugar de 10/30 observaciones. Sin embargo, después de algunas simulaciones, estaría de acuerdo en que sí.

• Un problema que tenía en mente (lo he encontrado con frecuencia, y recientemente había un artículo en Métodos en Ecología y Evolución sobre eso, aunque no pude encontrar la referencia) es que puede obtener casos degenerados con GLM en datos pequeños situaciones, donde el MLE está lejos de la verdad, o incluso en - / + infinito (debido al enlace no lineal, supongo). No está claro para mí cómo se deben tratar estos casos en la estimación del sesgo, pero según mis simulaciones, diría que parecen ser clave para el sesgo de eventos raros. Mi intuición sería eliminarlos, pero no está claro qué tan lejos deben estar para eliminarlos. Tal vez algo a tener en cuenta para la corrección de sesgos.

• Además, estos casos degenerados parecen propensos a causar problemas numéricos (por lo tanto, he aumentado el maxit en la función glm, pero también se podría pensar en aumentar el épsilon para asegurarme de que uno realmente informe el verdadero MLE).

De todos modos, aquí hay un código que calcula la diferencia entre las estimaciones y la verdad para la intercepción, la pendiente y las predicciones en una regresión logística, primero para una situación de bajo tamaño de muestra / incidencia moderada:

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)

El sesgo resultante y los errores estándar de intercepción, pendiente y predicción son

-0.120429315  0.296453122 -0.001619793
 0.016105833  0.032835468  0.002040664

Llegaría a la conclusión de que hay bastante buena evidencia de un ligero sesgo negativo en la intersección, y un ligero sesgo positivo en la pendiente, aunque un vistazo a los resultados trazados muestra que el sesgo es pequeño en comparación con la varianza de los valores estimados.

Si estoy configurando los parámetros a una situación de evento raro

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

Estoy obteniendo un sesgo mayor para la intercepción, pero aún NINGUNO en la predicción

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

En el histograma de los valores estimados, vemos el fenómeno de las estimaciones de parámetros degenerados (si los llamamos así)

Eliminemos todas las filas para las cuales las estimaciones de intercepción son <20

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

El sesgo disminuye y las cosas se vuelven un poco más claras en las figuras: las estimaciones de los parámetros claramente no se distribuyen normalmente. Me pregunto que eso significa para la validez de los IC que se informan.

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

Llegaría a la conclusión de que el sesgo de eventos raros en la intersección está impulsado por eventos raros en sí, es decir, esas estimaciones raras y extremadamente pequeñas. No estoy seguro de si queremos eliminarlos o no, no estoy seguro de cuál sería el límite.

Sin embargo, una cosa importante a tener en cuenta es que, de cualquier manera, parece que no hay sesgo en las predicciones en la escala de respuesta: la función de enlace simplemente absorbe estos valores extremadamente pequeños.

El sesgo de eventos raros solo ocurre cuando hay regresores. No ocurrirá en un modelo de solo intercepción como el simulado aquí. Vea esta publicación para más detalles: http://statisticalhorizons.com/linear-vs-logistic#comment-276108

La Figura 7 en el documento parece abordar más directamente la cuestión del sesgo en las predicciones. No entiendo completamente la figura (específicamente, la interpretación "las probabilidades de evento estimadas son demasiado pequeñas" parece una simplificación excesiva) pero logré reproducir algo similar en base a su breve descripción de su simulación en la Sección 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')

La primera gráfica es mi réplica de su figura 7, con la adición de curvas discontinuas que representan el rango completo de resultados en 10 ensayos.

Según el artículo, xaquí hay una variable predictora en la regresión extraída de una normal estándar. Por lo tanto, como se ilustra en la segunda gráfica, la frecuencia relativa de observaciones para x > 3(donde ocurre el sesgo más visible en la primera gráfica) es decrecientemente pequeña.

El tercer gráfico muestra las probabilidades de simulación "verdaderas" en el proceso de generación en función de x. Parece que el mayor sesgo ocurre donde xes raro o inexistente.

Tomados en conjunto, estos sugieren que el OP malinterpretó por completo la afirmación central del documento al confundir "evento raro" (es decir x > 3) con "evento para el cual P(Y = 1)es muy pequeño". Presumiblemente, el documento se refiere al primero en lugar del segundo.