INTRODUÇÃO
O banco de dados para análise geoestatística com observações da precipitação no estado da Paraíba.
Vamos aos cálculos
Para da início as atividades, se necessário fazer a istalação dos pacotes (geoR) e (MASS), e a requerição dos mesmos.
if(!require(geoR)){install.packages("geoR")}
if(!require(MASS)){install.packages("MASS")}
if(!require(moments)){install.packages("moments")}
if(!require(akima)){install.packages("akima")}
if(!require(car)){install.packages("car")}
setwd("F:/Estatistica Espacial")
dados <- read.table("Km_PB.txt", head=T)
verificando parte dos dados e caracteristicas do conjunto de dados
# head(dados)
# dim(dados)
# colnames(dados)
# rownames(dados)
# dimnames(dados)
Bd <-as.geodata(dados)
# summary(Bd)
| min |
13.55211 |
17.54344 |
| max |
436.42574 |
235.04460 |
| Precipitação |
241 |
37.4000 |
74.3500 |
106.9000 |
108.1937 |
142.2000 |
194.0000 |
borda=read.table("bord_Km.txt",h=T)
Bd$borders <- with(borda, cbind(x,y))
plot(Bd,lowess = T)
points(Bd, pt.div="quartiles", xlab="Leste", ylab="Norte")
verificando e efeito de tendência por meio dos gráficos.
plot(Bd,lowess=T,trend="1st")
plot(Bd,lowess=T,trend="1st",lam=0.5)
plot(Bd,lowess=T,trend=~coords[,2])
plot(Bd,lowess=T,trend=~coords[,2],lam=0.5)
plot(Bd,lowess=T,trend=~coords[,1])
|Efeito de Tendência na Coordenada x com \(\lambda\) = 0.5| |——————————————————-|
plot(Bd,lowess=T,trend=~coords[,1],lam=0.5)
plot(Bd,lowess=T,trend=~coords[,1]+~coords[,2])
plot(Bd,lowess=T,trend=~coords[,1]+~coords[,2],lam=0.5)
plot(Bd,lowess=T,trend="2nd")
plot(Bd,lowess=T,trend="2nd",lam=0.5) # Melhor Trasformação!
plot(Bd,lowess=T,trend=~coords[,1]+~coords[,2])
plot(Bd,lowess=T,trend=~coords[,1]+~coords[,2],lam=0.5)
Aparentemente a tendêndia quadrática de Segunda Ordem com \(\lambda\) = 0.5 obteve melhores resultados, dessa forma será utilizada para an-alises sub-sequentes.
Bdv1=variog(Bd,trend="2nd",lam=0.5, max.dist = 140)
## variog: computing omnidirectional variogram
plot(Bdv1)
A distância máxima utilizada no Semivariograma foi de 140 Km.
var.1 = variog(Bd,low=T,trend="2nd",lam=0,max.dist = 140);plot(var.1)
## variog: computing omnidirectional variogram
var.2 = variog(Bd,low=T,trend="2nd",lam=0,max.dist = 140,option="cloud");plot(var.2)
## variog: computing omnidirectional variogram
var.3 = variog(Bd,low=T,trend="2nd",lam=0,max.dist = 140,uvec = seq(0,140,l=20));plot(var.3)
## variog: computing omnidirectional variogram
env.var=variog.mc.env(Bd,obj.variog = var.3,nsim=999)
## variog.env: generating 999 simulations by permutating data values
## variog.env: computing the empirical variogram for the 999 simulations
## variog.env: computing the envelops
plot(var.3,env=env.var)
Utilizando como base de interpretação o P-Valor de 0,005 com (999) simulações, observa-se que constatou-se a dependência espacial dos dados por meio das observações que se alocam fora do envelope.
Modelagem das Estimativas de máxima verossimilhança para as Funções de Correlação:
mod0<-likfit(Bd,trend="1st",kappa = 0.5,ini=c(421.43,69.13),nug=202.75 ,cov.model="matern",lam=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.
# summary(mod0)
| 17.7184 |
0.0014 |
0.0178 |
1.2592 |
15.5944 |
181.3035 |
mod1=likfit(Bd,kappa = 0.5,ini=c(421.43,69.13 ),nug=202.75 ,cov.model="matern",lam=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.
# summary(mod1)
mod2=likfit(Bd,trend="2nd",kappa = 0.5,ini=c(421.43,69.13),nug=202.75 ,cov.model="matern",lam=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.
# summary(mod2)
| 23.4502 |
-0.0830 |
0.0193 |
0.0002 |
1.2613 |
5.9594 |
69.1027 |
| \(\beta4\) |
\(\beta5\) |
|
| ——– |
——– |
|
| 0.0000 |
0.0000 |
|
mod3=likfit(Bd,trend=~coords[,2],kappa = 0.5,ini=c(421.43,69.13),nug=202.75 ,cov.model="matern",lam=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.
# summary(mod3)
| 17.9928 |
0.0178 |
1.2573 |
15.2462 |
176.7026 |
mod4=likfit(Bd,trend=~coords[,1],kappa = 0.5,ini=c(421.43,69.13),nug=202.75 ,cov.model="matern",lam=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.
# summary(mod4)
| 20.219 |
0.002 |
1.252 |
18.798 |
217.107 |
mod5=likfit(Bd,trend=~coords[,1]+~coords[,2],kappa = 0.5,ini=c(421.43,69.13 ),nug=202.75 ,cov.model="matern",lam=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.
# summary(mod5)
| 17.7184 |
0.0014 |
0.0178 |
1.2592 |
15.5944 |
181.3035 |
Exp0<-likfit(Bd,trend="1st",ini=c(594,121.79),nug=220.75,lam=0.5,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.
# summary(Exp0)
| 17.6891 |
-0.0010 |
0.0184 |
1.2253 |
11.0947 |
121.7495 |
Exp1=likfit(Bd,ini=c(594,121.79),nug=220.75,lam=0.5,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.
# summary(Exp1)
| 19.936 |
1.203 |
11.444 |
121.753 |
Exp2=likfit(Bd,trend="2nd",ini=c(594,121.79),nug=220.75,lam=0.5,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.
# summary(Exp2)
| 22.8423 |
-0.0801 |
0.0162 |
0.0002 |
1.3312 |
9.2951 |
121.7534 |
| \(\beta4\) |
\(\beta5\) |
|
| ——– |
——– |
|
| 0.0000 |
0.0000 |
|
Exp3=likfit(Bd,trend="2nd",ini=c(594,121.79),nug=220.75,lam=0.5,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.
# summary(Exp3)
| 22.8423 |
-0.0801 |
0.0162 |
0.0002 |
1.3312 |
9.2951 |
121.7534 |
| \(\beta4\) |
\(\beta5\) |
|
| ——– |
——– |
|
| 0.0000 |
0.0000 |
|
Exp4=likfit(Bd,trend=~coords[,1],ini=c(594,121.79),nug=220.75,lam=0.5,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.
# summary(Exp4)
| 20.2335 |
-0.0014 |
0.0178 |
1.2592 |
15.5944 |
Tomando como base para o critério de escolha do melhor modelo, aquele que apresenta menor AIC, dessa forma, o que apresentou esta característica Foi o Modelo Matern “mod2” com AIC = 2076, isso implica dizer que para a precipitação média dos dados apresentaram tendência de segunda ordem.
Definindo a Grid
par(mfrow=c(1,1))
mod2=likfit(Bd,trend="2nd",kappa = 0.5,ini=c(421.43,69.13),nug=202.75 ,cov.model="matern",lam=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.
gr=pred_grid(Bd$borders,by=3)
points(Bd)
points(gr,col=2,pch=19,cex=.3)
gr0=gr[.geoR_inout(gr,Bd$borders),]
points(gr0,col=3,pch=19,cex=0.3)
Processo de Krigagem
O processo de Krigagem permite construir mapas de incerteza espacial, no qual fornece uma visão global da congruência dos resultados.
kc <- krige.conv(Bd, loc=gr, bor=borda, krige=krige.control(obj=mod2))
## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with mean given by a 2nd order polynomial on the coordinates
## krige.conv: performing the Box-Cox data transformation
## krige.conv: back-transforming the predicted mean and variance
## krige.conv: Kriging performed using global neighbourhood
# names(kc)
kc$pred[1:10]
## [1] 60.05314 59.26807 58.45540 57.62099 61.31288 60.57698 59.79977
## [8] 58.96019 58.07335 57.17339
image(kc,main="Mapa de Interpolação",x.leg=c(10,150),y.leg=c(0,15),xlab="Sul" ,ylab="Oeste")
As cores mais claras indicam a região que recebe maior precipiutação, neste caso no Sertão e no litoral.
image(kc, main="Mapa de Interpolação",col=gray(seq(1,0.2,len=10)),x.leg=c(10,150),y.leg=c(0,15) ,xlab="Sul" ,ylab="Oeste")
Este gráfico indica sa mesmas caracteristicas do anterior, mudadando apenas para a cor na escala de cinza.
kc1 <- krige.conv(Bd, loc=gr, bor=borda, krige=krige.control(obj=mod2))
## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with mean given by a 2nd order polynomial on the coordinates
## 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,main="Mapa de Interpolação",col= terrain.colors(200),x.leg=c(10,150),y.leg=c(0,15), xlab="Sul",ylab="Oeste")
contour(kc1,main="Mapa de Interpolação",c=T,col= terrain.colors(200), nlev = 200)
numlevel=200
contour(kc1,main="Mapa de Interpolação",f=T,col= terrain.colors(200), nlev = 150)
numlevel=50
contour(kc1,main="Mapa de Interpolação",c=T,col= adjustcolor(4, .4), nlev = 200)
numlevel=200
