Introducción

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: procesos geoestadísticos¹.

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

Por ejemplo podemos utilizar la función (genérica) spoints() para representar los datos (genera un gráfico de dispersión de las posiciones con una escala de color para la respuesta). En este trabajo utilizaremos como ejemplo dos conjuntos de datos suministrado con el paquete npsp:

• aquifer: nivel del agua subterránea (pies sobre el nivel del mar) en 85 localizaciones del acuífero Wolfcamp.

  x <- aquifer[,1:2]      # coordenadas espaciales
  y <- aquifer$head/100   # respuesta (en cientos de pies)
  spoints(x$lon, x$lat, y, 
          main = "Datos del acuífero Wolfcamp",
          xlab = "Longitud", ylab = "Latitud", 
          legend.lab = "nivel piezométrico (cientos de pies)")

  

• precipitation: total de precipitaciones (raíz cuadrada de pulgadas de lluvia) durante marzo de 2016 registradas en 1053 localidades de la parte continental de Estados Unidos (objeto de clase sp::SpatialGridDataFrame).

  library(sp)
  spoints(precipitation, asp = 1)

  

  attr <- attributes(precipitation)
  border <- attr$border 
  interior <- attr$interior 

scattersplot(x, y, main = "Wolfcamp aquifer data", 
             xlab = 'Longitude', ylab ='Latitude', 
             zlab = 'piezometric-head', col = jet.colors(128))

scattersplot(precipitation)

Estimación de la tendencia

Como ya se comentó, se emplea regresión polinómica local para la estimación de la tendencia (e.g. Fan and Gijbels, 1996). El estimador \(\hat{\mu}_{\mathbf{H}}(\mathbf{x})\) se obtiene aplicando un suavizado polinómico a los datos \(\{(\mathbf{x}_{i},Y(\mathbf{x}_{i})):i=1,\ldots,n\}\).

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)\).

• La estimación en \(x_{0}\) es \(\hat{\mu}_{h}(x_{0})=\hat{\beta}_{0}\).

• Podemos obtener también estimaciones de las derivadas: \(\widehat{\mu_{h}^{r)}}(x_{0}) = r!\hat{\beta}_{r}\).

Habitualmente se considera:

• \(p=0\): Estimador Nadaraya-Watson.

• \(p=1\): Estimador lineal local.

La ventana \(h\) controla el grado de suavizado (el tamaño del vecindario local empleado en los ajustes). Puede ser de utilidad esta aplicación shiny para ilustrar este método (ver p.e. García-Portugués, 2021, Sección 4.1.2; Fernández-Casal y Costa, 2021, Sección 7.1.2).

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:

• \(\mathbf{H}\) es la matriz de ventanas \(d\times d\) (simétrica no singular; controla la forma y el tamaño del vecindario).

• \(K_{\mathbf{H}}(\mathbf{u})=\left\vert \mathbf{H}\right\vert ^{-1}K(\mathbf{H}^{-1}\mathbf{u})\) y \(K\) núcleo multivariante.

De esta forma, las estimaciones en las posiciones de los datos se pueden expresar como una transformación lineal de las observaciones: \[\hat{\boldsymbol{\mu}} = \mathbf{SY},\] siendo \(\mathbf{S}\) la denominada matriz de suavizado (este tipo de métodos también se denominan suavizadores lineales). Esta matriz resulta de utilidad, entre otras cosas, para la selección de la ventana y la estimación de la dependencia.

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

• Emplea una matriz de ventanas completa y un núcleo triweight para el cálculo de los pesos.

• Permite calcular eficientemente la matriz de suavizado \(\mathbf{S}\).

• Código optimizado para minimizar el tiempo de computación y los requerimientos de memoria (especialmente para validación cruzada, estimación del variograma, …). Los cálculos se realizan en FORTRAN empleando la librería LAPACK.

Por ejemplo, en el caso de los datos de precipitación, podemos obtener la estimación de la tendencia mediante los comandos:

x <- sp::coordinates(precipitation)
y <- precipitation$y
lp <- locpol(x, y, nbin = c(120, 120), h = diag(c(5, 5)))

Podemos representar las estimaciones de la tendencia con la función simage() (o spersp()). Además, en este caso enmascaramos previamente las estimaciones fuera de la parte continental de EE.UU..

lp <- mask(lp, window = border, warn = FALSE)
# Representar
slim <- range(y)
col <- jet.colors(256)
simage(lp, slim = slim, main = "Trend estimates",
       col = col, xlab = "Longitude", ylab = "Latitude", asp = 1)
plot(border, border = "darkgray", lwd = 2, add = TRUE)
plot(interior, border = "lightgray", lwd = 1, add = TRUE)

Para acelerar los cálculos se emplea binning lineal (e.g. Wand and Jones, 1995), el resultado final es una rejilla. Para obtener estimaciones en otras posiciones se puede emplear las funciones predict() e interp() para interpolar².

Por ejemplo, en lugar del código anterior podríamos emplear (seleccionando de paso una ventana más pequeña para comparar el resultado):

bin <- binning(x, y, nbin = c(120, 120), window = border)
# lp <- locpol(bin, h = diag(c(5, 5)))
lp2 <- locpol(bin, h = diag(c(2, 2)))
simage(lp2, slim = slim, main = "Trend estimates",
       col = col, xlab = "Longitude", ylab = "Latitude", asp = 1)
plot(border, border = "darkgray", lwd = 2, add = TRUE)
plot(interior, border = "lightgray", lwd = 1, add = TRUE)

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.

bin <- binning(x, y, window = border)
lp0.h <- h.cv(bin)$h
lp0 <- locpol(bin, h = lp0.h, hat.bin = TRUE) 

Estimación del variograma

Al igual que en la geoestadística paramétrica tradicional, 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.

Si la media se supone constante (proceso estacionario), se puede obtener una estimación piloto lineal local \(\hat\gamma(\mathbf{u})\) mediante el suavizado de las semivarianzas \(\left(\mathbf{x}_{i}-\mathbf{x}_{j},(Y(\mathbf{x}_i)- Y(\mathbf{x}_j))^2/2 \right)\). La correspondiente ventana puede seleccionarse por ejemplo minimizando el error cuadrático relativo de validación cruzada.

En el caso general, cuando la tendencia no es constante, la aproximación habitual consiste en eliminar la tendencia y estimar el variograma a partir de los residuos \(\mathbf{r}=\mathbf{Y}-\mathbf{SY}\). Sin embargo, la variabilidad de los residuos puede ser muy diferente a la de los verdaderos errores (ver, por ejemplo, Cressie, 1993, sección 3.4.3, para el caso de una tendencia lineal). Como el sesgo debido al uso directo de los residuos en la estimación del variograma puede tener un impacto significativo en la inferencia, se puede considerar una aproximación similar a la descrita en Fernández-Casal y Francisco-Fernández (2014).

En la práctica, las estimaciones piloto del variograma pueden calcularse con las funciones np.svar() y np.svariso.corr(). La función genérica h.cv() puede utilizarse para seleccionar la ventana correspondiente⁴. Actualmente en el paquete npsp solo está implementado el caso isotrópico.

svar.bin <- svariso(x, residuals(lp0), nlags = 60, maxlag = 20)  
svar.h <- h.cv(svar.bin)$h

svar.np <- np.svar(svar.bin, h = svar.h)
svar.np2 <- np.svariso.corr(lp0, nlags = 60, maxlag = 20, 
                            h = svar.h, plot = FALSE) 

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 puede ajustar un modelo no paramétrico de Shapiro-Botha (Shapiro and Botha, 1991) empleando mínimos cuadrados ponderados (mediante programación cuadrática, lo que evita problemas de mínimos locales).

En la practica se pueden ajustar este tipo de modelos empleando la función fitsvar.sb.iso()⁵.

svm0 <- fitsvar.sb.iso(svar.np, dk = 0) 
svm1 <- fitsvar.sb.iso(svar.np2, dk = 0) 
# plot...
plot(svm1, main = "Nonparametric bias-corrected semivariogram\nand fitted models", 
     legend = FALSE, xlim = c(0,max(coords(svar.np2))), 
     ylim = c(0,max(svar.np2$biny, na.rm = TRUE)))
plot(svm0, lwd = c(1, 1), add = TRUE)
plot(svar.np, type = "p", pch = 2, add = TRUE)
# abline(h = c(svm1$nugget, svm1$sill), lty = 3)
# abline(v = 0, lty = 3)
legend("bottomright", legend = c("corrected", 'biased'),
            lty = c(1, 1), pch = c(1, 2), lwd = c(1, 1))

Empleando el modelo ajustado de variograma se podría seleccionar una nueva ventana “óptima” para la estimación de la tendencia. Aunque si cambia la estimación de la tendencia 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 procedimiento se podría repetir iterativamente hasta convergencia (sería equivalente al método de estimación paramétrica tradicional propuesto por Neuman y Jacobson, 1984).

Modelado automático

La función np.fitgeo() implementa un procedimiento automático para la estimación conjunta de la tendencia y el semivariograma. El algoritmo es el siguiente:

1. Seleccionar una ventana inicial para la estimación de la tendencia empleando el criterio MCV (función h.cv()).

2. Calcular la estimación lineal local de la tendencia (función locpol()).

3. Obtener una estimación piloto del semivariograma con corrección de sesgo
   a partir de los residuos correspondientes (función np.svariso.corr()).

4. Ajustar un modelo válido de Shapiro-Botha a las estimaciones piloto del
   semivariograma (función fitsvar.sb.iso()).

5. Seleccionar una nueva matriz de ventanas para la estimación de la tendencia
   con el criterio CGCV, utilizando el modelo ajustado de variograma
   para aproximar la matriz de covarianzas de los datos (función h.cv()
   con objective = "GCV" y cov.bin = modelo de variograma ajustado).

6. Repetir los pasos 2-5 hasta convergencia (por ejemplo, cuando no hay
   cambios significativos en la estimación del variograma, ya que el criterio de selección de la ventana proporcionaría un suavizado similar). En la práctica, normalmente dos iteraciones⁶ suelen ser suficientes (el valor por defecto, iter = 2).

geomod <- np.fitgeo(x, y, nbin = c(30, 30), maxlag = 20, 
                    svm.resid = TRUE, window = border)
plot(geomod)

Establecer corr.svar = TRUE puede ser muy lento (y requerir mucha memoria) cuando el número de datos es muy grande (hay que tener en cuenta también que el sesgo del variograma residual disminuye a medida que aumenta el tamaño de la muestra).

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).

krig.grid <- np.kriging(geomod, ngrid = c(120, 120))
# Plot kriging predictions and kriging standard deviations
old.par <- par(mfrow = c(1,2), omd = c(0.0, 0.98, 0.0, 1))
krig.grid2 <- krig.grid
simage(krig.grid2, 'kpred', main = 'Kriging predictions', slim = slim,
      xlab = "Longitude", ylab = "Latitude" , col = col, asp = 1)
plot(border, border = "darkgray", lwd = 2, add = TRUE)
plot(interior, border = "lightgray", lwd = 1, add = TRUE)
simage(krig.grid2, 'ksd', main = 'Kriging sd',
      xlab = "Longitude", ylab = "Latitude" , col = hot.colors(256), asp = 1)
plot(border, border = "darkgray", lwd = 2, add = TRUE)
plot(interior, border = "lightgray", lwd = 1, add = TRUE)

par(old.par)

Actualmente solo está implementado kriging simple y kriging residual con vecindario global (para el cálculo de las predicciones se usa la factorización de Cholesky de la matriz de covarianzas, implementada en el paquete spam para matrices dispersas). Alternativamente se puede emplear gstat::krige() (o gstat::krige.cv()) junto con la función as.vgm() para kriging local (ver por ejemplo la vignette Kriging with gstat).

Hay que tener en cuenta que las predicciones kriging son robustas y no dependen demasiado del modelo de tendencia (la variabilidad no explicada por la tendencia es capturada por la variabilidad de pequeña escala). Lo importante es la estimación del variograma cerca del origen.

Sin embargo, si se emplean las varianzas kriging (por ejemplo para el cálculo de intervalos de predicción) o otros métodos de inferencia (como contrastes de hipótesis o bootstrap) puede haber grandes diferencias entre los resultados obtenidos con distintos modelos.

Futuras implementaciones y temas abiertos

Bootstrap (condicional e incondicional):

• Francisco-Fernandez M., Quintela-del-Rio A. and Fernandez-Casal R. (2012). Nonparametric methods for spatial regression. An application to seismic events, Environmetrics, 23, 85-93.

• Fernández-Casal R., Castillo-Páez S. and Francisco-Fernández M. (2018). Nonparametric geostatistical risk mapping, Stochastic Environmental Research and Risk Assessment, 32, 675–684.

• Castillo-Páez S., Fernández-Casal R. and García-Soidán P. (2019). A nonparametric bootstrap method for spatial data, Computational Statistics and Data Analysis, 137, 1-15.

Modelos anisotrópicos de semivariograma (caso espacio-temporal):

• Fernández-Casal, R. (2003). Geoestadística espacio-temporal: modelos flexibles de variogramas anisotrópicos no separables. Tesis doctoral (USC). Director: Wenceslao González-Manteiga.

• Fernandez-Casal R., Gonzalez-Manteiga W. and Febrero-Bande M. (2003). Flexible Spatio-Temporal Stationary Variogram Models, Statistics and Computing, 13, 127-136.

Procesos espaciales heterocedásticos:

• Fernández-Casal R., Castillo-Páez S. and García-Soidán P. (2017). Nonparametric estimation of the small-scale variability of heteroscedastic spatial processes, Spatial Statistics, 22, 358-370.

• Castillo-Páez S., Fernández-Casal R. and García-Soidán P. (2020). Nonparametric bootstrap approach for unconditional risk mapping under heteroscedasticity, Spatial Statistics, 40, 100389.

Otros temas abiertos (entre muchos otros):

• Modelos aditivos y parcialmente lineales.

• Kriging local.

• Métodos geoestadísticos no paramétricos para datos masivos.

  • Estimación no paramétrica del variograma mediante FFT en datos reticulares (raster).

  • Estimación del variograma mediante FFT de datos discretizados (binning).

  • Predicción kriging empleando binning.

• Modelos dinámicos (inclusión de variables explicativas cuyo efecto varía con el espacio o el tiempo).

Agradecimientos

Este trabajo ha sido financiado por el MINECO a través del proyecto MTM2017-82724-R, y por la Xunta de Galicia (Grupos de Referencia Competitiva ED431C-2020-14 y Centro Investigación del Sistema universitario de Galicia ED431G 2019/01). Todas estas ayudas fueron cofinanciadas por el FEDER.

Referencias

• Castillo-Páez, S. (2017). Aportaciones a la Geoestadística no paramétrica. Tesis doctoral (UVigo). Directores: Pilar García-Soidán, Rubén Fernández-Casal.

• Chu, C.K. and Marron, J.S. (1991). Comparison of Two Bandwidth Selectors with Dependent Errors. The Annals of Statistics, 19, 1906-1918.

• Cressie, N. (1993). Statistics for Spatial Data. Wiley, New York.

• Fan, J. and Gijbels, I. (1996). Local polynomial modelling and its applications. Chapman & Hall, London.

• Fernandez-Casal R. (2021). npsp: Nonparametric Spatial Statistics. R package version 0.7.8.

• Fernández-Casal R. y Costa J. (2021). Aprendizaje Estadístico. https://rubenfcasal.github.io/aprendizaje_estadistico.

• Fernandez-Casal R. and Francisco-Fernandez M. (2014) Nonparametric bias-corrected variogram estimation under non-constant trend, Stochastic Environmental Research and Risk Assessment, 28, 1247-1259.

• Fernandez-Casal R., Gonzalez-Manteiga W. and Febrero-Bande M. (2003). Space‐time dependency modeling using general classes of flexible stationary variogram models, Journal of Geophysical Research-Atmospheres, 108, 8779, D24.

• Francisco-Fernandez M. and Opsomer J.D. (2005) Smoothing parameter selection methods for nonparametric regression with spatially correlated errors, Canadian Journal of Statistics, 33, 539-558.

• García-Portugués E. (2021). Notes for Nonparametric Statistics. https://bookdown.org/egarpor/NP-UC3M/.

• Opsomer, J. D., Wang, Y. and Yang, Y. (2001). Nonparametric regression with correlated errors. Statistical Science 16, 134–153.

• Rupert D. and Wand M.P. (1994). Multivariate locally weighted least squares regression. The Annals of Statistics, 22, 1346-1370.

• Shapiro A. and Botha J.D. (1991). Variogram fitting with a general class of conditionally non-negative definite functions. Computational Statistics and Data Analysis, 11, 87-96.

• Wand M.P. (1994). Fast Computation of Multivariate Kernel Estimators. Journal of Computational and Graphical Statistics, 3, 433-445.

• Wand M.P. and Jones M.C. (1995). Kernel Smoothing. Chapman and Hall, London.

Ver también referencias en la Sección Futuras implementaciones.