Geoestadistica no parametrica con el paquete *npsp*'

knitr::opts_chunk$set(comment=NA, fig.height=5, fig.width=7, fig.align = 'center', cache = FALSE)
# chunk options: results='hold', fig.show = 'hold'
library(npsp)

Geoestadística no paramétrica

Nos centraremos en el caso de procesos espaciales $\left\lbrace Y(\mathbf{x}), \mathbf{x} \in D\subset \mathbb{R}^{d} \right\rbrace$, con dominio $D$ continuo fijo: proceso geoestadístico.

Normalmente sólo se observa un conjunto de valores $\left{ y(\mathbf{x}{1}),\ldots,y(\mathbf{x}{n})\right}$ (realización parcial).

Ejemplo

Se utilizará como ejemplo el conjunto de datos aquifer: nivel del agua subterránea en el acuífero Wolfcamp.

str(aquifer)
#   Scatter plot with a color scale
with(aquifer, spoints(lon, lat, head, main = "Wolfcamp aquifer data"))

Modelo general

Se supone que el proceso se descompone en variabilidad de gran escala y variabilidad de pequeña escala: $$Y(\mathbf{x})=\mu(\mathbf{x})+\varepsilon(\mathbf{x}),$$ donde:

La dependencia espacial se modela normalmente a través del variograma: $$2\gamma(\mathbf{u}) = Var(\varepsilon \left( \mathbf{x}\right) -\varepsilon \left(\mathbf{x}+\mathbf{u}\right) ),$$ verificando $\gamma(\mathbf{u}) = C(0) - C(\mathbf{u})$.

Se empleará un modelo no paramétrico: no se supone ninguna forma concreta para $\mu(\mathbf{\cdot})$ y $\gamma(\mathbf{\cdot})$ (funciones suaves).

Ejemplo

Con el conjunto de datos de ejemplo puede ser preferible emplear métodos paramétricos. El tamaño muestral es pequeño y un modelo lineal para la tendencia parece ser adecuado.

x <- aquifer[,1:2]      # coordenadas espaciales
y <- aquifer$head/100   # respuesta (en cientos de pies)
scattersplot(x, y, main = "Wolfcamp aquifer data", xlab = 'lon', ylab ='lat', 
             zlab = 'piezometric-head', col = jet.colors(128))

Inferencia (no paramétrica) en procesos geoestadísticos

A partir de los valores observados $\mathbf{Y}=(Y(\mathbf{x}_1),\ldots,Y(\mathbf{x}_n))^{t}$ puede interesar:

Aunque hay una gran cantidad de paquetes de R con herramientas de utilidad en geoestadística (geoR, gstat, fields, ...) o en estadística no paramétrica (mgcv, KernSmooth, sm, ...), la mayoría presentan una funcionalidad muy limitada para estadística espacial no paramétrica o resulta difícil emplearlos para implementar nuevos métodos (además de otras limitaciones, como su eficiencia computacional o el número de dimensiones que soportan). Por este motivo se decidió desarrollar el paquete npsp:

La idea sería:

Binning lineal

Para acelerar los cálculos se emplea binning (WARPing; e.g. Wand and Jones, 1995):

Normalmente el resultado final del análisis es una rejilla y el último paso sólo interesa en cálculos intermedios (e.g. obtención de residuos).

Discretización

En la estimación polinómica local se emplea binning lineal.

Las principales funciones son:

Normalmente se emplearan dos rejillas, una de baja resolución para acelerar los cálculos durante el proceso de modelado (selección de la ventana, cálculos intermedios, ...) y otra con la resolución deseada para obtener los resultados finales.

cpu.time(reset=TRUE)
bin <- binning(x, y, nbin = c(41,41), set.NA = TRUE) 
simage(bin, main = 'Binning averages')
# Data points
points(bin$data$x, col = 'darkgray')
# Binning grid
coordvs <- coordvalues(bin)
abline(v = coordvs[[1]], lty = 3)
abline(h = coordvs[[2]], lty = 3)
cpu.time(total = FALSE)

Interpolación

En el paquete se emplea interpolación lineal (algoritmo análogo al binning lineal). Métodos:

nbin.hd <- c(128, 128)
bin.hd <- binning(x, y, nbin = nbin.hd)
summary(abs( y - interp(bin.hd, newx = x)$y )/ y)
cpu.time(total = FALSE)

Estimación de la densidad

La estimación de la densidad está implementada como un caso particular de regresión (se suavizan los pesos binning reescalados):

Por defecto emplea el estimador lineal local (establecer degree = 0 para Nadaraya-Watson).

den <- np.den(bin, h = diag(60, 2), degree = 0) 
# alternatively, set h = h.cv(as.bin.den(bin))$h
plot(den, main = 'Estimated log(density)')

Se puede emplear por ejemplo para establecer una máscara que filtre las posiciones alejadas de los datos.

bin <- mask(bin, mask = log(den$est) > -14)
str(bin)
cpu.time(total = FALSE)

Estimación de la tendencia

Regresión polinómica local

En el caso univariante, para cada $x_{0}$ se ajusta un polinomio: $$\beta_{0}+\beta_{1}\left(x - x_{0}\right) + \cdots + \beta_{p}\left( x-x_{0}\right)^{p}$$ por mínimos cuadrados ponderados, con pesos $w_{i} = \frac{1}{h}K\left(\frac{x-x_{0}}{h}\right)$.

Habitualmente se considera:

El caso multivariante es análogo. La estimación lineal local multivariante $\hat{\mu}{\mathbf{H}}(\mathbf{x})=\hat{\beta}{0}$ se obtiene al minimizar: $$\min_{\beta_{0},\boldsymbol{\beta}{1}}\sum{i=1}^{n} \left( Y(\mathbf{x}{i})-\beta{0}-{\boldsymbol{\beta}}{1}^{t} (\mathbf{x}{i}-\mathbf{x})\right)^{2} K_{\mathbf{H}}(\mathbf{x}_{i}-\mathbf{x}),$$ donde:

Explícitamente: $$\hat{\mu}{\mathbf{H}}(\mathbf{x}) = \mathbf{e}{1}^{t} \left( \mathbf{X}{\mathbf{x}}^{t} {\mathbf{W}}{\mathbf{x}} \mathbf{X}{\mathbf{x}} \right)^{-1} \mathbf{X}{\mathbf{x}}^{t} {\mathbf{W}}{\mathbf{x}}\mathbf{Y} \equiv {s}{\mathbf{x}}^{t}\mathbf{Y},$$ donde $\mathbf{e}{1} = \left( 1, \cdots, 0\right)^{t}$, $\mathbf{X}{\mathbf{x}}$ es la matriz con $(1,(\mathbf{x}{i}-\mathbf{x})^{t})$ en fila $i$, y $\mathbf{W}{\mathbf{x}} = \mathtt{diag} \left( K_{\mathbf{H}}(\mathbf{x}{1} - \mathbf{x}), ..., K{\mathbf{H}}(\mathbf{x}_{n}-\mathbf{x}) \right)$ es la matriz de pesos.

Se puede pensar que se obtiene aplicando un suavizado lineal a $(\mathbf{x}i, Y(\mathbf{x}_i))$: $$\hat{\boldsymbol{\mu}} = \mathbf{SY},$$ siendo $\mathbf{S}$ la matriz de suavizado con $\mathbf{s}{\mathbf{x}_{i}}^{t}$ en la fila $i$.

Implementación en el paquete npsp

La estimación polinómica local está implementada en la función genérica locpol:

Ejemplo

lp <- locpol(bin, h = diag(100, 2), hat.bin = TRUE)  
# Perspective plot with a color scale
spersp(lp, main = 'Trend estimates', zlab = 'piezometric-head', theta = 120)   
cpu.time(total = FALSE)

Estimación del variograma

Al igual que en la geoestadística paramétrica tradicional, el modelado de la dependencia se realiza a partir de los residuos: $$e(\mathbf{x}{i})=Y(\mathbf{x}{i})-\hat{\mu}(\mathbf{x}{i}).$$ Si la media se supone cte. (proceso estacionario) se trabaja directamente con las observaciones (se sustituiría $e(\mathbf{x}{i})$ por $Y(\mathbf{x}_{i})$ en las expresiones).

La estimación se realiza en dos pasos, en primer lugar se obtiene una estimación piloto del variograma y posteriormente se ajusta un modelo válido.

Estimación piloto del variograma

En los métodos tradicionales se trata como un caso particular de regresión: $$\gamma\left( \mathbf{x}{i}-\mathbf{x}{j}\right) =\frac{1}{2}E\left( \varepsilon(\mathbf{x}{i})-\varepsilon(\mathbf{x}{j})\right)^{2}$$ con $N=\frac{n(n-1)}{2}$ observaciones: $${\left(\mathbf{x}{i}-\mathbf{x}{j}, (e(\mathbf{x}{i})-e(\mathbf{x}{j}))^2 / 2 \right)},$$ asumiendo que la variabilidad de los residuos es (similar a) la de los errores.

La estimación polinómico local $\hat{\gamma}{\mathbf{G}}(\mathbf{u}) = \hat{\beta}{0}$ se obtendría minimizando: $$\begin{aligned} \min_{\beta_{0}, \boldsymbol{\beta}{1}, \cdots} \sum{i=1}^{n}\left( \frac{1}{2} \left( e(\mathbf{x}{i})-e(\mathbf{x}{j}) \right)^{2} - \beta_{0} - {\boldsymbol{\beta}}{1}^{t} (\mathbf{x}{i} - \mathbf{x}{j} - \mathbf{u}) - \cdots\right) ^{2}\times & \ K{\mathbf{G}}(\mathbf{x}{i}-\mathbf{x}{j}-\mathbf{u}) & \end{aligned},$$ siendo $\mathbf{G}$ la correspondiente matriz de ventanas.

Actualmente en el paquete npsp solo está implementado el caso isotrópico:

Estas funciones son adecuadas para procesos estacionarios y se podrían emplear en el caso no estacionario para estimar el variograma residual (sesgado) calculando previamente los residuos.

lp.resid <- residuals(lp)
esvar0 <- np.svariso(x, lp.resid, nlags = 50, h = 50)
plot(esvar0, main = "Nonparametric (residual) pilot semivariogram")
cpu.time(total = FALSE)

Estimación piloto del variograma con corrección de sesgo

El uso directo de los residuos $\boldsymbol{\hat{\varepsilon}} = \mathbf{Y} - \hat{\boldsymbol{\mu}} = (\mathbf{I - S})\mathbf{Y}$ introduce un sesgo en la estimación del variograma, ya que su variabilidad es distinta de la de los errores: $$Var(\boldsymbol{\hat{\varepsilon}}) = {\boldsymbol{\Sigma}}_{\hat{\varepsilon}} = \boldsymbol{\Sigma} + \mathbf{B},$$ siendo $\mathbf{B} = \mathbf{S} \boldsymbol{\Sigma} \mathbf{S}^{t} - \boldsymbol{\Sigma} \mathbf{S}^{t} - \mathbf{S} \boldsymbol{\Sigma}$ la matriz de sesgos.

En términos del variograma: $$Var\left(\hat{\varepsilon}(\mathbf{x}i) - \hat{\varepsilon}(\mathbf{x}_j) \right) = Var\left(\varepsilon(\mathbf{x}_i) - \varepsilon(\mathbf{x}_j) \right) + b{ii} + b_{jj} - 2 b_{ij}.$$ Este sesgo es normalmente negativo y mayor en saltos grandes (ver p.e. Cressie, 1993, sección 3.4.3). En Fernandez-Casal y Francisco-Fernandez (2013) se propone un método iterativo para su corrección. La función np.svariso.corr implementa un algoritmo similar (totalmente no paramétrico). Comenzando con el estimador lineal local residual anterior, en cada iteración $k$ se actualizan las diferencias $(\hat{\varepsilon}(\mathbf{x}i)- \hat{\varepsilon}(\mathbf{x}_j))^2$ reemplazándolas por: $$(\hat{\varepsilon}(\mathbf{x}_i)- \hat{\varepsilon}(\mathbf{x}_j))^2 - \hat{b}{ii}^{(k-1)} - \hat{b}{jj}^{(k-1)} + 2 \hat{b}{ij}^{(k-1)},$$ siendo $\hat{\mathbf{B}}^{(k-1)}$ la aproximación del sesgo obtenida en la iteración anterior.

Ejemplo:

esvar <- np.svariso.corr(lp, nlags = 50, h = 50, plot = TRUE)
cpu.time(total = FALSE)

Ajuste de un modelo no paramétrico

Los estimadores piloto no verifican las propiedades de un variograma (condicionalmente semidefinidos negativos) y no pueden ser usados en predicción espacial (kriging). Para resolver este problema se ajusta un modelo válido.

La función fitsvar.sb.iso permite ajustar un modelo no paramétrico de Shapiro-Botha. En el caso isotrópico son de la forma: $$\gamma(\left\Vert \mathbf{u} \right\Vert ) = \nu_{0} - \sum\limits_{k=1}^{K}\kappa_{d}(x_{k}\left\Vert \mathbf{u}\right\Vert ) z_{k},$$ siendo:

El ajuste por WLS a un conjunto de estimaciones piloto se puede realizar fácilmente mediante programación cuadrática (modificación de la función solve.QP del paquete quadprog para obtener la solución de problemas no estrictamente convexos).

Estos modelos son muy flexibles por lo que es habitual considerar una dimensión $d$ mayor que la de los datos para obtener ajustes más suaves (empleando el parámetro dk). En el límite obtendríamos $\kappa_{\infty}(x)\equiv e^{-x^{2}}$, que se correspondería con un modelo válido en cualquier dimensión (se selecciona estableciendodk = 0).

Ejemplo:

svm <- fitsvar.sb.iso(esvar, dk = 0)
plot(svm, main = "Nonparametric semivariogram and fitted model")
svm0 <- fitsvar.sb.iso(esvar0, dk = 0) 
with(svm0$fit, lines(u, fitted.sv, lty = 3, lwd = 2))
cpu.time(total = FALSE)

Estos modelos son extensibles al caso anisotrópico (Fernandez-Casal et al., 2003), adecuado para procesos espacio-temporales.

Selección de las ventanas

La función genérica h.cv permite seleccionar la ventana de una estimación polinómico local (de la tendencia, densidad o variograma) usando criterios de validación cruzada (CV), validación cruzada generalizada (GCV) o MASE (estándar y modificados).

Estimación de la tendencia

Con los criterios de selección de la ventana diseñados para datos independientes se tiende a infrasuavizar la estimación (e.g. Opsomer et al, 2001):

Por ello se han propuesto distintos criterios alternativos para el caso de datos dependientes:

Adicionalmente, reescribiendo el criterio MCV de la forma: $$MCV(\mathbf{H}) = \frac{1}{n} \left( \mathbf{Y}-\mathbf{S}{-N} \mathbf{Y}\right) ^{t} \left( \mathbf{Y}-\mathbf{S}{-N} \mathbf{Y} \right),$$ donde $N = \left\lbrace N(1),\ldots, N(n) \right\rbrace$, puede verse que: $$\mathrm{E}\left( MCV(\mathbf{H})\right) \simeq MASE(\mathbf{H}) + \sigma^{2} - \frac{2}{n}tr\left(\mathbf{S}_{-N}\boldsymbol{\Sigma}\right),$$ siendo: $$MASE(\mathbf{H})= \frac{1}{n}\mathbb{E}\left( \left( \mathbf{SY}-\boldsymbol{\mu}\right)^{t} \left( \mathbf{SY}-\boldsymbol{\mu}\right) \right) = \ = \frac{1}{n}\left(\mathbf{S}\boldsymbol{\mu}-\boldsymbol{\mu}\right)^{t} \left( \mathbf{S}\boldsymbol{\mu}-\boldsymbol{\mu}\right) + \frac{1}{n}tr\left( \mathbf{S}\boldsymbol{\Sigma }\mathbf{S}^{t}\right).$$

#   Example:
bin2 <- binning(x, y)     # sets nbin automatically
h.cv(bin2, h.start = c(100, 100))                # MCV (ncv = 2)
# cov.bin <-  varcov(svm2, coords = coords(bin2))
lp.h <- h.cv(bin2, cov.bin = svm, h.start = c(100, 100)) # CMCV (ncv = 2)
lp.h
cpu.time(total = FALSE)

Hay que tener cuidado con los algoritmos de optimización automáticos. Por defecto se emplea el método "L-BFGS-B" de la función optim. Como se trata de una minimización no lineal multidimensional puede haber problemas de mínimos locales. Para intentar evitar este problema se puede establecer unos valores iniciales para los parámetros de la ventana mediante el argumento h.start. Alternativamente se puede establecer el parámetro DEalgorithm = TRUE para utilizar el algoritmo genético de optimización del paquete DEoptim (aunque puede aumentar considerablemente el tiempo de computación).

Empleando la ventana seleccionada se podría reestimar de nuevo la tendencia, aunque también habría que volver a estimar el variograma (si se modifica la variabilidad de gran escala, habría que actualizar la variabilidad de pequeña escala). Este procedimento se podría repetir iterativamente. Se tratará más adelante en la sección "Estimación conjunta (automática)".

Estimación del variograma

La ventana para la estimación del variograma puede seleccionarse por ejemplo minimizando el correspondiente error cuadrático de validación cruzada: $$\sum\limits_{i=1}^{n-1}\sum\limits_{j=i+1}^{n} \left( \left( \hat{\varepsilon}(\mathbf{x}{i}) - \hat{\varepsilon}(\mathbf{x}{j})\right)^{2} - 2\hat{\gamma}{-(i,j)}\left( \left\Vert \mathbf{x}{i} - \mathbf{x}{j}\right\Vert \right) \right) ^{2},$$ siendo $\hat{\gamma}{-(i,j)}(\cdot )$ la estimación obtenida descartando $\left( \hat{\varepsilon}(\mathbf{x}{i}) - \hat{\varepsilon}(\mathbf{x}{j})\right) ^{2}$. Sin embargo, para tener un mejor ajuste cerca del origen, suele ser preferible emplear el error cuadrático relativo de validación cruzada: $$\sum\limits_{i=1}^{n-1}\sum\limits_{j=i+1}^{n}\left( \frac{\left( \hat{ \varepsilon}(\mathbf{x}{i})-\hat{\varepsilon}(\mathbf{x}{j})\right) ^{2} }{2\hat{\gamma}{-(i,j)}\left( \left\Vert \mathbf{x}{i}-\mathbf{x} _{j}\right\Vert \right) }-1\right) ^{2}.$$ Adicionalmente, si hay datos atípicos, sería recomendable emplear un error absoluto en lugar de cuadrático.

Ninguno de estos criterios tiene en cuenta la dependencia entre las semivarianzas, por lo que las ventanas seleccionadas tenderán a infrasuavizar (especialmente en saltos grandes). Si se emplea un procedimiento automático, puede ser recomendable incrementar la ventana obtenida con estos criterios (en el algoritmo automático implementado en np.fitgeo descrito a continuación, si no se especifica la ventana para la estimación del variograma, se incrementa por defecto la ventana obtenida con h.cv en un 50%).

NOTA: Actualmente sólo está implementada la estimación con ventana global. En este caso, como las semivarianzas no son homocedásticas, sería preferible emplear ventanas locales.

Estimación conjunta (automática)

Como se comentó en la sección anterior, para poder seleccionar una ventana "óptima" para la estimación de la tendencia es necesario disponer de una estimación del variograma. Sin embargo, la estimación de la dependencia requiere de una estimación de la tendencia. Para resolver este problema circular se puede emplear un algoritmo iterativo (similar al método de estimación paramétrica tradicional propuesto por Neuman y Jacobson, 1984):

  1. Seleccionar una ventana inicial $\mathbf{H}^{(0)}$ para $\hat{\mu}$ (usando por ejemplo $MCV$).

  2. Para $k\geq 1$, usando la ventana $\mathbf{H}^{(k-1)}$ calcular $\hat{\boldsymbol{\mu}} = \mathbf{S}\mathbf{Y}$ y los correspondientes residuos $\hat{\boldsymbol{\varepsilon}}$.

  3. Empleando el algoritmo para la corrección de sesgo, obtener una estimación piloto $\tilde{\gamma}(\cdot)$, y ajustar un modelo S-B.

  4. Construir $\hat{\Sigma}^{(k)}$ y seleccionar una nueva ventana $\mathbf{H}^{(k)}$ (empleando $CMCV$ o $CGCV$).

  5. Repetir los pasos 2 al 4 hasta obtener convergencia.

Este algoritmo está implementado en la función np.fitgeo y normalmente dos iteraciones (el valor por defecto) son suficientes. Por ejemplo, se podría hacer un modelado automático ejecutando:

geom <- np.fitgeo(x, y, svm.resid = TRUE, h.start = c(100, 100))
cpu.time(total = FALSE)

El parámetro iter controla el número máximo de iteraciones del algoritmo completo. Estableciendo iter = 0 se emplea el variograma residual:

geom0 <- np.fitgeo(x, y, iter = 0, h.start = c(100, 100))

La rutina devuelve el modelo no paramétrico ajustado, un objeto np.geo que contiene las estimaciones de la tendencia y del variograma. Al representarlo gráficamente:

plot(geom)

Además de las estimaciones finales de la tendencia y del variograma, se muestra también el variograma ajustado obtenido en la iteración anterior (y si svm.resid = TRUE el variograma residual). Si no hay grandes diferencias no habría necesidad de seguir iterando, ya que el criterio de selección de la ventana proporcionaría un suavizado similar.

Predicción kriging

La función genérica np.kriging permite obtener predicciones mediante kriging residual (combinando la estimación de la tendencia con el kriging simple de los residuos). Actualmente solo está implementado kriging simple y kriging residual con vecindario global. Para la resolución del sistema de ecuaciones kriging, a partir de la factorización de Choleski de la matriz de covarianzas, se emplea el paquete spam para matrices dispersas.

Como ejemplo, compararemos las predicciones kriging obtenidas con el ajuste automático con las correspondientes a la estimación con el variograma residual.

krig.grid0 <- np.kriging(geom0, ngrid = c(96, 96)) # 9216 predicciones
krig.grid <- np.kriging(geom, ngrid = c(96, 96))
cpu.time(total = FALSE)
old.par <- par(mfrow = c(1,2), omd = c(0.01, 0.9, 0.05, 0.95))
scale.range <- range(krig.grid0$trend, krig.grid$trend, finite = TRUE)
simage( krig.grid0, 'trend', slim = scale.range, 
        main = 'Trend estimates (residuals)', col = jet.colors(256), legend = FALSE)
simage( krig.grid, 'trend', slim = scale.range, 
        main = 'Trend estimates (iteration 2)', col = jet.colors(256), legend = FALSE)
par(old.par)
splot(slim = scale.range, col = jet.colors(256), add = TRUE)


old.par <- par(mfrow = c(1,2), omd = c(0.01, 0.9, 0.05, 0.95))
scale.range <- range(krig.grid0$kpred, krig.grid$kpred, finite = TRUE)
simage( krig.grid0, 'kpred', slim = scale.range, 
       main = 'Kriging predictions (residuals)', col = jet.colors(256), legend = FALSE)
simage( krig.grid, 'kpred', slim = scale.range, 
       main = 'Kriging predictions (iteration 2)', col = jet.colors(256), legend = FALSE)
par(old.par)
splot(slim = scale.range, col = jet.colors(256), add = TRUE)

Hay que tener en cuenta que se trata de un método robusto y no depende demasiado del modelo de tendencia (la variabilidad no explicada por la tendencia es capturada por la variabilidad de pequeña escala). En este caso es de especial importancia la estimación del variograma cerca del origen.

Sin embargo, si se comparan las varianzas kriging (o se emplean otros métodos de inferencia, como por ejemplo bootstrap) puede haber grandes diferencias entre los resultados obtenidos con distintos modelos.

old.par <- par(mfrow = c(1,2), omd = c(0.01, 0.9, 0.05, 0.95))
scale.range <- range(krig.grid0$ksd, krig.grid$ksd, finite = TRUE)
simage( krig.grid0, 'ksd', slim = scale.range, 
        main = 'Kriging sd (residuals)', col = hot.colors(256), legend = FALSE)
with(aquifer, points(lon, lat, cex = 0.75))
simage( krig.grid, 'ksd', slim = scale.range, 
        main = 'Kriging sd (iteration 2)', col = hot.colors(256), legend = FALSE)
with(aquifer, points(lon, lat, cex = 0.75))
par(old.par)
splot(slim = scale.range, col = hot.colors(256), add = TRUE)

cpu.time()

Futuras implementaciones

Referencias



Try the npsp package in your browser

Any scripts or data that you put into this service are public.

npsp documentation built on July 2, 2019, 9:08 a.m.