Cargar los paquetes necesarios

library(raster)
library(gstat)
library(sp)
library(tidyverse)
library(sf)
library(rgdal)
library(leaflet)
library(sm)
library(variosig)
library(ape)

Cargar datos

Para poder continuar con el análisis de los datos geoestadísticos, necesitamos primero contar con dos componentes indispensables. En primer lugar hablamos de el archivo .shp que contiene toda la información recogida sobre precipitación en los puntos medidos de Ethiopia. En segundo y último lugar tenemos el contorno geográfico del área sobre la que estamos interesados y, por tanto, en el que estaremos trabajando.

Luego de esto, también debemos obtener las coordenadas del contorno de Ethiopia con el fin de, posteriormente, formar un polygono válido en la cual predecir con kriging.

eth = st_read('eth_rainst.shp')
## Reading layer `eth_rainst' from data source `D:\Estudio\Universidad\Estadistica Espacial\Trabajos\DatosReales\eth_rainst.shp' using driver `ESRI Shapefile'
## Simple feature collection with 146 features and 26 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -22353.4 ymin: 392184.7 xmax: 1083995 ymax: 1578402
## CRS:            NA
frontera = shapefile('FronteraEtiopia.shp')
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
## = dumpSRS, : Discarded ellps Clarke 1880 (RGS) in Proj4 definition: +proj=utm
## +zone=37 +a=6378249.145 +rf=293.465 +towgs84=-166,-15,204,0,0,0,0 +units=m
## +no_defs
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
## = dumpSRS, : Discarded datum Adindan in Proj4 definition: +proj=utm +zone=37
## +a=6378249.145 +rf=293.465 +towgs84=-166,-15,204,0,0,0,0 +units=m +no_defs
## Warning in showSRID(wkt2, "PROJ"): Discarded ellps Clarke 1880 (RGS) in Proj4
## definition: +proj=utm +zone=37 +a=6378249.145 +rf=293.465 +units=m +no_defs
## +type=crs
## Warning in showSRID(wkt2, "PROJ"): Discarded datum Adindan in Proj4 definition
eth.coords <- data.frame(x=frontera@polygons[[1]]@Polygons[[1]]@coords[,1],y=frontera@polygons[[1]]@Polygons[[1]]@coords[,2])

set.seed(1234)

Asimismo, ahora definimos una variable con las coordenadas de los puntos medidos. Y además, específicamos el sistema de referencia de coordenadas manualmente pues no estaba descrito en los datos al ser cargados.

coords.eth = st_coordinates(eth)

st_crs(eth) <- 32637 


ethDF = data.frame(coords.eth)
ethDF = cbind(ethDF,eth$RAIN_SUM)
names(ethDF)[3] = 'Rain'

Localizaciones de las medidas

ethWGS = st_transform(eth,4326)
ethWGS %>% leaflet()%>% addTiles()%>%
  addCircles(radius = 100,col='blue')

El gráfico anterior nos muestra los puntos geográficos en los que se tomaron medidas de precipitación; es claro que la intención es entender el fenómeno en todo el territorio ethiope, pues como vemos hay observaciones a lo largo y ancho de todo el país.

Estadísticos descriptivos

Localizaciones según el nivel de precipitación

rbPal = colorRampPalette(c('green','yellow','red'))
color = rbPal(10)[as.numeric(cut(eth$RAIN_SUM,breaks=10))]

plot(frontera,axes=T)
grid()
points(coords.eth,col=color,pch=19)

Este primer gráfico descriptivo nos muestra una vez más las posiciones de los puntos medidos en Ethiopia, sumado a esto podemos, en este caso, diferenciar además entre el nivel de la variable de interés para cada punto, siendo los puntos verdes aquellos donde la cantidad de precipitación anual fue menor y los rojos aquellos puntos en los que el total de lluvia en el año fue mayor.

De aquí podemos obtener un primer indicio de autocorrelación especial; dado que parece que los puntos geograficamente cercanos comparten valores similares de la variable de estudio. Específicamente podemos ver como las mayores cantidades de lluvia se midieron en la parte izquierda del país.

Gráficos de tendencia

lox = loess(Rain~X,data=ethDF)
loy = loess(Y~Rain,data=ethDF)

j = order(ethDF$X)
i = order(ethDF$Rain)

par(mfrow = c(1,2))
plot(Rain~X,data=ethDF)
lines(ethDF$X[j],lox$fitted[j],col = 'red')
grid()
plot(Y~Rain,data = ethDF)
lines(ethDF$Rain[i],loy$fitted[i],col='red')
grid()

La idea de estos gráficos de puntos es ser capaces de identificar tendencia en la variable de interés según su posición (coordenada X y Y), en este caso, se podría observar cierto tipo de tendencía aunque no muy marcada, una vez más nos lleva a darnos cuenta de que la parte donde se presentó mayor precipitación está localizada en la zona centro-izquierda del país. Sin embargo, no parece que hayan muchos indicios para pensar que el fenómeno no es estacionario.

Aún así, tendremos en cuenta que la variable podría no tener una media constante y luego compararemos que modelo resulta más efectivo para nuestra tarea.

Prueba de normalidad

QQ-plot

qqnorm(eth$RAIN_SUM)
qqline(eth$RAIN_SUM)

El qq-plot nos muestra como los cuantiles muestrales parecen ajustarse relativamente bien a los cuantiles teóricos de una distribución normal, esto significa un buen indicio para pensar que la precipitación en Ethiopia se distribuye normalmente. No obstante, es necesario confirmar esto con el uso de una prueba formal de normalidad.

Test de Shapiro-Wilk

shapiro.test(eth$RAIN_SUM)
## 
##  Shapiro-Wilk normality test
## 
## data:  eth$RAIN_SUM
## W = 0.98434, p-value = 0.09569

El test de shapiro-wilk para normalidad nos devuelve un valor p de 0.09569 que, si tomamos un alpha típico de 0.05, no nos lleva a rechazar que nuestra variable de interés siga una distribución normal, por lo que asumiremos que se ajusta suficientemente bien a esta.

Este resultado es muy importante para cumplir los supuestos del Kriging y que este arroje resultados significativos y de valor en sus predicciones.

Variograma cloud

coordinates(ethDF) = ~X+Y

cloudEth = variogram(Rain~1,data=ethDF,cloud = T)

plot(cloudEth)

El variograma cloud nos muestra que existen algunos valores atípicos. Aun así no parecen ser muchos y por tanto no será necesario tomar en cuenta su influencia a la hora de estimar el variograma. En caso de que hubiera muchos valores atípicos deberiamos usar un estimador robusto para tomar en cuenta el efecto que estos outliers tienen sobre el fenómeno de interés.

Variograma bin

set.seed(1)

binEth = variogram(Rain~1,data=ethDF,locations = coords.eth)

plot(binEth)

El variograma de tipo bin nos muestra como, a medida que aumenta la distancia, la semivarianza en las observaciones aumenta, hasta apróximadamente un valor de 200000 donde esta parece que se establiza sin importar que tan grande se haga la distancia entre los puntos. Esto es un buen indicio de autocorrelación espacial también pues parece que el variograma podría ajustarse bien a un modelo teórico como veremos posteriormente.

Test de autocorrelación espacial

No obstante, primero probaremos que en efecto los datos cumplen con el supuesto de autocorrelación espacial.

Test de envelopes

test.ace = envelope(binEth,data=ethDF,formula = Rain~1,cluster = T,n.cluster = 10)

envplot(test.ace)

## [1] "There are 6 out of 15 variogram estimates outside the 95% envelope."

El test de envelopes parece concordar con nuestra hipótesis de no aleatoriedad espacial, pues, aunque no por mucho, podemos observar que el variograma en ciertos puntos se encuentra fuera de las bandas de confianza aleatorias propiciadas por el test.

Sin embargo podriamos estar reacios a aceptar la autocorrelación de los datos, pues una gran cantidad de puntos se mantiene dentro de las bandas. Por ello buscaremos usar otros tests que puedan confirmar nuestra creencia inicial.

Intervalos de confianza y pruebas basadas en suavizamiento del variograma empírico (paquete sm)

sm.variogram(ethDF@coords,ethDF$Rain,model ='independent')
## Test of spatial independence: p =  0.014

El test basado en suavizamiento del variograma nos arroja un valor p pequeño, lo que nos lleva a rechazar la hipótesis nula que, en este caso, es de aleatoriedad lo que significa que aceptariamos que existe autocorrelación espacial entre los datos.

Test de Mantel

dist_Rain <- as.matrix(dist(ethDF$Rain),diag=T,upper=T)
distancia <- as.matrix(dist(ethDF@coords,diag=T,upper = T))


mantel1 <- mantel.test(distancia, dist_Rain, graph = TRUE,
                       main = "Mantel test: Rainfall",
                       xlab = "M-statistic", ylab = "Density",
                       sub = "The vertical line shows the observed M-statistic",cex.lab=0.8,
                       cex.axis=0.8,cex.main=1,cex.sub=0.6,lwd=2,xlim=c(2.05e+12,2.40e+12))

mantel1
## $z.stat
## [1] 2.372262e+12
## 
## $p
## [1] 0.001
## 
## $alternative
## [1] "two.sided"

Finalmente, como se puede observar por el estadístico en la gráfica y el reducido valor p que nos devuelve el test de mantel, rechazamos una vez más que exista aleatoriedad entre los puntos medidos.

Viendo entonces los resultados que obtuvimos a partir de los diferentes tests de autorrelación espacial, podemos conluir que, en efecto, los puntos medidos estan relacionados espacialmente entre sí.

Ajuste de un modelo teórico para el variograma

Pesos clásicos

model.wls = fit.variogram(binEth,model = vgm(c('Exp','Sph','Gau','Mat')))

model.wls
##   model    psill    range
## 1   Nug  27820.8      0.0
## 2   Exp 208465.0 263034.2
plot(binEth,model.wls)

Utilizando pesos clásicos obtenemos un variograma que sigue un modelo teórico exponencial, este parece ajustarse de manera correcta, específicamente obtenemos un valor de sill de 208465 que representa la asíntota horizontal y un rango de 263034. Sin embargo es necesario comprobar que tan bueno es el ajuste con alguna medida, en este caso la suma de residuos cuadrados.

Mínimos cuadrados ordinarios

model.ols = fit.variogram(binEth,model = vgm(c('Exp','Sph','Gau','Mat')),fit.method = 6)

model.ols
##   model     psill    range
## 1   Nug  35953.25      0.0
## 2   Sph 160305.12 436459.8
plot(binEth,model.ols)

Cuando usamos lo mínimos cuadrados ordinarios obtenemos un modelo esférico con psill igual a 160305.12 y rango 436459, el cual parece también ajustarse correctamente al variograma empírico que construimos.

Es importante notar que, a simple vista, parece que cualquiera de los modelos calculados podría ser usado para ajustar el variograma que calculamos con los datos. Sin embargo, esto es sólo un efecto visual, para selecionar el modelo correcto para nuestro caso, compararemos las sumas de residuos cuadrados que obtenemos con cada uno de los modelos propuestos.

Elección del mejor modelo segun la suma de cuadrados de los residuos

attr(model.wls,'SSErr')
## [1] 18.89199
attr(model.ols,'SSErr')
## [1] 1030879686

A pesar de que ambos modelos parecían ajustarse muy adecuadamente al variograma construido, es claro a partir de la suma de cuadrados de los residuos que el modelo que hace uso de los pesos clásicos es más correcto para ajustarse al comportamiento de la varianza de nuestro fenómeno que el modelo que hace uso de los mínimos cuadrados ordinarios, razón por la que es este primero el que seleccionaremos a partir de ahora.

Como mencionamos antes queremos tambien comparar que tan bueno sería el modelo si específicamos una tendencia en los datos.

binEthTrend = variogram(Rain~X+Y,data=ethDF,locations = coords.eth)

model.wlsTrend = fit.variogram(binEthTrend,model = vgm(c('Exp','Gau','Sph','Mat')))

attr(model.wlsTrend,'SSErr')
## [1] 21.7641

Efectivamente como se pensaba, el modelo que toma a la precipitación como una variable estacionaria, nos devuelve un menor valor de la suma de cuadrados, lo que significa que este es más efectivo para la predicción que el fenómeno que se asume con media no constante. Sin embargo lo volveremos a corroborar con la validación cruzada correspondiente.

Validación cruzada

xv.rain = krige.cv(Rain~1, ethDF,model.wls)

ECM

mean(xv.rain$residual^2)
## [1] 79217.04

Modelo lineal entre valores predichos y observados

lm.xv.rain = lm(observed~var1.pred,data=xv.rain)


plot(xv.rain$var1.pred,xv.rain$observed,main = 'Valor predicho por el Kriging Vs Valor observado',xlab = 'Valor predicho',
     ylab = 'Valor observado',pch = 19)
abline(lm.xv.rain,col='red')
text(600,1700,labels = paste0('Corr = ', round(cor(xv.rain$var1.pred,xv.rain$observed),3)))
grid()

La Gráfica anterior nos muestra una gran correlación lineal entre los valores predichos en la validación cruzada y los verdaderos valores de la variable precipitación, de hecho observamos un valor para la correlación de pearson de 0.763 lo cual demuestra una gran similitud entre los valores observados y los predichos. Esto nos da razones para creer que la utilización de Kriging para predecir valores desconocidos, aunque podría ser mejor, será suficientemente correcta.

Además como mencionamos antes queremos ver que tan bien funcionaría el modelo si en vez de media constante asumieramos una no constante.

xv.rainTrend= krige.cv(Rain~X+Y, ethDF,model.wls)

mean(xv.rainTrend$residual^2)
## [1] 80113.08
lm.xv.rainTrend = lm(observed~var1.pred,data=xv.rainTrend)


plot(xv.rainTrend$var1.pred,xv.rainTrend$observed,main = 'Valor predicho por el Kriging Vs Valor observado',xlab = 'Valor predicho',
     ylab = 'Valor observado',pch = 19)
abline(lm.xv.rainTrend,col='red')
text(600,1700,labels = paste0('Corr = ', round(cor(xv.rainTrend$var1.pred,xv.rainTrend$observed),3)))
grid()

Lo resultados son muy similares. Sin embargo, observamos que el ECM de la validación cruzada de la variable cuando se asume estacionaria es menor y su coeficiente de correlación es mayor que cuando asumimos que la variable sigue una tendencia. Esto añadido a la suma de residuos cuadrados vista anteriormente, así como por facilidad, tomaremos finalmente el modelo que utiliza la variable como un fenómeno con media constante para realizar la predicción final con Kriging.

Crear grid para el kriging

Primero debemos generar una malla de puntos para predecir que no se salga del área de interés, en este caso Ethiopía.

matrizGrid = as.matrix(eth.coords)
gridPol = st_polygon(list(matrizGrid))

grid = expand.grid(x=seq(min(eth.coords[,1]),max(eth.coords[,1]),l=100),y=seq(min(eth.coords[,2]),max(eth.coords[,2]),l=100))

grid = st_as_sf(grid,coords = c('x','y'))

datos.grid = st_intersection(grid,gridPol)

datos.gridF = data.frame(st_coordinates(datos.grid))
coordinates(datos.gridF) = ~X+Y
gridded(datos.gridF) = T

Kriging

krigeRain = krige(Rain~1,ethDF,datos.gridF,model.wls)
## [using ordinary kriging]
spplot(krigeRain['var1.pred'],contour = T,xlab='X',ylab='Y',main ='Predicción de precipitación en Ethiopia usando Kriging')

La gráfica anterior nos muestra la distribución predicha del nivel de precipitación en un año en el país africano de Etiopia. Vemos como, la mayor cantidad de lluvia anual se concentra en la parte del centro y de la izquierda del país alcanzando valores de hasta 2000 mm, a partir de ahi el nivel va reduciendose hasta que finalmente llega a un valor mínimo cercano o incluso menor a 500 mm. No obstante la predicción en las zonas más alejadas no es tan precisa pues no contamos con muchas observaciones en este área por lo que se debe ser cuidadoso al concluir acerca de esta parte del país con los datos con los que contamos.

Es importante tener claras estas zonas, pues por ejemplo en el área de la agricultura se es bien sabido que muchos de los cultivos más predominantes del país etiope como los son el café y los cereales se presentan en condiciones óptimas en lugares con precipitación anual relativamente más alta de lo normal, alrededor de los 1800 mm por año para el caso del café del país.

Varianza de la predicción

spplot(krigeRain['var1.var'],contour = T,xlab='X',ylab='Y',main =' Varianza de la predicción de precipitación \n en Ethiopia usando Kriging')

Respaldando lo que mencionamos anteriormente, la parte más alejada hacia la derecha es donde la varianza toma calores máximos, esto pues como ya dijimos, en esta zona no contamos con muchas observaciones, mientras que en la parte del centro y de la izquierda la varianza es mucho más reducida pues ahi se concentran las mayoria de las observaciones. De hecho, podemos observar los puntos exactos donde se midió la precipitación pues son las parte más oscuras de la gráfica y, por tanto, los de menor varianza.