Pasó que corría Maxent en R con todas las variables bioclim que encontraba, esperaba un mapa con manchas de bosque donde correspondía, y me salía el país entero pintado de rojo. Una tarde de otoño en Valdivia, llueve que llueve, cuando intentaba mi primer modelado de nicho con el bosque nativo de por aquí como excusa. El ventilador del portátil empezó a soplar fuerte, como si le hubiera pedido una maratón, y todo ese esfuerzo terminó en una mancha uniforme que no distinguía un bosque de una vereda. Ahí me quedé pensando si esto del R espacial y la teledetección para principiantes era para mí, o si estaba haciendo algo básico mal sin enterarme.
La respuesta corta, por si vienes con apuro: casi nunca hacen falta las diecinueve variables. Se eligen unas pocas que tengan sentido biológico para la especie, se mira cómo se llevan entre sí con una matriz de correlación y se saca cualquier variable que pase de un VIF de 10. Todo lo demás es afinar (y equivocarse un rato, que a mí me tomó varios fines de semana).
Un aviso antes de seguir: algunos enlaces de esta bitácora son de afiliado. Si compras algo por ahí me llega una comisión chiquitita y a ti no te cuesta nada más (con lo que cuesta llegar a fin de mes, mejor dejarlo claro). Solo aparecen clases que abrí yo misma cuando me trabé con algo.
Diecinueve variables bioclim en la licuadora
Mi primer error fue creer que mientras más datos le diera al modelo, mejor iba a salir. Me bajé las variables bioclimáticas de WorldClim, las que casi todo el mundo usa para esto, vi que venían diecinueve y metí las diecinueve al código. Ahí empezó el desastre. Varias de ellas dicen casi lo mismo (la temperatura media anual y la temperatura del trimestre más cálido, por ejemplo), y al pobre algoritmo eso lo marea muchísimo.
Casi todo el sábado lo pasé buscando al culpable, porque mis puntos de presencia estaban bien (que por cierto me costó un mundo dejarlos así, y lo conté en mi experiencia limpiando puntos de presencia para modelado de nicho r). Buscando ayuda seguí un tutorial de Python que daba por sentado que uno ya sabía programar, cosa que yo no, y entre funciones que no reconocía y errores que no entendía terminé más perdida que al principio. El mapa final seguía pareciendo una pintura abstracta sin pies ni cabeza.
La multicolinealidad y por qué importa el VIF
A mitad de invierno, con el frío metido en los huesos, me puse a estudiar en serio qué estaba pasando. Resulta que existe algo llamado multicolinealidad, que es básicamente cuando tus variables están tan emparentadas que el modelo no sabe cuál es la que de verdad importa. Es como si dos personas te gritaran lo mismo al mismo tiempo: al final no entiendes nada.
Boris, el vecino del edificio que trabaja en TI y no sabe nada de mapas, me pilló en el pasillo con el portátil abierto y me preguntó por qué el programa no podía elegir solo la variable que sirve. Tuve que explicárselo en voz alta, con las manos, y me di cuenta de que la mitad de mi explicación tenía huecos. Esa conversación me obligó a ordenar el asunto: primero ver quién se parece a quién, después decidir a quién sacar.
Para sacar a los redundantes se usa el factor de inflación de la varianza. El valor que más se ve es un umbral de Factor de Inflación de la Varianza (VIF) de 10: si una variable pasa de ahí, es redundante y hay que sacarla del baile. Me pasé varias tardes de domingo probando combinaciones, quitando la Bio1, luego la Bio12, y mirando cómo cambiaba todo. La lluvia de afuera no paró en toda la tarde, y esperar a que R terminara de calcular la matriz de correlación de mis rasters se volvió la banda sonora de mi invierno.
Denise Alarcón, del foro de R en español, me mandó por esos días el gráfico de calor de correlación que armó con ggplot2, y quedó tan limpio que me dieron ganas de rehacer el mío entero. Ella domina ese lado, pero se peleaba con el paso del raster al data frame, y yo con el VIF, así que terminamos intercambiando trucos por audio.
Recortar el área de estudio antes de modelar
Dejé de mirar los rasters como fotos bonitas cuando empecé a prepararlos: recortar exacto al área de estudio (cropping) para que los cálculos de Maxent no se pierdan en el infinito. El recorte con un shapefile da para su propia entrada de esta bitácora, así que acá va solo lo mínimo.
Fue en esa etapa, atascada entre el recorte y las capas que no calzaban, que abrí la primera clase del curso de Teledetección Satelital con R y R-Studio. Lo estuve pensando varias semanas porque el precio asusta cuando uno aprende por su cuenta, y las primeras clases se me hicieron densas por mi poca base de R. Aun así, ahí vi por qué la resolución espacial de Landsat 8/9 de 30 metros es tan distinta a la de los datos climáticos globales, y cómo hacer que ambos conversen sin que R lance uno de esos errores crípticos que te arruinan la tarde.
Por el lado del modelado, también abrí la clase de Modelado de Nicho Ecológico con Maxent en R, que es más corta y me sirvió como un punto de partida más liviano cuando la parte satelital me mareaba. Eso sí, es otra rama: no continúa lo de teledetección, y como tiene pocas reseñas fui con las expectativas bien medidas.
Un domingo lluvioso de la quinta semana con el proyecto, el script de R corrió de arriba a abajo y la consola quedó en puro texto negro, sin una sola línea roja; me quedé mirándola un rato, medio desconfiada, esperando que el error apareciera tarde.
Borrar variables también tiene su costo
Mientras limpiaba variables leí algo que me dejó pensando: eliminar variables muy correlacionadas puede reducir artificialmente la capacidad predictiva del modelo en escenarios de cambio climático. ¿Por qué? Porque en la naturaleza la colinealidad es una propiedad intrínseca, y si el clima cambia, esa relación entre la lluvia y la temperatura que hoy damos por sentada podría romperse.
Entonces, para un modelo de hoy el VIF ayuda a limpiar el ruido, pero para proyectar al futuro conviene ir con más cuidado y no borrar variables por puro reflejo. Es un equilibrio delicado. Un sábado en la mañana, caminando por el Mercado Fluvial de Valdivia con el piso mojado y los charcos entre los puestos, lo pensaba así: hay que saber dónde poner el pie para no resbalarse. Y mis scripts, la verdad, son un intento bien humilde de entender algo mucho más grande que un mapa.
Mi rutina para elegir variables en un modelado de nicho
Primero, no uso las diecinueve: elijo las que tienen sentido biológico para la especie (si es un anfibio, la humedad manda; si es un árbol resistente, quizás las heladas). Después armo la matriz de correlación antes de correr Maxent y miro cómo se llevan mis rasters entre sí; hay paquetes de R que dibujan un gráfico de calor precioso para eso. Cuando dos variables se parecen demasiado, me quedo con la que puedo explicarle a Boris sin dibujar nada, y si alguna pasa del VIF de 10, sale.
El área de estudio también cuenta: no modelo un continente si la especie vive en una sola cuenca, y si R se cae al cargar capas pesadas, te dejo la solución a errores comunes al cargar rasters pesados en R que a mí me sirvió. Y si el mapa vuelve a salir uniforme, lo primero que reviso es cuántas variables entraron y qué tan grande dejé el área, antes de tocar cualquier otra cosa. Suena a poco, pero esa revisión me habría ahorrado ese sábado entero.
No te va a salir a la primera, y está bien. Yo pasé fines de semana enteros de invierno probando y fallando hasta que los números empezaron a cuadrar; si me quedo con una sola cosa, es que pocas variables bien elegidas, que puedes explicar en voz alta, le ganan a un modelo gigante que nadie entiende.
Al final, el alivio de ver un mapa de idoneidad de hábitat que coincide con el bosque que conozco vale todos esos sábados perdidos. Es como si la computadora y el bosque se pusieran a hablar el mismo idioma. Las clases que fui abriendo en el camino me ayudaron a entender cómo crear mapas de idoneidad de hábitat en R para especies locales de una forma que hace sentido.
Cualquier duda que tengas, ando por aquí, probablemente con otro mate y esperando que termine de correr algún proceso pesado en RStudio. ¡Ánimo con esos mapas!