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

Question

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

En particular, King y Zeng afirman

"... en los datos de sucesos raros los sesgos en las probabilidades pueden ser significativos con muestras de miles de ejemplares y están en una dirección predecible: las probabilidades estimadas de los sucesos 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 verdadero (y evidentemente muy pequeño) sesgo? ¿Requiere el sesgo que se incluya algún tipo de covariable (más que el intercepto)?

Actualización 1: King y Zeng incluyen una aproximación para el sesgo de $\beta_0$ en la ecuación 12 de su documento. Teniendo en cuenta la N en el denominador, reduje drásticamente N ser 5 y volví a ejecutar la simulación, pero seguía sin apreciarse ningún sesgo en las probabilidades de sucesos estimadas. (Lo he utilizado sólo como inspiración. Tenga en cuenta que mi pregunta anterior se refiere a las probabilidades de sucesos estimadas, no a las probabilidades de sucesos estimadas. $\hat \beta_0$ .)

Actualización 2: A raíz de una sugerencia en los comentarios, he incluido una variable independiente en la regresión, lo que ha dado lugar 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: He utilizado p como variable independiente, donde p es un vector con repeticiones de un valor pequeño (0,01) y un valor mayor (0,2). Al final, sim almacena únicamente las probabilidades estimadas correspondientes a $p = 0.01$ y no hay indicios de parcialidad.

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 por el recordatorio ). R no lanza una advertencia porque su " tolerancia de convergencia positiva se satisface realmente". En términos más liberales, la MLE existe y es menos infinito, lo que corresponde a $p=0$ de ahí la actualización de mi función. La única otra cosa coherente que se me ocurre hacer es descartar aquellas ejecuciones de la simulación en las que y es idénticamente cero, pero eso conduciría claramente a resultados aún más contrarios a la afirmación inicial de que "las probabilidades de sucesos estimadas son demasiado pequeñas".

Answer 1

Se trata de una cuestión interesante. He realizado algunas simulaciones que expongo a continuación con la esperanza de que susciten un debate más profundo.

En primer lugar, algunos comentarios generales:

• El artículo que citas trata sobre el sesgo de los sucesos raros. Lo que no me quedaba claro antes (también con respecto a los comentarios que se han hecho más arriba) es si hay algo especial en los casos en los que se tienen 10/10000 frente a 10/30 observaciones. Sin embargo, después de algunas simulaciones, estoy de acuerdo en que sí lo hay.

• Un problema que tenía en mente (me he encontrado con esto a menudo, y hubo recientemente un artículo en Methods in Ecology and Evolution sobre eso, no pude encontrar la referencia sin embargo) es que usted puede conseguir casos degenerados con GLMs en situaciones de datos pequeños, donde el MLE es FAAAR lejos de la verdad, o incluso en - / + infinito (debido a la relación no lineal supongo). No me queda claro cómo se deberían tratar estos casos en la estimación del sesgo, pero a partir de mis simulaciones diría que parecen clave para el sesgo de eventos raros. Mi intuición sería eliminarlos, pero entonces no está muy claro qué tan lejos tienen que estar para ser eliminados. Tal vez algo a tener en cuenta para la corrección del sesgo.

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

De todas formas, aquí hay algo de código que calcula la diferencia entre las estimaciones y la verdad para el intercepto, la pendiente y las predicciones en una regresión logística, primero para una situación de tamaño de muestra bajo / 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 y los errores estándar resultantes para el intercepto, la pendiente y la predicción son

-0.120429315  0.296453122 -0.001619793
 0.016105833  0.032835468  0.002040664

Yo concluiría que hay pruebas bastante sólidas de un ligero sesgo negativo en el intercepto y un ligero sesgo positivo en la pendiente, aunque si se observan los resultados trazados se ve que el sesgo es pequeño en comparación con la varianza de los valores estimados.

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

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

Obtengo un mayor sesgo para la intercepción, pero todavía 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, observamos el fenómeno de las estimaciones degeneradas de los parámetros (si es que debemos llamarlas así)

Vamos a eliminar todas las filas para las que 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 cifras: las estimaciones de los parámetros claramente no están distribuidas normalmente. Me pregunto qué significa esto para la validez de los IC que se presentan.

-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 el intercepto se debe a los propios eventos raros, es decir, a 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 punto de corte.

Sin embargo, es importante señalar que, en cualquier caso, no parece haber ningún sesgo en las predicciones a escala de respuesta: la función de enlace simplemente absorbe estos valores extremadamente pequeños.

Answer 2

La figura 7 del documento parece abordar más directamente la cuestión del sesgo en las predicciones. No entiendo del todo la figura (en concreto, la interpretación "las probabilidades de sucesos estimadas son demasiado pequeñas" parece una simplificación excesiva), pero he conseguido reproducir algo parecido basándome en la escueta 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')

El primer gráfico es mi réplica de su figura 7, con la adición de curvas discontinuas que representan toda la gama de resultados a lo largo de 10 ensayos.

Según el documento, x aquí es una variable predictora en la regresión extraída de una normal estándar. Por lo tanto, como se ilustra en el segundo gráfico, la frecuencia relativa de las observaciones de x > 3 (donde se produce el sesgo más visible en el primer gráfico) es cada vez menor.

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 se produce cuando x es rara o inexistente.

En conjunto, esto sugiere que el PO malinterpretó por completo la afirmación central del artículo al confundir "acontecimiento raro" (es decir. x > 3 ) con "acontecimiento para el que P(Y = 1) es muy pequeño". Es de suponer que el documento se refiere a lo primero y no a lo segundo.

Answer 3

El sesgo de eventos raros sólo se produce cuando hay regresores. No ocurrirá en un modelo de sólo intercepción como el simulado aquí. Véase esta entrada para d http://statisticalhorizons.com/linear-vs-logistic#comment-276108