Modelos geoestadísticos aplicados a datos de lluvia en Guatemala
Olinda, R. A.
Noviembre, 2021
Facultad de Ingeniería de la Universidad de San Carlos de Guatemala
Unidad de Modelación Matemática e Investigación
UMMI
Use R!
##R Markdown
Sobre Programación básica en R, puede ser consultada Aqui
Paquete que utilizan datos de la estadística espacial en software R Analysis of Spatial Data
Organización de datos espaciales
La arquitectura y los Sistemas de Información Geográfica; Geometría: Punto2D, Muestra, Polígono y Representación Geométrica de Rejilla Regular
La Infraestructura de Datos Espaciales de Guatemala (IDEG) desde el año 2010, puede ser consultada Aqui
##Tipo de datos Espaciales
El problema en el análisis espacial considera tres tipos de datos:
- Análisis de patrones de puntos : Este tipo de análisis permite caracterizar la estructura espacial de un conjunto de puntos en función de parâmetros como la densidad o las distancias entre puntos y su configuraci?n en el espacio.
Ejemplo:
La investigación sobre fracturas maxilofaciales resultante de violencia física y accidente de trânsito en la ciudad de Campina Grande, estado de Paraíba, Brasil.
La investigación sobre fracturas maxilofaciales resultante de violencia física y accidente de trânsito en la ciudad de Campina Grande, estado de Para?ba, Brasil.
- Autocorrelación espacial : La autocorrelación espacial es un procedimiento intrínsecamente geográfico que nos puede decir mucho acerca del comportamiento de la información georreferenciada a diferentes escalas, en particular el tipo de asociación existente entre unidades espaciales vecinas.
Ejemplo:
La base de datos posee informaciones sobre el índice de Desarrollo Humano (IDH) de 2006, 2011 y 2014, de los departamentos en Guatemala
- Datos de Superficies Continuas : Son fenômenos que se distribuyen continuamente en una región. Por ejemplo: medidas de concentración de un elemento químico en el suelo.
##Una aplicación con datos de superficies continuas (Geoestadística)
Comienza cargando las librerías y preparando los datos a utilizar.
#El análisis de los datos Geoestadísticos
library(geoR) ## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.8-1 (built on 2020-02-08) is now loaded
## --------------------------------------------------------------#Estadística descriptiva
library(MASS)
#Estadística descriptiva
library(moments)
#Estadística descriptiva
library(fBasics)## Carregando pacotes exigidos: timeDate##
## Attaching package: 'timeDate'## The following objects are masked from 'package:moments':
##
## kurtosis, skewness## Carregando pacotes exigidos: timeSeries#Obtener gráficos de dispersión 3d
library(scatterplot3d)
# leitura de arquivos em Excel
library(readxl) Para realizar los análisis exploratorios a seguir, utilizaremos el paquete geoR e inicialmente necesitaremos crear un objeto del tipo geodata, como sigue.
Importacion del datos para el R
# Extensión en .xlsx texto separado por espacios
lluvia <- "lluvia_mensual_decadas.xlsx"
Datos<- read_excel(path = lluvia, sheet = 1)
Datos## # A tibble: 192 x 5
## ID LONG LAT ALT DEC
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Labor_ovalle 391. 1645. 2385 4.3
## 2 Alameda_icta 467. 1618. 1759 2.8
## 3 Asunción_Mita 586. 1585. 469 3.2
## 4 Cahabón 574. 1726. 242 183.
## 5 Camantulul 441. 1584. 284 37.1
## 6 Camotán 621. 1639. 453 5.4
## 7 Catarina 330. 1643. 238 31
## 8 Chinique 444. 1664. 1918 12.7
## 9 Chuitinamit 436. 1691. 1207 7
## 10 Cobán 510. 1710. 1337 141.
## # ... with 182 more rowsDatos_lluvia <- as.data.frame(Datos)
# dimensión de los datos
dim(Datos_lluvia)## [1] 192 5# Leyendo las seis primeras líneas
head(Datos_lluvia)## ID LONG LAT ALT DEC
## 1 Labor_ovalle 390.9058 1644.524 2385 4.3
## 2 Alameda_icta 467.3279 1618.179 1759 2.8
## 3 Asunción_Mita 585.6604 1585.174 469 3.2
## 4 Cahabón 573.9183 1726.073 242 183.2
## 5 Camantulul 440.5905 1584.058 284 37.1
## 6 Camotán 621.3140 1639.117 453 5.4# Mostrar de forma compacta la estructura de un objeto R arbitrario
str(Datos_lluvia) ## 'data.frame': 192 obs. of 5 variables:
## $ ID : chr "Labor_ovalle" "Alameda_icta" "Asunción_Mita" "Cahabón" ...
## $ LONG: num 391 467 586 574 441 ...
## $ LAT : num 1645 1618 1585 1726 1584 ...
## $ ALT : num 2385 1759 469 242 284 ...
## $ DEC : num 4.3 2.8 3.2 183.2 37.1 ...
##Descripción
Este conjunto de datos contiene el contenido de lluvia medido en cuatro decadas….
Guardar la hoja de datos en diferentes extensiones.
write.table(Datos_lluvia, sep=";", row.names = FALSE, file="lluvia.csv")
##Estadística descriptiva espacial
Análise descritiva del datos de lluvia medidos en estaciones meteorológicas en Guatemala
Eligiendo una variable (lluviaen Deciembre) para análisis y una como covariable (altitude)
geolluvia <- as.geodata(Datos_lluvia, coords.col=c(2, 3), data.col=5,covar.col=4)
class(geolluvia)## [1] "geodata"attach(geolluvia)Argumentos utilizados por la función as.geodata() anterior:
obj - datos de lluvia medidos en estaciones meteorológicas en Guatemala (base de datos);
coords.col - 2:3 (columnas con la longitud y la latitud en UTM);
data.col - 2 (variable de interés).
Medidas de resumen son facilmente obtenidas, utilizando la función summary().
summary(geolluvia)## Number of data points: 192
##
## Coordinates summary
## LONG LAT
## min 286.1198 1526.924
## max 775.8430 2024.941
##
## Distance summary
## min max
## 0.6614913 574.0199120
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.70000 6.90000 18.30000 51.64167 70.37500 362.20000
##
## Covariates summary
## ALT
## Min. : 4.0
## 1st Qu.: 116.8
## Median : 545.5
## Mean : 804.0
## 3rd Qu.:1381.0
## Max. :3144.0A seguir son mostrados dos gráficos que ayudan a visualizar los datos.
plot(geolluvia, low=TRUE)
scatterplot3d(Datos_lluvia[,2],Datos_lluvia[,3],Datos_lluvia[,5],pch=1,xlab="Longitud",ylab="Latitud",zlab="lluvia",color="red", main="Gráfico de Dispersión")
Generamos la grilla de trabajo se puede usar otra funcion expand.grid, en help pueden ver las opciones
bord = read.table('bord.txt')
geolluvia$borders <- with(bord, cbind(V1,V2))
##Cómo crear un borde en la región de estudios?
Definición de un borde (arbitraria) en la región de estudio haciendo clic en puntos
#bor <- locator(type="p", pch=21)
#polygon(bor)
#geolluvia$borders <- with(bor, cbind(x,y))
#points(geolluvia, pt.div="quintile", xlab="este", ylab="norte")
#summary(geolluvia)Otras formas de visualizar los datos, pueden ser obtenida de la siguinte manera.
par(mfrow=c(1,2))
points(geolluvia,xlab="Longitud",ylab="Latitud")
points(geolluvia,xlab="Longitud",ylab="Latitud",col=gray(seq(1,0.1,l=90)))
Es posible calcular todas las estadísticas descriptivas (para todas las variables) a través del paquete fBasics, utilizando la funciónbasicStats
round(basicStats(Datos_lluvia[,c(-1,-2,-3)], ci = 0.95),2)## ALT DEC
## nobs 192.00 192.00
## NAs 0.00 0.00
## Minimum 4.00 0.70
## Maximum 3144.00 362.20
## 1. Quartile 116.75 6.90
## 3. Quartile 1381.00 70.38
## Mean 803.98 51.64
## Median 545.50 18.30
## Sum 154364.00 9915.20
## SE Mean 56.45 5.00
## LCL Mean 692.63 41.77
## UCL Mean 915.33 61.51
## Variance 611900.62 4808.23
## Stdev 782.24 69.34
## Skewness 0.76 2.12
## Kurtosis -0.59 4.92
##Otras funciones descriptivas importantes
par(mfrow=c(2,2))
histPlot(as.timeSeries(Datos_lluvia$DEC))
densityPlot(as.timeSeries(Datos_lluvia$DEC))
qqnormPlot(as.timeSeries(Datos_lluvia$DEC))
Con el fin de hacer constante la variabilidad en los datos, se ha evaludado la posibilidad de transformarlos usando la familia de transformaciones potenciales de Box-Cox.
boxcox(geolluvia)
abline(v=1, col="red")
Comprobando posible relación con la covariable altitude
with(Datos_lluvia, plot(DEC ~ ALT));
with(Datos_lluvia, lines(lowess(DEC ~ ALT)))
plot(geolluvia, low=TRUE) #Comprobación de los supuestos
Comprobando gráficamente la dependencia espacial con 100 simulaciones
par(mfrow=c(2,2))
var1 = variog(geolluvia, max.dist=350)## variog: computing omnidirectional variogramplot (var1)
var2 = variog(geolluvia, option="cloud")## variog: computing omnidirectional variogramplot (var2)
var3 = variog(geolluvia,uvec=seq(0,350,l=50))## variog: computing omnidirectional variogramplot(var3)
env.var = variog.mc.env(geolluvia, obj.v=var3, nsim=40) ## variog.env: generating 40 simulations by permutating data values
## variog.env: computing the empirical variogram for the 40 simulations
## variog.env: computing the envelopsplot(var3, env=env.var)
Si al menos un punto está fuera de “sobre” del simulación hay evidencia de que existe una dependencia espacial.
De acuerdo con lo que discutimos es necesario estudiar la superficie de tendencia.
plot(geolluvia,lowess=TRUE)#Superficie de tendencia constante.
plot(geolluvia,lowess=TRUE,lam=0) #Superficie de tendencia constante con los datos transformados (log10)
plot(geolluvia,lowess=TRUE,trend="1st") #Superficie de tendencia de 1er orden.
plot(geolluvia,lowess=TRUE,trend="1st",lam=0) #Superficie de tendencia de 1er orden con los datos transformados.
plot(geolluvia,lowess=TRUE,trend="2nd") #Superficie de tendencia de 2er orden.
plot(geolluvia,lowess=TRUE,trend=~ALT) #Superficie de tendencia con altitude.
De acuerdo con los gráficos anteriores, con cual superficie de tendencia trabajaríamos?
##Ajuste de un modelo de variograma
Los estimadores empíricos no pueden ser empleados en la práctica (no verifican necesariamente las propiedades de un variograma válido), por lo que se suele recurrir en la práctica al ajuste de un modelo válido. Con el paquete geoR podemos realizar el ajuste:
“A ojo” : representando diferentes modelos sobre un variograma empírico (usando la función lines.variomodel o la función eyefit);
Por mínimos cuadrados : ajustando por mínimos cuadrados ordinarios (OSL) o ponderados (WLS) al variograma empírico (usando la función variofit);
Por máxima verosimilitud : estimando por máxima verosimilitud (ML) o máxima verosimilitud restringida (REML) los parâmetros a partir de los datos (utilizando la función likfit);
Métodos bayesianos (utilizando la función krige.bayes).
Además de la superficie de tendencia, debemos entender la estructura de covarianza presente en los datos (esto si ella existe). A seguir son obtenidos los variogramas por medio de la función variog().
varc=variog(geolluvia, max.dist=350)## variog: computing omnidirectional variogramvarct=variog(geolluvia, lam=0, max.dist=350)## variog: computing omnidirectional variogramvar1st=variog(geolluvia, trend="1st", max.dist=350)## variog: computing omnidirectional variogramvar1stt=variog(geolluvia, lam=0,trend=~ALT, max.dist=350)## variog: computing omnidirectional variogramArgumentos utilizados por la función variog() anterior:
geodata - geolluvia (un objeto geolluvia);
trend - 1st (tipo de tendencia);
max.dist - 350 (distancia máxima).
Abajo se encuentran los comandos para plotar los variogramas.
par(mfrow=c(2,2))
plot(varc,main="Tendencia constante",xlab="Distancia en metros",ylab="Semivarianza")
plot(varct,main="Tendencia constante/transformación",xlab="Distancia en metros",ylab="Semivarianza")
plot(var1st,main="Tendencia de primer orden",xlab="Distancia en metros",ylab="Semivarianza")
plot(var1stt,main="Tendencia de altitude/transformacion",xlab="Distancia en metros",ylab="Semivarianza")
Para concluir este tópico mostramos una función auxiliar que facilita la construcción de variogramas empíricos en direcciones diferentes que pueden ser utilizados para investigar anisotropía.
Para estudiar si hay anisotropía, se pueden cálcular de forma rápida variogramas direccionales con la función variog4. Por defecto calcula cuatro variogramas direccionales, correspondientes a los ângulos 0, 45, 90 y 135 grados:
vario4c=variog4(geolluvia, max.dist=350)## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogramvario4ct=variog4(geolluvia, lam=0, max.dist=350)## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogramvario41st=variog4(geolluvia, trend="1st", max.dist=350)## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogramvario41stt=variog4(geolluvia, lam=0,trend=~ALT, max.dist=350)## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogrampar(mfrow=c(1,2))
plot(vario4c, omni=TRUE,legend=FALSE) #Incluyendo el variograma omnidireccional
legend("bottomright", c("omnid.", "0°","45°","90°","135°"), lty=2, col=c("black","red","green","blue"),xpd=TRUE,cex=0.8)
plot(vario4ct, omni=TRUE,legend=FALSE)
legend("bottomright", c("omnid.", "0°","45°","90°","135°"), lty=2, col=c("black","red","green","blue"),xpd=TRUE,cex=0.8)
# Gráfico con la firma, la fecha y la versión
mtext(side=3, at=par("usr")[3], adj=0.7,
cex=0.8, col="gray40", line=1,
text=paste("Ricardo Olinda --",
format(Sys.time(), "%d/%m/%Y %H:%M:%S --"),
R.version.string))
Variograma con diferentes gráficos en comparación con el efecto omnidireccional con Tendencia constante
plot(vario4c, omni=TRUE, same=FALSE)
Variograma con diferentes gráficos en comparación con el efecto omnidireccional con Tendencia de primer orden
plot(vario41st, omni=TRUE, same=FALSE)
Tendrás que estudiar adecuadamente la cuestión de la anisotropía en su datos.
Esto no le exime de la anisotropóa efecto combinado + tendencia.
Sólo para probar las posibilidades, vamos a hacer el variograma hacia 45 grados.
va.45 <- variog(geolluvia, lam=0,max.dist=300, dir=pi/4,estimator="classical")## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)plot(va.45)
Variograma visual sólo para dar “disparos” iniciales para los parâmetros semivariagrama
v1 = variog(geolluvia, max.dist=350,dir=pi/4,estimator="classical") ## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)plot(v1)
## Ajuste (visual) del modelo para el variograma
#ef1 = eyefit(v1)
#summary(ef1)Estimación de los parâmetros (máxima verosimilitud) los valores iniciales “ef1” sugerido visualmente al ver la función variograma del modelo ajustado por máxima verosimilitud (ML).
Las estimativas de máxima verosimilitud son obtenidas ajustando el modelo a los datos. Los modelos pueden ser comparados por medio de medidas de comparación de modelos, como, por ejemplo, el AIC y BIC.
ml1 <- likfit(geolluvia, ini=c(4,225), nug=1,
cov.model= "matern",kappa=0.5,lam=0,trend=~ALT)## ---------------------------------------------------------------
## 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.ml1## likfit: estimated model parameters:
## beta0 beta1 tausq sigmasq phi
## " 3.2491" " 0.0002" " 0.1781" " 3.8170" "224.9779"
## Practical Range with cor=0.05 for asymptotic range: 673.9735
##
## likfit: maximised log-likelihood = -825.4summary(ml1)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1
## 3.2491 0.0002
##
## Parameters of the spatial component:
## correlation function: exponential
## (estimated) variance parameter sigmasq (partial sill) = 3.817
## (estimated) cor. fct. parameter phi (range parameter) = 225
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.1781
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 0 (log-transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 673.9735
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-825.4" "5" "1661" "1677"
##
## non spatial model:
## log.L n.params AIC BIC
## "-922.7" "3" "1851" "1861"
##
## Call:
## likfit(geodata = geolluvia, trend = ~ALT, ini.cov.pars = c(4,
## 225), nugget = 1, kappa = 0.5, lambda = 0, cov.model = "matern")# información disponible en la función likfit
names(ml1)## [1] "cov.model" "nugget"
## [3] "cov.pars" "sigmasq"
## [5] "phi" "kappa"
## [7] "beta" "beta.var"
## [9] "lambda" "aniso.pars"
## [11] "tausq" "practicalRange"
## [13] "method.lik" "trend"
## [15] "loglik" "npars"
## [17] "AIC" "BIC"
## [19] "parameters.summary" "info.minimisation.function"
## [21] "max.dist" "trend"
## [23] "trend.matrix" "transform.info"
## [25] "nospatial" "model.components"
## [27] "call"ml2 <- likfit(geolluvia, ini=c(4,225), nug=1,
cov.model= "matern",kappa=1.5,lam=0)## ---------------------------------------------------------------
## 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.ml2## likfit: estimated model parameters:
## beta tausq sigmasq phi
## " 1.9192" " 0.5698" " 30.2784" "224.9458"
## Practical Range with cor=0.05 for asymptotic range: 1067.112
##
## likfit: maximised log-likelihood = -846.5summary(ml1)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1
## 3.2491 0.0002
##
## Parameters of the spatial component:
## correlation function: exponential
## (estimated) variance parameter sigmasq (partial sill) = 3.817
## (estimated) cor. fct. parameter phi (range parameter) = 225
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.1781
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 0 (log-transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 673.9735
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-825.4" "5" "1661" "1677"
##
## non spatial model:
## log.L n.params AIC BIC
## "-922.7" "3" "1851" "1861"
##
## Call:
## likfit(geodata = geolluvia, trend = ~ALT, ini.cov.pars = c(4,
## 225), nugget = 1, kappa = 0.5, lambda = 0, cov.model = "matern")ml3 <- likfit(geolluvia, ini=c(4,225), nug=1,
cov.model= "gaussian",lam=0)## kappa not used for the gaussian 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.ml3## likfit: estimated model parameters:
## beta tausq sigmasq phi
## " 3.5381" " 0.5733" " 1.5716" "84.6174"
## Practical Range with cor=0.05 for asymptotic range: 146.4574
##
## likfit: maximised log-likelihood = -846.5summary(ml3)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta
## 3.5381
##
## Parameters of the spatial component:
## correlation function: gaussian
## (estimated) variance parameter sigmasq (partial sill) = 1.572
## (estimated) cor. fct. parameter phi (range parameter) = 84.62
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.5733
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 0 (log-transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 146.4574
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-846.5" "4" "1701" "1714"
##
## non spatial model:
## log.L n.params AIC BIC
## "-927" "2" "1858" "1865"
##
## Call:
## likfit(geodata = geolluvia, ini.cov.pars = c(4, 225), nugget = 1,
## lambda = 0, cov.model = "gaussian")ml4 <- likfit(geolluvia,ini=c(4,225), nug=1,
cov.model="spherical",lam=0)## 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.ml4## likfit: estimated model parameters:
## beta tausq sigmasq phi
## " 3.5553" " 0.1888" " 2.6952" "243.3253"
## Practical Range with cor=0.05 for asymptotic range: 243.3253
##
## likfit: maximised log-likelihood = -827.4summary(ml4)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta
## 3.5553
##
## Parameters of the spatial component:
## correlation function: spherical
## (estimated) variance parameter sigmasq (partial sill) = 2.695
## (estimated) cor. fct. parameter phi (range parameter) = 243.3
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.1888
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 0 (log-transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 243.3253
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-827.4" "4" "1663" "1676"
##
## non spatial model:
## log.L n.params AIC BIC
## "-927" "2" "1858" "1865"
##
## Call:
## likfit(geodata = geolluvia, ini.cov.pars = c(4, 225), nugget = 1,
## lambda = 0, cov.model = "spherical")ml5 <- likfit(geolluvia, ini=c(4,225), nug=1,
cov.model="circular",lam=0)## kappa not used for the circular 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.ml5## likfit: estimated model parameters:
## beta tausq sigmasq phi
## " 3.5526" " 0.1866" " 3.0847" "233.9366"
## Practical Range with cor=0.05 for asymptotic range: 233.9366
##
## likfit: maximised log-likelihood = -826.5summary(ml5)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta
## 3.5526
##
## Parameters of the spatial component:
## correlation function: circular
## (estimated) variance parameter sigmasq (partial sill) = 3.085
## (estimated) cor. fct. parameter phi (range parameter) = 233.9
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.1866
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 0 (log-transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 233.9366
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-826.5" "4" "1661" "1674"
##
## non spatial model:
## log.L n.params AIC BIC
## "-927" "2" "1858" "1865"
##
## Call:
## likfit(geodata = geolluvia, ini.cov.pars = c(4, 225), nugget = 1,
## lambda = 0, cov.model = "circular")Comparando los ajustes por el criterio de información de Akaike (AIC). El criterio de AIC es una medida de la calidad relativa de un modelo estadístico, para un conjunto dado de datos. Como tal, el AIC proporciona un medio para la selección del modelo.
AIC_modelo=data.frame(ml1$AIC,ml2$AIC,ml3$AIC,ml4$AIC,
ml5$AIC);AIC_modelo## ml1.AIC ml2.AIC ml3.AIC ml4.AIC ml5.AIC
## 1 1660.828 1701.046 1700.982 1662.71 1661.051plot(variog(geolluvia, max.dist=350, lam=0,trend=~ALT,estimator="classical"))## variog: computing omnidirectional variogramlines.variomodel(ml1,col="yellow")
lines.variomodel(ml2,col="red")
lines.variomodel(ml3,col="green")
lines.variomodel(ml4,col="black")
lines.variomodel(ml5,col="blue")
legend("bottomright", c("matérn k=0.5", "matérn k=1.5","gaussiano","esférico",
"circular"), lty=1, col=c("yellow","red","green","black","blue"))
#http://listas.inf.ufpr.br/pipermail/r-br/2014-October/014618.html
title("Ajuste de variograma con diferentes funciones de correlación")
Los modelos suponiendo una función de correlación espacial mattern tendrian un mejor ajuste comparado com los ajustes suponiendo una función de correlación matérn
plot(variog(geolluvia, max.dist=350, lam=0,estimator="classical"))## variog: computing omnidirectional variogramlines.variomodel(ml1,col="black")
legend("bottomright", c("matérn k=0.5"), lty=1, col=c("black"))
title("Ajuste del variograma con función de correlación Matérn")
##Predicción espacial (kriging)
El paquete geoR dispone de opciones para los métodos kriging tradicionales, que dependiendo de las suposiciones acerca de la función de tendencia se clasifican en:
Kriging simple (KS) : media conocida;
Kriging ordinario (KO) : se supone que la media es constante y desconocida;
Kriging universal (KU) : también denominado kriging con modelo de tendencia, se supone que la media es una combinación lineal (desconocida) de las coordenadas o de otras variables explicativas.
Existen también opciones adicionales para kriging trans-normal (con transformaciones Box-Cox para aproximarse a la normalidad y transformación de nuevo de resultados a la escala original manteniendo insesgadez). También admite modelos de variograma geométricamente anisotrópicos.
Malla de puntos de prediccón
gr <- pred_grid(geolluvia$borders, by=10)
points(geolluvia)
points(gr, col=2, pch=19, cex=0.3)
gr0 <- gr[.geoR_inout(gr, geolluvia$borders),]
points(gr0, col=3, pch=19, cex=0.3)
#title("Malla de puntos de predicción")
dim(gr)## [1] 3127 2dim(gr0)## [1] 2192 2Para ver todas las opciones de kriging disponibles ejecutar ?krige.control. Para kriging con vecindario local (archivos de datos grandes) se puede utilizar la función ksline.
Para representar las superficies se puede utilizar la función image:
Viendo el mapa del kriging
par(mfrow=c(1,2))
KC <- krige.control(type="OK", obj.model=ml1)
kc1 <- krige.conv(geolluvia, loc=gr, krige=KC)## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: performing the Box-Cox data transformation
## krige.conv: back-transforming the predicted mean and variance
## krige.conv: Kriging performed using global neighbourhoodimage(kc1,main="Mapa de interpolación")
points(geolluvia$coords, pch=20) #añadir posiciones datos
contour(kc1,add=TRUE) #añadir gráfico de contorno
image(kc1, col=terrain.colors(200),x.leg=c(270, 600),
y.leg=c(1430, 1450), main="Mapa de interpolación")
contour(kc1,add=TRUE) #añadir gráfico de contorno
### Simulación!!
gr = pred_grid(geolluvia$borders, by=10)
s.out = output.control(n.predictive = 1000, n.post=1000, quant=0.95, thres=51.6)
geolluvia.kc = krige.conv(geolluvia, loc=gr, krige=krige.control(obj=ml1), output = s.out)## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: performing the Box-Cox data transformation
## krige.conv: sampling from the predictive distribution (conditional simulations)
## krige.conv: back-transforming the simulated values
## krige.conv: back-transforming the predicted mean and variance
## krige.conv: Kriging performed using global neighbourhoodpar(mfrow=c(1,2))
image(geolluvia.kc, border=geolluvia$borders, loc=gr,
main="Mapa de los cuantiles",
col=terrain.colors(200),val=geolluvia.kc$quan,
x.leg=c(270, 750),
y.leg=c(1420, 1440))
image(geolluvia.kc, col = terrain.colors(200), main="Mapa de la P(Y > 51.6)", val=(1-geolluvia.kc$probabilities.simulations),
x.leg=c(270, 600),
y.leg=c(1430, 1450))
geolluvia.kc1 <- krige.conv(geolluvia, loc=gr, krige=krige.control(obj=ml1), out=output.control(n.pred=1000, threshold =51.6, quant=c(0.10, 0.5, 0.9)))## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: performing the Box-Cox data transformation
## krige.conv: sampling from the predictive distribution (conditional simulations)
## krige.conv: back-transforming the simulated values
## krige.conv: back-transforming the predicted mean and variance
## krige.conv: Kriging performed using global neighbourhooddim(geolluvia.kc1$simul)## [1] 2192 1000image(geolluvia.kc1,border=geolluvia$borders,loc=gr,main="Mapa del cuantil 0.10",col=gray(seq(1,0,l=21)), val=geolluvia.kc1$quan[,1],
x.leg=c(270, 600),
y.leg=c(1430, 1450))
##Referencia
R Core Team (2018). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
Paulo J. Ribeiro Jr and Peter J. Diggle (2016). geoR: Analysis of Geostatistical Data. R package version 1.7-5.2. https://CRAN.R-project.org/package=geoR
Gracias!