Geoestadística aplicada a suelos, Agronomía y ciencias afines
Olinda, R. A.
Octubre, 2018
Universidad de San Carlos de Guatemala
Centro Universitario de Oriente
CUNORI
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.
Ejemplo
Precipitación pluviométrica en el estado de Paraíba, Brasil
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.7-5.2 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------#Estadística descriptiva
library(MASS)
#Estadística descriptiva
library(moments)
#Estadística descriptiva
library(fBasics)## Loading required package: timeDate##
## Attaching package: 'timeDate'## The following objects are masked from 'package:moments':
##
## kurtosis, skewness## Loading required package: timeSeries## ## Rmetrics Package fBasics## Analysing Markets and calculating Basic Statistics## Copyright (C) 2005-2014 Rmetrics Association Zurich## Educational Software for Financial Engineering and Computational Science## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.## https://www.rmetrics.org --- Mail to: info@rmetrics.org#Obtener gráficos de dispersión 3d
library(scatterplot3d) 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 .txt texto separado por tabulaciones
suelo = read.table("regiones.txt", header=TRUE)
# Leyendo las seis primeras líneas
head(suelo)## este norte altitud ca20 mg20 ctc20 ca40 mg40 ctc40 muestra region
## 1 5710 4829 6.10 52 18 106.0 40 16 86.3 a1 R3
## 2 5727 4875 6.05 57 20 131.1 49 21 123.1 a2 R3
## 3 5745 4922 6.30 72 22 114.6 63 22 101.7 a3 R3
## 4 5764 4969 6.60 74 11 114.4 74 12 95.6 a4 R3
## 5 5781 5015 6.60 68 34 124.4 44 36 106.3 a5 R3
## 6 5799 5062 5.75 45 27 132.9 27 22 105.3 a6 R3# Mostrar de forma compacta la estructura de un objeto R arbitrario
str(suelo) ## 'data.frame': 179 obs. of 11 variables:
## $ este : int 5710 5727 5745 5764 5781 5799 5817 5837 5873 5907 ...
## $ norte : int 4829 4875 4922 4969 5015 5062 5109 5161 5254 5347 ...
## $ altitud: num 6.1 6.05 6.3 6.6 6.6 5.75 5.35 5.1 5 5.2 ...
## $ ca20 : int 52 57 72 74 68 45 47 49 38 53 ...
## $ mg20 : int 18 20 22 11 34 27 27 30 25 26 ...
## $ ctc20 : num 106 131 115 114 124 ...
## $ ca40 : int 40 49 63 74 44 27 35 35 29 23 ...
## $ mg40 : int 16 21 22 12 36 22 25 25 28 18 ...
## $ ctc40 : num 86.3 123.1 101.7 95.6 106.3 ...
## $ muestra: Factor w/ 179 levels "a1","a10","a12",..: 1 7 8 9 10 11 12 13 2 3 ...
## $ region : Factor w/ 3 levels "R1","R2","R3": 3 3 3 3 3 3 3 3 2 2 ...Es posible importar en otras extensiones?
# Extensión en .csv texto separado por espacios
suelo1 = read.csv("regiones.csv",sep=";", header=TRUE)
# Leyendo las seis primeras líneas
head(suelo1)## ï..este norte altitud ca20 mg20 ctc20 ca40 mg40 ctc40 muestra region
## 1 5710 4829 6.10 52 18 106.0 40 16 86.3 a1 R3
## 2 5727 4875 6.05 57 20 131.1 49 21 123.1 a2 R3
## 3 5745 4922 6.30 72 22 114.6 63 22 101.7 a3 R3
## 4 5764 4969 6.60 74 11 114.4 74 12 95.6 a4 R3
## 5 5781 5015 6.60 68 34 124.4 44 36 106.3 a5 R3
## 6 5799 5062 5.75 45 27 132.9 27 22 105.3 a6 R3
Descripción
Este conjunto de datos contiene el contenido de calcio medido en muestras de suelo tomadas de la capa de 0-20 cm en 178 ubicaciones dentro de una determinada área de estudio dividida en tres subáreas. La elevación en cada ubicación también se registró.
La primera región suele estar inundada durante la temporada de lluvias y no se usa como área experimental. Los niveles de calcio representarían el contenido natural en la región. La segunda región ha recibido fertilizantes hace un tiempo y está típicamente ocupada por campos de arroz. La tercera región ha recibido fertilizantes recientemente y se usa con frecuencia como área experimental.
Estos datos están disponibles en el paquete geoR (Paulo J. Ribeiro Jr and Peter J. Diggle (2016)).
Guardar la hoja de datos en diferentes extensiones.
write.table(suelo, sep=" ", row.names = FALSE, file="regiones1.txt")
Estadística descriptiva espacial
Análise descritiva del contenido de calcio medido en muestras de suelo
geosuelo <- as.geodata(suelo, coords.col=c(1, 2), data.col=4)
class(geosuelo)## [1] "geodata"attach(geosuelo)Argumentos utilizados por la función as.geodata() anterior:
obj - Contenido de calcio medido en muestras de suelo (base de datos);
coords.col - 1:2 (columnas con la longitud y la latitud en UTM);
data.col - 8 (variable de interés).
Medidas de resumen son facilmente obtenidas, utilizando la función summary().
summary(geosuelo)## Number of data points: 179
##
## Coordinates summary
## este norte
## min 4957 4829
## max 5961 5720
##
## Distance summary
## min max
## 43.01163 1138.11774
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 21.00000 43.00000 50.00000 50.67039 58.00000 78.00000A seguir son mostrados dos gráficos que ayudan a visualizar los datos.
plot(geosuelo, low=TRUE)
scatterplot3d(suelo[,1],suelo[,2],suelo[,4],pch=1,xlab="Longitud",ylab="Latitud",zlab="ca20",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')
geosuelo$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)
#geosuelo$borders <- with(bor, cbind(x,y))
#points(geosuelo, pt.div="quintile", xlab="este", ylab="norte")
#summary(geosuelo)Otras formas de visualizar los datos, pueden ser obtenida de la siguinte manera.
par(mfrow=c(1,2))
points(geosuelo,xlab="Longitud",ylab="Latitud")
points(geosuelo,xlab="Longitud",ylab="Latitud",col=gray(seq(1,0.1,l=90)))
Estadística descriptiva de la variable ca20
Es posible calcular todas las estadísticas descriptivas (para todas las variables) a través del paquete fBasics, utilizando la funciónbasicStats
round(basicStats(suelo[,c(-1,-2,-10,-11)], ci = 0.95),2)## altitud ca20 mg20 ctc20 ca40 mg40 ctc40
## nobs 179.00 179.00 179.00 179.00 179.00 179.00 179.00
## NAs 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## Minimum 3.30 21.00 11.00 72.00 21.00 10.00 73.10
## Maximum 6.60 78.00 46.00 186.20 79.00 54.00 212.40
## 1. Quartile 5.20 43.00 23.00 122.30 34.00 21.00 114.50
## 3. Quartile 6.00 58.00 32.00 141.00 56.00 30.00 139.10
## Mean 5.51 50.67 27.31 132.23 45.01 25.74 127.26
## Median 5.65 50.00 27.00 132.30 44.00 26.00 128.00
## Sum 986.74 9070.00 4889.00 23669.60 8056.00 4608.00 22779.00
## SE Mean 0.05 0.83 0.47 1.37 1.01 0.52 1.67
## LCL Mean 5.41 49.04 26.39 129.53 43.01 24.72 123.97
## UCL Mean 5.61 52.30 28.24 134.93 47.00 26.76 130.55
## Variance 0.47 122.10 39.28 335.64 182.22 47.88 497.31
## Stdev 0.68 11.05 6.27 18.32 13.50 6.92 22.30
## Skewness -1.12 -0.09 0.04 -0.11 0.10 0.41 0.19
## Kurtosis 0.94 -0.37 -0.28 1.10 -0.80 0.75 0.91
Otras funciones descriptivas importantes
par(mfrow=c(2,2))
histPlot(as.timeSeries(suelo$ca20))
#savePlot("figdbc1.jpg", type="jpg")
densityPlot(as.timeSeries(suelo$ca20))
#savePlot("figdbc2.jpg", type="jpg")
qqnormPlot(as.timeSeries(suelo$ca20))
#savePlot("figdbc3.jpg", type="jpg")
newca20=data.frame(suelo$ca20,suelo$region)
plot.design(newca20, xlab="Factores", ylab="Media de ca20")
#savePlot("figdbc4.jpg",type="jpg")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(geosuelo)
abline(v=1,col="red")
plot(geosuelo, low=TRUE) #Comprobación de los supuestos
Comprobando gráficamente la dependencia espacial con 100 simulaciones
par(mfrow=c(2,2))
var1 = variog(geosuelo, max.dist=700)## variog: computing omnidirectional variogramplot (var1)
var2 = variog(geosuelo, option="cloud")## variog: computing omnidirectional variogramplot (var2)
var3 = variog(geosuelo,uvec=seq(0,700,l=50))## variog: computing omnidirectional variogramplot(var3)
env.var = variog.mc.env(geosuelo, 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 “envelope” del simulación hay evidencia de que existe una dependencia espacial.
De acuerdo con lo que discutimos por la mañana, es necesario estudiar la superficie de tendencia.
plot(geosuelo,lowess=TRUE)#Superficie de tendencia constante.
plot(geosuelo,lowess=TRUE,lam=0) #Superficie de tendencia constante con los datos transformados
plot(geosuelo,lowess=TRUE,trend="1st") #Superficie de tendencia de 1er orden.
plot(geosuelo,lowess=TRUE,trend="1st",lam=0) #Superficie de tendencia de 1er orden con los datos transformados.
plot(geosuelo,lowess=TRUE,trend="2nd") #Superficie de tendencia de 2er orden.
plot(geosuelo,lowess=TRUE,trend=~coords[,2]) #Superficie de tendencia en la coordenada y.
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(geosuelo, max.dist=700)## variog: computing omnidirectional variogramvarct=variog(geosuelo, lam=0, max.dist=700)## variog: computing omnidirectional variogramvar1st=variog(geosuelo, trend="1st", max.dist=700)## variog: computing omnidirectional variogramvar1stt=variog(geosuelo, lam=0,trend="1st", max.dist=700)## variog: computing omnidirectional variogramArgumentos utilizados por la función variog() anterior:
geodata - geosuelo (un objeto geodata);
trend - 1st (tipo de tendencia);
max.dist - 700 (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 primer orden/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(geosuelo, max.dist=700)## 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(geosuelo, lam=0, max.dist=700)## 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(geosuelo, trend="1st", max.dist=700)## 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(geosuelo, lam=0,trend="1st", max.dist=700)## 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 variogramVariograma 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(geosuelo, trend="1st",max.dist=700, 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(geosuelo, max.dist=700,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(geosuelo, ini=c(150,225), nug=25,
cov.model= "matern",kappa=0.5,trend="1st")## ---------------------------------------------------------------
## 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 beta2 tausq sigmasq phi
## "179.2046" " 0.0036" " -0.0285" " 0.0000" "101.7030" " 69.4838"
## Practical Range with cor=0.05 for asymptotic range: 208.1548
##
## likfit: maximised log-likelihood = -632.5summary(ml1)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2
## 179.2046 0.0036 -0.0285
##
## Parameters of the spatial component:
## correlation function: exponential
## (estimated) variance parameter sigmasq (partial sill) = 101.7
## (estimated) cor. fct. parameter phi (range parameter) = 69.48
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 208.1548
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-632.5" "6" "1277" "1296"
##
## non spatial model:
## log.L n.params AIC BIC
## "-654.8" "4" "1318" "1330"
##
## Call:
## likfit(geodata = geosuelo, trend = "1st", ini.cov.pars = c(150,
## 225), nugget = 25, kappa = 0.5, 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(geosuelo, ini=c(150,225), nug=25,
cov.model= "matern",kappa=1.5,trend="1st")## ---------------------------------------------------------------
## 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:
## beta0 beta1 beta2 tausq sigmasq phi
## "184.0651" " 0.0041" " -0.0298" " 20.0389" " 78.9739" " 35.1854"
## Practical Range with cor=0.05 for asymptotic range: 166.9149
##
## likfit: maximised log-likelihood = -632.9summary(ml1)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2
## 179.2046 0.0036 -0.0285
##
## Parameters of the spatial component:
## correlation function: exponential
## (estimated) variance parameter sigmasq (partial sill) = 101.7
## (estimated) cor. fct. parameter phi (range parameter) = 69.48
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 208.1548
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-632.5" "6" "1277" "1296"
##
## non spatial model:
## log.L n.params AIC BIC
## "-654.8" "4" "1318" "1330"
##
## Call:
## likfit(geodata = geosuelo, trend = "1st", ini.cov.pars = c(150,
## 225), nugget = 25, kappa = 0.5, cov.model = "matern")ml3 <- likfit(geosuelo, ini=c(150,225), nug=25,
cov.model= "gaussian",trend="1st")## 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:
## beta0 beta1 beta2 tausq sigmasq phi
## "181.4353" " 0.0043" " -0.0296" " 41.6387" " 56.1442" " 97.6096"
## Practical Range with cor=0.05 for asymptotic range: 168.9445
##
## likfit: maximised log-likelihood = -634.1summary(ml3)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2
## 181.4353 0.0043 -0.0296
##
## Parameters of the spatial component:
## correlation function: gaussian
## (estimated) variance parameter sigmasq (partial sill) = 56.14
## (estimated) cor. fct. parameter phi (range parameter) = 97.61
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 41.64
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 168.9445
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-634.1" "6" "1280" "1299"
##
## non spatial model:
## log.L n.params AIC BIC
## "-654.8" "4" "1318" "1330"
##
## Call:
## likfit(geodata = geosuelo, trend = "1st", ini.cov.pars = c(150,
## 225), nugget = 25, cov.model = "gaussian")ml4 <- likfit(geosuelo, ini=c(150,225), nug=25,
cov.model="spherical",trend="1st")## 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:
## beta0 beta1 beta2 tausq sigmasq phi
## "180.2111" " 0.0045" " -0.0296" " 26.9864" " 73.8827" "199.0425"
## Practical Range with cor=0.05 for asymptotic range: 199.0425
##
## likfit: maximised log-likelihood = -633.6summary(ml4)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2
## 180.2111 0.0045 -0.0296
##
## Parameters of the spatial component:
## correlation function: spherical
## (estimated) variance parameter sigmasq (partial sill) = 73.88
## (estimated) cor. fct. parameter phi (range parameter) = 199
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 26.99
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 199.0425
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-633.6" "6" "1279" "1298"
##
## non spatial model:
## log.L n.params AIC BIC
## "-654.8" "4" "1318" "1330"
##
## Call:
## likfit(geodata = geosuelo, trend = "1st", ini.cov.pars = c(150,
## 225), nugget = 25, cov.model = "spherical")ml5 <- likfit(geosuelo, ini=c(150,225), nug=25,
cov.model="powered.exponential",trend="1st")## ---------------------------------------------------------------
## 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:
## beta0 beta1 beta2 tausq sigmasq phi
## "165.1189" " 0.0023" " -0.0242" " 0.0000" "134.0219" "180.9213"
## Practical Range with cor=0.05 for asymptotic range: 1623.662
##
## likfit: maximised log-likelihood = -635.3summary(ml5)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2
## 165.1189 0.0023 -0.0242
##
## Parameters of the spatial component:
## correlation function: powered.exponential
## (estimated) variance parameter sigmasq (partial sill) = 134
## (estimated) cor. fct. parameter phi (range parameter) = 180.9
## (fixed) extra parameter kappa = 0.5
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 1623.662
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-635.3" "6" "1283" "1302"
##
## non spatial model:
## log.L n.params AIC BIC
## "-654.8" "4" "1318" "1330"
##
## Call:
## likfit(geodata = geosuelo, trend = "1st", ini.cov.pars = c(150,
## 225), nugget = 25, cov.model = "powered.exponential")ml6 <- likfit(geosuelo, ini=c(150,225), nug=25,
cov.model="circular",trend="1st")## 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.ml6## likfit: estimated model parameters:
## beta0 beta1 beta2 tausq sigmasq phi
## "178.2657" " 0.0047" " -0.0294" " 32.0082" " 66.3511" "178.6892"
## Practical Range with cor=0.05 for asymptotic range: 178.6892
##
## likfit: maximised log-likelihood = -634summary(ml6)## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2
## 178.2657 0.0047 -0.0294
##
## Parameters of the spatial component:
## correlation function: circular
## (estimated) variance parameter sigmasq (partial sill) = 66.35
## (estimated) cor. fct. parameter phi (range parameter) = 178.7
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 32.01
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 178.6892
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-634" "6" "1280" "1299"
##
## non spatial model:
## log.L n.params AIC BIC
## "-654.8" "4" "1318" "1330"
##
## Call:
## likfit(geodata = geosuelo, trend = "1st", ini.cov.pars = c(150,
## 225), nugget = 25, 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,ml6$AIC);AIC_modelo## ml1.AIC ml2.AIC ml3.AIC ml4.AIC ml5.AIC ml6.AIC
## 1 1276.991 1277.761 1280.218 1279.191 1282.605 1279.966plot(variog(geosuelo, max.dist=700, trend="1st",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("red","yellow","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 esférica?
plot(variog(geosuelo, max.dist=700, trend="1st",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 Matern")
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 predicción
gr <- pred_grid(geosuelo$borders, by=20)
points(geosuelo)
points(gr, col=2, pch=19, cex=0.3)
gr0 <- gr[.geoR_inout(gr, geosuelo$borders),]
points(gr0, col=3, pch=19, cex=0.3)
#title("Malla de puntos de predicción")
dim(gr)## [1] 2496 2dim(gr0)## [1] 1419 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(geosuelo, loc=gr, krige=KC)## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with mean given by a 1st order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhoodimage(kc1,main="Mapa de interpolación")
points(geosuelo$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(4950, 5450),
y.leg=c(4800, 4860), main="Mapa de interpolación")
contour(kc1,add=TRUE) #añadir gráfico de contorno
### Simulación!!
gr = pred_grid(geosuelo$borders, by=20)
s.out = output.control(n.predictive = 1000, n.post=1000, quant=0.95, thres=50)
geosuelo.kc = krige.conv(geosuelo, 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 mean given by a 1st order polynomial on the coordinates
## krige.conv: sampling from the predictive distribution (conditional simulations)
## krige.conv: Kriging performed using global neighbourhoodpar(mfrow=c(1,2))
image(geosuelo.kc, border=geosuelo$borders, loc=gr,
main="Mapa de los cuantiles",
col=terrain.colors(200),val=geosuelo.kc$quan,
x.leg=c(4950, 5450), y.leg=c(4770, 4830))
image(geosuelo.kc, col = terrain.colors(200), main="Mapa de la P(Y > 50)", val=(1-geosuelo.kc$probabilities.simulations),
x.leg=c(4950, 5500), y.leg=c(4800, 4860))
geosuelo.kc1 <- krige.conv(geosuelo, loc=gr, krige=krige.control(obj=ml1), out=output.control(n.pred=1000, thre=50, 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 mean given by a 1st order polynomial on the coordinates
## krige.conv: sampling from the predictive distribution (conditional simulations)
## krige.conv: Kriging performed using global neighbourhooddim(geosuelo.kc1$simul)## [1] 1419 1000hist(geosuelo.kc1$simul,probability = TRUE)
par(mfrow=c(1,2))
image(geosuelo.kc1, border=geosuelo$borders,loc=gr,main="Mapa del cuantil 0.10",col=gray(seq(1,0,l=21)), val=geosuelo.kc1$quan[,1],
x.leg=c(4950, 5450),
y.leg=c(4800, 4860))
image(geosuelo.kc1, border=geosuelo$borders,loc=gr,main="Mapa de la mediana",col=gray(seq(1,0,l=21)), val=geosuelo.kc1$quan[,2],
x.leg=c(4950, 5450),
y.leg=c(4800, 4860))
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!