Modelos binomiales en R

Muchos problemas estadísticos implican el análisis de una respuesta binaria. Por ejemplo, muchas veces se clasifican los individuos muestrales como:

• Vivos o muertos,
• Presencia o ausencia de especies,
• Hombre o mujer,
• Sano o enfermo,
• Deforestado o no deforestado.

En estos y otros casos, los modelos binarios suponen una herramienta de análisis muy interesante. Para ello es conveniente que, al menos, una de las variables explicativas sea continua. En el caso de la deforestación (o cualquier otro cambio de uso binario), los datos consistirán, por ejemplo, en una lista de variables (ej. distancia a caminos, elevación, pendiente) medidas en píxeles que han sido deforestados en un periodo determinado, y una lista similar para píxeles que no han sido deforestados en ese mismo periodo. La cuestión es tratar de averiguar cuál de estas variables, si es que hay alguna, aumentan la probabilidad de que un píxel individual sea deforestado para ese periodo.

La variable respuesta contiene sólo 0s y 1s. Por ejemplo, 0 representaría píxeles que permanecen como bosque y 1 píxeles que han sido deforestados. Por tanto, habría una única columna para la variable respuesta. Los pasos a seguir para modelar una variable respuesta de tipo binomial son los siguientes:

1. Crear un vector de 0s y 1s como variable respuesta.
2. Ajustar un modelo lineal generalizado (funcion glm) con una familia de distribución de errores de tipo binomial (argumento family = binomial).
3. Reducir el modelo eliminando los términos que no sean significativos y comparando el cambio en devianza con respecto al modelo total.
4. Interpretación de los resultados.
   

La función de vínculo en los glm hace referencia a la transformación de la(s) variable(s) antes de ajustar el modelo. En el caso de los modelos binomiales, la función de vínculo que se usa habitualmente es la logística (logit), que transforma con logaritmos naturales la variable respuesta. Por defecto, no hace falta especificar ésto en los modelos.

A la hora de comparar modelos se utiliza la devianza como un indicador de la bonanza del modelo. Para aquellos que no estén muy familiarizados con la jerga estadística, la devianza es una medida del ajuste del modelo en relación al número de parámetros del mismo. Si un modelo explica lo mismo que otro, pero es más sencillo (menor número de párametros) es, por regla general, mejor modelo. Y esto es lo que va a indicar la devianza.

Un ejemplo: El caso de la deforestación en Chile Central

En este ejemplo, la variable respuesta va a ser la ocurrencia de deforestación en Chile Central para el periodo 1975-1985, en dónde 1 = deforestación y 0 = no deforestación. Para obtener estos datos, se llevó a cabo un muestreo sistemático cada 1 km en el área de estudio y se registró la cobertura en 1975 y 1985 para todos aquellos píxeles que eran bosque en 1975. Se midieron también toda una serie de variables explicativas que se detallan a continuación:

• Elevación (m),
• Pendiente (grados),
• Orientación (grados),
• Distancia a ríos (m),
• Distancia a ciudades grandes (m),
• Distancia a pueblos (m),
• Distancia a carreteras principales (m),
• Distancia a carreteras secundarias (m),
• Distancia a suelo agrícola (m),
  
• Distancia al borde del bosque (m).

Los datos originales están en Excel. Para leerlos directamente en R suele ser más rápido pasarlos a *.txt, aunque R puede leer también datos en Excel. Los datos están accesibles en una página web, por lo que se pueden leer desde R utilizando el comando url(). Basta con copiar el código que se ve en azul y pegarlo en R. Funcionará.

fnf75_85 <- read.table(url("http://archivos-para-subir.googlegroups.com/web/F-NF1975-85.txt?gda=f500bkEAAACkOfwqHAVd4YqgfIB09GDRCcOUO-z8uvHRarwVOPM5xOEDwyLfg5d1mmtcx5AUJ3JTCT_pCLcFTwcI3Sro5jAzlXFeCn-cdYleF-vtiGpWAA"), header = T, dec = ".", sep = "\t")

Ya se ha creado un objeto en R llamado fnf75_85, que contiene los datos de las distintas variables que se utilizarán en el análisis. Con el argumento header = T se especifica que la primera línea contiene los nombres de las variables. Con el argumento dec = "." se dice cómo vienen especificados los decimales (alternativamente podría ser ","). Por último, el argumento sep = "\t" especifica que los campos están separados por tabulaciones. Con la función str() se obtiene la estructura del objeto fnf75_85 e información sobre la naturaleza (integer, factor, character, etc) de cada una de las variables.

str(fnf75_85)

'data.frame': 1980 obs. of 15 variables:
$ ID_ALL : int 10 20 77 80 81 118 158 208 255 262 ...
$ ID_FOREST : int 1 2 3 4 5 6 7 8 9 10 ...
$ Class1975 : int 1 1 1 1 1 1 1 1 1 1 ...
$ Class1985 : Factor w/ 8 levels "1","2","3","5",..: 4 3 2 2 2 3 1 1 2 1 ...
$ Change : Factor w/ 3 levels "0","1","NoData": 2 2 2 2 2 2 1 1 2 1 ...
$ Elevation : num 435 136 299 300 300 ...
$ Slope : num 19.6 0.4 0.2 0 0 0 0.2 12.3 0.3 3.8 ...
$ Aspect : num 149 290 357 -1 -1 ...
$ Dist_river : num 4536 1025 751 2005 2672 ...
$ Dist_City.20T: num 27750 52254 28529 25994 25180 ...
$ Dist_village : num 2271 2620 1142 1084 742 ...
$ Dist_road_P : num 4601 2313 4621 1973 1442 ...
$ Dist_road_S : num 752 395 323 552 0 ...
$ Dist_agri : num 94.9 60.0 30.0 30.0 30.0 ...
$ Dist_f_forest: num 84.9 524.8 67.1 161.6 67.1 ...

La variable respuesta es $Change. Se observa que es un factor con 3 niveles: "0", "1" y "NoData". No interesa tener muestras con "NoData", por lo que se deben de eliminar estos datos. Se puede hacer sobreescribiendo el archivo fnf75_85 (! indica exclusión de casos).

fnf75_85 <- fnf75_85[!fnf75_85$Change == "NoData", ]
fnf75_85$Change <- as.integer(fnf75_85$Change)
fnf75_85$Change <- ifelse(fnf75_85$Change == 1, 0, 1)

Y ahora un resumen de la variable respuesta.

table(fnf75_85$Change)

0 1
774 1184

Así que hay 774 píxeles que continúan siendo bosque en 1985 y 1184 píxeles que han sido deforestados entre 1975 y 1985. El diseño es claramente no balanceado (distinto número de casos en la variable respuesta), aunque esto no tiene por qué influir necesariamente en el ajuste del modelo. El ajuste del modelo completo (excluyendo la variable Dist_f_forest, que ha sido mal calculada) se haría de la siguiente forma:

model.full <- glm(Change ~ Elevation + Slope +Aspect + Dist_river + Dist_City.20T + Dist_village + Dist_road_P + Dist_road_S + Dist_agri, data = fnf75_85, family = binomial)

Para explorar los resultados se puede usar la función summary().

summary(model.full)

Call:
glm(formula = Change ~ Elevation + Slope + Aspect + Dist_river +
Dist_City.20T + Dist_village + Dist_road_P + Dist_road_S +
Dist_agri, family = binomial, data = fnf75_85)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.5265 -1.3419 0.9214 1.0039 1.3176

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.177e-01 1.650e-01 3.742 0.000182 ***
Elevation 3.529e-04 2.209e-04 1.597 0.110201
Slope 3.913e-03 4.902e-03 0.798 0.424748
Aspect -8.098e-04 4.641e-04 -1.745 0.081044 .
Dist_river 1.350e-05 6.302e-05 0.214 0.830348
Dist_City.20T -3.136e-06 4.025e-06 -0.779 0.435960
Dist_village 2.353e-05 3.188e-05 0.738 0.460615
Dist_road_P -2.579e-05 9.662e-06 -2.669 0.007608 **
Dist_road_S -8.984e-05 6.592e-05 -1.363 0.172909
Dist_agri -4.340e-05 5.029e-05 -0.863 0.388227
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2627.9 on 1957 degrees of freedom
Residual deviance: 2607.7 on 1948 degrees of freedom
AIC: 2627.7

Number of Fisher Scoring iterations: 4

Cómo se observa, las variables que influyen significativamente sobre la deforestación en el periodo 1975-1985 son la distancia a carreteras principales (Dist_road_P) y, de manera más marginal (p = 0.081), la orientación (Aspect). La primera variable influye negativamente sobre la probabilidad de deforestacion (Estimate = -2.579 e-05), es decir, que a medida que nos alejamos de las carreteras principales la probabilidad de deforestación disminuye. La segunda es más difícil de intepretar, puesto que la orientación va de 0º (Norte) a 360º (Norte). En este caso, sería más interpretable utilizar una variable como insolación, que está directamente relacionada con la orientación y la pendiente.

Puesto que el resto de variables no son significativas se puede reducir el modelo.

model.red1 <- glm(Change ~ Aspect + Dist_road_P, data = fnf75_85, family = binomial)
summary(model.red1)

Call: glm(formula = Change ~ Aspect + Dist_road_P, family = binomial, data = fnf75_85) Deviance Residuals: Min 1Q Median 3Q Max -1.5055 -1.3468 0.9385 1.0064 1.2628 Coefficients: Estimate Std. Error z value Pr(>|z|)

(Intercept) 7.514e-01 1.092e-01 6.881 5.92e-12 ***
Aspect -7.502e-04 4.593e-04 -1.633 0.10236
Dist_road_P -2.496e-05 8.126e-06 -3.072 0.00213 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2627.9 on 1957 degrees of freedom
Residual deviance: 2615.3 on 1955 degrees of freedom
AIC: 2621.3

Number of Fisher Scoring iterations: 4

Y se comparan los dos modelos con la función anova().

anova(model.full, model.red1, test = "Chi")

Analysis of Deviance Table

Model 1: Change ~ Elevation + Slope + Aspect + Dist_river + Dist_City.20T +
Dist_village + Dist_road_P + Dist_road_S + Dist_agri
Model 2: Change ~ Aspect + Dist_road_P
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 1948 2607.66
2 1955 2615.32 -7 -7.66 0.36

Lo que indica que el modelo más simple no es significativamente peor que el modelo completo (p = 0.36). Así que aceptamos este modelo por el momento. Sin embargo, cuando analizamos los resultados se observa que la variable orientación no es significativa (p = 0.102), por lo que probamos a reducir el modelo un poco más.

model.red2 <- glm(Change ~ Dist_road_P, data = fnf75_85, family = binomial)
summary(model.red2)

Call: glm(formula = Change ~ Dist_road_P, family = binomial, data = fnf75_85) Deviance Residuals: Min 1Q Median 3Q Max -1.4540 -1.3500 0.9432 1.0034 1.2042 Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.333e-01 8.129e-02 7.791 6.66e-15 ***
Dist_road_P -2.551e-05 8.115e-06 -3.144 0.00167 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2627.9 on 1957 degrees of freedom
Residual deviance: 2618.0 on 1956 degrees of freedom
AIC: 2622

Number of Fisher Scoring iterations: 4

Y volvemos a comparar los dos modelos:

anova(model.red1, model.red2, test = "Chi")

Analysis of Deviance Table

Model 1: Change ~ Aspect + Dist_road_P
Model 2: Change ~ Dist_road_P
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 1955 2615.32
2 1956 2617.99 -1 -2.67 0.10

Así que de nuevo nos encontramos que el modelo reducido (model.red2) no es significativamente peor que el modelo anterior (model.red1). Aceptamos este último como válido, por lo que concluimos que la probabilidad de deforestación en el periodo 1975-1985 estuvo influida fundamentalmente por la presencia de carreteras principales (¿se debió la deforestación a un uso maderero por parte de empresas forestales?).

Por último, se representan gráficamente los resultados.

xv<-seq(0,30000,10)
yv<-predict(model.red2, list(Dist_road_P = xv), type = "response")
plot(as.integer(fnf75_85$Change) ~ fnf75_85$Dist_road_P, type = "p", ylab = "Probability of forest loss", xlab = "Distance to main roads (m)")
lines(xv,yv)

¿Qué nos dice este gráfico? Básicamente, que la probabilidad de deforestación disminuye con la distancia a carretera principales, pero disminuye poco. Así, cuando uno se aleja 25 kms de las carreteras principales (una distancia bastante grande), la probabilidad de deforestación disminuye del 65% a apenas el 50%.

Cómo hay muchos puntos en el gráfico, no se aprecia realmente si hay más 0s que 1s en distancias próximas a las carreteras principales y a la inversa... sería bueno intentar ajustar funciones de densidad de los 0s y 1s en el mismo gráfico.

¿Qué ocurre si la variable respuesta es una proporción?

Puede ser que la variable respuesta no sean 0s y 1s sino que sea una proporción o porcentaje, ej. porcentaje de semillas que sobreviven en un experimento, sex ratios, proporción de píxeles deforestados en un determinado área, etc. En este caso, tambińe podemos llevar a cabo un análisis glm con distribución de errrores de tipo binomial. Las únicas diferencias con respecto a los pasos descritos anteriormente serían:

• Para definir la variable respuesta creamos un objeto con dos columnas, utilizando la función cbind() para juntar dos vectores que contengan el número de éxitos y fracasos en cada uno de los casos.

n.fracasos <- denominador - n.exitos
var.respuesta <- cbind(n.exitos, n.fracasos)

• Comprobar que no hay sobredispersión (la devianza residual no es mayor que los grados de libertad residuales) y corregirlo en caso necesario utilizando una familia de tipo cuasibinomial (family = quasibinomial).
  

Bibliografía recomendada:

- Crawley, M.J. 2007. The R Book. Wiley, London.