1.2 Geoestadística

La geoestadística (Matheron 1962) surgió como una mezcla de varias disciplinas: ingeniería de minas, geología, matemáticas y estadística, para dar respuesta a problemas como, por ejemplo, el de la estimación de los recursos de una explotación minera (se desarrolló principalmente a partir de los años 80). La diferencia (ventaja) respecto a otras aproximaciones es que, además de tener en cuenta la tendencia espacial (variación de gran escala), también tiene en cuenta la correlación espacial (variación de pequeña escala). Otros métodos sin embargo, sólo incluyen la variación de larga escala y suponen que los errores son independientes (Sección 1.2.1). Hoy en día podemos decir que la geoestadística es la rama de la estadística espacial que estudia los procesos con índice espacial continuo.

Uno de los problemas iniciales más importantes de la geoestadística fue la predicción de la riqueza de un bloque minero a partir de una muestra observada. A este proceso Matheron (1963) lo denominó kriging1, y también predicción espacial lineal óptima (estos métodos de predicción se muestran en el Capítulo 4).

El modelo general habitualmente considerado en geoestadística considera que el proceso se descompone en variabilidad de gran escala y variabilidad de pequeña escala: \[\begin{equation} Z(\mathbf{s}) = \mu(\mathbf{s}) + \varepsilon(\mathbf{s}), \tag{1.1} \end{equation}\] donde:

  • \(\mu(\mathbf{s}) = E \left( Z(\mathbf{s}) \right)\) es la tendencia (función de regresión, determinística).

  • \(\varepsilon(\mathbf{s})\) es un proceso de error de media cero que incorpora la dependencia espacial.

Como en condiciones normales únicamente se dispone de una realización parcial del proceso, se suelen asumir hipótesis adicionales de estacionariedad sobre el proceso de error \(\varepsilon(\mathbf{s})\) para hacer posible la inferencia. En la Sección 1.3 se definen los principales tipos de estacionariedad habitualmente considerados en geoestadística y se introducen dos funciones relacionadas con procesos estacionarios, el covariograma y el variograma. Algunas propiedades de estas funciones, que podríamos decir que son las herramientas fundamentales de la geoestadística, se muestran en la Sección 1.3.

1.2.1 Modelos clásicos y modelos espaciales

En general, cuando se considera que la componente espacial (o espacio-temporal) puede ser importante en el modelado y el análisis de los datos es necesaria una aproximación estadística distinta a la tradicionalmente usada.

Uno de los modelos más utilizados en estadística para el caso de datos no homogéneos es el conocido modelo clásico de regresión lineal. Si \(\left\{ Z(\mathbf{s}):\mathbf{s}\in D\subset \mathbb{R}^{d} \right\}\) es un proceso espacial, podemos suponer que: \[\begin{equation} Z(\mathbf{s})=\sum\limits_{j=0}^{p}X_{j}(\mathbf{s})\beta_{j} +\varepsilon(\mathbf{s}),\ \mathbf{s}\in D, \tag{1.2} \end{equation}\] (un caso particular del modelo general (1.1)), donde \(\boldsymbol{\beta }=(\beta_{0}, \ldots,\beta_{p})^{\top}\in \mathbb{R}^{p+1}\) es un vector desconocido, \(\left\{ X_{j} (\cdot):j=0, \ldots,p\right\}\) un conjunto de variables explicativas (típicamente \(X_0(\cdot)=1\)) y \(\varepsilon(\cdot)\) un proceso de media cero incorrelado (i.e. \(Cov(\varepsilon (\mathbf{u}),\varepsilon (\mathbf{v}))=0\) si \(\mathbf{u}\neq \mathbf{v}\)) con \(Var(\varepsilon (\mathbf{s}))=\sigma^{2}\).

Supongamos por el momento que el objetivo es la estimación eficiente de la tendencia, o lo que es lo mismo la estimación óptima de los parámetros de la variación de gran escala \(\boldsymbol{\beta}\), a partir de los datos observados en un conjunto de posiciones espaciales \(\left\{ \mathbf{s}_{1}, \ldots,\mathbf{s}_{n} \right\}\). Bajo las hipótesis anteriores: \[\mathbf{Z} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon},\] siendo \(\mathbf{Z}=\left( Z(\mathbf{s}_{1}), \ldots,Z(\mathbf{s}_{n}) \right)^{\top}\), \(\mathbf{X}\) una matriz \(n\times (p+1)\) con \(\mathbf{X}_{ij}=X_{j-1}(\mathbf{s}_{i})\) y \(\boldsymbol{\varepsilon} = \left( \varepsilon (\mathbf{s}_{1}), \ldots,\varepsilon (\mathbf{s}_{n})\right)^{\top}\); y el estimador lineal insesgado de \(\boldsymbol{\beta}\) más eficiente resulta ser el estimador de mínimos cuadrados ordinarios (OLS, ordinary least squares): \[\begin{equation} \hat{\boldsymbol{\beta}}_{ols}=(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Z}, \tag{1.3} \end{equation}\] con \[Var(\hat{\boldsymbol{\beta}}_{ols})=\sigma^{2}(\mathbf{X}^{\top}\mathbf{X})^{-1}.\]

Sin embargo la suposición de que los errores son independientes e idénticamente distribuidos influye crucialmente en la inferencia. En el modelo anterior, en lugar de errores incorrelados, si suponemos que: \[Var\left( \boldsymbol{\varepsilon} \right) =\boldsymbol{\Sigma},\] obtenemos el modelo lineal de regresión generalizado y en este caso el estimador lineal óptimo de \(\boldsymbol{\beta}\) es el estimador de mínimos cuadrados generalizados (GLS, generalized least squares): \[\begin{equation} \hat{\boldsymbol{\beta}}_{gls} =(\mathbf{X}^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{X})^{-1} \mathbf{X}^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{Z}. \tag{1.4} \end{equation}\]

Si \(\boldsymbol{\Sigma}=\sigma^{2} \mathbf{I}_{n}\), siendo \(\mathbf{I}_{n}\) la matriz identidad \(n\times n\), los estimadores (1.3) y (1.4) coinciden; pero en caso contrario las estimaciones basadas en el modelo anterior pueden llegar a ser altamente ineficientes. Puede verse fácilmente que en el caso general: \[Var\left( \hat{\boldsymbol{\beta}}_{gls} \right)=(\mathbf{X}^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{X})^{-1}, \\ Var\left( \hat{\boldsymbol{\beta}}_{ols} \right) =(\mathbf{X}^{\top}\mathbf{X})^{-1} (\mathbf{X}^{\top}\boldsymbol{\Sigma}\mathbf{X})(\mathbf{X}^{\top}\mathbf{X})^{-1},\] resultando además que \(Var( \hat{\boldsymbol{\beta}}_{ols}) - Var( \hat{\boldsymbol{\beta}}_{gls} )\) es una matriz semidefinida positiva (e.g. Searle, 1971, Sección 3.3).

En muchos casos el objetivo final es la predicción del proceso en una posición espacial \(\mathbf{s}_{0}\): \[Z(\mathbf{s}_{0} )=\mathbf{x}_0^{\top}\boldsymbol{\beta}+\varepsilon (\mathbf{s}_{0} ),\] donde \(\mathbf{x}_0=\left( X_{0} (\mathbf{s}_{0} ), \ldots,X_{p} (\mathbf{s}_{0})\right)^{\top}\). Bajo las hipótesis del modelo clásico el predictor óptimo sería la estimación de la tendencia \(\hat{\mu}(\mathbf{s}_{0} ) = \mathbf{x}_0^\top \hat{\boldsymbol{\beta}}_{ols}\) (el predictor óptimo de un error independiente sería cero). En el caso general, siguiendo esta aproximación, podríamos pensar en utilizar como predictor el estimador más eficiente de la tendencia: \[\hat{Z} (\mathbf{s}_{0})=\mathbf{x}_0^\top \hat{\boldsymbol{\beta}}_{gls},\] sin embargo no es el predictor lineal óptimo. Puede verse (e.g. Goldberger, 1962; Sección 4.3.2) que en este caso el mejor predictor lineal insesgado es: \[\begin{equation} \tilde{Z}(\mathbf{s}_{0}) = \mathbf{x}_0^{\top}\hat{\boldsymbol{\beta}}_{gls} + \mathbf{c}^{\top} \boldsymbol{\Sigma}^{-1} \left( \mathbf{Z} - \mathbf{X}\hat{\boldsymbol{\beta}}_{gls} \right), \tag{1.5} \end{equation}\] (el denominado predictor del kriging universal, Sección 4.3), siendo \[\mathbf{c} = \left( Cov\left( \varepsilon (\mathbf{s}_{1} ),\varepsilon (\mathbf{s}_{0} )\right), \ldots, Cov\left( \varepsilon (\mathbf{s}_{n} ),\varepsilon (\mathbf{s}_{0} )\right) \right)^{\top},\] y la diferencia \(Var( \hat{Z} (\mathbf{s}_{0} ) ) - Var( \tilde{Z} (\mathbf{s}_{0} ) ) \ge 0\) puede ser significativamente mayor que cero2 (si la dependencia espacial no es muy débil). Naturalmente, si se ignora por completo la dependencia y se emplea únicamente el estimador \(\hat{\boldsymbol{\beta}}_{ols}\) disminuye aún más la eficiencia de las predicciones.

Teniendo en cuenta los resultados anteriores podemos afirmar que al explotar la dependencia presente en los datos el incremento en eficiencia puede ser importante. Sin embargo el principal inconveniente es que en la práctica normalmente la matriz \(\boldsymbol{\Sigma}\) y el vector \(\mathbf{c}\) son desconocidos. El procedimiento habitual, para evitar la estimación de \(n+n(n+1)/2\) parámetros adicionales a partir del conjunto de \(n\) observaciones, suele ser la elección de un modelo paramétrico adecuado (ver Sección 3.2.1): \[C(\mathbf{u},\mathbf{v}\left| \boldsymbol{\theta}\right. )\equiv Cov\left( \varepsilon (\mathbf{u}),\varepsilon (\mathbf{v})\right),\] i.e. suponer que \(\boldsymbol{\Sigma} \equiv \boldsymbol{\Sigma}\left( \boldsymbol{\theta}\right)\) y \(\mathbf{c}\equiv \mathbf{c}\left( \boldsymbol{\theta}\right)\). Una hipótesis natural es suponer que los datos cercanos en el espacio o en el tiempo están correlados y que la correlación disminuye al aumentar la separación entre ellos; por tanto es normal pensar en errores espacialmente correlados. Por ejemplo, podemos considerar: \[C(\mathbf{u},\mathbf{v}\left| \boldsymbol{\theta}\right. )=\sigma^{2} \rho^{\left\| \mathbf{u}-\mathbf{v}\right\| },\] con \(\sigma^{2} \geq 0\) y \(0<\rho <1\). De esta forma, si \(\hat{\boldsymbol{\theta}}\) es un estimador de \(\boldsymbol{\theta}\) (ver Sección 3.3), podemos obtener por ejemplo una aproximación del predictor óptimo de \(Z(\mathbf{s}_{0} )\) sustituyendo en (1.5) \(\boldsymbol{\Sigma}(\boldsymbol{\theta})\) por \(\boldsymbol{\Sigma}(\hat{\boldsymbol{\theta}} )\) y \(\mathbf{c}(\boldsymbol{\theta})\) por \(\mathbf{c}(\hat{\boldsymbol{\theta}} )\).

1.2.2 Ventajas de la aproximación espacial (y espacio-temporal)

Algunos de los beneficios de utilizar modelos espaciales para caracterizar y explotar la dependencia espacial de un conjunto de datos son los siguientes:

  • Modelos más generales: en la mayoría de los casos, los modelos clásicos no espaciales son un caso particular de un modelo espacial.

  • Estimaciones más eficientes: de la tendencia, de los efectos de variables explicativas, de promedios regionales…

  • Mejora de las predicciones: más eficientes, con propiedades de extrapolación más estables…

  • La variación espacial no explicada en la estructura de la media debe ser absorbida por la estructura del error, por lo que un modelo que incorpore la dependencia espacial puede decirse que esta protegido frente a una mala especificación de este tipo. Esto en muchos casos tiene como resultado una simplificación en la especificación de la tendencia; en general los modelos con dependencia espacial suelen tener una descripción más parsimoniosa (en ocasiones con muchos menos parámetros) que los clásicos modelos de superficie de tendencia.

  1. D. G. Krige fue un ingeniero de minas de Sudáfrica que desarrolló en los años 50 métodos empíricos para determinar la distribución de la riqueza de un mineral a partir de valores observados. Sin embargo la formulación de la predicción espacial lineal óptima no procede del trabajo de Krige. Al mismo tiempo que la geoestadística se desarrollaba en la ingeniería de minas por Matheron en Francia, la misma idea se desarrollaba en la meteorología por L.S. Gandin en la antigua Unión Soviética. El nombre que Gandin le dio a esta aproximación fue análisis objetivo y utilizó la terminología de interpolación óptima en lugar de kriging. Para más detalles sobre el origen del kriging ver p.e. Cressie (1990).↩︎

  2. Por ejemplo, para un caso particular, Goldberger (1962, pp.374-375) observó que la mejora en la predicción puede llegar a ser del 50%. En Cressie (1993, Sección 1.3) se muestran también otros ejemplos del efecto de la presencia de correlación en la estimación.↩︎