Creo que todos sabéis la historia de las admisiones de la Universidad de Berkeley y la paradoja de Simpson. Con palabras, muchas palabras, está contado, por ejemplo, aquí. Y si buscáis ubc admissions simpson en Google la encontraréis también en modo --verbose
en muchos más sitios.
En R puede resumirse en
Echad un vistazo a los coeficientes de Gender
en ambos modelos y veréis.
Lo que se me ha ocurrido es lo siguiente. Imaginad que hombres y mujeres solicitan admisiones a esos departamentos sin que el sexo induzca sesgo alguno en las preferencias. Por su parte, los departamentos aceptan o rechazan candidatos en proporciones fijas (e iguales a las del caso que nos ocupa).
De otra manera, un candidato tiene una probabilidad de ser hombre, una probabilidad
de solicitar la admisión en el departamento
y una probabilidad
de ser admitido en él.
La idea es comparar las regresiones logísticas anteriores: una que incluya Dept
y la otra no. Lo que yo esperaba es que los coeficientes de Gender
en muchas simulaciones:
- Estuviesen centrados en cero en ambos casos: por construcción, el efecto del sexo es nulo.
- Que la varianza del del modelo que incluye
Dept
fuese menor. Esto está inspirado por el ejemplo anterior. En los casos en que por efecto del muestreo haya un desequilibrio por sexos en las solicitudes a los distintos departamentos —como ocurrió en el ejemplo con el que abro la entrada—, controlar porGender
ayudaría a corregir el sesgo. Reducir el sesgo en cada simulación conlleva reducir la varianza global, ¿no?
Pues parece que algo se me ha escapado porque, salvo error u omisión,
accept.rate <- dcast(raw, Dept ~ Admit, fun.aggregate = sum, value.var = "Freq") accept.rate$size <- accept.rate$Admitted + accept.rate$Rejected accept.rate$rate <- accept.rate$Admitted / accept.rate$size gender.rate <- melt(dat) gender.rate <- tapply(gender.rate$value, gender.rate$Gender, sum) extract.coefs <- function(accept.rate, gender.rate){ Gender <- sample(names(gender.rate), sum(gender.rate), prob = gender.rate / sum(gender.rate), replace = T) Dept <- sample(accept.rate$Dept, sum(gender.rate), prob = accept.rate$rate, replace = T) Admitted <- accept.rate$rate[match(Dept, accept.rate$Dept)] Admitted <- runif(length(Admitted)) < Admitted dat <- data.frame(Gender = Gender, Dept = Dept, Admitted = Admitted) dat <- ddply(dat, .(Gender, Dept), summarize, Rejected = sum(!Admitted), Admitted = sum(Admitted)) mod.0 <- glm(cbind(Admitted, Rejected) ~ Gender, data = dat, family = binomial) mod.1 <- glm(cbind(Admitted, Rejected) ~ Gender + Dept, data = dat, family = binomial) c(coefficients(mod.0)[2], coefficients(mod.1)[2]) } res <- replicate(1000, extract.coefs(accept.rate, gender.rate)) res <- as.data.frame(t(res)) colnames(res) <- c("mod0", "mod1") boxplot(res, col = "gray")
produce
justo al contrario de lo esperado. ¡Carajo!
P.D.: ¿Igual n
es demasiado grande como para que se manifiesten las ventajas de controlar por Dept
?
Curioso. No sé si va por aquí, pero lo enlazo de todos modos http://scott.fortmann-roe.com/docs/BiasVariance.html