1.1 Pacotes Utilizados

library(geoR)
library(MASS)
library(moments)
library(fBasics)
library(akima)
library(readxl)
library(skimr)
library(terra)
library(readxl)
library(knitr)
library(classInt)

1.2 Importação e vizualização dos dados

# Importação
dados_ph = read_xlsx("dados_ph.xlsx")
dados_geral = read_xlsx("Dados_Apresenta.xlsx")
borda = read.table('bord.txt')

# Ajustes
df_base = cbind(dados_geral[,1:5],dados_ph)
#polygon(borda)

geosolo <- as.geodata(df_base, coords.col=c(2,3), data.col=6)
geosolo$borders <- with(borda, cbind(V1,V2))
attach(geosolo)

# Vizualização
kable(head(df_base,10), align = "l",
      caption = "10 primeiras observações")

2.1 Resumo

summary(geosolo)

## Number of data points: 243 
## 
## Coordinates summary
##        X...2   Y...3
## min 420797.9 8806076
## max 985985.9 9686916
## 
## Distance summary
##         min         max 
##    120.1225 968602.3272 
## 
## Borders summary
##            V1      V2
## min  390646.6 8791849
## max 1012045.0 9696707
## 
## Data summary
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 3.110000 3.792500 4.550000 4.715905 5.412500 8.565000

     Verificamos assim, que a média dos dados de pH foi igual 4.715905, esta medida descritiva iremos utilizar no processo de krigagem mais adiante. Ainda podemos verificar que em relação as coordenadas temos do mínimo para o máximo:

• X: [420797.9 : 985985.9]
• Y: [8806076 : 9686916]

     A distância mínima encontrada foi de 120.1225 e a máxima de 968602.3272, como iremos utilizar 1/3 da distância para a construção do variograma, iremos ter utilizar a distância no valor 322867.4.

2.2 Descritiva Espacial

plot(geosolo, low=T,lam=0)

     Podemos obervar que nos gráficos 2 e 3, não existe uma evidência clara de efeitos de tendência sobre as coordenadas Y e as coordenadas X, respectivamente. No primeiro gráfico é apresentado a relação da variável pH com base nas coordenadas X e Y, enquanto que pelo quarto gráfico vemos uma distribição quase uniforme nos dados da variável entre 1.2 e 1.7.

     Podemos ainda verificar por uma mapa onde podemos verificar as maiores e menores níveis de pH no estado do Piauí, com base em um intervalo mais aproximado da escala de pH:

points(geosolo, pt.div=c(3,4,6,8,9), col = c("red","yellow","green3","blue3"),
       xlab="Coord X", ylab="Coord Y",
       x.leg=250000,
       y.leg=9650000)

     Verificamos que pelo mapa de pontos, na parte sudeste do estado há uma maior concentração de níveis de pH considerados neutro com uma escala entre 6 e 8, assim também como no norte do estado que possuem alguns pontos nesta escala. Ainda região sudeste, verificamos que é a única região do estado a possuir pontos evidenciando níveis mais altos de pH, ou seja, águas mais alcalinas.

     Na região norte e centro do estado vemos uma maior concentração de níveis que variam entre 4 a 6, ou seja, razoavelmente ácidas, possivelmente uma maior concetração de águas advindas pela chuva. Já na região sul do estado vemos uma maior distribuição de pontos com níveis mais báixos, ou seja, uma região que possui águas moderadamente alcalinas.

     Outro mapa que podemos mostrar é o com base nos quintis, portanto:

points(geosolo, pt.div="quintiles",
       xlab="Coord X", ylab="Coord Y",
       x.leg=250000,
       y.leg=9650000)

plot(geosolo, low = T, trend = "1st")

plot(geosolo, low = T, trend = "2nd")

2.3 Variogramas

par(mfrow=c(2,2))
var1 = variog(geosolo, coords = coords , max.dist=322867.4)

## variog: computing omnidirectional variogram

var2 = variog(geosolo, coords = coords , option="cloud")

## variog: computing omnidirectional variogram

var3 = variog(geosolo, coords = coords ,uvec=seq(0,322867.4,l=50))

## variog: computing omnidirectional variogram

env.var = variog.mc.env(geosolo, 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 envelops

plot(var1)
plot(var2)
plot(var3)
plot(var3, env=env.var)

     Podemos notar que em determinados intervalos de distância, há evidencias para uma dependencia espacial, como nos intervalos de 0 a mais ou menos 70000, de 80000 a uns 250000 e assim por diante.Ainda assim, pelo último plot da simulação de envelope, podemos considerar que existe evidências de dependência espacial.

3.1 Estimação dos parâmetros

     Utilizando os cov.models:

• spherical
• matern
• wave
• exponential
• gaussian

onde no último modelo wave, será realizado alguns ajustes extras afim de verificar um melhor ajuste. Portanto, teremos os seguintes ajustes:

# Estimação dos Parâmetros
pH.ml0 <- list()
pH.ml0$matern <- likfit(geosolo,  ini=c(1.31, 83679.35), nug=0.24,cov.model = "matern")

## ---------------------------------------------------------------
## 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.

pH.ml0$exp <- likfit(geosolo,  ini=c(1.31, 83679.35), nug=0.24,cov.model = "exponential")

## kappa not used for the exponential 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.

pH.ml0$gaus <- likfit(geosolo,  ini=c(1.31, 83679.35), nug=0.24,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.

pH.ml0$spher <- likfit(geosolo,  ini=c(1.31, 83679.35), nug=0.24,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.

pH.ml0$wave <- likfit(geosolo,  ini=c(1.31, 83679.35), nug=0.24,cov.model = "wave")

## kappa not used for the wave 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.

summary(pH.ml0$matern)

## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##   beta 
## 4.7988 
## 
## Parameters of the spatial component:
##    correlation function: exponential
##       (estimated) variance parameter sigmasq (partial sill) =  0.5305
##       (estimated) cor. fct. parameter phi (range parameter)  =  83679
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (estimated) nugget =  0.251
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 1 (no transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 250680.9
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
## "-254.9"      "4"  "517.8"  "531.8" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
## "-363.7"      "2"  "731.4"  "738.4" 
## 
## Call:
## likfit(geodata = geosolo, ini.cov.pars = c(1.31, 83679.35), nugget = 0.24, 
##     cov.model = "matern")

pH.ml0$matern$nospatial$AIC.ns

## [1] 731.412

AIC = (sapply(pH.ml0, AIC))
BIC = c(pH.ml0$matern$BIC,pH.ml0$exp$BIC,pH.ml0$gaus$BIC,pH.ml0$spher$BIC,pH.ml0$wave$BIC)

compar = rbind(AIC,BIC)

kable(compar, caption = "Comparação das métricas entre os modelos espaciais")

     Podemos verificar que dois modelos registraram os melhores AIC’s e BIC’s entre os demais e sendo iguais entre si, portanto, por critério de escolha iremos continuar as análises utilizando o modelo com base na função de correlação matern.

AIC_NS = c(pH.ml0$matern$nospatial$AIC.ns, pH.ml0$exp$nospatial$AIC.ns,
        pH.ml0$gaus$nospatial$AIC.ns,
        pH.ml0$spher$nospatial$AIC.ns,
        pH.ml0$wave$nospatial$AIC.ns)

BIC_NS = c(pH.ml0$matern$nospatial$BIC.ns, pH.ml0$exp$nospatial$BIC.ns,
        pH.ml0$gaus$nospatial$BIC.ns,
        pH.ml0$spher$nospatial$BIC.ns,
        pH.ml0$wave$nospatial$BIC.ns)

compar_ns = rbind(AIC_NS,BIC_NS)

colnames(compar_ns) = c("matern","exp","gaus","spher","wave")

kable(compar_ns, caption = "Comparação das métricas entre os modelos não espaciais")

     Com base na tabela acima, podemos verificar que as métricas AIC e BIC para os modelos não espcaciais são idênticas.

4.1 Grid para predição espacial

GR <- pred_grid(geosolo$borders, by=10000)
gr0 <- GR[.geoR_inout(GR, geosolo$borders),]
points(geosolo)
points(GR, pch=".",cex=1.5)

4.2 Predição espacial

KC <- krige.control(type="SK",obj.model = pH.ml0$matern)
geosolo.kc <- krige.conv(geosolo, loc=gr0, krige=KC)

## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

image(geosolo.kc, GR,
      geosolo$borders,
      col=terrain.colors(200),
      x.leg=c(858233, 1200000), y.leg=c(8859000, 8910972),cex=0.7)

     Portanto, conforme o que já vinha sido apresentado anteriormente, vemos que na região sudeste existe a maior intensidade de níveis mais elevados de pH com valores de 6 a mais.

4.3 Predição espacial com base na raiz de variância

image(geosolo.kc, loc=GR,
      bor=geosolo$borders, col=terrain.colors(21),
      val=sqrt(geosolo.kc$krige.var),
      x.leg=c(858233, 1250000), y.leg=c(8859000, 8910972),cex=0.6)

     Com base neste mapa verificamos uma maior concentração na parte sul, nordeste do estado.

4.4 Simulação

     Utilizando a média = 4.715905:

OC <- output.control(thres=4.715905,
                     quantile=c(0.1, 0.9))

geosolo.kc <- krige.conv(geosolo, loc=gr0, krige=KC, out=OC)

## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: sampling from the predictive distribution (conditional simulations)
## krige.conv: Kriging performed using global neighbourhood

image(geosolo.kc, loc=GR, bor=geosolo$borders, col=terrain.colors(21), 
      val=1-geosolo.kc$prob,x.leg=c(858233, 1200000), y.leg=c(8859000, 8910972),cex=0.7)

image(geosolo.kc, loc=GR, bor=geosolo$borders, col=terrain.colors(21), 
      val=geosolo.kc$quantile$q90,x.leg=c(858233, 1200000), y.leg=c(8859000, 8910972),cex=0.7)
title("90º percentil")

image(geosolo.kc, loc=GR, bor=geosolo$borders, col=terrain.colors(21), 
      val=geosolo.kc$quantile$q10,x.leg=c(858233, 1200000), y.leg=c(8859000, 8910972),cex=0.7)
title("10º percentil")