Averigüémoslo.
Para empezar, ¿qué ocurre con los conjuntos de datos equilibrados?
Este es un gráfico de dispersión de un conjunto de datos de $200$ observaciones, de las cuales $50\%$ son ceros y el resto son unos. En él he representado gráficamente las probabilidades subyacentes ("verdaderas") y las probabilidades ajustadas con la regresión logística. Los dos gráficos coinciden estrechamente, lo que indica que la regresión logística hizo un buen trabajo en este caso.
Para entenderlo mejor, he mantenido el mismo $x$ pero regeneró los valores de $y$ valores al azar $500$ tiempos. Cada ajuste produjo su estimación de la intercepción y la pendiente ( $\hat\beta$ ) en esta regresión logística. Aquí hay un gráfico de dispersión de esas estimaciones.
El triángulo rojo central representa los coeficientes reales $(0, -3).$ Las elipses son las aproximaciones de segundo orden a esta nube de puntos: una pretende encerrar aproximadamente la mitad de los puntos y la otra pretende encerrar aproximadamente el 95% de los puntos. El hecho de que lo hagan indica que dan una indicación sólida de la incertidumbre que puede tener cualquier estimación de ese conjunto de datos: el intercepto podría estar desviado en unos $\pm 0.45$ (la anchura de la elipse exterior) y la pendiente podría estar desviada en aproximadamente $\pm 1$ (la altura de la elipse exterior). Estos son márgenes de error.
¿Qué ocurre con los conjuntos de datos desequilibrados?
Aquí hay una configuración similar pero con sólo $5\%$ de los puntos de una clase (más o menos puntos, en función de la aleatoriedad de estas observaciones):
( $5\%$ es realmente pequeño: nos dice que esperemos ver sólo $10$ de una clase con la otra $190$ en la otra clase).
El ajuste ahora se aleja visiblemente del gráfico real, pero ¿es esto una prueba de que la regresión logística no es "robusta"? De nuevo, podemos averiguarlo repitiendo el proceso de generación de datos aleatorios y estimando el ajuste muchas veces. Aquí está el gráfico de dispersión de $500$ estimaciones.
En general, las estimaciones se acercan al valor real cerca de $(-4,-3).$ En este sentido, la regresión logística parece "robusta". (He mantenido la misma pendiente de $-3$ como antes y ajustamos el intercepto para reducir la tasa de $+1$ observaciones).
Los márgenes de error han cambiado: el semieje de la elipse exterior que (más o menos) describe la incertidumbre en el intercepto ha pasado de $0.45$ a más de $4$ mientras que el otro semieje se ha reducido un poco de $1$ a cerca de $0.8;$ y todo el panorama se ha inclinado.
Las elipses ya no describen la nube de puntos tan bien como antes: ahora, hay cierta tendencia a que la regresión logística estime pendientes e interceptos extremadamente negativos. La inclinación indica una correlación notable entre las estimaciones: los interceptos bajos (negativos) tienden a asociarse con pendientes negativas bajas (que compensan los interceptos pequeños al predecir algunos $1$ valores cercanos a $x=-1.$ ) Pero esta correlación es de esperar: se parece a la regresión por mínimos cuadrados ordinarios siempre que el punto de las medias de los datos no esté cerca del eje vertical.
¿Qué muestran estos experimentos?
Para conjuntos de datos de este tamaño (o mayores), al menos:
-
La regresión logística suele funcionar bien y dar valores razonablemente cercanos a los parámetros correctos incluso cuando los resultados están desequilibrados.
-
Las descripciones de segundo orden de la correlación entre las estimaciones de los parámetros (que son resultados rutinarios de la regresión logística) no captan del todo la posibilidad de que las estimaciones puedan estar simultáneamente muy lejos de la verdad.
Una meta-conclusión
Se puede evaluar la "robustez" (o, de forma más general, las propiedades estadísticas más destacadas) de cualquier procedimiento, como la regresión logística, ejecutándolo repetidamente con datos generados según un modelo realista conocido y haciendo un seguimiento de los resultados que son importantes para usted.
Esta es la R
código que produjo las cifras. Para las dos primeras figuras, la primera línea se modificó como p <- 50/100
. Retire el set.seed
para generar ejemplos aleatorios adicionales.
Experimentar con simulaciones como ésta (ampliada a más variables explicativas) podría persuadirle de la utilidad de un regla general:
Deje que el número de observaciones en la clase más pequeña guíe la complejidad del modelo.
Mientras que en la regresión por mínimos cuadrados ordinarios podríamos estar cómodos teniendo diez observaciones (en total) para cada variable explicativa, para la regresión logística querremos tener diez observaciones en la clase más pequeña para cada variable explicativa.
p <- 5/100 # Proportion of one class
n <- 200 # Dataset size
x <- seq(-1, 1, length.out=n) # Explanatory variable
beta <- -3 # Slope
#
# Find an intercept that yields `p` as the true proportion for these `x`.
#
logistic <- function(z) 1 - 1/(1 + exp(z))
alpha <- uniroot(function(a) mean(logistic(a + beta*x)) - p, c(-5,5))$root
#
# Create and plot a dataset with an expected value of `p`.
#
set.seed(17)
y <- rbinom(n, 1, logistic(alpha + beta*x))
plot(range(x), c(-0.015, 1.015), type="n", bty="n", xlab="x", ylab="y",
main="Data with True (Solid) and Fitted (Dashed) Probabilities")
curve(logistic(alpha + beta*x), add=TRUE, col="Gray", lwd=2)
rug(x[y==0], side=1, col="Red")
rug(x[y==1], side=3, col="Red")
points(x, y, pch=21, bg="#00000020")
#
# Fit a logistic model.
#
X <- data.frame(x=x, y=y)
fit <- glm(y ~ x, data=X, family="binomial")
summary(fit)
#
# Sketch the fit.
#
b <- coefficients(fit)
curve(logistic(b[1] + b[2]*x), add=TRUE, col="Black", lty=3, lwd=2)
#
# Evaluate the robustness of the fit.
#
sim <- replicate(500, {
X$y.new <- with(X, rbinom(n, 1, logistic(alpha + beta*x)))
coefficients(glm(y.new ~ x, data=X, family="binomial"))
})
plot(t(sim), main="Estimated Coefficients", ylab="")
mtext(expression(hat(beta)), side=2, line=2.5, las=2, cex=1.25)
points(alpha, beta, pch=24, bg="#ff0000c0", cex=1.6)
#
# Plot second moment ellipses.
#
V <- cov(t(sim))
obj <- eigen(V)
a <- seq(0, 2*pi, length.out=361)
for (level in c(.50, .95)) {
lambda <- sqrt(obj$values) * sqrt(qchisq(level, 2))
st <- obj$vectors %*% (rbind(cos(a), sin(a)) * lambda) + c(alpha, beta)
polygon(st[1,], st[2,], col="#ffff0010")
}