La pregunta central de este capítulo es: ¿cómo expresamos en una función la estructura de dependencia espacial presente en la realización observada? Esto se conoce como análisis estructural, y es clave para la predicción óptima (kriging), porque el éxito del kriging depende de las funciones que dan información sobre esa dependencia.

Esas funciones son los covariogramas y los semivariogramas. Como en la práctica solo tenemos la realización observada, se calcula primero una versión empírica de estas funciones y luego se le ajusta un modelo teórico válido.

El covariograma es útil como fundamento teórico, pero en la práctica geoestadística el protagonista es el semivariograma: cubre un espectro más amplio de variables regionalizadas y no exige conocer la media de antemano. Por eso este documento repasa el covariograma con brevedad y le da el peso principal al semivariograma.

Como se vio en la Definición 2.2.5, la función de covarianza de una fa espacial es:

Bajo estacionariedad de segundo orden, depende solo del vector h que une las ubicaciones, no de las ubicaciones mismas:

Definición 3.2.1 y 3.2.2Si \(C(h)\) depende solo de la distancia \(|h|\) (no de la dirección), se llama isotrópico. Si depende también de la dirección, es anisotrópico —justo el tema que exploramos en la Sección 3.5.

Solo un puñado de funciones han sido demostradas válidas como covariograma en \(\mathbb{R}^d\). Las tres más usadas en la práctica:

El esférico (el más usado en la práctica) llega exactamente a cero en el rango \(a\); es continuo pero no diferenciable. El exponencial y el Gaussiano nunca llegan exactamente a cero: se usa el rango práctico (la distancia al 5% del valor en el origen): \(3a\) para el exponencial, \(a\sqrt{3}\) para el Gaussiano. El Gaussiano es infinitamente diferenciable (fa muy suaves, poco comunes en la práctica).

Este ya lo conocesEl modelo esférico es exactamente el que usamos en la Parte 2 (rango \(a=110\)) para mostrar por qué ignorar la correlación espacial dispara la varianza real de la media muestral.

En la práctica solo tenemos la realización observada. En el marco de estacionariedad de segundo orden, el estimador de momentos (MoM) más usado es:

con \(\bar Z\) la media muestral (estimador de \(\mu\)) y \(\#N(h)\) el número de pares de ubicaciones separados por la distancia \(h\).

Como se usa \(\bar Z\) en vez del verdadero \(\mu\), \(\hat C(h)\) es un estimador sesgado de \(C(h)\): entre más pequeño \(\#N(h)\), mayor el sesgo. Por eso el covariograma empírico no debe usarse directamente para predicción —puede no ser definido positivo, y solo está definido en unas pocas distancias.

Para calcular un covariograma con sentido bajo estacionariedad de segundo orden, necesitamos datos sin tendencia —el transecto A (MO creciente) no sirve aquí, porque su media no es constante. Usamos el transecto B: 24 puntos cada 10 m, oscilando alrededor de una media aproximadamente constante (\(\bar Z=2.971\)):

Nota importanteA partir de h=60 m, \(\hat C(h)\) se vuelve negativa. Esto es válido para un covariograma (a diferencia del semivariograma, que nunca puede serlo): significa que a esa distancia, cuando un punto está por encima de la media, el otro tiende a estar por debajo.

Bajo estacionariedad de segundo orden o intrínseca (sin deriva):

En el marco de segundo orden, semivariograma y covariograma son equivalentes: \(C(h)=C(0)-\gamma(h)\) (3.13). Pero el semivariograma se prefiere porque no requiere conocer la media de la fa (que en la práctica hay que estimar, introduciendo sesgo) y porque cubre un espectro más amplio de fenómenos: sirve incluso cuando la covarianza ni siquiera está definida.

El semivariograma tiene varias propiedades matemáticas que, más que un ejercicio formal, existen porque garantizan que el kriging funcione correctamente. Vale la pena conocerlas por su utilidad práctica, no por la demostración en sí:

Por definición \(\gamma(0)=0\): la disimilitud de un punto consigo mismo es nula. En la práctica, casi siempre aparece una discontinuidad ahí de todas formas —el efecto pepita que veremos en la Sección 3.4.4—, y esa discontinuidad es justamente la señal de variabilidad que no se explica por la posición.

Es una función par (\(\gamma(h)=\gamma(-h)\)): la disimilitud entre dos puntos no depende de en qué orden los mires, solo de la distancia y dirección que los separa. En la práctica, esto significa que basta con calcular y graficar el semivariograma para \(h\ge0\) —no hay que preocuparse por el lado "negativo".

Nunca toma valores negativos, a diferencia del covariograma (que sí puede serlo, como vimos con el transecto B). Su utilidad práctica: cualquier curva o modelo ajustado que cruce por debajo de cero se puede descartar de entrada como inválido, sin necesidad de revisar nada más.

Es condicionalmente definida negativa —una condición técnica cuya utilidad práctica es garantizar que la varianza de cualquier predicción de kriging construida con ese semivariograma nunca salga negativa. Sin esta propiedad, un modelo ajustado "a mano" o de forma ingenua podría producir errores de predicción sin sentido matemático.

Por último, el semivariograma de una fa estacionaria de segundo orden es siempre finito (alcanza una meseta); el de una fa intrínseca sin deriva puede crecer sin límite, pero nunca más rápido que \(|h|^2\). Su utilidad práctica: sirve como diagnóstico rápido de si un ajuste a grandes distancias tiene sentido, o si el modelo se está extrapolando de forma poco realista.

El semivariograma es, en general, monótono no decreciente. Sube desde el origen y se aproxima a su valor límite, la meseta (o sill, que coincide con \(C(0)\)). La distancia a la que se alcanza se llama rango: el umbral de dependencia espacial. Más allá del rango, dos observaciones ya no aportan información una sobre la otra —están, para efectos prácticos, descorrelacionadas.

Este es el diagrama que resume, en una sola imagen, casi todo lo que hemos definido hasta aquí: los puntos rojos son el semivariograma experimental (la Ec. 3.12 calculada de los datos, con su dispersión natural); la curva negra es el modelo que se le ajusta. El nugget es la discontinuidad en el origen; el scale (meseta parcial) es lo que sube desde el nugget hasta la meseta total; el sill es la meseta completa (nugget + scale); y el length (rango) marca la distancia a partir de la cual los datos pasan de estar correlacionados a estar no correlacionados.

Lectura, con el vocabulario de la imagen clave\(\gamma(h)\) sube con claridad hasta \(h\approx90\text{–}100\) m (el length/rango) y luego se estabiliza alrededor de 0.63–0.64 (el sill). Como partimos de \(\gamma(10)=0.051\) y no de 0 exactamente, este transecto también insinúa un pequeño nugget.

El transecto A (MO creciente, del documento anterior) produce un semivariograma que nunca se aplana. Eso ocurre porque tiene deriva (no es estacionario): semivariogramas sin meseta son comunes con fa no estacionarias, con fa intrínsecas sin meseta finita, o cuando el rango excede la distancia máxima estimable.

Aclaración honesta antes de ver los paneles: la estacionariedad estricta y la de segundo orden producen variogramas indistinguibles a simple vista —ambos suben y se aplanan en una meseta. La diferencia entre ellas no vive en el variograma, sino en toda la distribución conjunta. Por eso el primer panel representa a las dos al tiempo.

El comportamiento del semivariograma justo cerca del origen —no su forma general, sino los primeros metros— dice mucho sobre qué tan regular o irregular es el fenómeno que se está modelando, y tiene una utilidad práctica muy directa: ayuda a decidir qué modelo teórico probar primero.

Si la curva sube de forma lineal cerca de \(h=0\) (como el modelo esférico o el exponencial), el fenómeno es continuo pero no especialmente suave —es el comportamiento más común en variables naturales como la materia orgánica del suelo.

Si sube de forma parabólica (como el modelo Gaussiano), el fenómeno es extremadamente regular, casi sin cambios bruscos entre puntos vecinos. Esto es poco común en la práctica, y cuando aparece conviene sospechar: puede ser indicio de que los datos ya venían suavizados antes de llegar a ti (por ejemplo, por una interpolación previa) en vez de una propiedad genuina del fenómeno.

Y si la curva ni siquiera vuelve a cero cuando \(h\to0\) (una discontinuidad clara), eso es exactamente el efecto pepita de la siguiente sección.

Utilidad prácticaEsta lectura rápida te dice, antes de ajustar nada formalmente, qué familia de modelo probar primero: esférico o exponencial para el caso típico (lineal cerca del origen), Gaussiano solo si tienes evidencia genuina de un fenómeno muy suave —y si no se cumple ninguno de los dos, hay que sospechar de un efecto pepita real.

Aunque \(\gamma(0)=0\) por definición, en la práctica el semivariograma empírico muchas veces no vuelve a cero cerca del origen —es justo el nugget de la imagen clave. Suele indicar que la variable es muy irregular, quizás discontinua a escalas menores que la distancia mínima muestreada.

Causas típicas: una microestructura con rango más corto que la distancia mínima entre muestras, o errores de medición. Si \(Z^*(s)=Z(s)+\varepsilon(s)\), con \(\varepsilon(s)\) el error de medición, independiente y sin correlación espacial:

es decir, los datos adquieren un efecto pepita adicional al que ya pudiera tener la fa. El caso límite es el efecto pepita puro: el semivariograma es constante para cualquier distancia —ausencia total de correlación espacial.

EjemploSi tu equipo de laboratorio para medir materia orgánica tiene un error de medición con desviación estándar de 0.1% (independiente de dónde se tomó la muestra), ese \(0.1\%^2=0.01\) se sumará como nugget a cualquier semivariograma que calcules, incluso si el suelo mismo no tuviera ninguna microestructura de corto alcance.

Los semivariogramas teóricos son funciones con una expresión analítica simple que cumplen la condición de definición negativa condicional (Sección 3.4.1). Se usan para representar a los semivariogramas empíricos y son esenciales para el kriging. No se deducen de ninguna hipótesis especial: son simplemente herramientas matemáticas válidas para modelar dependencia espacial.

Se dividen en tres grupos: (1) con meseta (o de transición), (2) con meseta y efecto de hoyo, y (3) sin meseta.

Se asocian con la hipótesis estacionaria de segundo orden. La distancia a la que se alcanza la meseta se llama rango, y marca la transición de correlación espacial a ausencia de ella.

Válido en \(\mathbb{R}^1,\mathbb{R}^2,\mathbb{R}^3\). Comportamiento lineal cerca del origen (continuo pero algo irregular), con pendiente \(1.5m/a\). La tangente en el origen cruza la meseta en \(|h|=2a/3\). Es el modelo más usado en la práctica.

Válido en \(\mathbb{R}^d, d\ge1\). Lineal cerca del origen, con pendiente \(m/a\) (1.5 veces menor que el esférico). Alcanza la meseta solo asintóticamente; se usa el rango práctico \(a'\approx3a\) (donde \(\gamma(a')=0.95m\)).

Comportamiento parabólico cerca del origen (pendiente nula) —implica regularidad infinita, poco realista salvo para fenómenos extremadamente suaves. Rango práctico \(a'=a\sqrt3\).

Válido en \(\mathbb{R}^1,\mathbb{R}^2,\mathbb{R}^3\). Comportamiento parabólico cerca del origen, similar al Gaussiano, pero finitamente diferenciable (más realista). Alcanza una meseta plana exactamente en \(h=a\).

Generaliza al exponencial (\(\alpha=1\)) y al Gaussiano (\(\alpha=2\)) en un solo parámetro de forma \(\alpha\). Rango práctico \(a'=a\cdot3^{1/\alpha}\).

Ejemplo con α=1.5 (entre exponencial y Gaussiano)

Parabólico cerca del origen. Para \(\alpha\) pequeño (\(<2\)) alcanza la meseta muy lentamente —útil para modelar dependencia espacial de largo alcance.

Ejemplo con α=0.5 (aproximación lenta a la meseta)

Donde \(K_\alpha\) es la función de Bessel modificada de segunda especie. Es el más flexible de la familia: según \(\alpha\) puede comportarse de forma lineal, parabólica o cualquier punto intermedio cerca del origen. Para \(\alpha=1/2\) se reduce exactamente al modelo exponencial.

Ejemplo con α=1.5

Estos modelos no son funciones monótonas: oscilan, permitiendo dependencia espacial negativa. Reflejan componentes subyacentes con una periodicidad física real (por ejemplo, patrones de siembra en hileras, o estratos geológicos repetitivos).

Donde \(J_\alpha\) es la función de Bessel de primera especie. Válido en \(\mathbb{R}^d\), \(d\le2(\alpha+1)\).

Ejemplo con α=1: sube, pasa la meseta, y oscila amortiguándose

Es la particularización del J-Bessel para \(\alpha=1/2\). Modelo pseudo-periódico, con comportamiento parabólico cerca del origen, asociado a estructuras muy continuas. Rango práctico \(\approx20.37a\); la amplitud del efecto de hoyo es 1.2 veces la meseta.

Van más allá de la hipótesis estacionaria de segundo orden: corresponden a funciones aleatorias intrínsecas pero no estacionarias de segundo orden. Ya los vimos en acción en el documento anterior, con el transecto de materia orgánica creciente.

Válido en \(\mathbb{R}^d, d\ge1\). Gran variedad de comportamientos posibles cerca del origen según \(\alpha\). Para \(\alpha\ge2\) deja de cumplir la hipótesis intrínseca; \(\alpha=0\) corresponde a un efecto pepita puro.

Ejemplo con α=0.7

Caso particular del modelo potencial con \(\alpha=1\). Puede indicar una fa que en realidad sí es estacionaria de segundo orden, pero cuya meseta simplemente no se alcanza dentro de las distancias observadas —el mismo fenómeno que vimos con el transecto A (MO creciente) en el documento anterior.

Modelo problemático: no está definido en el origen (\(\gamma\to-\infty\) cuando \(h\to0\)) y tampoco tiene meseta (\(\gamma\to\infty\) cuando \(h\to\infty\)). Se usa en contextos específicos con variables regularizadas —no como primera opción de modelado.

En la práctica, los semivariogramas experimentales rara vez se ajustan a un solo modelo simple. Se pueden combinar (superponer, "anidar") para crear modelos más complejos, sumando estructuras a distintas escalas.

Ejemplo: efecto pepita puro + dos esféricos anidados

Cada escala de variación integra la variabilidad de todas las escalas menores. La suma de semivariogramas corresponde, conceptualmente, a la suma de procesos espaciales independientes actuando a distintas distancias.

m₁=0.2 (nugget), m₂=0.3 con a₂=3 (corto alcance), m₃=0.5 con a₃=10 (largo alcance)

Todo lo anterior asume implícitamente que \(\gamma(h)\) depende solo de la distancia \(|h|\), no de la dirección (isotropía). En la práctica, muchos fenómenos tienen una dirección preferente —por ejemplo, la textura del suelo puede estar más correlacionada a lo largo de las curvas de nivel que perpendicular a ellas. Calcular el variograma por separado en distintas direcciones permite detectar esto.

Para comparar objetivamente qué tan bien ajusta cada dirección, se usa el mismo criterio que emplea fit.variogram() internamente —el criterio de Cressie: una suma de cuadrados ponderada por el número de pares \(N(h)\) en cada distancia, comparando los puntos empíricos de esa dirección contra una curva teórica de referencia:

Menor \(SSR_\theta\) significa que esa dirección se parece más al modelo de referencia.

Que 90° tenga el SSR más bajo no basta para concluir que hay anisotropía real: con datos finitos, siempre habrá alguna dirección con el valor más bajo, aunque el campo sea perfectamente isotrópico. Para obtener una conclusión con respaldo estadístico, se usa una prueba de Montecarlo:

1. Se calcula la dispersión observada entre direcciones: \( \text{estadístico} = \max_\theta(SSR_\theta) - \min_\theta(SSR_\theta) \).
2. Se simulan muchos campos nuevos bajo el modelo isotrópico ajustado —asumiendo que \(H_0\) (no hay dirección preferente) es cierta.
3. En cada campo simulado se repite el cálculo, obteniendo una distribución nula de cuánta dispersión produce el puro azar de muestreo.
4. El p-valor es la proporción de campos simulados cuya dispersión iguala o supera la observada en los datos reales.

Interpretación: p-valor pequeño (<0.05) ⇒ se rechaza la isotropía, hay evidencia real de una dirección preferente. P-valor grande ⇒ la diferencia observada es indistinguible de lo que produciría un campo isotrópico por simple variabilidad de muestreo.

Resultado con nuestros datosEstadístico observado = 233.3; p-valor Montecarlo = 0.97. No se rechaza la isotropía —coherente con que el campo se generó sin ninguna dirección preferente real.

Hasta ahora hemos usado el variograma para predecir (kriging). Pero tiene otra aplicación igual de importante: corregir un modelo de regresión cuando sus residuos están correlacionados espacialmente. La ecuación general es:

donde \(Y(s)\) es la variable respuesta en la ubicación \(s\), \(X(s)\) son las covariables, \(\beta\) los coeficientes de regresión, y \(\delta(s)\) un error espacialmente correlacionado, cuya estructura de covarianza se define a partir del variograma.

La clave: el variograma no se ajusta "aparte" y ya —se estima junto con los coeficientes \(\beta\) mediante máxima verosimilitud (o mínimos cuadrados generalizados, GLS). El variograma define la matriz de covarianza \(\Omega\) que corrige tanto los coeficientes como sus errores estándar.

El error no es ruido blanco: su covarianza está definida por el variograma ajustado a los residuos,

y la matriz de covarianza completa se construye como \(\Omega_{ij}=\sigma^2\cdot\rho(|s_i-s_j|)\), con \(\rho(h)\) la función de correlación derivada del variograma. Esa matriz \(\Omega\) es la que "corrige" el modelo.

Precisión importante: es un error común (y bastante extendido) decir que OLS da coeficientes "sesgados" cuando hay correlación espacial. No es así: \(\hat\beta_{OLS}\) sigue siendo insesgado —en promedio, sobre muchas muestras, acierta igual que GLS. Lo que falla en OLS son los errores estándar, casi siempre subestimados (exactamente el mismo fenómeno que vimos en la Parte 2 de este curso con la varianza de la media muestral). Eso es lo que produce valores p artificialmente pequeños e intervalos de confianza artificialmente angostos —no coeficientes equivocados en promedio, sino una falsa sensación de precisión.

Antes de mostrarte los resultados, hay un bug real en la versión original del script que vale la pena señalar: usa corSpatial = TRUE dentro de corSpher(), un argumento que no existe en esa función de nlme —el modelo GLS ni siquiera llega a ajustarse, R se detiene con Error: unused argument (corSpatial = TRUE). También conviene darle a corSpher() un valor inicial explícito (rango y proporción de nugget tomados del variograma de los residuos OLS); sin eso, aunque el código corra sin error, el optimizador puede quedarse atascado y el GLS termina siendo numéricamente idéntico al OLS —lo comprobé aquí mismo antes de corregirlo.

Cómo leer esto honestamenteEn esta muestra concreta, los coeficientes de OLS no quedaron sistemáticamente peor que los de GLS —de hecho \(\beta_1\) y \(\beta_2\) de OLS cayeron un poco más cerca de la verdad que los de GLS. Eso es normal: ambos estimadores son insesgados, así que en una sola muestra cualquiera de los dos puede tocarle "ganar" por azar. Lo que sí mejora de forma clara e inequívoca es: (1) los errores estándar de GLS son más chicos y correctos (más eficiencia real, no un artificio), y (2) el AIC cae más de 100 puntos —el modelo GLS describe genuinamente mejor los datos, porque no le está pidiendo a las covariables \(X_1,X_2\) que expliquen una estructura espacial que en realidad viene del error.

1. Para diagnosticar: se ajusta un variograma a los residuos del modelo OLS. Si muestra estructura real (no es solo pepita), hay correlación espacial que las covariables no capturaron.
2. Para corregir: los parámetros de ese variograma (nugget, sill, rango) definen la estructura de correlación \(\Omega\) del modelo GLS —lo que corrige los errores estándar, y con ellos, los valores p y los intervalos de confianza.

Resumen: ¿qué demuestra este modelo?

1. El variograma no es solo para kriging —también define la matriz de covarianza en modelos de regresión espacial.

2. Ignorar la correlación espacial no sesga los coeficientes, pero sí invalida los errores estándar, los valores p y los intervalos de confianza —una falsa sensación de precisión, el mismo fantasma de la Parte 2.

3. El GLS con estructura espacial corrige eso, y las medidas de ajuste (AIC, verosimilitud) mejoran porque el modelo refleja mejor la verdadera estructura de los datos.

Este es, en esencia, el tipo de modelo que describen Hoeting et al. (2006) —Model selection for geostatistical models, Ecological Applications—: un modelo de regresión donde el error tiene una estructura espacial definida por un variograma o una función de correlación tipo Matérn, y donde la selección de variables se hace con criterios como el AIC que ya incorporan esa estructura, en vez de ignorarla.