Sección 4 Predicción con múltiples capas de datos

La predicción del valor de una variable en un sitio no observado, realizada por interpolación kriging local, se alimenta de datos de la misma variable que se quiere predecir, pero obtenidos en sitios del vecindario de aquel para el que se requiere la predicción. Aún cuando la predicción pueda ser óptima, resulta informativo conocer el impacto de otras variables (variables secundarias, auxiliares o explicativas) sobre la predicción.

El método conocido como regresión kriging (Hengl, Heuvelink, and Rossiter 2007) permite predecir o interpolar los valores de una variable en sitios no muestreados que se encuentran entre los sitios con observaciones a partir de capas de datos secundarios que actúan como variables predictoras. La incorporación de covariables de sitio puede mejor sustancialmente la predicción de una variable respuesta espacializada. Para tal fin es necesario el ajuste de un modelo de regresión múltiple que describe la relación entre la variable observada y las capas de variables predictoras. Luego de ajustar la regresión, se obtienen los residuos (diferencia entre valor observado y valor predicho por el modelo a partir de las variables secundarias) y se realiza una interpolación kriging sobre los residuos para contemplar la espacialidad que no se encuentra relacionada con las capas de información usadas en el ajuste previo. Finalmente, los resultados del ajuste de regresión y de la interpolación kriging son combinados para producir la predicción.

Matemáticamente, esta interpolación es equivalente a la que hemos llamado kriging universal, ya que allí se realiza un ajuste de la variable de interés y las coordenadas espaciales (un tipo de variable auxiliar) cuyo impacto es descontado antes de realizar la interpolación espacial. Kriging desde modelo de regresión es también equivalente matemáticamente al método de interpolación conocido como kriging con deriva externa donde múltiples variables auxiliares son usadas como predictoras (covariables de sitio que son externas o distintas a las coordenadas espaciales) e intervienen directamente para resolver los pesos de la interpolación espacial.

La predicción sitio-especifica de una variable respuesta en función de covariables de sitio puede realizarse también desde el marco teórico de los modelos lineales mixtos de covarianza residual (Balzarini, Macchiavelli, and Casanoves 2004). Es decir, ajustando un único modelo, donde la componente determinística es un modelo de regresión lineal múltiple que relaciona la respuesta con las covariables de sitio y las componentes sistemáticas (patrón espacial) y estocástica (variación sin estructura) de la variación espacial son expresadas conjuntamente en la matriz de varianza y covarianza de los términos de error (matriz de covarianza residual). Bajo este marco teórico, se utiliza máximo verosimilitud restringida o REML (Patterson and Thompson 1971) para estimar el modelo predictivo. El modelo relaciona la variable de interés para el estudio de variación espacial con las variables explicativas también distribuidas espacialmente y contemple la dependencia espacial en las observaciones de la variable a predecir incorporando en la matriz de covarianza residual covariaciones derivadas de una función basada en la distancia que separa dichas observaciones en el espacio.

La construcción de un modelo predictivo para datos espaciales también puede realizarse desde un marco teórico Bayesiano, donde sistemática relacionada con la media de la respuesta sigue siendo el modelo de regresión lineal múltiple, pero se incorporan parámetros aleatorios de sitio para modelar la componente sistemática de los datos espaciales, los cuales se predicen independientemente de la componente de error. Trabajar el efecto de sitio como aleatorio, junto a la suposición de procesos espaciales locales, ha permitido eficientizar la estimación en un único paso de modelos para datos espaciales. Bajo el marco teórico bayesiano y el supuesto de campo aleatorio gaussiano markoviano para el proceso espacial, la estimación por aproximación de Laplace anidada conocida como INLA (del término en inglés Integrated Nested Laplace Aproximation) es uno de los métodos elegidos para ajustar el modelo espacial (Rue, Martino, and Chopin 2009).

Alternativamente, la construcción del modelo predictivo espacial puede realizarse mediante métodos de aprendizaje automático, también bajo el concepto de ajustar un modelo para la componente determinística, obtener los residuos y modelar el proceso espacial en los residuos para finalmente combinar ambos resultados (Li et al. 2011). Los modelos basados en árboles de regresión (Breiman et al. 2017) son particularmente útiles para abordar la primera etapa, aquella donde la variable respuesta es relacionado con las múltiples covariables que definen capas de información espacial.

Cualquiera sea la aproximación usada en la construcción del modelo para realizar la predicción espacial y posterior mapeo de la variable respuesta, para lograr buenos resultados se requiere hacer uso del criterio experto en la disciplina, realizar un adecuado preprocesamiento de los datos, identificar los predictores influyentes y validar el modelo estimado (Kuhn and Johnson 2013).

La validación cruzada es una técnica de amplio uso para evaluar el desempeño de modelos predictivos con múltiples capas de información. La forma más adecuada de muestrear para dividir el conjunto de datos en estudio y los criterios de validación es discutida (Brenning 2012; Efron and Tibshirani 1997), pero como regla general se puede establecer que en cuanto menor es la proporción de datos que conforma el subconjunto de validación, menor capacidad de extrapolación tendrá el modelo ajustado y que éste no puede crecer indefinidamente ya que también se requiere buen tamaño del subconjunto de calibración para ajustar un buen modelo. El compromiso se relaja a medida que la base de datos observados es de mayor tamaño.

Una medida comúnmente usada en la evaluación de modelos de regresión lineal estimados por mínimos cuadrados es el \(R^2\) o coeficiente de determinación. Éste expresa la bondad del ajuste del modelo respecto a los datos observados más que la capacidad predictiva que interesa evaluar cuando el modelo es usado como predictor. Otras medidas de bondad de ajuste, como los criterios de información de Akaike y BIC usados para evaluar modelos lineales de covarianza residual estimados por métodos basados en la verosimilitud o DIC usado en contextos de ajustes de modelos bayesianos, permite calificar el ajuste del modelo a los datos observados, pero no su desempeño para predecir observaciones no realizadas. La validación cruzada es, por el contrario, una estrategia transversal a distintas aproximaciones metodológicas o marcos teóricos que provee información de la capacidad predictiva del modelo, y por tanto una medida de cómo podría funcionar en la interpolación espacial. No obstante, es de resaltar que un buen modelo predictivo debería representar también un buen ajuste de los datos de entrenamiento, y por ello métricas comunes de bondad de ajuste suelen ser usadas para una primera selección de modelos predictivos.

4.1 Regresión con errores correlacionados espacialmente vía REML

Se asume una relación lineal determinística entre la variable respuesta y las covariables espacializadas que se modela a través de coeficientes de regresión, y el proceso espacial subyacente en los datos se modela sobre los términos de error del modelo. Éstos, en lugar de considerarse independientes como en aproximaciones estadísticas clásicas, se suponen espacialmente correlacionados, i.e. con correlaciones entre pares de términos de error expresadas como función de la distancia que separa los sitios desde los que se obtuvieron las observaciones. La varianza residual del modelo mide la variabilidad no estructurada espacialmente y caracteriza la componente estocástica del modelo. Se ajusta un modelo directamente a los datos, y no a las semivarianzas como en la geoestadística clásica. La estructura de covariación espacial se define en función de la distancia entre las observaciones al igual que en la geoestadística, pero se estima simultáneamente con los parámetros del modelo de regresión, i.e. los coeficientes de regresión que interpretan como pendientes o cambios en la respuesta por unidad de cambio de la covariable de sitio.

El modelo de regresión lineal múltiple con errores correlacionados espacialmente para una variable \(Y\) asume la siguiente distribución para la i-ésima observación:

\[Y_i\sim N \big(E(Y),V(Y) \big)\]

\[Y_i=\beta_0+\sum_{j=1}^{p}{x_{ij}\beta_j}+\varepsilon_i\]

\[\varepsilon(s_i)\sim N(0,\sigma^2)\]

\[Cov(\varepsilon_i,\varepsilon_j)=\left\{\ \begin{array}{lcc} \sigma_s^2+\sigma_e^2\ &\ si\ &\ i=j \\\sigma_e^2\ f(d_{ij})\ &\ si\ &\ i\ \neq\ j \end{array} \right.\]

donde \(\beta_0\) es el intercepto; \(\beta_j\) es el vector de efectos fijos para las variables regresoras \(x_j\); \(x_{ij}\) la valuación de la covariable \(x_j\) en el sitio \(i\) y \(\varepsilon_i\) es el término de error que se asume distribuido Normal con media 0 y varianza \(\sigma^2\). \(Cov(\varepsilon_i,\varepsilon_j)\) es la covarianza de los errores de los sitios \(i\) y \(j\) determinada a partir de una función de covarianza espacial dependiente de la distancia \(d_{ij}\) entre observaciones.

La estimación REML es la elegida ya que ha demostrado que reduce el sesgo en las estimaciones de los parámetros de covarianza (Morrell 1998). Si el proceso de ajuste se realiza en etapas, se recomienda postular un modelo saturado o con máximo número de covariables en primera instancia y sobre los residuos de este modelo, vía REML, estimar las varianzas y covarianzas en los datos. Seleccionado el modelo para la matriz de covarianza residual, se intentará reducir la componente determinística y en este momento, cuando se evalúan los coeficientes de regresión con un modelo de varianza-covarianza adecuado, se pude utilizar estimaciones máximum likelihood clásica (ML). El método de REML, ajusta los grados de libertad de los efectos fijos (estructura de medias) antes de estimar los componentes de varianza y por ello es preferido a ML al momento de identificar la matriz de covarianza residual más apropiada para modelar los términos de error. La correlación espacial en los errores podría explicarse con una función de correlación espacial exponencial, gaussiana, esférica o lineal, entre otras. Ellas expresan como la correlación entre dos observaciones decae con la interdistancia (usualmente Euclídea) entre los sitios desde los cuales se obtienen y consecuentemente la selección de un modelo de correlación espacial es equivalente a la selección de un modelo teórico para el ajuste de un semivariograma experimental de la geoestadística clásica. Los modelos de correlación espacial pueden contener o no el efecto nugget, i.e. una estimación de la varianza a una escala más fina que la de la grilla entre observaciones. Finalmente, el modelo para la matriz de covarianza residual puede ser homocedástico (i.e. con varianza residual única) o heterocedástico (i.e. con varianza residual diferente para distintos subgrupos de datos). Comúnmente se evalúan varios modelos alternativos para la matriz de covarianza residual y se selecciona aquel para cual el ajuste del modelo tenga un menor valor para los criterios de información penalizada como AIC y BIC.

4.2 Regresión con efectos aleatorios de sitio vía INLA

Los modelos lineales de covarianza residual pueden no ser eficientes para modelar procesos espaciales continuos que involucran matriz de covarianzas grandes y densas, por ejemplo, el modelado de correlaciones espaciales en un conjunto de datos con más de 10000 observaciones es computacionalmente intensivo y podría no ser logrado en computadores de escritorio actuales. Nuevas implementaciones de la regresión lineal múltiple para datos espaciales, desarrolladas en el marco teórico bayesiano, facilitan esta estimación

En estadística bayesiana se considera que los parámetros del modelo son variables aleatorias y se calculan distribuciones de probabilidad para los parámetros de las cuales se deriva medidas de incertidumbre (Correa Morales, Causil, and Javier 2018). La información previa sobre los parámetros debe resumirse en distribuciones de probabilidad denominadas distribuciones a priori, a partir de las cuales se estima la distribución de probabilidad a posteriori dadas las observaciones. Estimaciones puntuales de los parámetros de interés se pueden obtener calculando medidas resumen de la distribución a posteriori, como la media o el modo, y se informan juntos a intervalos de credibilidad calculados desde percentiles de la distribución a posteriori. La credibilidad se interpreta como la probabilidad de que el valor estimado para el parámetro pertenezca al intervalo reportado, dado los datos observados.

Los métodos de simulación por cadenas de Markov Monte Carlo (MCMC), han permitido resolver modelos complejos sin la necesidad de imponer estructuras que lo simplifiquen. Éstos han sido usados para la estimación de modelos con datos espaciales. Sin embargo, el método MCMC también conlleva desafíos computacionales. Rue, Martino, and Chopin (2009) propusieron una alternativa para obtener la distribución a posteriori en contextos de datos espaciales a partir de aproximaciones basadas en el algoritmo INLA que bajo el supuesto de que la variación espacial subyacente se describe como un campo aleatorio gaussiano Markoviano simplifica las estimaciones de covarianzas espaciales en procesos espaciales continuos con una gran cantidad de observaciones. Sobre la base de las aproximaciones por INLA y la implementación de la alternativa en el lenguaje de programación R (R-INLA) se han popularizado las aplicaciones de la regresión bayesiana a datos espacial y espaciotemporales (Cameletti et al. 2013). En este contexto es posible obtener la matriz de precisión de los efectos aleatorio de sitio utilizando aproximaciones por ecuaciones diferenciales parciales estocásticas (spde) (Lindgren and Rue 2015). Bajo este enfoque se construye una malla de predicción con unidades de celda triangulares que cubre todo el dominio espacial para el cual se requiere la predicción, cada vértice de los triángulos representa un nodo sobre el que se predice la variable respuesta por interpolación (Blangiardo and Cameletti 2015). Es posible trabajar con dominios espaciales de límites y bordes complejos (Bakka et al. 2018) y asignar medidas de incertidumbre de cada predicción puntual ya que lo que se obtiene del modelo bayesiano es la distribución a posteriori de los valores predichos para cada sitio más que un único valor de predicción.

Bajo este enfoque el modelo de regresión bayesiana para una variable \(Y\) asume la siguiente distribución para la i-ésima observación:

\[Y_i\sim N(\eta_i,\sigma_e^2)\]

\[\eta_i = \beta_0 + \sum_{j-1}^{p} {x_{ij}\beta_j}+ \xi(s_i)\]

donde \(\beta_0\) es el intercepto; \(\beta_j\) es el vector de efectos fijos de las variables explicativas \(x_j\); \(x_{ij}\) la valuación de la covariable \(x_j\) en el sitio \(i\) y \(\xi(s_i)\) el efecto aleatorio de sitio que se asume una realización de un proceso gausiano markoviano latente \(\xi(s_i) \sim MVN(0,\Sigma)\), siendo \(\Sigma\) la matriz de varianza y covarianza de los efectos de sitio definidos por la función de covariación espacial de Matérn (Matérn, 1986). En R-INLA la estimación de la inversa de \(\Sigma\) (matriz de precisión) se resuelve eficientemente usando spde.

Para estimar el modelo, es necesario obtener una representación de la estructura de dependencia desde una estructura de vecindario para datos continuos (malla). La malla se obtiene mediante triangulación de Delaunay restringida comenzando sobre la estructura base correspondiente a los vértices iniciales de las observaciones, luego se agregan o eliminan vértices adicionales para satisfacer las restricciones de la triangulación que se encuentran definidas por los siguientes parámetros: 1) offset define hasta qué punto se debe extender la malla hacia lo interno (es decir, dentro del área donde se pretende predecir) y hacia el exterior (es decir, fuera del área donde se pretende predecir); 2) cutoff define la distancia mínima entre vértices permitida; 3) maximumedge que refiere a la longitud máxima del borde de cada triángulo; 4) minimumangle o ángulo mínimo interior de cada triángulo. Construida la malla es necesario seleccionar un modelo de correlación espacial. En R-INLA con SPDE esta función se define parametrizando la función de correlación espacial de Matérn definiendo su parámetro alpha (variando entre 0 y 2). El valor por defecto de alpha es 2 y corresponde a una función de correlación espacial del tipo exponencial. La matriz de covarianza de los efectos aleatorios es una matriz rala dado el proceso espacial que se supone, y sus valores son aproximados por suavizado vía spde.

4.3 Regresión vía modelos basados en árbol

El término aprendizaje de máquina o aprendizaje automático corresponde a una rama de la inteligencia artificial que hace referencia a algoritmos o procedimientos de cálculo basados en intenso proceso computacional que “aprenden” de los datos intentando minimizar la intervención humana. Algunos se basan en particiones binarias recurrentes de los datos y evaluaciones de éstas hasta identificar la mejor para explicar variabilidad en la variable respuesta. Los árboles de clasificación y regresión (algoritmos CART) (Breiman 2001) se utilizan con fines predictivos y son particularmente útiles para interpretar relaciones (no necesariamente lineales) en contextos de regresión múltiple con variables explicativas correlacionadas. Estos algoritmos pueden ser empoderados mediante métodos de remuestreo que se usan para obtener muestras aleatorias a partir de los datos observados o muestra original, derivar modelos para cada muestra y ensamblar los resultados para optimizar la predicción. Otros algoritmos de aprendizaje automático basados en árboles son la regresiones por bosques aleatorios (RF por el término en inglés Random Forest) y los árboles de regresión generalizados (GBR por el término en inglés Generalized Boosting Regression) (Efron and Hastie 2016). Los algoritmos basados en árboles engloban así a un conjunto de técnicas supervisadas no paramétricas (i.e. sin supuestos distribucionales) para segmentar el espacio de los predictores en regiones simples con máxima diferencia en la variable respuesta. Es necesario tener cuidados particulares al momento de estimar el árbol que será usados como modelo predictivo ya que puede existir sobreajuste, i.e. construirse un modelo sólo útil para los datos disponibles cuyas predicciones pueden cambiar con pequeños cambios en el conjunto de datos observados.

Si bien estos algoritmos se han utilizados para datos espaciales (Kanevski et al. 2009), es raro que se modele la estructura espacial. Una propuesta para incorporar la correlación espacial en los datos es utilizar las coordenadas geográficas o las matrices de distancias entre observaciones covariables en la construcción del modelo (Pejović et al. 2018). Otra propuesta, es modelar el residuo remanente del ajuste del algoritmo de aprendizaje automático con una función de autocorrelación espacial (Li et al. 2011) y finalmente combinar los resultados de la predicción determinística dada por el árbol y la predicción espacial obtenida mediante kriging de los residuos.

4.3.1 Bosques aleatorios

El método de bosques aleatorios o Random Forest (RF) es una modificación del proceso de ensamblaje de varios árboles (Bagging) donde se ajustan múltiples árboles desde cada muestra obtenida por remuestreo formando un “bosque”. En cada nueva predicción, todos los árboles que forman el “bosque” participan aportando su predicción. Como valor final, se toma la media de todas las predicciones en el caso de variables respuesta continuas. El método RF a diferencia de Bagging realiza una selección aleatoria de \(m\) predictores antes de evaluar cada división. Si \(m=p\) los resultados de RF y Bagging son equivalentes. Este método trabaja bien con grandes bases de datos presentando mayor facilidad en la implementación y baja tendencia al sobreajuste. Para implementar RF es necesario optimizar el parámetro \(m\), no obstante, existe la recomendación de usar \(m=\frac{p}{3}\) por defecto para regresión.

Para contemplar la estructura espacial es posible combinar los resultados de RF con una interpolación geoestadística basada en kriging (Li et al. 2011). La predicción de los residuos se complementa con la predicción de RF de manera aditiva. Es recomendable que el ajuste del modelo espacial para realizar kriging se logre con un subconjunto de datos de entrenamiento, diferente al grupo de validación, para evitar sobreajustes.

4.3.2 Árboles de regresión generalizados

El método Boosting es otro método de ensamblaje que consiste en ajustar secuencialmente modelos sencillos, de manera que cada modelo aprende de los errores del anterior. Los algoritmos de Boosting trabajan minimizando una función de pérdida (deviance) para maximizar la proporción de varianza que explica el modelo. Los árboles de regresión generalizados, conocido en inglés como Generalized Boosted Regression Trees (GBR) particularmente utiliza Boosting para ensamblar los árboles obtenidos de múltiples muestras obtenida mediante remuestreo de la muestra original. El algoritmo ajusta árboles de regresión a los datos de entrenamiento de manera iterativa. El primer árbol que se ajusta es aquel que, según la complejidad del árbol seleccionada, minimiza la deviance. El siguiente árbol se ajusta a los residuos del primer árbol. Luego, se vuelven a realizar predicciones para las observaciones que tienen en cuenta los dos árboles. En cada uno de los pasos siguientes se ajusta un nuevo árbol sobre los residuos de la combinación de los árboles anteriores. El procedimiento es parametrizado por varias constantes que es necesario identificar probando numerosas o todas las combinaciones posibles de parámetros ya que estos son dependientes entre sí.

Una vez encontrada la combinación optima de parámetros se ajusta GBR en el grupo de entrenamiento. Luego, a partir de los residuos del modelo se ajusta un kriging y se guardan los parámetros de la función de semivarianza ajustados. Finalmente se utiliza el modelo GBR construido para predecir los datos en el grupo de validación adicionando a los resultados obtenidos desde el árbol la predicción de los residuos de cada sitio.