Estadística Espacio Temporal

Interpolación espacial por regresión Kriging




logo

Maestría en Estadística Aplicada

Roberto Trespalacios

Temas

  • Interpolación espacial por regresión Kriging
    • Kriging ordinario.
    • Kriging simple.
    • Kriging universal.
    • Kriging por bloques.
    • Kriging indicador.

Introducción

Kriging (Georges Matheron 1963) debe su nombre a Daniel G. Krige en 1945, un ingeniero de minas de Sudáfrica y primero se aplicó en datos de minería. Kriging asume un campo aleatorio expresado a través de un variograma o función de covarianza.

Georges Matheron

Daniel G. Krige

Interpolación espacial por regresión Kriging

  • Es un método muy popular para resolver el problema de predicción espacial.
  • Sea \( Z =(Z(s_1), Z(s_2),\dots, Z(s_n)) \) el vector de los datos observados en las ubicaciones espaciales \( s_1, s_2,\dots, s_n \).
  • El objetivo es estimar el valor no observado \( Z(s_0) \) (el valor de \( Z \) en la ubicación \( s_0 \)).

  • Modelo: La suposición del modelo es: \[ Z(s) = \mu + \delta(s) \] donde \( \delta(s) \) es un término aleatorio con media cero con variograma \( 2\gamma(·) \) (varianza).

  • Sistema Kriging: La suposición del predictor es \[ \hat{Z}(s_0) = \sum_{i=1}^n w_iZ(s_i) \]

es decir, es un promedio ponderado de los valores de muestra, y \( \sum_{i=1}^nw_i=1 \) para asegurar imparcialidad. Los \( w_i \), son los pesos que se estimarán.

  • Kriging minimiza el error cuadratico medio de la predicción; es decir,

\[ \min \sigma_{\varepsilon}^2 = E[Z(s_0) - \hat{Z}(s_0)]^2 \] o \[ \min \sigma_{\varepsilon}^2 = E\left[Z(s_0) - \sum_{i=1}^nw_iZ(s_i)\right]^2 \]

Error cuadratico medio para predicción Kriging

Proposición

Para procesos estacionarios de primer orden, el error cuadratico medio para predicción Kriging:

\[ \sigma_{\varepsilon}^2 = E\left[Z(s_0) - \sum_{i=1}^nw_iZ(s_i)\right]^2 \]

puede ser escrita como:

\[ \sigma_{\varepsilon}^2 = 2 \sum_{i=1}^nw_i\gamma(s_0-s_i)-\sum_{i=1}^n\sum_{j=1}^nw_iw_j\gamma(s_i-s_j) \]

Contaminación del aire en California

Utilizamos los datos de la contaminación del aire de California para ilustrar la interpolación geostatistcal (Kriging Ordinario).

Preparación de datos

Descargue “airqual.csv” para interpolar los niveles de ozono de California (promedios de 1980-2009). Use la variable OZDLYAV (la unidad es partes por billón).

Lectua de los datos

Para obtener números más fáciles de leer, multiplico OZDLYAV con 1000

x <- read_csv("~/Desktop/curso_glm_utb/geo-estadistica/pres.geo/sec_5_geo/airqual.csv")
x$OZDLYAV <- x$OZDLYAV * 1000

Creemos un SpatialPointsDataFrame y transformemoslo coordenadas Teale Albers. Tenga en cuenta las unidades = km que se necesitaba para adaptarse al variograma.

library(sp)
coordinates(x) <- ~LONGITUDE + LATITUDE
proj4string(x) <- CRS('+proj=longlat +datum=NAD83')
TA <- CRS("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=km +ellps=GRS80")
library(rgdal)
aq <- spTransform(x, TA)

Raster plantilla para interpolar.

Cree un ráster de plantilla para interpolar. Por ejemplo, dado un SpatialPolygonsDataFrame de California, lo volvemos un objeto 'SpatialGrid' (una representación diferente, pero más útil para nustro obgetivo)

cageo <- readRDS("~/Desktop/curso_glm_utb/geo-estadistica/pres.geo/sec_4_geo/counties.rds")
ca <- spTransform(cageo, TA)
r <- raster(ca)
res(r) <- 10  # 10 km si queremos resolucion en factor de 10 km
g <- as(r, 'SpatialGrid')

Ajustar un variograma

Use gstat para crear un variograma empírico.

library(gstat)
gs <- gstat(formula=OZDLYAV~1, locations=aq)
v <- variogram(gs, width=20)
head(v)

plot(v)

Ahora, ajuste un modelo de variograma teorico.

fve <- fit.variogram(v, vgm(85, "Exp", 75, 20))
fve

plot(variogramLine(fve, 400), type='l', ylim=c(0,120))
points(v[,2:3], pch=20, col='red')

Ajustemos otro tipo diferente de variograma

Observemos como se ajusta el variograma esférico en lugar de exponencial.

fvs <- fit.variogram(v, vgm(85, "Sph", 75, 20))
fvs

plot(variogramLine(fvs, 400), type='l', ylim=c(0,120) ,col='blue', lwd=2)
points(v[,2:3], pch=20, col='red')

Gráfica de Kriging ordinario (iterpolación)

Usa el variograma fve en una interpolación kriging.

k <- gstat(formula=OZDLYAV~1, locations=aq, model=fve)

kp <- predict(k, g)

spplot(kp)

Usemos la función plot para graficar.

# varianza
ok <- brick(kp)
ok <- mask(ok, ca)
names(ok) <- c('prediccción', 'varianza')
plot(ok)

Kriging universal

  • El Kriging Ordinario (OK) que se discutió anteriormente se basa en el modelo de media constante dada por \[ Z(S) = \mu + \delta(S) \] donde \( \delta(·) \) es un error y tiene una media cero y un variograma \( 2\gamma(·) \).
  • Muchas veces este es un modelo demasiado simple de usar.
  • La media puede ser una función de las coordenadas \( X, Y \), en algunos casos lineal, cuadrático o superior formar.
  • Por ejemplo, el valor de \( Z \) en la ubicación \( S \) se puede expresar ahora como lineal o cuadrático, \[ Z(S_i) = \beta_0 + \beta_1 X_i + \beta_2 Y_i + \delta (S_i) \] \[ Z(S_i) = \beta_0 + \beta_1 X_i + \beta_2 Y_i + \beta_3 X_i^2 + \beta_4 X_iY_i + \beta_5Y_i^2 + \delta(S_i) \]
  • Si este es el caso, decimos que hay una tendencia del tipo polinomial.
  • Necesitamos tomar esto en cuenta cuando encontramos los pesos kriging.
  • El valor predicho \( Z(S_0) \) en la ubicación \( S_0 \), será nuevamente una combinación lineal de los valores observados de \( Z(S_i) \), para \( i = 1, \dots, n \): \[ \hat{Z}(s_i) = w_1Z(S_1) + w_1Z(S_1) + \cdots + w_nZ(S_n)=\sum_{i=1}^nw_iZ(S_i) \],

Kriging universal

Supongamos ahora que está presente una tendencia de la forma lineal. Entonces el valor \( Z(S_0) \) se puede expresar como

\[ \hat{Z}(S_0) = \sum_{i=1}^nw_iZ(S_i)= \sum_{i=1}^nw_i\beta_0 + \sum_{i=1}^nw_i\beta_1 X_i + \sum_{i=1}^nw_i\beta_2 Y_i + \sum_{i=1}^nw_i\delta (S_i) \]

o

\[ \hat{Z}(S_0) = \sum_{i=1}^nw_iZ(S_i)= \beta_0 + \beta_1 \sum_{i=1}^nw_iX_i + \beta_2 \sum_{i=1}^nw_iY_i + \sum_{i=1}^nw_i\delta (S_i) \qquad (1) \]

Pero también se puede expresar el valor de \( Z(S_0) \) (basado en la tendencia lineal) como

\[ Z(S_0) = \beta_0 + \beta_1 X_0 + \beta_2Y_0 + \delta(S_0) \qquad (2) \]

Kriging universal

De la diapositiva anterior, compare (1) y (2).

\[ \hat{Z}(S_0) = \sum_{i=1}^nw_iZ(S_i)= \beta_0 + \beta_1 \sum_{i=1}^nw_iX_i + \beta_2 \sum_{i=1}^nw_iY_i + \sum_{i=1}^nw_i\delta (S_i) \qquad (1) \]

\[ Z(S_0) = \beta_0 + \beta_1 X_0 + \beta_2Y_0 + \delta(S_0) \qquad (2) \]

Para asegurarnos de tener un estimador insesgado, necesitaremos las siguientes condiciones:

\[ \sum_{i=1}^nX_i=X_0, \qquad \sum_{i=1}^nY_i=Y_0 \qquad y \qquad \sum_{i=1}^nw_i=1 \]

Al igual que con kriging ordinario, para encontrar los pesos cuando una tendencia está presente, necesitamos minimizar el error cuadrático medio (MSE)

\[ \min \left[Z(S_0) - \sum_{i=1}^nw_iZ(S_i)\right]^2 \]

Ejemplo universal kriging

El siguiente ejemplo usa los datos de meuse para los llacimientos de zinc

library(gstat)
# universal kriging
v <- variogram(log(zinc)~sqrt(dist), meuse)
m <- fit.variogram(v, vgm(1, "Exp", 300, 1))
plot(v, model = m)

u.kr <- krige(log(zinc)~sqrt(dist), meuse, meuse.grid, model = m)
spplot(u.kr["var1.pred"], main = "prediccion universal kriging of log-zinc")
spplot(u.kr["var1.var"], main = " varianza universal kriging")

Kriging por bloque

  • La modificación de las ecuaciones de kriging a estimar un valor promedio \( Z(A) \) de la variable \( Z \) sobre un bloque del área \( A \).

El valor promedio de \( Z(A) \) sobre el bloque A es dada por

\[ Z(A) =\int_A \frac{Z(s)}{area(A)}ds \]

estimada por:

\[ \hat{Z}(A)=\sum_{i=1}^nw_iZ(S_i) \]

donde $ \sum_{i=1}nw_i=1$

La varianza mínima, está dada por: \( \sigma^2(A) = \sum_{i=1}^nw_i\gamma(S_i,A) + \delta(B) - \gamma(A,A) \)

Ejemplo bloque kriging

library(gstat)
vg<- variogram(log(zinc)~1, meuse)
vg.fit <- fit.variogram(vg, model=vgm(psill=1, model="Sph", range=900, nugget=1))

plot(vg, vg.fit)

Krigin Ordinario

#Krigin Ordinario
ok.riging <- krige(log(zinc)~1, meuse, meuse.grid, model = vg.fit )
spplot(ok.riging["var1.pred"], main = "O.kriging predicciones")
spplot(ok.riging["var1.var"],  main = "O.kriging varianza")

#definamos el bloque para el kriging
b.k <- krige(log(zinc)~1, meuse, meuse.grid, model = vg.fit, block=c(500,500))
spplot(b.k["var1.pred"], main = "Predicciones del bloque kriging")

Indicador kriging

El indicador kriguin, es uno de los más comunes y usados.Se define de la siguiente forma:

Pas un punto \( S \), del cual no tenemos valores, se estima un valor de acuerdo a un rango definido.

\[ I(S)= \begin{cases} 1, & Si \quad Z(S_i) \leqslant valor \\ 0, & Si \quad Z(S_i) > valor \end{cases} \]

Usualmente, se definen los valores de corte, segmentando el rango de los valores \( Z(S_i) \) y de esta forma, las estimaciones son más suaves. 

#Variograma del indicadior kriging 
i.vg <- variogram(I(log(zinc)>6)~1, meuse)
i.vg.fit <- fit.variogram(i.vg, model=vgm(psill=0.25, model="Sph", range=900, nugget=0.1))
plot(i.vg, i.vg.fit)

# estimaciones indicadior kriging 
i.kriging <- krige(I(log(zinc)>6.21)~1, meuse, meuse.grid, i.vg.fit)

spplot(i.kriging["var1.pred"], main = " predicciones del indicador kriging")
spplot(i.kriging["var1.var"],  main = "varianza del indicador kriging")