¿Has estudiado estadística espacial?

La estadística espacial tiene aplicaciones en áreas como climatología, ecología, salud y bienes racíces, esto se debe a que en estas áreas es común contar con información espacial que deseamos modelar y es necesario incorporar el efecto de la autocorrelación espacial en la estimación y predicción.

La autocorrelación espacial se puede expresar informalmente como (primera ley de la geografía de Tobler):

Todo esta relacionado con todo, pero las cosas cercana
están más relacionadas que las lejanas.

Usualmente la estadística espacial se divide en tres ramas:

Análisis de datos de variación continua comprende procesos estocásticos cuyos valores pueden ser conocidos en todos los puntos del área de estudio.

El proceso \(Y(s)\) es un vector aleatorio en una ubicación \(s \in \mathbb{R}^r\) donde \(s\) varía de manera continua sobre el dominio \(\mathcal{D} \subset \mathbb{R}^r\) (usualmente r= 2).

¿Qué preguntas interesantes puedes responder?

Análisis de datos de variación discreta (datos sobre áreas): Se refiere a la distribución de eventos cuya localización se asocia a zonas delimitadas por polígonos.

\(\mathcal{D}\) es nuevamente un subconjunto fijo de \(\mathbb{R}^r\), pero en este caso \(\mathcal{D}\) esta particionado en un número finito de unidades de area con fronteras bien definidas.

## Regions defined for each Polygons

Análisis de procesos espaciales de puntos: Conjunto de puntos distribuidos en un área fija cuya ubicación fue generada por un mecanismo estocástico.

\(\mathcal{D}\) es aleatorio y sus índices indican la ubicación de eventos aleatorios que son el proceso puntual.

Supongamos \(Z_1,...,Z_n\) variables aleatorias i.i.d. Gaussianas, con media \(\mu\) y varianza conocida \(\sigma_0^2\).

Entonces: \[\hat{\mu}=\frac{1}{n}\sum_{i=1}^nZ_i\] es un estimador insesgado para \(\mu\) y \[Var(\hat{\mu})=\frac{\sigma_0^2}{n}\]
Supongamos ahora que \(Z_1,...,Z_n\) son v.a. Gaussianas con media \(\mu\) desconocida y varianza \(\sigma_0^2\), supongamos además que las variables están correlacionadas positivamente: \[Cov(Z_i, Z_j)=\sigma_0^2\rho^{|i-j|}\] donde \(0<\rho<1\).

Sea \(\hat{\mu}=\frac{1}{n}\sum_{i=1}^nZ_i\). ¿\(\hat{\mu}\) es un estimador insesgado para \(\mu\)?

Sí
No
No se puede saber

\(Var(\hat{\mu})\) es

Menor a \(\frac{\sigma_0^2}{n}\).
Mayor a \(\frac{\sigma_0^2}{n}\).
Igual a \(\frac{\sigma_0^2}{n}\).
No se puede saber.

Análisis de datos de variación continua

\(Y(s)\) es un vector aleatorio en una ubicación \(s \in \mathbb{R}^r\) donde \(s\) varía de manera continua sobre el dominio \(\mathcal{D} \subset \mathbb{R}^r\) (usualmente r= 2).

Si consideramos un número finito de ubicaciones \(\{s_1,...,s_n\} \in \mathcal{D}\) entonces \((Y(s_1),...,Y(s_n))\) es un vector aleatorio de dimensión \(n\) cuya distribución debe reflejar la dependencia espacial entre las variables. Si observamos datos \(y_1,...,y_n\) en las ubicaciones \(s_1,...,s_n\) entonces \(y_1,...,y_n\) es una realización del vector aleatorio \((Y(s_1),...,Y(s_n))\).

Ahora, para especificar la distribución del proceso espacial utilizaremos las distribuciones finito dimensionales: \[p_{s_1,...,s_n}(y_1,...,y_n)\] para toda \(n \geq 1\) y para todo \(s_1,...,s_n \in \mathcal{D}\).

Las distribuciones finito-dimensionales definen una distribución válida sobre el proceso estocástico espacial \(\{Y(s):s\in \mathcal{D}\}\) si satisfacen las siguientes condiciones (compatibilidad de Kolmogorov):

Son invariantes ante permutaciones.
Son consistentes bajo marginalización.

Un ejemplo de un proceso espacial que satisface las propiedades anteriores es el proceso Gaussiano.

\(\{Y(s):s\in \mathcal{D}\}\) es un proceso Gaussiano si para toda \(n \geq 1\) y \(s_1,...,s_n \in \mathcal{D}\) la distribución de \((Y(s_1),...,Y(s_n))\) es normal multivariada.

Para especificar un proceso Gausiano sólo hace falta especificar las medias y las matrices de covarianza para las distribuciones finito dimensionales. Más aún en el caso del proceso Gaussiano la condición de consistencia de Kolmogorov se reduce a que la matriz de covarianzas sea semidefinida positiva.

En general para especificar un proceso espacial \(Y(s)\) lo descomponemos en dos partes (primera ley de la geoestadística):

\[Y(s) = \mu(s) + \eta(s)\]

donde

\(\mu(s)\) es la media de \(Y(s)\). Es una función determinística de \(s\) y se le llama la tendencia espacial. Explica la variación a gran escala en el proceso espacial \(Y(s)\).
\(\eta(s)\) explica la varianza a pequeña escala del proceso \(Y(s)\). Tiene media cero en cada \(s\) y explica la dependencia espacial en \(Y(s)\) a través de la función de covarianza.

Estimación de \(\mu(s)\)

Suponemos que \(\mu(s)\) es una función suave, continua que varía espacialmente. Podemos modelar \(\mu(s)\) usando polinomios locales, splines o cualquier función de covariables y parámetros.

Las covariables pueden incluir las coordenadas geográficas y variables que expliquen la variabilidad espacial (por ejemplo, si analizamos contaminación del aire podemos usar temperatura o dirección del aire.)

Adicionalmente se pueden añadir otras covariables que contribuyan a explicar la variación en el proceso (por ejemplo, si estamos analizando concentración de contaminantes en el aire podemos incluir covariables como temperatura o dirección del viento).

Veamos un ejemplo donde se midió el pH del agua en 250 sitios.

soil <- read.table("data/Soil_data.txt", header=T)
head(soil)

##   x  y elevation sand silt clay pH_water pH_KCl calcium magnesium
## 1 0  0       578    9   26   43      5.8    4.9     3.6       0.8
## 2 0  5       578    9   26   42      5.9    4.9     3.4       0.8
## 3 0 10       578    9   25   41      5.9    4.9     3.7       0.9
## 4 0 15       579   11   28   40      5.8    4.9     3.7       0.8
## 5 0 20       579    9   27   41      5.8    5.0     4.2       0.9
## 6 0 25       579    8   27   43      6.1    5.1     4.1       0.9
##   potassium aluminum hydrogen carbon nitrogen cation_exchange sulfur
## 1      0.50        0      3.1    1.2     0.12             8.0    4.9
## 2      0.44        0      3.2    1.1     0.12             7.8    4.6
## 3      0.59        0      2.5    1.2     0.13             7.7    5.2
## 4      0.52        0      3.5    1.3     0.12             8.5    5.0
## 5      0.56        0      3.4    1.4     0.13             9.1    5.7
## 6      0.56        0      3.1    1.2     0.13             8.7    5.6
##   volume mass   NC CEC CN
## 1     61    0 1.05 4.9 10
## 2     59    0 1.27 4.6  9
## 3     67    0 0.29 5.2  9
## 4     59    0 1.42 5.0 10
## 5     62    0 1.02 5.7 10
## 6     64    0 0.75 5.6  9

# Vamos a modelar el PH del agua
class.ph <- classIntervals(soil$pH_water, 9, style = "fixed", 
  fixedBreaks = seq(min(soil$pH_water), max(soil$pH_water), length = 10))
soil$ph <- cut(soil$pH_water, class.ph$brks, include.lowest = TRUE)

Comenzamos con un análisis exploratorio:

Gráficas de los datos y de curvas de nivel nos pueden ayudar en la detección de atípicos espaciales.

ggplot(soil, aes(x = x, y = y, color = ph)) + 
  geom_point(size = 2.3) +
  scale_colour_brewer(palette="Blues")

library(akima)
# Interpolar los valores observados (akima)
surf.ph <- interp(soil$x, soil$y, soil$pH_water)

# hacer una gráfica
surf.ph.r <- cbind(expand.grid(surf.ph$x, surf.ph$y), c(surf.ph$z))
colnames(surf.ph.r) <- c("x", "y", "z")
surf.ph.r$ph <- cut(surf.ph.r$z, class.ph$brks)
ggplot(surf.ph.r, aes(x = x, y = y, z = z)) + 
  geom_tile(aes(fill = ph)) + scale_fill_brewer(palette = "Blues") + 
  stat_contour()

Podemos explorar la dependencia espacial de \(\mu(s)\) en las coordenadas geográficas graficando \(Y(s_i)\) contra las coordenadas, o si los datos están muestreados sobre una cuadrícula graficando la media y mediana de la fila contra el ínidce de la fila o columna.

# EDA: gráficas de la media/mediana del pH vs el número de fila y de columna
x.est <- soil %>%
  group_by(x) %>%
  summarise(
    media = mean(pH_water), 
    mediana = median(pH_water)
    ) %>%
  gather(variable, valor, -x)

y.est <- soil %>%
  group_by(y) %>%
  summarise(
    media = mean(pH_water), 
    mediana = median(pH_water)
    ) %>%
  gather(variable, valor, -y)

x <- ggplot(x.est, aes(x = x, y = valor, color = variable)) + 
  geom_point() + geom_line()

y <- ggplot(y.est, aes(x = y, y = valor, color = variable)) + 
  geom_point() + geom_line()

grid.arrange(x, y, ncol = 2, main = "Gráficas exploratorio")

Podemos modelar \(\mu(s)\) como una función lineal de las covariables. Por ejemplo, si modelamos \(\mu(s)\) como una función lineal únicamente de las coordenadas tenemos:

\[\mu(s) = \beta_0 + \beta_1 s_x + \beta_2 s_y + \beta_3 s_x^2 + \beta_4 s_y^2 + \beta_5 s_x s_y\]

Cuando se utiliza un polinomio de las coordenadas geográficas es conveniente utilizar un polinomio completo, esto es, incluir todos los monomios de grado menor o igual al grado del polinomio. Esto garantiza que la superficie que ajustamos es invariante a la elección del origen y la orientación de las coordenadas (transformaciones lineales).
Es común centrar las coordenadas geográficas restando la media.

Supongamos que queremos modelar la tendencia espacial utilizando únicamente información de las coordenadas, sea \(X(s)\) un vector de tamaño \(p \times 1\) que contiene la información geográfica. Entonces, para obtener una estimación provisional de la tendencia espacial de \(Y(s)\) podemos ajustar un modelo de regresión: \[Y(s) = \mu(s) + \epsilon(s) = X(s)^T\beta + \epsilon(s)\] donde \(\epsilon(s)\) son iid. Notemos que si queremos incorporar otras covariables (además de información de las coordenadas), simplemente las añadimos a la matriz \(X(s)\).

Volvamos al ejemplo del ph del agua.

#### Estimación de la tendencia espacial
mod.1 <- lm(pH_water ~ x + y, data = soil)

# La tendencia espacial estimada está dada por los valores ajustados y podemos 
# obtener una superficie utilizando interpolación bilineal o predict
surf.ph.r$z <- predict(mod.1, surf.ph.r[, 1:2])

# Discretizamos pH utilzando los mismos cortes que el plot exploratorio
surf.ph.r$z <- predict(mod.1, surf.ph.r[, 1:2])
surf.ph.r$ph <- cut(surf.ph.r$z, class.ph$brks)
ggplot(surf.ph.r, aes(x = x, y = y, z = z)) + 
  geom_tile(aes(fill = ph)) + 
  scale_fill_brewer(palette = "Blues", drop = FALSE)

# Polinomio de segundo orden
mod.1 <- lm(pH_water ~ 1 + x + y + I(x * y) + I(x ^ 2) + I(y ^ 2), data = soil)
surf.ph.r$z <- predict(mod.1, surf.ph.r)
surf.ph.r$ph <- cut(surf.ph.r$z, class.ph$brks)
ggplot(surf.ph.r, aes(x = x, y = y, z = z)) + 
  geom_tile(aes(fill = ph)) + 
  scale_fill_brewer(palette = "Blues", drop = FALSE)

# Polinomio de tercer orden
mod.1 <- lm(pH_water ~ 1 + x + y + I(x * y) + I(x ^ 2) + I(y ^ 2) + I(x ^2 * y) +
  I(x * y ^ 2) + I(x ^ 3) + I(y ^ 3), data = soil)
surf.ph.r$z <- predict(mod.1, surf.ph.r)
surf.ph.r$ph <- cut(surf.ph.r$z, class.ph$brks)
ggplot(surf.ph.r, aes(x = x, y = y, z = z)) + 
  geom_tile(aes(fill = ph)) + 
  scale_fill_brewer(palette = "Blues", drop = FALSE)

Variación de pequeña escala \(\eta(s)\)

Recordemos que dividimos el proceso espacial en: \(Y(s) = \mu(s)+\eta(s)\) donde \(\mu(s)\) era la tendencia espacial (función determinística de covariables) y \(\eta(s)\) es un proceso espacial en \(\mathcal{D} \subset \mathbb{R}^r\).

Suponemos que \(\eta(s)\) es un proceso intrínisicamente estacionario con media cero, esto es:

\(E(\eta(s)) = 0\) para toda \(s\) en \(\mathcal{D}\)
\(\frac{1}{2}Var(\eta(s) - \eta(s'))=\gamma(s-s')\)

Un proceso intríniscamente estacionario es un proceso que satisface las dos condiciones enumeradas, con media constante (pero no necesariamente cero). Notemos que podemos escribir la segunda condición como \[\frac{1}{2}Var(\eta(s) - \eta(s+h))=\gamma(h)\] y vemos que el lado derecho de la igualdad depende únicamnte de \(h\) (donde \(h \in \mathbb{R}^r\)) y no de la elección particular de \(s\).

La función \(2\gamma(h) = Var(\eta(s) - \eta(s+h))\) donde \(h \in \mathbb{R}^r\) se conoce como variograma, mientras que a \(\gamma(h)\) se le denota semivariograma.

¿Cómo se relaciona el variograma con la función de covarianzas \(C(h)\) (también conocida como covariograma)?

\[ \begin{align} \nonumber 2\gamma(h) &= Var(\eta(s+h) - \eta(s))\\ \nonumber &=Var(\eta(s+h)) + Var(\eta(s)) - 2Cov(\eta(s+h), \eta(s))\\ \nonumber &= C(0) + C(0) -2Cov(\eta(s+h), \eta(s))\\ \nonumber &=2[C(0) - C(h)] \end{align} \]

Vemos entonces que \(\gamma(h) = C(0)-C(h)\) y por tanto si conocemos \(C\) podemos recuperar \(\gamma\); sin embargo, para ir de \(\gamma\) a \(C\) nesesitamos un supuesto adicional: si el proceso espacial cumple \(C(h) \rightarrow 0\) cuando \(\|h\| \to \infty\) (ergódico), tomando limites en la igualdad \(\gamma(h) = C(0)-C(h)\) vemos que \[\lim_{\|h\| \to \infty} \gamma(h) = C(0)\] y podemos recuperar \(C\) de \(\gamma\) utilizando la variable auxiliar \(u\) de la siguiente manera, \[C(h) = C(0) - \gamma(h) = \lim_{\|u\| \to \infty} \gamma(u) - \gamma(h)\] El límite del lado derecho no existe siempre (si existe se decieque el proceso es estacionario de segundo orden), estudiaremos casos en los que existe.

Veamos que el supuesto de ergodicidad es sensible (intuitivamente) pues nos dice que la covarianza entre dos puntos desaparece conforme los puntos se separan en el espacio.

Utilizaremos estos conceptos en la estimación de \(\eta(s)\), primero aproximamos los residuales como \[\hat{\eta}(s_i) = Y(s_i) - \hat{\mu}(s_i)\] para \(i = 1,...,n\) donde \(\hat{\mu}(s)\) es la tendencia estimada.

Semivariograma empírico

Como parte del análisis exploratorio de datos es útil graficar el semivariograma empírico \(\hat{\gamma}(h)\) definido como (si las observaciones están en una cuadrilla regular):

\[\hat{\gamma}(h) = \frac{1}{2N(h)}\sum_{i,j:s_i-s_j = h}(\hat{\eta}(s_i)- \hat{\eta}(s_j))^2\]

donde \(N(h)\) es el número de pares \((s_i, s_j)\) tales que \(s_i-s_j = h\).

Si las observaciones no están en una cuadrilla regular definimos regiones de tolerancia \(T(h)\) (vecindades alrededor de h) para cada \(h\):

\[\hat{\gamma}(h) = \frac{1}{2\#T(h)}\sum_{i,j:s_i-s_j \in T(h)}(\hat{\eta}(s_i)- \hat{\eta}(s_j))^2\]

En la práctica para construir el semivariograma empírico seguimos los siguientes pasos:

Particionamos el espacio de todas las diferencias \(\mathcal{H}=\{s-s': s, s' \in \mathcal{D} \}\) en clases (bins) \(B_m\).
Calculamos las diferencias \(s_i-s_j\) para toda \(i\neq j\) y las asignamos a una clase \(B_m\).
Determinamos \(h_m\) el vector diferencia representativo de la clase \(B_m\), \(h_m\) puede ser el promedio de todas las diferencias \(s_i-s_j\) que caen en la clase \(m\).
Determinamos el número de pares ordenados \((s_i,s_j)\) tales que \(s_i-s_j \in B_m\) y calculamos \(N(h_m)=\#B_m\).
Para \(h_m\) claculamos:

\[\hat{\gamma}(h_m) = \frac{1}{2N(h_m)}\sum_{i,j:s_i-s_j = h}(\hat{\eta}(s_i)- \hat{\eta}(s_j))^2\]

Grafiquemos el semivariograma empírico para el ejmplo de ph.

library(gstat)

## 
## Attaching package: 'gstat'
## 
## The following object is masked from 'package:spatstat':
## 
##     idw

mod.1 <- lm(pH_water ~ x + y, data = soil)
soil$eta.hat <- soil$pH_water - mod.1$fitted.values
# formato para gstat
water.data <- data.frame(x = soil$x, y = soil$y, eta.hat = soil$eta.hat)
coordinates(water.data) = ~ x + y
emp.variog <- variogram(eta.hat ~ 1, water.data)
emp.variog

##      np dist gamma dir.hor dir.ver   id
## 1   465  5.0 0.012       0       0 var1
## 2   432  7.1 0.018       0       0 var1
## 3  1228 10.8 0.023       0       0 var1
## 4   368 14.1 0.032       0       0 var1
## 5  1127 15.5 0.029       0       0 var1
## 6   674 18.0 0.036       0       0 var1
## 7  1946 21.1 0.033       0       0 var1
## 8  1483 25.2 0.033       0       0 var1
## 9   802 27.4 0.036       0       0 var1
## 10 1324 29.8 0.033       0       0 var1
## 11 1380 32.4 0.036       0       0 var1
## 12 1745 35.7 0.032       0       0 var1
## 13  734 38.5 0.034       0       0 var1
## 14 1480 40.7 0.033       0       0 var1

# plot(emp.variog)
# graficar el semivariograma empírico (notemos que estamos asumiendo
# un proceso isotrópico)
# plot(emp.variog, col = "black", type = "p",pch=20)
ggplot(emp.variog, aes(x = dist, y = gamma)) +
  geom_point() + ylim(0, 0.038) +
  labs(title = expression("Semivariograma empírico"), 
       x = "distancia", y = "semivarianza")

Una definición más:

Un proceso intrínsicamente estacionario es isotrópico si su semivariograma \(\gamma(h)\) depende del vector de separación \(h\) únicamente a través de su norma \(\|h\|\).

Entonces, para un proceso istrópico \(\gamma(h)\) es una función univariada que se puede escribir como \(\gamma(\|h\|)\). Los procesos istrópicos son populares debido a su simplicidad y a que existen varias funciones paramétricas que se pueden utilizar para definir el semivariograma del proceso.

Nos podemos preguntar si en el ejemplo es razonable suponer que el proceso es isotrópico.

water.data <- data.frame(x = soil$x, y = soil$y, eta.hat = soil$eta.hat)
coordinates(water.data) = ~ x + y
dir.variog <- variogram(eta.hat ~ 1, water.data, alpha = c(0,45,90,135))
# plot(dir.variog)
# ggplot
dir.variog$direction <- factor(dir.variog$dir.hor, levels = c(0, 45, 90, 135),
  labels = c("E-O", "NE-SO", "N-S", "NO-SE"))
ggplot(dir.variog, aes(x = dist, y = gamma, colour = direction)) + 
  geom_point(size = 1.8) + 
  labs(title = expression("Semivariograma direccional"), 
       x = "distancia", y = "semivariance", colour = "dirección") + geom_line()

# cuando invesigamos anisotropías podemos tener celdas con pocos datos
dir.variog

##     np dist  gamma dir.hor dir.ver   id direction
## 1  240  5.0 0.0085       0       0 var1       E-O
## 2  230 10.0 0.0130       0       0 var1       E-O
## 3  616 15.5 0.0216       0       0 var1       E-O
## 4  588 20.4 0.0253       0       0 var1       E-O
## 5  560 25.3 0.0275       0       0 var1       E-O
## 6  320 26.9 0.0364       0       0 var1       E-O
## 7  532 30.3 0.0275       0       0 var1       E-O
## 8  304 31.6 0.0361       0       0 var1       E-O
## 9  792 35.7 0.0302       0       0 var1       E-O
## 10 748 40.6 0.0314       0       0 var1       E-O
## 11 216  7.1 0.0190      45       0 var1     NE-SO
## 12 399 11.2 0.0241      45       0 var1     NE-SO
## 13 184 14.1 0.0295      45       0 var1     NE-SO
## 14 337 18.0 0.0327      45       0 var1     NE-SO
## 15 460 22.0 0.0317      45       0 var1     NE-SO
## 16 279 25.0 0.0317      45       0 var1     NE-SO
## 17 126 28.3 0.0323      45       0 var1     NE-SO
## 18 250 29.2 0.0307      45       0 var1     NE-SO
## 19 446 32.8 0.0305      45       0 var1     NE-SO
## 20 298 35.8 0.0311      45       0 var1     NE-SO
## 21 367 38.5 0.0299      45       0 var1     NE-SO
## 22 247 41.0 0.0301      45       0 var1     NE-SO
## 23 225  5.0 0.0159      90       0 var1       N-S
## 24 200 10.0 0.0274      90       0 var1       N-S
## 25 511 15.5 0.0374      90       0 var1       N-S
## 26 438 20.4 0.0361      90       0 var1       N-S
## 27 365 25.3 0.0316      90       0 var1       N-S
## 28 230 26.9 0.0333      90       0 var1       N-S
## 29 292 30.3 0.0330      90       0 var1       N-S
## 30 184 31.6 0.0339      90       0 var1       N-S
## 31 357 35.7 0.0325      90       0 var1       N-S
## 32 238 40.6 0.0368      90       0 var1       N-S
## 33 216  7.1 0.0176     135       0 var1     NO-SE
## 34 399 11.2 0.0261     135       0 var1     NO-SE
## 35 184 14.1 0.0336     135       0 var1     NO-SE
## 36 337 18.0 0.0403     135       0 var1     NO-SE
## 37 460 22.0 0.0423     135       0 var1     NO-SE
## 38 279 25.0 0.0462     135       0 var1     NO-SE
## 39 126 28.3 0.0443     135       0 var1     NO-SE
## 40 250 29.2 0.0452     135       0 var1     NO-SE
## 41 446 32.8 0.0436     135       0 var1     NO-SE
## 42 298 35.8 0.0395     135       0 var1     NO-SE
## 43 367 38.5 0.0384     135       0 var1     NO-SE
## 44 247 41.0 0.0345     135       0 var1     NO-SE

Observaciones

¿Cómo elegimos el tamaño de los bins \(B_m\)? ¿Cuántos bins?
Más bins puede ser más informativo; sin embargo, bins con pocas observaciones pueden ser muy ruidosos. Una regla de dedo es definirlos para que contengan al menos 30 observaciones.

set.seed(1362)
soil.2 <- soil[sample(1:nrow(soil), size = 90), ]
water.data <- data.frame(x = soil.2$x, y = soil.2$y, eta.hat = soil.2$eta.hat)
coordinates(water.data) = ~ x + y
emp.variog.1 <- variogram(eta.hat ~ 1, water.data)
emp.variog.1$variogram <- "default"
emp.variog.2 <- variogram(eta.hat ~ 1, water.data, width = 0.1)
emp.variog.2$variogram <- "pocas obs."
emp.variog.3 <- variogram(eta.hat ~ 1, water.data, width = 9)
emp.variog.3$variogram <- "pocos bins"
emp.variogs <- rbind(emp.variog.1, emp.variog.2, emp.variog.3)

ggplot(emp.variogs, aes(x = dist, y = gamma, color = variogram, group = variogram)) +
  geom_point() + 
  geom_line() +
  ylim(0, 0.06) +
  labs(title = expression("Semivariograma empírico"), 
       x = "distancia", y = "semivarianza")

La variabilidad en el valor de \(\hat{\gamma}(h)\) aumenta conforme \(\|h\|\) aumenta, por lo que \(\hat{\gamma}(h)\) es poco confiable para valores grandes de \(\|h\|\). Es aconsejable establecer una cota superior a la magnitud de \(\|h\|\), una regla de dedo es usar 1/2 de la máxima distancia.

Corrobora el punto 2 en el ejemplo de ph.

Semivariogramas paramétricos

Modelamos el semivariograma empírico utilizando semivariogramas paramétricos, hay varias razones por las cuales esto es conveniente:

El semivariograma empírico no esta disponible en distancias que no aparecen en los datos.
El semivariograma empírico puede ser muy ruidoso por lo que una verisón suavizada puede ser más confiable.
El semivariograma empírico puede no cumplir los supuestos necesarios para estimar el proceso.

Hay varias características que debe cumplir el semivariogreama paramétrico, por ejemplo, es una función no-decreciente de la distancia. Es común que el semivariograma empírico se aplane a partir de un cierto punto, indicando que la dependencia espacial disminuye con la distancia, si es el caso el semivariograma paramétrico debe reflejar este escenario. Veamos algunos modelos de semivariogramas paramétricos:

Esférico

\[\gamma(h;\theta) = \left\{ \begin{array}{lr} 0 & : \|h\| = 0 \\ \tau^2 +\sigma^2\big(\frac{3\|h\|}{2\phi}-\frac{\|h\|^3}{2\phi^3}\big) & : 0 \le \|h\| \le \phi\\ \tau^2+\sigma^2 & : \|h\| > 0 \end{array} \right. \] donde \(\sigma^2, \phi>0\). Esta modelo del semivariograma tiene asociada la siguiente función de covarianza:

\[C(h;\theta) = \left\{ \begin{array}{lr} \tau^2 + \sigma^2 & : \|h\| = 0 \\ \sigma^2\left(1-\big(\frac{3\|h\|}{2\phi}-\frac{\|h\|^3}{2\phi^3}\big)\right) & : 0 < \|h\| \le \phi\\ 0 & : \|h\| > \phi \end{array} \right. \]

sph.variog <- function(sigma2, phi, tau2, dist){
  n <- length(dist)
  sph.variog.vec <- rep(0,n)
  for(i in 1:n){
    if(dist[i] < phi){
      sph.variog.vec[i] <- tau2 + (sigma2 * (((3 * dist[i]) / (2 * phi)) - 
          ((dist[i] ^ 3)/(2 * (phi ^ 3)))))
    }
    if(dist[i] >= phi){
      sph.variog.vec[i] <- tau2 + sigma2  
    }
  }
  return(sph.variog.vec)
}

dat <- data.frame(x = seq(0, 10, 1), y = seq(0, 1, 0.1))
ggplot(dat, aes(x = x, y = y)) + 
  labs(title = expression("Semivariograma esférico"), 
       x = "distancia", y = "semivarianza") +
  stat_function(fun = sph.variog, args = list(sigma2 = 1, phi = 1, tau2 = 0), 
    colour = "green3") +
  stat_function(fun = sph.variog, args = list(sigma2 = 1, phi = 8, tau2 = 0), 
    colour = "orange") +
  stat_function(fun = sph.variog, args = list(sigma2 = 0.5, phi = 1, tau2 = 0), 
    colour = "steelblue") +
  stat_function(fun = sph.variog, args = list(sigma2 = 1, phi = 1, tau2 = 0.2), 
    colour = "darkgray")

Un semivariograma paramétrico tiene las siguientes características:

Por definición \(\gamma(0) = 0\). Sin embargo, en ocasiones vemos que el valor en cero del semivariograma empírico es \(\tau^2\neq 0\). En geostadística \[\tau^2=\gamma(0^+)=\lim_{\|h\|\to 0^+} \gamma(t)\] se conoce como el efecto nugget o pepita. El efecto nugget equivale a descomponer \(\eta(s)\) como \[\eta(s)=\eta'(s)+\epsilon\] donde \(\epsilon\) es ruido blanco. La siguiente figura muestra el efecto del nugget en un proceso simulado, en el lado izquierdo \(\tau^2 = 0\) y en el del lado derecho \(\tau^2 = 0.5\).

\(lim_{\|h\|\to \infty} \gamma(\|h\|) = \tau^2 + \sigma^2\), este valor (asintótico normalmente) del semivariograma se conoce como el silo, el silo menos el nugget se conoce como silo parcial y es simplemente \(\sigma^2\). El silo no siempre exite.
El rango es el menor valor de \(\|h\|\) en el que el semivariograma alcanza el silo (\(\phi\)).

En el caso del semivariograma gris del ejemplo de arriba, ¿Cuál es el nugget, silo parcial y el rango?

0.2, 1.2, 1
0.2, 1, 1
0, 1, 1
0, 1.2, 1
Ninguna de las anteriores

Gaussiano

\[\gamma(h;\theta) = \left\{ \begin{array}{lr} 0 & : \|h\| = 0 \\ \tau^2 +\sigma^2\big(1-exp\big(-\frac{\|h\|^2}{\phi^2}\big)\big) & : \|h\| > 0 \end{array} \right. \] donde \(\sigma^2, \phi>0\). Con función de covarianza:

\[C(h;\theta) = \left\{ \begin{array}{lr} \tau^2 + \sigma^2 & : \|h\| = 0 \\ \sigma^2\left(1- \big(1-exp\big(-\frac{\|h\|^2}{\phi^2}\big)\big)\right) & : \|h\| > 0 \end{array} \right. \]

# existen funciones en los paquetes pero escribirlas puede ser ilustrativo
gau.variog <- function(sigma2,phi,tau2,dist){
  n <- length(dist)
  gau.variog.vec <- rep(0,n)
  for(i in 1:n){
    gau.variog.vec[i] <- tau2 + sigma2 * (1-exp(-dist[i]^2/phi^2))
  }
  gau.variog.vec[0] <- 0 
  return(gau.variog.vec)
}
dat <- data.frame(x = seq(0, 10, 1), y = seq(0, 1, 0.1))
ggplot(dat, aes(x = x, y = y)) + 
  labs(title = expression("Semivariograma gaussiano"), 
       x = "distancia", y = "semivarianza") +
  stat_function(fun = gau.variog, args = list(sigma2 = 1, phi = 1, tau2 = 0), 
    colour = "green3") +
  stat_function(fun = gau.variog, args = list(sigma2 = 1, phi = 8, tau2 = 0), 
    colour = "orange") +
  stat_function(fun = gau.variog, args = list(sigma2 = 0.5, phi = 1, tau2 = 0), 
    colour = "steelblue") +
  stat_function(fun = gau.variog, args = list(sigma2 = 1, phi = 1, tau2 = 0.2), 
    colour = "darkgray") + xlim(0, 20)

Matérn

\[\gamma(h;\theta) = \left\{ \begin{array}{lr} 0 & : \|h\| = 0 \\ \tau^2 +\sigma^2 \left(1-\frac{\big(2\sqrt{\nu} \frac{\|h\|}{\phi}\big)^{\nu} \mathcal{K}_\nu\big(2\sqrt{\nu}\frac{\|h\|}{\phi}\big)} {2^{\nu-1}\Gamma(\nu)}\right) & : \|h\| > 0 \end{array} \right. \] donde \(\sigma^2, \phi>0\).

Notemos que esta función tiene el parámetro adicional \(\nu\), este controla el suavizamiento del proceso espacial. El semivariograma exponencial es un caso especial del Matern con \(\nu = 1/2\), el Gaussiano se puede derivar del Matern haciendo \(\nu \to \infty\).

library(geoR)
mat.variog <- function(sigma2, phi, tau2, nu = 1, dist){
  #phi <- 2*phi
  n <- length(dist)
  mat.variog.vec <- rep(0,n)
  for(i in 1:n){ # in fun matern u = 2*sqrt(nu)*dist
    mat.variog.vec[i] <- tau2 + sigma2 * (1 - geoR::matern(2*sqrt(nu)*dist[i], 
      phi, nu)) 
  }
  return(mat.variog.vec)
}
dat <- data.frame(x = seq(0, 10, 1), y = seq(0, 1, 0.1))
mat_1 <- ggplot(dat, aes(x = x, y = y)) + 
  labs(x = "distancia", y = "semivarianza") +
  stat_function(fun = mat.variog, args = list(sigma2 = 1, phi = 1, tau2 = 0, 
    nu = 1/4), colour = "green3") +
  stat_function(fun = mat.variog, args = list(sigma2 = 1, phi = 8, tau2 = 0, 
    nu = 1/4), colour = "orange") +
  stat_function(fun = mat.variog, args = list(sigma2 = 0.5, phi = 1, tau2 = 0, 
    nu = 1/4), colour = "steelblue") +
  stat_function(fun = mat.variog, args = list(sigma2 = 1, phi = 1, tau2 = 0.2, 
    nu = 1/4), colour = "darkgray") 
mat_2 <- ggplot(dat, aes(x = x, y = y)) + 
  labs(x = "distancia", y = "semivarianza") +
  stat_function(fun = mat.variog, args = list(sigma2 = 1, phi = 1, tau2 = 0, 
    nu = 10), colour = "green3") +
  stat_function(fun = mat.variog, args = list(sigma2 = 1, phi = 8, tau2 = 0, 
    nu = 10), colour = "orange") +
  stat_function(fun = mat.variog, args = list(sigma2 = 0.5, phi = 1, tau2 = 0, 
    nu = 10), colour = "steelblue") +
  stat_function(fun = mat.variog, args = list(sigma2 = 1, phi = 1, tau2 = 0.2, 
    nu = 10), colour = "darkgray") + xlim(0, 20)

grid.arrange(mat_1, mat_2, ncol = 2, main = "Semivariograma Matern")

Vimos distintos modelos de variogramas paramétricos, ¿Cuál es la diferencia práctica entre ellos? Producen distintos tipos de dependencia espacial y distintos procesos espaciales.

grid.xy <- expand.grid(x = seq(0, 2, 0.02), y = seq(0, 2, 0.02))

# Esférico
esf_sin <- gstat(formula = z ~ 1, locations = ~ x + y, beta = 1, dummy = TRUE, 
  model = vgm(psill = 1, model = "Sph", range = 1), nmax = 20)
# make 4 simulaciones
set.seed(23)
esf <- predict(esf_sin, newdata = grid.xy, nsim = 4)

## [using unconditional Gaussian simulation]

esf_df <- melt(esf, id.vars = c("x", "y"))
cuts <- seq(min(esf_df$value) - 1, max(esf_df$value) + 1, length = 10)
esf_df$data_cat <- cut(esf_df$value, breaks = cuts, include.lowest = TRUE)
esf_df$tipo <- "Sin nugget"

# Esférico con nugget
esf_con <- gstat(formula = z ~ 1, locations = ~ x + y, beta = 1, dummy = TRUE, 
  model = vgm(psill = 1, model = "Sph", range = 1, nugget = 0.2), nmax = 20)
# hacemos 4 simulaciones
set.seed(23)
esf <- predict(esf_con, newdata = grid.xy, nsim = 4)

## [using unconditional Gaussian simulation]

esf_df_2 <- melt(esf, id.vars = c("x", "y"))
esf_df_2$data_cat <- cut(esf_df_2$value, breaks = cuts, include.lowest = TRUE)
esf_df_2$tipo <- "Con nugget"
esf_df_3 <- rbind(esf_df, esf_df_2)

ggplot(esf_df_3, aes(x = x, y = y, color = data_cat)) +
  geom_point(size = 1.8) +
  scale_colour_brewer(palette = "Blues", drop = FALSE) + 
  facet_grid(tipo ~ variable) +
  labs(title = "Esférico")

# Gaussiano
grid.xy <- expand.grid(x = seq(0, 2, 0.02), y = seq(0, 2, 0.02))
gau_sin <- gstat(formula = z ~ 1, locations = ~ x + y, beta = 1, dummy = TRUE, 
  model = vgm(psill = 1, model = "Gau", range = 1, nugget = 0.0001), nmax = 20)
# make 4 simulations
set.seed(23)
gau <- predict(gau_sin, newdata = grid.xy, nsim = 4)

## [using unconditional Gaussian simulation]

gau_df <- melt(gau, id.vars = c("x", "y"))
cuts <- seq(min(gau_df$value)-0.5, max(gau_df$value)+0.5, length = 10)
gau_df$data_cat <- cut(gau_df$value, breaks = cuts, include.lowest = TRUE)
gau_df$tipo <- "Sin nugget"

# Gaussiano con nugget
gau_con <- gstat(formula = z ~ 1, locations = ~ x + y, beta = 1, dummy = TRUE, 
  model = vgm(psill = 1, model = "Gau", range = 1, nugget = 0.05), nmax = 20)
set.seed(23)
gau <- predict(gau_con, newdata = grid.xy, nsim = 4)

## [using unconditional Gaussian simulation]

gau_df_2 <- melt(gau, id.vars = c("x", "y"))
gau_df_2$data_cat <- cut(gau_df_2$value, breaks = cuts, include.lowest = TRUE)
gau_df_2$tipo <- "Con nugget"

gau_df_3 <- rbind(gau_df, gau_df_2)

ggplot(gau_df_3, aes(x = x, y = y, color = data_cat)) +
  geom_point(size = 1.8) +
  scale_colour_brewer(palette = "Blues", drop = FALSE) + 
  facet_grid(tipo ~ variable) +
  labs(title = "Gaussiano")

Estimación de \(\eta(s)\)

Veremos la estimación de \(\eta(s)\) utilizando técnicas de geosestadística clásica y de máxima verosimilitud.

Geoestadística clásica

Una vez que elegimos el modelo paramétrico que vamos a usar (esférico, gaussiano, etc.) etsimamos los valores de los parámetros usando mínimos cuadrados ponderados (donde se da un mayor peso a observaciones cercanas pues la varianza suele ser menor). \[WSS(\theta) = \sum_{m=1}^k \frac{N(h_m)}{(\gamma(h_m:\theta))^2}(\hat{\gamma}(h_m) -\gamma(h_m:\theta))^2\] El estimador de \(\theta\) se encuentra minimizando de manera iterativa la suma de cuadrados ponderada (\(WSS(\theta)\)).

Ajustemos los semivariogramas paramétricos a nuestros datos usando WLS.

#### Variograma esférico
sph.variog.w <- fit.variogram(emp.variog, vgm(psill= 0.02, model="Sph",
  range = 20, nugget = 0.012),fit.method = 2)
sph.variog.w

##   model   psill range
## 1   Nug 0.00091     0
## 2   Sph 0.03274    21

print(plot(emp.variog, sph.variog.w))

# Estimaciones
sigma2.sph <- sph.variog.w$psill[2]
phi.sph <- sph.variog.w$range[2]
tau2.sph <- sph.variog.w$psill[1]
dist.vec <- emp.variog$dist # vector de distancias
sph.variog.2 <- sph.variog(sigma2.sph,phi.sph,tau2.sph,dist.vec)

# Gaussiano
gau.variog.w <- fit.variogram(emp.variog, vgm(psill= 0.02, model="Gau",
  range = 20, nugget = 0.012), fit.method = 2)
gau.variog.w

##   model  psill range
## 1   Nug 0.0069   0.0
## 2   Gau 0.0270   9.7

print(plot(emp.variog, gau.variog.w))

sigma2.gau <- gau.variog.w$psill[2]
phi.gau <- gau.variog.w$range[2]
tau2.gau <- gau.variog.w$psill[1]
gau.variog.2 <- gau.variog(sigma2.gau, phi.gau, tau2.gau, dist.vec)

# Matern
mat.variog.w <- fit.variogram(emp.variog, vgm(psill = 0.02, model = "Mat",
  range = 20, nugget = 0.012, kappa = 1), fit.method = 2)
mat.variog.w

##   model psill range kappa
## 1   Nug 0.000   0.0     0
## 2   Mat 0.035   5.6     1

print(plot(emp.variog, mat.variog.w))

sigma2.mat <- mat.variog.w$psill[2]
phi.mat <- 2*mat.variog.w$range[2]
tau2.mat <- mat.variog.w$psill[1]
nu.mat <- mat.variog.w$kappa[2]
mat.variog.2 <- mat.variog(sigma2.mat,phi.mat,tau2.mat,nu.mat,dist.vec)

# Evaluación de los modelos
### Calculamos la suma de cuadrados ponderada (WSS)
weights.sph <- emp.variog$np/((sph.variog.2)^2)
squared.diff.sph <- (emp.variog$gamma - sph.variog.2)^2
wss.sph <- sum(weights.sph * squared.diff.sph)
wss.sph

## [1] 40

weights.gau <- emp.variog$np/((gau.variog.2)^2)
squared.diff.gau <- (emp.variog$gamma-gau.variog.2)^2
wss.gau <- sum(weights.gau * squared.diff.gau)
wss.gau

## [1] 56

weights.mat <- emp.variog$np/((mat.variog.2)^2)
squared.diff.mat <- (emp.variog$gamma-mat.variog.2)^2
wss.mat <- sum(weights.mat * squared.diff.mat)
wss.mat

## [1] 53

ggplot(emp.variog, aes(x = dist, y = gamma)) + 
  geom_point() +
  labs(title = expression("Semi-variogramas"), 
       x = "distancia", y = "semivarianza") +
  stat_function(fun = gau.variog, args = list(sigma2 = sigma2.gau, phi = phi.gau, tau2 = tau2.gau), colour = "green3") +
  stat_function(fun = sph.variog, args = list(sigma2 = sigma2.sph, phi = phi.sph, tau2 = tau2.sph), colour = "orange") +
  stat_function(fun = mat.variog, args = list(sigma2 = sigma2.mat, phi = phi.mat, 
    tau2 = tau2.mat, nu = nu.mat), colour = "steelblue") +
  ylim(0, 0.038) + xlim(0, 45)

Una vez que tenemos la estimación de los parámetros del semivariograma debemos recalcular la tendencia, par esto usaremos Mínimos cuadrados generalizados (GLS): \[ Y = X\beta+ \epsilon\] donde \(\epsilon\sim N_n(0, \Sigma)\), si conocemos \(\Sigma\) el estimador de \(\beta\) es: \(\hat{\beta}_{GLS} = (X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}Y\).

De esta manera, el procedimiento de geostadística clásica es como sigue:

Obtenemos un estimador provisional de \(\beta\) usando OLS.
Calculamos los residuales usando \(\hat{\beta}_{OLS}\) y ajustamos un semivariograma paramétrico usando WLS.
Construimos la matriz de covarianzas utilizando la función de covarianzas \(C(h:\theta)\) al semivariograma paramétrico que elegimos \(\gamma(h:\theta)\), \(\hat{\Sigma}_{ij} = C(\|s_i-s_j\|:\hat{\theta})\).
Volvemos a estimar \(\beta\) usando la matriz estimada \(\hat{\Sigma}\): \(\hat{\beta}_{GLS} = (X'\hat{\Sigma}^{-1}X)^{-1}X'\hat{\Sigma}^{-1}Y\)
Se puede proceder iterativamente y recalcular los residuales: \[\hat{\eta}(s_i)=Y(s_i)-X(s_i)'\hat{\beta}_{GLS}\] para después reètir 3 y 4.

Métodos de máxima verosimilitud

Los métodos de máxima verosimilitud implican supuestos ditrsibucionales, lo más común en estadística espacial es suponer que \(Y(s)\) es un proceso Gaussiano. Escribimos nuevamente:

\[ \begin{align} \nonumber Y(s) = \mu(s) + \eta(s)\\ \nonumber \mu(s) = X(s)'\beta \end{align} \]

de manera que \(\eta(s)\) se modela como un proceso Gaussiano de media cero y covarianza \(Cov(\eta(s'),\eta(s)) = C(s',s:\theta)\). Este supuesto implica \(Y \sim MVN(X\beta, \Sigma)\), y podemos escribir la verosimilitud: \[\mathcal{L} = \frac{1}{(2\pi)^{n/2}}\frac{1}{(\| \Sigma(\theta) \|^{1/2}}exp\left(-\frac{1}{2}(Y-X\beta)'\Sigma(\theta)^{-1}(Y-X\beta)\right)\] donde \(\Sigma(\theta)_{ij}=C(s_i,s_j:\theta)\).

De tal manera que buscamos estimar \(\beta\) y \(\theta\), esto no se puede hacer de manera directa por lo que utilizamos métodos iterativos (como Newton-Raphson, o el método de puntajes de Fisher).

Algunos aspectos importantes:

Podemos comparar la verosimilitud (o el AIC, BIC) de los modelos y tener así un criterio adicional para seleccionar el modelo paramétrico.
En el proceso iterativo de estimación podemos utilizar como valores iniciales \(\theta^{(0)}\) las estimaciones de geostadística clásica (WLS).
El MLE puede no existir, la función podría ser multimodal.
Es conveniente utilizar varios valores iniciales para evaluar la convergencia.
Las estimaciones de los parámetros de covarianza \(\theta\) son sesgados, y esto podría ser problemático en muestras chicas.

Una solución al problema del sesgo de \(\hat{\theta}\) es utilizar máxima verosimilitud restringida (REML), esta consiste en realizar la máxima verosimilitud en contrastes de \((Y_1,...,Y_n)\) en lugar de en los observados \((Y_1,...,Y_n)\). En la práctica se han observado mejoras en el sesgo; sin embargo, REML puede conllevar mayor error cuadrático medio comparado a máxima verosimilitud usual.

water.geo <- as.geodata(soil, coords.col = c(1,2), data.col = 7)
water.ml <- likfit(water.geo, trend = "1st", cov.model="spherical", 
  ini=c(sigma2.sph, phi.sph), nugget = tau2.sph, fix.nugget = FALSE, 
  lik.met="ML")

## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

water.reml <- likfit(water.geo, trend = "1st", cov.model="spherical", 
  ini=c(sigma2.sph, phi.sph), nugget = tau2.sph, fix.nugget = FALSE, 
  lik.met="REML")

## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

data.frame(method = c("Classical", "ML", "REML"), 
  sigma2 = c(sigma2.sph, water.ml$sigmasq, water.reml$sigmasq), 
  phi = c(phi.sph, water.ml$phi, water.reml$phi), 
  tau2 = c(tau2.sph, water.ml$tausq, water.reml$tausq), 
  AIC = c(NA, water.ml$AIC, water.reml$AIC), 
  BIC = c(NA, water.ml$BIC, water.reml$BIC))

##      method sigma2 phi    tau2  AIC  BIC
## 1 Classical  0.033  21 0.00091   NA   NA
## 2        ML  0.032  21 0.00087 -339 -318
## 3      REML  0.034  21 0.00041 -340 -319

Compara el AIC y BIC del modelo con polinomio de segundo grado para la tendencia contra el modelo con polinomio de tercer grado.

Kriging: Predicción espacial

Dado un proceso espacial \(\{Y(s),s\in\mathbb{R}^r\}\) para el cual hemos recolectado datos en \(n\) ubicaciones \(s_1,...,s_n\) buscamos predecir una función \(g(Y(s))\) condicional a los datos observados. Iniciemos con el caso simple \(g(Y(s_0)) = Y(s_0)\), es decir, queremos predecir el valor del proceso en una nueva ubicación \(s_0\). Para hacer las estimaciones de \(Y(s_0)\) utilizaremos un método de interpolación conocido como kriging.

Las características principales de kriging son:

Las ecuaciones de kriging reflejan que las observaciones se deben ponderar distinto: las observaciones cercanas a la predicción deben contribuir más.
La dependencia espacial está incluida en las ecuaciones de kriging por medio de los variogramas (o función de covarianzas).
Los valores interpolados se modean con un proceso Gaussiano.

Vamos a plantear el problema desde un enfoque de pérdida. Sea \(\hat{Y}(s_0)\) la predicción de \(Y(s_0)\) y sea: \[L[Y(s_0), \hat{Y}(s_0)]\] la pérdida incurrida cuando predecimos \({Y}(s_0)\) usando \(\hat{Y}(s_0)\), entonces un predictor óptimo es el que minimiza la pérdida esperada también conocida como riesgo de Bayes: \[E(L[Y(s_0), \hat{Y}(s_0)]|Y)\]

donde \(Y=(Y(s_1),...,Y(s_n))'\) es el vector de observaciones de un proceso espacial \(Y(s)\).

Ahora, si escogemos como función de pérdida la pérdida cuadrática \[L[Y(s_0), \hat{Y}(s_0)] = (Y(s_0) - \hat{Y}(s_0))^2\] entonces el predictor que minimiza el riesgo de Bayes es: \[\hat{Y}(s_0) = E[Y(s_0)|Y]\] En general, no es fácil estimar \(E[Y(s_0)|Y]\), pero si \(\{Y(s),s\in\mathbb{R}^r\}\) es un proceso gaussiano \(E[Y(s_0)|Y]\) es una función lineal de \(Y\) y kriging consiste en encontrar el parámetro \(\lambda\) tal que \[\hat{Y}(s_0)=\sum_{i=1}^n \lambda_i Y(s_i)\] minimizando \(E[(Y(s_0) - \hat{Y}(s_0))^2]\).

Hay tres tipos de kriging dependiendo de los supuestos que realicemos de la tendencia \(\mu(s)\).

Kriging simple. Suponemos que la tendencia espacial \(\mu(s)\) es conocida, suponemos además que conocemos la función de covarianza (incluyendo parámetros). Bajo estos supuestos, \(\hat{Y}(s_0)\) tiene la forma \[\hat{Y}(s_0)=\sum_{i=1}^n \lambda_i Y(s_i) + k\] y \(E[(Y(s_0) - \hat{Y}(s_0))^2]\) se minimiza con: \[\lambda = (C(s_0,s_1),...,C(s_0,s_n))'\Sigma^{-1}\] \[k = \mu(s_0)-\sum_{i=1}^n\lambda_i\mu(s_i)\] por lo tanto, el predictor lineal optimo de kriging simple es: \[\hat{Y}(s_0) = \mu(s_0) + c'\Sigma^{-1}(Y-\mu)\] donde \(\mu=(\mu(s_1),...,\mu(s_n))'\) y \(c = (C(s_0,s_1),...,C(s_0,s_n))'\).

El error cuadrático medio de predicción (minimizado) es: \[\sigma_{sk}^2(s_0) = E((Y(s_0) - \hat{Y}(s_0))^2) =E((Y(s_0) - \mu(s_0) - c'\Sigma^{-1}(Y-\mu)\] \[=Var(Y(s_0)) - c'\Sigma^{-1}c\] \(\sigma_{sk}^2(s_0)\) se interpreta como la varianza de kriging y se utiliza para construir intervalos de confinaza de la predicción.

# supongamos que queremos predecir en x0, y0
x0 <- 13
y0 <- 13
# utilizaremos el objeto water.geo que creamos antes
# especificamos las opciones para el kriging, sk = simple kriging 
kc.sk.control <- krige.control(type.krige = "sk", trend.d = "cte", 
  trend.l="cte",beta = 5.8, cov.model="spherical", cov.pars = c(0.02, 20.0),
  nugget = 0.0001)
# Ponemos las coordenadas del punto en el que haremos kriging:
loc.sk <- matrix(c(x0,y0),nrow=1,ncol=2)
kc.sk.s0 <- krige.conv(water.geo,locations=loc.sk,krige=kc.sk.control)

## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

kc.sk.s0$predict

## data 
##  5.9

kc.sk.s0$krige.var

## [1] 0.0042

Kriging ordinario. Suponemos que la tendencia espacial \(\mu(s)\) es consatnte y desconocida, suponemos además que conocemos la función de covarianza. El predictor lineal optimo de kriging ordinario se obtiene utilizando multiplicadores de Langrange, buscamos: \[\hat{Y}(s_0) = \sum_{i=1}^{n} \lambda_i Y(s_i)\] sujeto a \(\sum_{i=1}^n\lambda_i = 1\), la restricción adicional viene de la media desconocida pero constante, en este obtenemos tenemos: \[\lambda'=(c + I(1-I\Sigma^{-1}c(I'\Sigma^{-1}I)^{-1})'\Sigma^{-1}\]

# Primero estimaremos los parámetros usand WLS y covarinza esférica
# Vamos a estimar una media constante (aunque desconocida) por lo que no hace
# falta estimarla por separado, es decir usar lm 
water.data <- data.frame(x = soil$x, y = soil$y, pH = soil$pH_water)
coordinates(water.data) = ~ x + y
emp.variog.water <- variogram(pH ~ 1, water.data)
plot(emp.variog.water)

# Ajustamos un variograma esférico usando WLS
sph.variog.water <- fit.variogram(emp.variog.water, vgm(psill=0.01, "Sph", 
  range=35, nugget = 0.01), fit.method = 2)
sph.variog.water

##   model   psill range
## 1   Nug 8.3e-05     0
## 2   Sph 6.1e-02    39

sigma2.wls <- sph.variog.water$psill[2]
phi.wls <- sph.variog.water$range[2]
tau2.wls <- sph.variog.water$psill[1]
# especificamos las opciones para el kriging, ok = ordinary kriging 
kc.ok.control <- krige.control(type.krige = "ok", trend.d = "cte", 
  trend.l = "cte", cov.model = "spherical", cov.pars=c(sigma2.wls, phi.wls),
  nugget = tau2.wls)
# Ponemos las coordenadas del punto en el que haremos kriging:
loc.ok <- matrix(c(x0,y0),nrow=1,ncol=2)
kc.ok.s0 <- krige.conv(water.geo, locations = loc.ok, krige = kc.ok.control)

## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

kc.ok.s0$predict

## data 
##  5.9

kc.ok.s0$krige.var

## [1] 0.0065

Kriging universal. Suponemos que la tendencia espacial es desconocida de la forma \(\mu(s)=X(s)\beta\) donde \(X(s)\) es un vector con \(p\) covariables y \(\beta\) es el vector de parámetros desconocidos. Suponemos además que conocemos la función de covarianza. Entonces, minimizando con multiplicadores de Lagrange: \[\hat{Y}(s_0) = \sum_{i=1}^{n} \lambda_i Y(s_i)\] sujeto a \(\lambda'X = X(s_0)\) y obtenemos: \[\lambda'=(c+X(X'\Sigma^{-1}X)(X(s_0)-X'\Sigma^{-1}c))'\Sigma^{-1}\]

# Primero ajustamos un modelo lineal a los datos para obtener los residuales eta.hat
mod.water.lin <- lm(pH_water ~ x + y, data = soil)
trend.water.lin <- mod.water.lin$fitted.values
res.water.lin <- soil$pH_water - trend.water.lin

# Derivamos la estimación WLS del semivariograma esférico
water.res.data <- data.frame(x = soil$x, y = soil$y, res.water.lin = res.water.lin)
coordinates(water.res.data) = ~ x + y
emp.variog.res.water <- variogram(res.water.lin ~ 1, water.res.data)
plot(emp.variog.res.water)

sph.variog.res.water <- fit.variogram(emp.variog.res.water, vgm(psill = 0.03, 
  "Sph", range = 15, nugget = 0.012), fit.method = 2)
sigma2.sph <- sph.variog.res.water$psill[2]
phi.sph <- sph.variog.res.water$range[2]
tau2.sph <- sph.variog.res.water$psill[1]

# Calculamos las estimaciones REML
water.reml.sph <- likfit(water.geo, trend = ~ soil$x + soil$y, 
  cov.model="spherical",ini=c(sigma2.sph, phi.sph), nugget=tau2.sph, 
  fix.nug = FALSE, lik.met="REML")

## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

water.reml.sph

## likfit: estimated model parameters:
##     beta0     beta1     beta2     tausq   sigmasq       phi 
## " 5.9083" "-0.0073" " 0.0000" " 0.0004" " 0.0339" "20.9611" 
## Practical Range with cor=0.05 for asymptotic range: 21
## 
## likfit: maximised log-likelihood = 176

sigma2.reml <- water.reml.sph$sigmasq
phi.reml <- water.reml.sph$phi
tau2.reml <- water.reml.sph$tausq

# Definimos los parámetros para universal kriging
# Como estamos asumiendo un polinomio lineal como tendencia cambiamos el 
# valor de trend.d trend.l
kc.uk.control <- krige.control(type.krige = "ok", trend.d = "1st", 
  trend.l = "1st", cov.model = "spherical",
cov.pars=c(sigma2.reml, phi.reml), nugget = tau2.reml)
loc.uk <- matrix(c(x0,y0), nrow=1, ncol=2)
kc.uk.s0 <- krige.conv(water.geo,locations=loc.uk,krige=kc.uk.control)

## krige.conv: model with mean given by a 1st order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhood

# estimaciones de beta:
kc.uk.s0$beta.est

##    beta0    beta1    beta2 
##  5.90831 -0.00734 -0.00002

# predicción
kc.uk.s0$predict

## data 
##  5.9

kc.uk.s0$krige.var

## [1] 0.0071

# Predecimos en un grid
x <- seq(0, 45, 1)
y <- seq(0, 120, 1)
grid.xy <- expand.grid(list(x, y))
head(grid.xy)

##   Var1 Var2
## 1    0    0
## 2    1    0
## 3    2    0
## 4    3    0
## 5    4    0
## 6    5    0

kc.uk.grid <- krige.conv(water.geo,locations=grid.xy,krige=kc.uk.control)

## krige.conv: model with mean given by a 1st order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhood

grid.xy$preds.uk <- kc.uk.grid$predict
grid.xy$preds.uk.c <- cut(grid.xy$preds.uk, class.ph$brks, include.lowest = TRUE)
ggplot(grid.xy, aes(x = Var1, y = Var2, colour = preds.uk.c)) +
  geom_point(size = 2.3) +
  scale_colour_brewer(palette = "Blues")

# varianzas
grid.xy$se.uk <- sqrt(kc.uk.grid$krige.var)
ggplot(grid.xy, aes(x = Var1, y = Var2, colour = se.uk)) +
  geom_point(size = 2.3)

Introducción a estadística espacial

Teresa Ortiz