Leandro Valter
leandro.vvalter@gmail.com
Universidade Estadual da Paraíba
Centro de Ciência e Tecnologia
Departamento de Estatística

Dados

setwd("C:/Users/Win7/Desktop/Estatística Espacial/Chuva")
dados=read.table("Km_PB.txt",h=T)
borda=read.table("bord_Km.txt")
attach(dados)

Pacotes

require(akima)
require(geoR)
require(MASS)
require(moments)
require(car)

Formatando para geodata+borda

chuva=as.geodata(dados, coords.col=c(1, 2), data.col=3)
chuva$borders <- with(borda, cbind(V1,V2))

Descritiva

Tabela 1:Estatísticas descritivas: média, mínimo, máximo, desvio padrão (D.P), coeficiente de variação (C.V.), curtose e assimetria para os dados de precipitação média.

Variável N Média Mediana Mínimo Máximo D.P Cs Ck Cv%
Precipitação 241 108,19 106,90 37,40 194,00 40,03 0,13 -1,13 36,99

O coeficiente de variação foi maior 20% representando uma baixa homogeneidade. A precipitação apresentou uma assimetria à esquerda (positiva). Verificando a curtose, a precipitação possue uma função de distribuição platicúrtica, ou seja, é mmais achatada do que a distribuição normal. Provavelmente será necessário fazer o uso de algum artifício para corrigir as baixas homogeneidades, as assimetrias positivas e o achatamento maior que a da distribuição normal. Será utilizada a transformação Box-Cox, caso seja necessário.

Descritiva Espacial

summary(chuva)
## Number of data points: 241 
## 
## Coordinates summary
##             x         y
## min  13.55211  17.54344
## max 436.42574 235.04460
## 
## Distance summary
##         min         max 
##   0.9745422 428.5703722 
## 
## Borders summary
##         V1     V2
## min   0.00   0.00
## max 439.69 251.95
## 
## Data summary
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  37.4000  74.3500 106.9000 108.1937 142.2000 194.0000

A distância máxima é de 428.57 km, é usual ter um terço da distância máxima nos semivariogramas, porém se tratando de quilometros será usado menos que um terço, pois o raio de visualização pode ser que se extenda mais que a região de estudo.

plot(chuva,low=T)

points(chuva, pt.div="quartiles", xlab="Leste", ylab="Norte")

Aparentemente tem-se tendência nas coordenadas, após a verificação se vai ser necessário usar alguma transformação, dado a baixa homogeneidade dos dados, será estudado qual o tipo de tendência.

Investigando uso de transformação

boxcox(chuva)

Talvez trasformação raiz quadradado com \(\lambda=0.5\) seja interessante.

Investigando a tendência

Linear

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

plot(chuva,low=T,trend="1st",lam=0.5)

Apenas na longitude

plot(chuva,low=T,trend=~x)

plot(chuva,low=T,trend=~x,lam=0.5)

Apenas na latitude

plot(chuva,low=T,trend=~y )

plot(chuva,low=T,trend=~y,lam=0.5)

Quadrática

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

plot(chuva,low=T,trend="2nd",lam=0.5) #melhor

Quadrática pra lon e lat

plot(chuva,low=T,trend=~ x + I(x^2),lam=0.5)

plot(chuva,low=T,trend=~ y + I(y^2),lam=0.5)

A tendência quadrática com transformação raiz quadrada foi a que melhor verificou-se aos dados, nas análises seguintes esta será utilizada.

Semivariograma Teórico

var.chuva=variog(chuva,low=T,trend="2nd",lam=0.5,max.dist = 90)
## variog: computing omnidirectional variogram
plot(var.chuva)

A distância máxima utilizada no semivariograma foi a de 90 (km).

Dependência Espacial

par(mfrow=c(2,2))
var1=variog(chuva,low=T,trend="2nd",lam=0.5,max.dist = 90);plot(var1)
## variog: computing omnidirectional variogram
var2=variog(chuva,low=T,trend="2nd",lam=0.5,max.dist = 90,option = "cloud");plot(var2)
## variog: computing omnidirectional variogram
var3=variog(chuva,low=T,trend="2nd",lam=0.5,max.dist = 90,uvec=seq(0,90,l=20))
## variog: computing omnidirectional variogram
plot(var3)
env.var=variog.mc.env(chuva,obj.variog = var3,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(var3,env=env.var)

Com um p-valor de 0,005 (999 simulações) constatou-se a dependência espacial dos dados pelas observações fora do envelope.

Semivariogramas Teóricos

par(mfrow=c(2,4))
#cauchy
var.chuva=variog(chuva,low=T,trend="2nd",lam=0.5,max.dist = 90)
## variog: computing omnidirectional variogram
plot(var.chuva,ylab=expression("Semivariância ("*Cauchy*")"),
     xlab="Distancia (Km)", main = "Semivariograma")
lines.variomodel(cov.model="cauchy", cov.pars=c(4.59,18.71),nug=2, max.dist=90,kappa =.5,lam=.5)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "lam" não é um
## parâmetro gráfico
# circular
plot(var.chuva,ylab=expression("Semivariância ("*Circular*")"),
     xlab="Distancia (Km)", main = "Semivariograma")
lines.variomodel(cov.model="circular",lam=.5, cov.pars=c(3.22,44.44), 
                 nug=1.84, max.dist=90)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "lam" não é um
## parâmetro gráfico
#cubic ######################
plot(var.chuva,ylab=expression("Semivariância ("*Cubico*")"),
     xlab="Distancia (Km)", main = "Semivariograma")
lines.variomodel(cov.model="cubic",lam=.5, cov.pars=c(2.91,51.45), 
                 nug=2, max.dist=90)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "lam" não é um
## parâmetro gráfico
#matern k=0.5  
plot(var.chuva,ylab=expression("Semivariância ("*Matern*"k=.5)"),
     xlab="Distancia (Km)", main = "Semivariograma")
lines.variomodel(cov.model="matern",lam=.5, cov.pars=c(3.22,25.73), 
                 nug=2, max.dist=90,kappa = 0.5)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "lam" não é um
## parâmetro gráfico
# gaussian
plot(var.chuva,ylab=expression("Semivariância ("*Gaussiano*")"),
     xlab="Distancia (Km)", main = "Semivariograma")
lines.variomodel(cov.model="gaussian", cov.pars=c(2.91,23.39), 
                 nug=2.14, max.dist=90,lam=.5)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "lam" não é um
## parâmetro gráfico
#matern k=1  
plot(var.chuva,ylab=expression("Semivariância ("*Matern*"k=1)"),
     xlab="Distancia (Km)", main = "Semivariograma")
lines.variomodel(cov.model="matern", cov.pars=c(3.52,18.71 ), 
                 nug=2, max.dist=90,kappa = 1,lam=.5)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "lam" não é um
## parâmetro gráfico
#powered.exponential   
plot(var.chuva,ylab=expression("Semivariância ("*Exp.Exponencial*")"),
     xlab="Distancia (Km)", main = "Semivariograma")
lines.variomodel(cov.model="powered.exponential", cov.pars=c(3.22,25.73), 
                 nug=2, max.dist=90,kappa = 1.49,lam=.5)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "lam" não é um
## parâmetro gráfico
# spherical   
plot(var.chuva,ylab=expression("Semivariância ("*Esferico*")"),
     xlab="Distancia (Km)", main = "Semivariograma")
lines.variomodel(cov.model="spherical", cov.pars=c(3.06,51.45 ), 
                 nug=2, max.dist=90,lam=.5)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "lam" não é um
## parâmetro gráfico

Pelo semivariograma teórico, visualmente as funções que melhor se ajustaram foram as Cauchy, Exponencial Exponenciada, Matern (kappa=1).

Modelagem

Tabela 2:Estimativas de máxima verossimilhança dos parâmetros associados ao modelo com diferentes funções de correlação e efeito de tendência.

Modelo \(\hat\beta_0\) \(\hat\beta_1\) \(\hat\beta_2\) \(\hat\beta_3\) \(\hat\tau^2\) \(\hat\sigma^2\) \(\hat\phi\) AIC BIC
Matern(Kappa=0.5) 22.5945 -0.0769 0.0002 0.0180 1.2689 6.3359 74.5726 2072.4140 2096.8080
Matern(Kappa=1) 23.4606 -0.0801 0.0002 0.0178 1.5722 5.3345 32.2864 2071.8180 2096.2110
Cauchy 23.1713 -0.0787 0.0002 0.0178 1.6755 5.8675 27.8870 2071.9760 2096.3700
Circular 24.6791 -0.0862 0.0002 0.0179 1.2102 4.3699 58.4410 2080.4410 2104.8350
Cúbico 24.5089 -0.0848 0.0002 -0.0848 1.8210 4.7453 92.6464 2078.2720 2102.665
Esférico 24.1202 -0.0825 0.0002 0.0183 1.2584 4.9001 85.0320 2078.4930 2102.8870
Exp. Exponenciada 23.8638 -0.0818 0.0002 0.0178 1.5594 5.1265 46.4417 2073.2400 2097.6340
Gaussiano 24.4831 -0.0846 0.0002 0.0177 1.8194 4.4007 36.8616 2076.9310 2101.3250

Para a precipitação média os dados apresentaram tendência quadrática na longitude e uma tendência linear na latitude. O melhor modelo pelos critérios de seleção foi o Matern (Kappa=1), apresentando um menor AIC e BIC.

Krigagem

Modelo escolhido pelo menor AIC/BIC

ml3.3=likfit(chuva,cov.model="matern",ini.cov.pars=c(3.52,18.71),nug=2,kappa = 1,lam=.5,trend =~x+I(x^2)+y)
## ---------------------------------------------------------------
## 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(chuva$borders,by=1)
points(chuva)
points(gr,col=2,pch=19,cex=.3)
gr0=gr[.geoR_inout(gr,chuva$borders),]
points(gr0,col=3,pch=19,cex=0.3)
title("Malha de pontos de predição")

dim(gr)
## [1] 110880      2
dim(gr0)
## [1] 56526     2
KC=krige.control(type="OK",obj.model=ml3.3)
kc1=krige.conv(chuva,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,main="Mapa de interpolação",col=terrain.colors(200),x.leg = c(10,150),
      y.leg =c(0,15),xlab = "Sul", ylab = "Oeste")

Gerando simulações

gr <- pred_grid(chuva$borders, by=3)
chuva.kc <- krige.conv(chuva, loc=gr, krige=krige.control(obj=ml3.3), out=output.control(n.pred=1000))
## 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

calculando probabilidades e quantis no grid de predição

chuva.loci1 <- krige.conv(chuva, loc=gr, krige=krige.control(obj=ml3.3), out=output.control(n.pred=1000, thre=106.9, quant=c(0.10, .5, .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(chuva.loci1$simul)
## [1] 6274 1000
hist(chuva.loci1$simul,probability = T)

Representando probabilidades no mapa

image(chuva.loci1, border=chuva$borders, loc=gr, col=gray(seq(1,0,l=21)), coords.dat=chuva$coords, val=1-chuva.loci1$prob,x.leg = c(10,150),
      y.leg =c(0,15))

A probabilidade de uma região do estado receber uma maior precipitação baseando-se na mediana (106.9) é maior no sertão e no litoral do estado.

## mapeando o quantil 0.1
image(chuva.loci1, border=chuva$borders, loc=gr, col=gray(seq(1,0,l=21)), coords.dat=chuva$coords, val=chuva.loci1$quan[,1], x.leg = c(10,150),
      y.leg =c(0,15))

## mapeando a mediana
image(chuva.loci1, border=chuva$borders, loc=gr, col=gray(seq(1,0,l=21)), coords.dat=chuva$coords, val=chuva.loci1$quan[,2], x.leg = c(10,150),
      y.leg =c(0,15))

## mapeando o quantil 0.9
image(chuva.loci1, border=chuva$borders, loc=gr, col=gray(seq(1,0,l=21)), coords.dat=chuva$coords, val=chuva.loci1$quan[,3], x.leg = c(10,150),
      y.leg =c(0,15))