4.5 Consideraciones acerca de los métodos kriging

Los métodos kriging descritos anteriormente permiten obtener el mejor predictor lineal insesgado (BLUP, best linear unbiased predictor). Como es bien sabido, en el caso de normalidad el predictor óptimo (tomando como función de pérdidas el error cuadrático) es lineal y va a coincidir con los predictores kriging. Pero si el proceso no es normal no tiene porque serlo, lo que ha motivado el desarrollo de métodos kriging no lineales (ver e.g. Rivoirard, 1994) y del kriging trans-normal (ver Sección 4.7.2).

En estos métodos se supone que el variograma (covariograma) es conocido, sin embargo en la práctica en realidad el variograma es estimado (kriging estimado). Puede verse (Yakowitz y Szidarovszky, 1985) que, bajo condiciones muy generales, el predictor kriging con variograma estimado converge al valor correcto si la densidad de datos tiende a infinito, incluso cuando la estimación del variograma no es muy buena. Además, Stein (1988) probó que para una predicción asintóticamente (de relleno) eficiente lo que se necesita generalmente es capturar la conducta del variograma cerca del origen. Por tanto, en cuanto a las predicciones, el factor más importante es que la aproximación al variograma verdadero cerca del origen no sea muy mala.

Al contrario que en el caso del predictor, la estimación del variograma afecta directamente a la varianza kriging y en general no es un estimador consistente del error en media cuadrática del predictor con variograma estimado. En general, por ejemplo si el proceso es normal, es de esperar que subestime la verdadera varianza de predicción y por tanto debería incrementarse para que tenga en cuenta el efecto de la estimación del variograma. Para más detalles sobre este problema, ver por ejemplo Cressie (1993, Sección 5.3), Christensen (1991, Sección 6.5), Stein (1999, Sección 6.8) o Chilès y Delfiner (2012, Sección 3.4.3). En la Sección 4.6 se sugiere una posible corrección de la varianza kriging a partir del método de validación cruzada.

4.5.1 Kriging como interpolador

Las ecuaciones kriging tienen en cuenta varios aspectos del problema de interpolación22:

  • La configuración espacial de los datos, a través de las matrices \(\boldsymbol{\Sigma}\) (o \(\boldsymbol{\Gamma}\)), donde el covariograma (o el variograma) actúa como una “distancia estadística” entre las observaciones y de forma que se tiene en cuenta la información redundante presente en los datos.
  • La situación de la posición de predicción respecto a los datos, a través de \(\mathbf{c}\) (o \(\boldsymbol{\gamma}\)).
  • La presencia de una función determinística de tendencia.

Adicionalmente también tienen en cuenta propiedades estadísticas del proceso \(Z(\cdot)\), a través del variograma o el covariograma (que como se comentó en la Sección 1.3, entre otras cosas, determina las propiedades de continuidad del proceso). Esta es la principal diferencia con otros métodos de interpolación que no tienen en cuenta la estructura de segundo orden del proceso.

Una propiedad importante de los predictores kriging es que son interpoladores exactos (suponiendo que no hay error de medida23), en el sentido de que \(p(\mathbf{Z},\mathbf{s}_{0})=Z(\mathbf{s}_{0})\) cuando \(\mathbf{s}_{0}=\mathbf{s}_{i}\) (la solución de los sistemas es \(\lambda_{i} =1\) y \(\lambda_{j} =0\), \(\forall j\neq i\)), y naturalmente en ese caso la estimación de la varianza es 0. Además, por lo general no son continuos en las posiciones de los datos, para que lo sean el efecto nugget debe ser nulo.

Otra característica de la interpolación kriging (y que no aparece en otros métodos como los que asignan pesos inversamente proporcionales a la distancia) es el denominado efecto pantalla. En general se observa (ver comentarios en la siguiente sección) que la influencia de un valor es menor si está oculto detrás de otro valor (e.g. Journel y Huijbregts, 1978, p. 346). Esto produce, por ejemplo en el caso del KO (incluso con un variograma isotrópico), que puntos situados a la misma distancia de la posición de predicción puedan tener distintos pesos y que los datos cercanos no apantallados reciban los mayores pesos, reduciéndose considerablemente (llegando a ser negativos) los pesos de los datos que quedan ocultos24.

La aparición de pesos negativos (o mayores que 1) en el KO como consecuencia del efecto pantalla, puede provocar (incluso suponiendo media constante) que el predictor kriging no esté necesariamente comprendido entre el valor mínimo y máximo de los datos. Esto que en principio puede ser una propiedad muy interesante puede conducir a resultados extraños en ciertas ocasiones, como por ejemplo dar lugar a predicciones negativas en casos en los que la variable considerada es necesariamente positiva. Para solucionar estos problemas se han propuesto numerosas alternativas, entre ellas la inclusión de restricciones adicionales sobre los pesos de forma que sean positivos (e.g. Chilès y Delfiner, 2012, Sección 3.9.1) (lo cual puede dar lugar a un incremento considerable del MSPE de predicción) o sobre el predictor (ver p.e. Goovaerts, 1997, Sección 7.4.2). Otra alternativa que puede ser preferible es la transformación del proceso \(Z(\cdot)\) a otra escala (de forma que se aproxime a la normalidad), realizar la predicción kriging del proceso transformado y volver a la escala original (pero asegurándose de que al hacer la transformación inversa el resultado tenga las propiedades de optimalidad deseadas); más detalles sobre este procedimiento se tienen en la Sección 4.7.2.

4.5.2 Efecto del variograma (covariograma) en el kriging

El variograma (o el covariograma) tiene un efecto determinante en la predicción espacial. Como ejemplo, a continuación se incluyen algunas observaciones acerca de la influencia en el kriging de las tres principales características de un variograma estacionario (normalmente tratadas como parámetros) definidas en la Sección 1.3.

Rango

Supongamos que la posición de predicción \(\mathbf{s}_{0}\) está a una distancia mayor que el rango de las posiciones de los datos \(\left\{ \mathbf{s}_{1}, \ldots,\mathbf{s}_{n} \right\}\) (i.e. la posición de predicción está fuera de la zona de influencia de los datos), entonces \(\mathbf{c}=\mathbf{0}\) en las ecuaciones kriging, obteniéndose que: \[p_{KS} (\mathbf{Z},\mathbf{s}_{0})=\mu(\mathbf{s}_{0}),\ \ p_{KU} (\mathbf{Z},\mathbf{s}_{0})=\mathbf{x}_0^\top\hat{\boldsymbol{\beta}}_{gls},\] por lo tanto la predicción kriging se reduce a la media25 (estimada en el caso del KU).

Nugget y umbral

Resulta claro que las estimaciones obtenidas de la varianza de los predictores kriging dependen en gran medida de estos parámetros. Es importante destacar que la escala del variograma (o del covariograma) no influye en las predicciones obtenidas, solamente en la varianza kriging. Si se multiplica el variograma (o el covariograma) por una constante, las ecuaciones de los predictores kriging quedan invariantes y consecuentemente los pesos kriging no cambian, aunque la varianza kriging resulta multiplicada por esa constante. Para estudiar su influencia en la predicción resulta de utilidad la proporción del efecto nugget en el umbral total \(c_{0} /\sigma^{2}\) (que como ya se comentó al final de la Sección 1.3, proporciona mucha información acerca del grado de dependencia espacial presente en los datos). Por ejemplo, en el caso en que toda la variabilidad es efecto nugget (i.e. el proceso \(Z(\cdot)\) es ruido blanco) entonces \(\boldsymbol{\Sigma}=\sigma^{2} \mathbf{I}_{n}\) y \(\mathbf{c}=\mathbf{0}\) (suponiendo que \(\mathbf{s}_{0}\neq \mathbf{s}_{i}, \forall i\)), y los predictores kriging se reducen a la estimación OLS de la tendencia: \[p_{KU} (\mathbf{Z},\mathbf{s}_{0}) = \mathbf{x}_0^\top(\mathbf{X }^\top\mathbf{X})^{-1} \mathbf{X}^\top\mathbf{Z} = \mathbf{x }_0^\top\hat{\boldsymbol{\beta}}_{ols}.\] En el caso del KO se obtiene la media muestral: \[p_{KO} (\mathbf{Z},\mathbf{s}_{0})=\bar{Z} =\dfrac{1}{n} \sum\nolimits_{i=1}^{n}Z(\mathbf{s}_{i}),\] predictor bien conocido que asigna igual peso a todos los datos (por tanto el efecto pantalla es nulo). Además, teniendo en cuenta los resultados de la Sección 4.3.2, como \(\sigma_{KS}^{2} (\mathbf{s}_{0})=\sigma^{2}\) (caso más desfavorable), entonces: \[\sigma_{KO}^{2} (\mathbf{s}_{0})=\sigma_{}^{2} +\dfrac{\sigma_{}^{2} }{n},\] es decir, la varianza del KO para el caso de procesos incorrelados es igual a la varianza del proceso más la varianza de la media muestral. En general (al menos en KS y KO), se puede ver que al aumentar el porcentaje de efecto nugget en el umbral total disminuye el efecto pantalla y aumenta la varianza kriging (ver e.g. Isaaks y Srivastava, 1989, pp. 305-306). Hay que tener en cuenta que cuando la media no es constante (e.g. KU), puede ocurrir incluso lo contrario del efecto pantalla, y observaciones alejadas pueden tener una gran influencia en la estimación (como es bien conocido en la estimación de la tendencia).

Teniendo en cuenta los comentarios anteriores, podemos afirmar que todos los parámetros (características) del variograma influyen en el kriging (aunque quizás el rango es el que tiene un menor efecto, ya que pequeñas variaciones en este parámetro producen resultados casi idénticos). Estas observaciones no contradicen el resultado de que asintóticamente lo importante es capturar el comportamiento del variograma cerca del origen (Stein, 1988). Es difícil determinar cuando los datos están suficientemente cerca como para tener sólo en cuenta el efecto nugget y la pendiente del variograma cerca del origen.

4.5.3 Elección del vecindario

Una práctica habitual en geoestadística, en lugar de considerar todas las observaciones disponibles, es incluir en las ecuaciones kriging únicamente los \(n(\mathbf{s}_{0})\) datos más “próximos” a la posición de predicción \(\mathbf{s}_{0}\) . Esto puede justificarse por varias razones:

  • Utilizar todos los datos puede dar lugar a un sistema de difícil solución debido a problemas numéricos. Por ejemplo, entre otros, el tiempo de computación (aproximadamente proporcional a \(n(\mathbf{s}_{0})^{3}\)) aumenta drásticamente al aumentar el numero de datos.
  • Las estimaciones del variograma son normalmente eficientes (o incluso el propio modelo geoestadístico válido) únicamente en pequeñas distancias.
  • El uso de vecindarios locales permite la relajación de las hipótesis del modelo (como la estacionariedad intrínseca en el caso del KO) o su simplificación (por ejemplo, en el caso del KU se puede suponer que locamente la estructura de la tendencia es más simple o incluso constante y utilizar en su lugar KO local).
  • En muchos casos los datos cercanos apantallan a los más alejados reduciendo su influencia (aunque no siempre de forma significativa; ver observaciones siguientes).

La selección “optima” de un vecindario local resulta no obstante un problema algo complejo. Por ejemplo, se acostumbra a pensar que el rango del variograma permite determinar por sí solo un criterio de vecindad, como incluir en las ecuaciones sólo los datos que estén dentro del rango de \(\mathbf{s}_{0}\), sin embargo esto puede no ser adecuado en muchos casos (aunque en la práctica el valor de este parámetro puede ser de gran utilidad como referencia inicial). Teniendo en cuentas las observaciones realizadas en secciones anteriores, cuando aumenta la proporción de efecto nugget disminuye el efecto pantalla, el predictor kriging se reduce a la media y observaciones a más distancia que el rango de la posición de predicción contribuyen (a veces de forma significativa) a la estimación de la tendencia. Hay que tener en cuenta que los pesos kriging dependen de la configuración espacial de todas las observaciones y observaciones fuera del rango pueden tener influencia en la predicción a través de su correlación con observaciones dentro del rango (esto es conocido en la literatura como efecto relay, que podríamos traducir por efecto transmisión).

Se han propuesto algunos criterios para la selección de un vecindario óptimo (e.g. Cressie, 1993, pp. 176-177), que son de utilidad cuando los datos están regularmente espaciados y el vecindario se puede fijar de antemano. Por ejemplo, se pueden ir incluyendo observaciones en el vecindario hasta que no se produzca una disminución “significativa” en la estimación de la varianza kriging.

En la práctica, la densidad de datos y su configuración espacial en torno a las posiciones de predicción pueden ser muy irregulares. Teniendo en cuenta esto se han desarrollado distintos algoritmos, algunos bastante sofisticados, para la selección de vecindarios (ver e.g. Isaaks y Srivastava, 1989, Capítulo 14; Deutsch y Journel, 1992, secciones II.4 y IV.6). Para la selección de los datos típicamente se fija un radio máximo de búsqueda y únicamente se consideran los datos dentro de una circunferencia (esfera) centrada en la posición de predicción. En el caso de anisotropía (Sección 2.2.2), se considera una elipse (elipsoide) con el radio mayor orientado en la dirección de máxima variación. Además suele interesar disponer de observaciones en todas direcciones, por lo que se divide la zona de búsqueda en sectores angulares (por ejemplo cuadrantes en el caso bidimensional o octantes en el caso tridimensional) y se selecciona un número mínimo de datos en todos o en la mayoría de esos sectores (esto evita también que clusters de datos tengan una excesiva influencia sobre predicciones en su entorno). Si se sospecha de la presencia de una tendencia en los datos (KU), puede ser deseable la inclusión de observaciones más alejadas de la posición de predicción para poder estimarla de forma más eficiente.

Utilizando un algoritmo de búsqueda que tenga en cuenta todas o alguna de las consideraciones anteriores, típicamente se selecciona un número pequeño de datos como vecindario (e.g. entre 10 y 20 observaciones) para cada posición de predicción. Sin embargo esto puede causar que las superficies de predicción presenten discontinuidades (especialmente cuando los datos están regularmente espaciados y se utiliza búsqueda por octantes). Una aproximación distinta sería la selección un único vecindario más grande (e.g. de 20 a 40 observaciones) para pequeños conjuntos de posiciones de predicción, de esta forma en condiciones normales los vecindarios correspondientes a conjuntos de predicción próximos se solapan y es de esperar que no aparezcan discontinuidades. Además el incremento en tiempo de computación debido a la inclusión de un número mayor de observaciones se compensa por el hecho de que sólo es necesario factorizar una vez la matriz sistema kriging para obtener las predicciones en el conjunto considerado26.

Otra forma de proceder que puede ser de interés en la práctica, es usar el semivariograma como distancia en lugar de la tradicional distancia euclidea; de esta forma los datos son seleccionados preferentemente en la dirección de máxima continuidad (y se evita tener que considerar elipsoides en el caso de anisotropía). Adicionalmente, cuando el número de datos es grande y sus posiciones irregulares, es recomendable utilizar alguna técnica de búsqueda, como el denominado kd-tree (e.g. ver paquete FNN), de forma que se puedan determinar eficientemente las observaciones más próximas a la posición de predicción \(\mathbf{s}_{0}\) (en lugar de comprobar todas las posiciones).

Independientemente del algoritmo de búsqueda que se vaya a utilizar, puede ser recomendable realizar antes algunas pruebas utilizando por ejemplo la técnica de validación cruzada descrita a continuación.

  1. Para un tratamiento más detallado ver por ejemplo Chilès y Delfiner (2012), secciones 3.3.2 y 3.4.2.↩︎

  2. En caso contrario interesaría predecir el proceso libre de ruido y habría que modificar ligeramente las ecuaciones kriging. Para más detalles, ver por ejemplo Cressie (1993, pp. 128-130) o Chilès y Delfiner (2012, Sección 3.7.1)↩︎

  3. En la literatura se muestran numerosos ejemplos sobre el comportamiento de los pesos kriging en distintos escenarios; una colección bastante completa se tiene en Wakernagel (1998, Capítulo 13).↩︎

  4. Puede verse fácilmente que la ecuaciones del kriging universal con \(\mathbf{c}=\mathbf{0}\), son las obtenidas en el denominado método kriging de estimación de la tendencia (e.g. Wackernagel, 1998, pp. 212-213; Chilès y Delfiner, 2012, Sección 3.4.5).↩︎

  5. Por ejemplo, si se pretenden obtener predicciones en una rejilla bidimensional, el tiempo de computación considerando un vecindario distinto con 16 observaciones para cada nodo resulta similar a considerar un vecindario con 25 (o incluso más) observaciones para grupos de 4 nodos (dos por fila y columna).↩︎