Ricardo Alves de Olinda


http://lattes.cnpq.br/7767223263366578

Universidad del Estado de Paraíba

http://departamentos.uepb.edu.br/estatistica/corpo-docente/



Universidad de San Carlos de Guatemala

Simposio Nacional de Estadística Aplicada 2024

en las ciencias agronómicas, ambientales y forestales. En conmemoración de los 25 años de fundación del CETE



Con el apoyo de:











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



Ubicación espacial de la ciudad donde vivo.



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 distribución espacial de los incendios, basada en la función Kernel, se puede ver en la Figura anterior, donde la interpretación dentro de los límites políticos del estado describe una regionalización de la información. Se observa que la mayor concentración de focos de incendio se da en las regiones Noreste y Noroeste del estado, en el dominio del Bioma Amazónico. El año 2018 presentó un patrón diferente a los demás años observados, con una concentración de focos de incendio en la región Norte.



    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:

Durante un evento de incendio, además de energía térmica, se liberan gases y material particulado (PM), que sufren dispersión ambiental a través del aire, influenciados por las condiciones meteorológicas. La biota (incluidos los humanos) está expuesta a estos materiales en diferentes concentraciones, que varían en el tiempo y el espacio. Esta variabilidad tiene una variación de escala espacial dependiendo de las proporciones de los incendios y puede modelarse utilizando técnicas de estadística espacial para datos de área.



    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) 
# Estadística descriptiva
library(MASS)  
# Estadística descriptiva
library(fBasics)
# 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
Carbon = read.table("arbol81.txt", header=TRUE) 
# Leyendo las seis primeras líneas
head(Carbon)
##       x y Carbon Nitrogenio  Fosforo N_pastos pastos tot_pastos
## 1  0.00 0  2.609      0.208 146.6285     0.00 306.64     306.64
## 2  2.62 0  2.080      0.186 153.0903     0.00 248.16     248.16
## 3  5.25 0  1.527      0.118 146.6285     5.04 166.28     171.32
## 4  7.87 0  2.338      0.204 190.6189     0.00  55.08      55.08
## 5 10.50 0  1.727      0.146 165.5170     0.00   9.66       9.66
## 6 13.12 0  2.874      0.270 212.9869     0.00  68.40      68.40
# Mostrar de forma compacta la estructura de un objeto R arbitrario
# str(Carbon)

    Es posible importar en otras extensiones?

# Extensión en .csv  texto separado por espacios
Carbon1 = read.csv("arbol81.csv",sep=";", header=TRUE) 

# Leyendo las seis primeras líneas
head(Carbon1)
##       x y Carbono Nitrogenio   Fosforo N_gramineas gramineas tot_gramineas
## 1  0.00 0   2.609      0.208 1.466.285        0.00    306.64        306.64
## 2  2.62 0   2.080      0.186 1.530.903        0.00    248.16        248.16
## 3  5.25 0   1.527      0.118 1.466.285        5.04    166.28        171.32
## 4  7.87 0   2.338      0.204 1.906.189        0.00     55.08         55.08
## 5 10.50 0   1.727      0.146   165.517        0.00      9.66          9.66
## 6 13.12 0   2.874      0.270 2.129.869        0.00     68.40          68.4



Descripción



    Los datos para esta investigación fueron recolectados en la finca ganadera productora de carne de la Agropecuária Jaçanã, ubicada en la ciudad de Custódia-PE (8º13.70’ Latitud Sur, 37º44.70’ Longitud Oeste y 540m al nivel del mar). La vegetación de la zona inicialmente estuvo compuesta por Caatinga.



    Bioma caatinga en época seca (sin lluvias) y en época de lluvias.



    Según Tiessen (2003), en un área de la finca de aproximadamente 3.000 ha, se plantó en 1985 la especie arbórea Algaroba (Prosopis juliflora Swartz DC.), a una distancia de 10 \(m^2\), intercalada con Pasto Buffel (Cenchrus ciliares L.) y en otra zona de la finca con alrededor de 2.000 ha, ya estaban plantados Umbuzeiro y Juazeiro, de la misma manera que se hizo con Algaroba, se incluyó Capim Buffel. El pasto nunca recibió fertilización, aplicación de pesticidas agrícolas ni se registraron incendios en el sitio, la carga de animales en los pastos de la finca rondaba, al momento del estudio, alrededor del 0,17% animales por ha.

    Se establecieron 81 puntos de muestreo distribuidos regularmente, donde la especie arbórea se encuentra siempre en el centro del área. Se recolectaron muestras de suelo a una profundidad de 0-0.15m, con las coordenadas X e Y apropiadas, se secaron al aire y se pasaron a través de un tamiz de 2 mm y se analizaron para determinar el carbono orgánico total (C). De igual manera se tomaron muestras de biomasa herbácea en cada posición de la rejilla (0.7cm × 0.7cm), toda la biomasa viva en pie se cortó a nivel del suelo y se colocó en bolsas de papel, las muestras se secaron en estufa por 48 horas a 60º C. (MENEZÉS, 1999).

Estadística descriptiva espacial



    Análise descritiva del contenido de Carcon medido en muestras de suelo

geocarbon <- as.geodata(Carbon, coords.col=c(1, 2), data.col=3)
class(geocarbon)
## [1] "geodata"
attach(geocarbon)

    Argumentos utilizados por la función as.geodata() anterior:

  1. obj - Contenido de Carbon medido en muestras de suelo (base de datos);

  2. coords.col - 1:2 (columnas con la longitud y la latitud en metros);

  3. data.col - 8 (variable de interés).

    Medidas de resumen son facilmente obtenidas, utilizando la función summary().

summary(geocarbon)
## Number of data points: 72 
## 
## Coordinates summary
##      x     y
## min  0  0.00
## max 21 18.37
## 
## Distance summary
##      min      max 
##  2.62000 27.90084 
## 
## Data summary
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.746000 1.353500 1.637500 1.826544 2.115000 3.947000

    A seguir son mostrados dos gráficos que ayudan a visualizar los datos.

plot(geocarbon, low=TRUE)

scatterplot3d(Carbon[,1],Carbon[,2],Carbon[,3],pch=1,xlab="Longitud",ylab="Latitud",zlab="Carbon",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')
geocarbon$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)
#geocarbon$borders <- with(bor, cbind(x,y))
#points(geocarbon, pt.div="quintile", xlab="este", ylab="norte")
#summary(geocarbon)

    Otras formas de visualizar los datos, pueden ser obtenida de la siguinte manera.

par(mfrow=c(1,2))
points(geocarbon,xlab="Longitud",ylab="Latitud")
points(geocarbon,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(Carbon[,c(-1,-2)], ci = 0.95),2)
##             Carbon Nitrogenio Fosforo N_pastos   pastos tot_pastos
## nobs         72.00      72.00   72.00    72.00    72.00      72.00
## NAs           0.00       0.00    0.00     0.00     0.00       0.00
## Minimum       0.75       0.00   31.06     0.00     0.00       0.00
## Maximum       3.95       0.33  213.98    44.88   585.70     390.60
## 1. Quartile   1.35       0.12  118.79     0.00    65.06      67.96
## 3. Quartile   2.12       0.19  154.77     0.00   207.68     207.68
## Mean          1.83       0.16  134.50     2.42   151.41     149.60
## Median        1.64       0.16  131.96     0.00   137.94     138.64
## Sum         131.51      11.74 9683.99   174.44 10901.30   10771.43
## SE Mean       0.08       0.01    4.32     0.96    13.08      11.68
## LCL Mean      1.67       0.15  125.89     0.51   125.32     126.32
## UCL Mean      1.99       0.18  143.11     4.34   177.50     172.89
## Variance      0.46       0.00 1341.96    66.27 12325.22    9816.55
## Stdev         0.68       0.05   36.63     8.14   111.02      99.08
## Skewness      1.15       0.44   -0.32     4.07     1.08       0.49
## Kurtosis      1.11       0.95    0.70    16.33     1.70      -0.64



Otras funciones descriptivas importantes



par(mfrow=c(2,2))

histPlot(as.timeSeries(Carbon$Carbon))
#savePlot("figdbc1.jpg", type="jpg")

densityPlot(as.timeSeries(Carbon$Carbon))
#savePlot("figdbc2.jpg", type="jpg")

qqnormPlot(as.timeSeries(Carbon$Carbon))
#savePlot("figdbc3.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(geocarbon)
abline(v=1,col="red")

plot(geocarbon, low=TRUE) #Comprobación de los supuestos

    Comprobando gráficamente la dependencia espacial con 99 simulaciones

par(mfrow=c(2,2))
var1 = variog(geocarbon, max.dist=18)
## variog: computing omnidirectional variogram
plot (var1)

var2 = variog(geocarbon, option="cloud")
## variog: computing omnidirectional variogram
plot (var2)

var3 = variog(geocarbon,uvec=seq(0,18,l=20))
## variog: computing omnidirectional variogram
plot(var3)

env.var = variog.mc.env(geocarbon, obj.v=var3, nsim=99) 
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
plot(var3, env=env.var)

    Si al menos un punto está fuera de “sobre” del simulación hay evidencia de que existe una dependencia espacial.

    Es necesario estudiar la superficie de tendencia.

plot(geocarbon,lowess=TRUE)#Superficie de tendencia constante.

plot(geocarbon,lowess=TRUE,lam=0) #Superficie de tendencia constante con los datos transformados

plot(geocarbon,lowess=TRUE,trend="1st") #Superficie de tendencia de 1er orden.

plot(geocarbon,lowess=TRUE,trend="1st",lam=0) #Superficie de tendencia de 1er orden con los datos transformados.

plot(geocarbon,lowess=TRUE,trend="2nd") #Superficie de tendencia de 2er orden.

plot(geocarbon,lowess=TRUE,trend=~coords[,2]) #Superficie de tendencia de 2er orden.

    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:

  1. “A ojo” : representando diferentes modelos sobre un variograma empírico (usando la función lines.variomodel o la función eyefit);

  2. Por mínimos cuadrados : ajustando por mínimos cuadrados ordinarios (OSL) o ponderados (WLS) al variograma empírico (usando la función variofit);

  3. 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);

  4. 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(geocarbon, max.dist=13)
## variog: computing omnidirectional variogram
varct=variog(geocarbon, lam=0, max.dist=13)
## variog: computing omnidirectional variogram
var1st=variog(geocarbon, trend="1st", max.dist=13)
## variog: computing omnidirectional variogram
var1stt=variog(geocarbon, lam=0,trend="1st", max.dist=13)
## variog: computing omnidirectional variogram

Argumentos utilizados por la función variog() anterior:

  1. geodata - geocarbon (un objeto geodata);

  2. trend - 1st (tipo de tendencia);

  3. max.dist - 13 (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/transformación",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(geocarbon, max.dist=27)
## 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 variogram
vario4ct=variog4(geocarbon, lam=0, max.dist=27)
## 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 variogram
vario41st=variog4(geocarbon, trend="1st", max.dist=27)
## 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 variogram
vario41stt=variog4(geocarbon, lam=0,trend="1st", max.dist=27)
## 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 variogram
par(mfrow=c(1,2))
plot(vario4c, omni=TRUE,legend=F) #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))

par(mfrow=c(1,2))
plot(vario41st, omni=TRUE,legend=F) #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(vario41stt, 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 visual sólo para dar “disparos” iniciales para los parâmetros semivariagrama

v1 = variog(geocarbon, max.dist=27,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(geocarbon, ini=c(0.074,6.41), nug=0.02,lam=0,
        cov.model= "matern",kappa=0.5)
## ---------------------------------------------------------------
## 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:
##     beta    tausq  sigmasq      phi 
## "0.5975" "0.0000" "0.1301" "4.4033" 
## Practical Range with cor=0.05 for asymptotic range: 13.19113
## 
## likfit: maximised log-likelihood = -48.44
summary(ml1)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##   beta 
## 0.5975 
## 
## Parameters of the spatial component:
##    correlation function: exponential
##       (estimated) variance parameter sigmasq (partial sill) =  0.1301
##       (estimated) cor. fct. parameter phi (range parameter)  =  4.403
##    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 = 0 (log-transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 13.19113
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
## "-48.44"      "4"  "104.9"    "114" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
## "-64.01"      "2"    "132"  "136.6" 
## 
## Call:
## likfit(geodata = geocarbon, ini.cov.pars = c(0.074, 6.41), nugget = 0.02, 
##     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(geocarbon, ini=c(0.074,6.41), nug=0.02,lam=0,
       cov.model= "matern",kappa=1.5)
## ---------------------------------------------------------------
## 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 
## "0.5897" "0.0241" "0.1018" "2.1594" 
## Practical Range with cor=0.05 for asymptotic range: 10.24375
## 
## likfit: maximised log-likelihood = -48.29
summary(ml2)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##   beta 
## 0.5897 
## 
## Parameters of the spatial component:
##    correlation function: matern
##       (estimated) variance parameter sigmasq (partial sill) =  0.1018
##       (estimated) cor. fct. parameter phi (range parameter)  =  2.159
##       (fixed) extra parameter kappa =  1.5
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (estimated) nugget =  0.0241
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 0 (log-transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 10.24375
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
## "-48.29"      "4"  "104.6"  "113.7" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
## "-64.01"      "2"    "132"  "136.6" 
## 
## Call:
## likfit(geodata = geocarbon, ini.cov.pars = c(0.074, 6.41), nugget = 0.02, 
##     kappa = 1.5, lambda = 0, cov.model = "matern")
ml3  <- likfit(geocarbon, ini=c(0.074,6.41), nug=0.02,lam=0,
       cov.model= "gaussian")
## 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 
## "0.5810" "0.0350" "0.0857" "4.4382" 
## Practical Range with cor=0.05 for asymptotic range: 7.681706
## 
## likfit: maximised log-likelihood = -48.18
summary(ml3)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##  beta 
## 0.581 
## 
## Parameters of the spatial component:
##    correlation function: gaussian
##       (estimated) variance parameter sigmasq (partial sill) =  0.0857
##       (estimated) cor. fct. parameter phi (range parameter)  =  4.438
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (estimated) nugget =  0.035
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 0 (log-transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 7.681706
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
## "-48.18"      "4"  "104.4"  "113.5" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
## "-64.01"      "2"    "132"  "136.6" 
## 
## Call:
## likfit(geodata = geocarbon, ini.cov.pars = c(0.074, 6.41), nugget = 0.02, 
##     lambda = 0, cov.model = "gaussian")
ml4  <- likfit(geocarbon, ini=c(0.074,6.41), nug=0.02,lam=0,
        cov.model="spherical")
## 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 
## "0.5673" "0.0056" "0.1069" "7.1170" 
## Practical Range with cor=0.05 for asymptotic range: 7.116973
## 
## likfit: maximised log-likelihood = -48.29
summary(ml4)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##   beta 
## 0.5673 
## 
## Parameters of the spatial component:
##    correlation function: spherical
##       (estimated) variance parameter sigmasq (partial sill) =  0.1069
##       (estimated) cor. fct. parameter phi (range parameter)  =  7.117
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (estimated) nugget =  0.0056
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 0 (log-transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 7.116973
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
## "-48.29"      "4"  "104.6"  "113.7" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
## "-64.01"      "2"    "132"  "136.6" 
## 
## Call:
## likfit(geodata = geocarbon, ini.cov.pars = c(0.074, 6.41), nugget = 0.02, 
##     lambda = 0, cov.model = "spherical")
ml5  <- likfit(geocarbon, ini=c(0.074,6.41), nug=0.02,lam=0,
        cov.model="cauchy")
## ---------------------------------------------------------------
## 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 
## "0.6042" "0.0234" "0.1261" "2.7708" 
## Practical Range with cor=0.05 for asymptotic range: 55.34761
## 
## likfit: maximised log-likelihood = -48.6
summary(ml5)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##   beta 
## 0.6042 
## 
## Parameters of the spatial component:
##    correlation function: cauchy
##       (estimated) variance parameter sigmasq (partial sill) =  0.1261
##       (estimated) cor. fct. parameter phi (range parameter)  =  2.771
##       (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.0234
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 0 (log-transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 55.34761
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
##  "-48.6"      "4"  "105.2"  "114.3" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
## "-64.01"      "2"    "132"  "136.6" 
## 
## Call:
## likfit(geodata = geocarbon, ini.cov.pars = c(0.074, 6.41), nugget = 0.02, 
##     lambda = 0, cov.model = "cauchy")

    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. El criterio de información de Akaike (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. Dado un conjunto de modelos candidatos para los datos, el modelo preferido es el que tiene el valor mínimo en el AIC.

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 104.8704 104.5733 104.3555 104.5877 105.1975
plot(variog(geocarbon, max.dist=13, lam=0,estimator="classical"))
## variog: computing omnidirectional variogram
lines.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","gaussian","spherical","cauchy"), 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 esférica.

plot(variog(geocarbon, max.dist=13, lam=0,estimator="classical"))
## variog: computing omnidirectional variogram
lines.variomodel(ml1,col="black")
legend("bottomright", c("gaussian"), lty=1, col=c("black"))

title("Ajuste del variograma con función de correlación Gaussian")



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(geocarbon$borders, by=0.5)
points(geocarbon)
points(gr, col=2, pch=19, cex=0.3)

gr0 <- gr[.geoR_inout(gr, geocarbon$borders),]
points(gr0, col=3, pch=19, cex=0.3)

#title("Malla de puntos de predicción")

dim(gr)
## [1] 1716    2
dim(gr0)
## [1] 1634    2

    Para 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

KC <- krige.control(type="OK", obj.model=ml3)
kc1 <- krige.conv(geocarbon, 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 neighbourhood
image(kc1,col=heat.colors(51),main="Mapa de interpolación", ylim=c(-3,20),xlim=c(-1,22),
      x.leg=c(5,21.5), y.leg=c(-2.5,-1.5))

image(kc1, col=terrain.colors(51), main="Mapa de interpolación",
      ylim=c(-3,20),xlim=c(-1,22),x.leg=c(5,21.5), y.leg=c(-2.5,-1.5))

### Simulación!!

gr = pred_grid(geocarbon$borders, by=0.5)

s.out = output.control(n.predictive = 1000, n.post=1000, quant=0.95,             thres=1.6)

geocarbon.kc = krige.conv(geocarbon, loc=gr, krige=krige.control(obj=ml3),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 neighbourhood
image(geocarbon.kc, border=geocarbon$borders, loc=gr, 
      main="Mapa de los cuantiles",
      col=terrain.colors(51),val=geocarbon.kc$quan,
      ylim=c(-3,20),xlim=c(-1,22),x.leg=c(5,21.5), y.leg=c(-2.5,-1.5))

image(geocarbon.kc, col=terrain.colors(21), main="Map of P(Carbon > 1.6)= 0.95", val=(1-geocarbon.kc$probabilities.simulations),
ylim=c(-3,20),xlim=c(-1,22),x.leg=c(5,21.5), y.leg=c(-2.5,-1.5))

geocarbon.kc1 <- krige.conv(geocarbon, loc=gr, krige=krige.control(obj=ml3), out=output.control(n.pred=1000, thre=1.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 neighbourhood
dim(geocarbon.kc1$simul)
## [1] 1634 1000
hist(geocarbon.kc1$simul,probability = TRUE)

image(geocarbon.kc1, border=geocarbon$borders,loc=gr,main="Mapa del cuantil 0.10",col=gray(seq(1,0,l=51)), val=geocarbon.kc1$quan[,1], 
ylim=c(-3,20),xlim=c(-1,22),x.leg=c(5,21.5), y.leg=c(-2.5,-1.5))

image(geocarbon.kc1, border=geocarbon$borders,loc=gr,main="Mapa de la mediana",col=gray(seq(1,0,l=51)), val=geocarbon.kc1$quan[,2], 
ylim=c(-3,20),xlim=c(-1,22), x.leg=c(5,21.5), y.leg=c(-2.5,-1.5))



library(remotes)
library(rgeoboundaries)
library(leaflet)



El pensamiento estadístico un día será tan necesario para la convivencia eficiente como la capacidad de leer y escribir Traducción de la famosa frase de H.G. Wells (1886-1946).



guate_boundary <- geoboundaries("Guatemala")

guate_belize_boundaries_adm1 <- geoboundaries(c("Guatemala", "Belize"), "adm1")
guate_belize_boundaries_adm1 %>%
  leaflet() %>%
  addTiles() %>%
  addPolygons(label = guate_belize_boundaries_adm1$shapeName)



Referencia



R Core Team (2024). 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!