sábado, enero 17, 2015

Dicotomizar variables numéricas en regresión, informar únicamente de resultados significativos: Malas ideas.

Imaginemos que queremos ver si hay relación entre el índice de masa corporal (kg / m^2) y la depresión, evaluada ésta a través de alguno de los muchos tests disponibles. Tendríamos, pues, dos variables representadas numéricamente. La técnica de análisis de elección sería, pues, una correlación o su prima hermana, la regresión lineal.

Otra opción habitualmente aplicada es dicotomizar la variable predictora. Podríamos cambiar el IMC, con sus infinitos valores posibles, por una nueva variable donde marcásemos si las personas presenta infranormopeso (IMC < 25) o sobrepeso y obsesidad (IMC > 25). Así, podríamos aplicar una t de Student.

Algunos prefieren este modo de analizar datos. Primero, por tradición. En algunos campos como psicología experimental es más común el aplicar análisis de variance o pruebas t que regresiones. Segundo, porque algunos investigadores estiman que es más comprensible la explicación de los resultados si éstos están basados en grupos. Puede ser más descriptivo hablar de "normopeso frente a sobrepeso" que hablar de "cambio esperado con incrementos de IMC".

Otro motivo habitual para dicotomizar es la presencia de más de dos variables predictoras. Si sospechamos que la relación entre IMC y depresión no es igual para hombres y mujeres, estamos suponiendo que la presencia de una interacción entre una variable numérica y otra categórica. Un modelo de regresión con un término de interacción es bien sencillo de calcular e interpretar, pero no suele verse en la carrera y muchos investigadores no se sienten cómodos con ese análisis. Desgraciadamente, el SPSS, programa de análisis estadístico de referencia en Psicología, no ofrece directamente estos modelos, sino que son necesarios algunos pasos como creación de variables o trabajar con la sintaxis. Sin embargo, los investigadores están mucho más acostumbrados a trabajar con interacciones en un ANOVA y el SPSS las calcula por defecto. Esto ayuda a explicar el recurso a la dicotomización.

El dicotimizar variables presenta amplios problemas. Éstos han sido extensamente descritos en la literatura científica. Demos una primera dos aproximaciones intuitivas:
- Al dicotomizar variables incrementamos algunas diferencias mínimas y minimizamos algunas diferencias enormes. En el ejemplo en el que nos encontramos, dos personas, una con un IMC de 24.9 y otra con IMC de 25.1, serán asignados a dos grupos separados, tan separados como puedan estar. Y dos personas con IMCs de 17 y 24 pertenecerán al mismo grupo.
- Imaginemos que queremos calcular la correlación entre nuestra altura calculada en metros y nuestra altura calculada en centímetros. Será de 1, porque conocida una variable podemos pronosticar sin error el valor en la otra, y la recta es creciente. Sin embargo, si intento pronosticar la altura en metros mediante una variable que marque a los participantes como bajos = 0 o altos = 1, será imposible encontrar una correlación perfecta. El saber si alguien ha sido marcado como alto no me permite determinar a la perfección su altura en metros.

Esto supone que la dicotomización del predictor ha de reducir el tamaño estimado de la asociación entre variables. Vamos a verlo con un pequeño estudio de simulación. La sintaxis en R la encontrarán al final del post.

Hemos trabajado con tamaños muestrales iguales a 50. Estos tamaños no son infrecuentes en psicología experimental o con muestras clínicas. Hemos simulado 10000 réplicas donde, conociendo la correlación poblacional real (valores desde 0 hasta 1 con incrementos de 0.01), vemos cuál es la correlación muestral estimada. Las dos variables siguen una distribución normal. La variable predictora se dicotomiza por la media. Tanto para el caso continuo como para el caso dicotomizado, la relación entre variables se evalúa con una correlación de Pearson.

Veamos los resultados. En la siguiente figura se representan los tres cuartiles de la correlación muestral (en negro, mediana; en rojo, percentiles 25 y 75) según valor de la correlación poblacional, tanto para la variable continua como para la dicotomizada. La línea de la diagonal sirve de referencia. En la medida que las líneas se acercan a esa diagonal, menor sesgo habrá.

Se observa que para la variable continua la mediana de la correlación recuparada coincide con el valor poblacional. La variabilidad en las correlaciones recuperadas (la separación entre líneas rojas) se va reduciendo a medida que aumenta la correlación poblacional. (Si tenemos una correlación poblacional igual a 0.95, es complicado tener una alta variabilidad).

Tal y como cabía esperar, la situación cambia para el predictor dicotimizado. Lo más reseñable es cómo este paso de variable numérica a variable categórica reduce el tamaño de estimado de la relación entre variables. Para el caso de una relación perfecta entre variables originales, la transformación reduce la correlación a 0.80. Esta infraestimación está presente para cualquier valor de la correlación poblacional (excepto cuando no hay relación) y se hace tanto mayor cuando más se acerca la relación a la unidad. Para correlaciones reales medias y altas, ni el error muestral lleva a alcanzar estimaciones cercanas a la correlación real (línea roja superior).


Esta infraestimación de los tamaños de los efectos tiene traducción en la potencia estadística de los contraste. En la siguiente figura se muestra la proporción de correlaciones marcadas como diferentes de 0 de forma estadísticamente significativa (alfa = 0.05). La línea continua es el caso de variable numérica; la línea discontinua, la variable dicotomizada.

Por definición de alfa, la probabilidad de marcar como estadísticamente significativa una correlación muestral proviniente de una población donde las dos variables son independientes es igual a 0.05. Podemos ver que, efectivamente, para ambos tipos de variables esa probabilidad al comienzo de las líneas se sitúa en el valor que corresponde.

En nuestro caso, la potencia estadística es la probabilidad de marcar como estadísticamente signitificativo un efecto poblacionalmente no nulo. Por tanto, las líneas representan, para cualquier valor de correlación poblacional distinto de 0, la potencia estadística. Vemos cómo ambas líneas siguen una sigmoide (arrancan más bien planas, se aceleran y, cuando se acercan al máximo, reduce la pendiente). Para valores poblacionales medios-altos o mayores, la potencia es básicamente igual a 1, esto es, siempre se marcarán las correlaciones como estadísticamente significativas. En la comparación que más nos interesa, para cualquier valor de correlación poblacional, la potencia siempre es menor cuando trabajamos con variables dicotomizadas.


Lamentablemente, es común que los resultados estadísticamente no significativos tengan problemas para ser publicados. Nos pondremos en el peor escenario: ps < 0.05 suponen publicación, ps > 0.5 suponen que el estudio se queda en el cajón. Esto equivale a publicar únicamente los estudios que caen por debajo de las líneas de la figura anterior y a malgastar los estudios que caen por encima.

En un caso así, los datos publicados no son una muestra representativa de los resultados posibles. Aquellos estudios con tamaños del efecto inflados por el error muestral tendrán más opciones de ser publicados. Esto lo vemos en la siguiente figura, donde se muestran el promedio de las correlaciones estadísticamente significativas. A medida que avanzamos por la correlación poblacional vamos aumentando desde un 5% de los estudios disponibles hacia el 100% de los estudios disponibles.

Cuando la correlación poblacional es igual a 0, la media de las correlaciones significativas es básicamente igual a 0, también: las correlaciones significativas lo son tanto con signo positivo como negativo y, así, se anulan. Sería un escenario en el que tendríamos correlaciones publicadas absolutamente contradictorias.

Para correlaciones menores a, aproximadamente, 0.4, el informar únicamente de los resultados significativos sesga positivamente y de forma marcada la correlación estimada. De un modo interesante, para un rango de valores poblacionales de la correlación (en torno a [0.10, 0.30]), apenas hay cambios en la correlación muestral. Valores muestrales claramente diferentes suponen valores recuperados casi iguales. Esto supone para variables tanto originales como dicotomizadas. Por último, para valores de correlación poblacional medios y mayores, dado que la potencia estadística se aproxima a 1, no hay diferencia entre lo esta figura y la anterior.



Tenemos, pues, que la dicotomización de variables sesga negativamente la distribución de las correlaciones muestrales y reduce la potencia estadística. El informar únicamente de resultados estadísticamente significativos infla de forma marca la estimación del tamaño de los efectos, salvo que la potencia estadística sea alta (y, así, pocos estudios queden en el cajón).

No nos hemos situado en el caso de dicotomización con mayor impacto. Si en lugar de dicotomizar por la media (o mediana), lo hacemos con un valor que haga los tamaños de los grupos muy desiguales, mayor será la reducción en la correlación recuperada. Si hacemos un análisis con grupos extremos (nos quedamos con el primer cuartil y el cuarto cuartil y tiramos el resto de los datos), nos metemos en un follón todavía peor.

Teniendo todo esto en cuenta, se desaconseja la dicotomización de variables y, ojalá, algún día los resultados no significativos pueden ser evaluados al margen de su valor-p. Hay multitud de estudios de simulación que señalan el problema que yo he ilustrado aquí. Los manuales de estadística indican el modo más adecuado de analizar datos.

La sintaxis permite ver el efecto para otros tamaños muestrales, distribuciones o modos de dicotomizar. Hay que tener en cuenta que el número de puntos con los que se evalúa el efecto de la correlación muestral (101) y el número de réplicas (10000) puede hacer lenta la ejecución del código. Para 'echar un vistazo', se pueden reducir estos números. Seguro que hay mil modos de acelarar la velocidad de ejecución del código. Soy un poco 'peatón' de R, asi que cualquier consejo será bien recibido.
rm(list=ls())
poblacion.cort1t2 <- seq(0, 1, 0.01)
res.cort1t2 <- matrix(nrow = length(poblacion.cort1t2), ncol = 3)
res.cord1t2 <- matrix(nrow = length(poblacion.cort1t2), ncol = 3)
res.sigt1t2 <- NULL
res.sigd1t2 <- NULL
res.corsigt1t2 <- NULL
res.corsigd1t2 <- NULL
muestra <- 50
replicas <- 10000
t1 <- replicate(replicas,rnorm(muestra))
et2 <- replicate(replicas,rnorm(muestra))
d1 <- t1>colMeans(t1)
d1[t1>colMeans(t1)] <- 1
d1[t1<=colMeans(t1)] <- 0
 
 
for (j in 1:length(poblacion.cort1t2)) {
  muestra.cort1t2 <- NULL
  muestra.pt1t2 <- NULL
  muestra.cord1t2 <- NULL
  muestra.pd1t2 <- NULL
  t2 <- t1*poblacion.cort1t2[j] + et2*(1-poblacion.cort1t2[j]^2)^.5
  for (i in 1:replicas) {
    cort1t2 <- cor.test(t1[,i],t2[,i])
    cord1t2 <- cor.test(d1[,i],t2[,i])
    muestra.cort1t2 <- c(muestra.cort1t2, cort1t2$estimate)
    muestra.pt1t2 <- c(muestra.pt1t2, cort1t2$p.value)
    muestra.cord1t2 <- c(muestra.cord1t2, cord1t2$estimate)
    muestra.pd1t2 <- c(muestra.pd1t2, cord1t2$p.value)
  }
  res.cort1t2[j,] <- quantile(muestra.cort1t2, c(.25, 0.50, .75))
  res.cord1t2[j,] <- quantile(muestra.cord1t2, c(.25, 0.50, .75))
  res.sigt1t2[j] <- mean(muestra.pt1t2<0.05)
  res.sigd1t2[j] <- mean(muestra.pd1t2<0.05)
  res.corsigt1t2[j] <- mean(muestra.cort1t2[muestra.pt1t2<0.05]) 
  res.corsigd1t2[j] <- mean(muestra.cord1t2[muestra.pd1t2<0.05])
}
 
par(mfrow=c(1,2))
plot(poblacion.cort1t2, poblacion.cort1t2, type="l", col="orange", lwd=5,
     xlab="Correlación poblacional", ylab="Correlación muestral", main="Variable continua")
lines(poblacion.cort1t2, res.cort1t2[,1], col="red", lwd=3)
lines(poblacion.cort1t2, res.cort1t2[,2], col="black", lwd=3)
lines(poblacion.cort1t2, res.cort1t2[,3], col="red", lwd=3)
plot(poblacion.cort1t2, poblacion.cort1t2, type="l", col="orange", lwd=5,
     xlab="Correlación poblacional", ylab="Correlación muestral", main="Variable dicotomizada")
lines(poblacion.cort1t2, res.cord1t2[,1], col="red", lwd=3, lty=2)
lines(poblacion.cort1t2, res.cord1t2[,2], col="black", lwd=3, lty=2)
lines(poblacion.cort1t2, res.cord1t2[,3], col="red", lwd=3, lty=2)
 
par(mfrow=c(1,1))
plot(poblacion.cort1t2, res.sigt1t2, type="l", col="black", lwd=3,
     xlab="Correlación poblacional", ylab="Proporción de ps < 0.05", main="Correlaciones estadisticamente significativas")
lines(poblacion.cort1t2, res.sigd1t2, col="black", lty=2, lwd=3)
abline(h=0.05, lty="dotted", col="blue", lwd=3)
abline(h=0.80, lty="dotted", col="blue", lwd=3)
 
par(mfrow=c(1,1))
plot(poblacion.cort1t2, poblacion.cort1t2, type="l", col="orange", lwd=5,
     xlab="Correlación poblacional", ylab="Correlación muestral | p < 0.05", main="Correlaciones muestrales de los resultados significativos")
lines(poblacion.cort1t2, res.corsigt1t2, col="black", lwd=3)
lines(poblacion.cort1t2, res.corsigd1t2, col="black", lwd=3, lty=2)
Created by Pretty R at inside-R.org

No hay comentarios:

Publicar un comentario