Geoestatística (parte 2)

1 Simulando um conjunto de dados

### Estatistica Espacial I
### Aula 09

library(geoR)

#Dados de precipitação no Estado do Paraná

#Análise exploratória
data(parana)
plot(parana)

points(parana, cex.min=1, cex.max=5, col=terrain.colors(5), ps="quintiles")

#Obtenção do variograma
variograma <- variog(parana)
## variog: computing omnidirectional variogram
plot(variograma)

#Definição da família e busca por valores iniciais
#eyefit(variograma)

#Estimando via MV
#opções para trend: cte (média constante), 1st (tendência linear) e 2st (tendência quadrática).

#Estimando o efeito pepita:
mv <- likfit(parana, ini=c(4000,100), trend="cte", cov.model="gaussian",
             nugget=0)
## 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.
#Fixando o efeito pepita em algum valor arbitrário:
mv.fix <-  likfit(parana, ini=c(4000,150) , trend="cte", cov.model="gaussian",
                  fix.nugget=T, nugget=400)
## 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.
plot(variograma)
lines(mv,lty=1,lwd=2,col="blue")
lines(mv.fix,lty=1,lwd=2,col="red")

########## Analisando os parâmetros estimados
summary(mv)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##     beta 
## 260.5917 
## 
## Parameters of the spatial component:
##    correlation function: gaussian
##       (estimated) variance parameter sigmasq (partial sill) =  6869
##       (estimated) cor. fct. parameter phi (range parameter)  =  336.2
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (estimated) nugget =  521.1
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 1 (no transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 581.861
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
## "-669.3"      "4"   "1347"   "1358" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
## "-781.5"      "2"   "1567"   "1573" 
## 
## Call:
## likfit(geodata = parana, trend = "cte", ini.cov.pars = c(4000, 
##     100), nugget = 0, cov.model = "gaussian")
summary(mv.fix)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##     beta 
## 266.3647 
## 
## Parameters of the spatial component:
##    correlation function: gaussian
##       (estimated) variance parameter sigmasq (partial sill) =  7260
##       (estimated) cor. fct. parameter phi (range parameter)  =  321.3
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (fixed) nugget = 400
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 1 (no transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 556.1328
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
## "-671.7"      "3"   "1349"   "1358" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
## "-781.5"      "2"   "1567"   "1573" 
## 
## Call:
## likfit(geodata = parana, trend = "cte", ini.cov.pars = c(4000, 
##     150), fix.nugget = T, nugget = 400, cov.model = "gaussian")
par(mfrow=c(1,3))
profile.sill<-proflik(mv,parana,sill.val=seq(6868,6869,l=50))
## proflik: computing profile likelihood for the sill
plot(profile.sill)

profile.phi<- proflik(mv,parana,range.val=seq(250,400,l=20))
## proflik: computing profile likelihood for the range
plot(profile.phi)

profile.nug<- proflik(mv,parana, nugget.val=seq(400,600,l=20))
## proflik: computing profile likelihood for the nugget
plot(profile.nug)

#envelopes de confiança
limites <- variog.model.env(parana, obj.var=variograma, model=mv,
                            nsim=30)
## variog.env: generating 30 simulations (with  143 points each) using the function grf
## variog.env: adding the mean or trend
## variog.env: computing the empirical variogram for the 30 simulations
## variog.env: computing the envelops
par(mfrow=c(1,1))
plot(variograma, envelope=limites, lwd=2)

### Interpolação ponderada pelo inverso da distância

library(gstat)
library(sp)

dados <- as.data.frame(cbind(parana$coords,parana$data))
coordinates(dados) <- c("east", "north")

#inverse distance weighted interpolation (idp is the weighting power)
summary(parana$coords)
##       east           north       
##  Min.   :150.1   Min.   : 70.36  
##  1st Qu.:263.9   1st Qu.:163.33  
##  Median :366.9   Median :223.70  
##  Mean   :399.2   Mean   :243.01  
##  3rd Qu.:512.7   3rd Qu.:319.76  
##  Max.   :768.5   Max.   :461.97
leste = seq(100,800,1)
norte = seq(0,600,1)
nx = length(leste)
ny = length(norte)
# criando o grid de pontos para interpolar
grid = as.data.frame(expand.grid(leste,norte))

# coordenadas dos pontos do grid
coordinates(grid) = ~Var1+Var2


# interpolacao
interpol <- idw(V3 ~ 1, dados, newdata = grid, idp = 2)
## [inverse distance weighted interpolation]
interpol.matriz = matrix(interpol$var1.pred,nx,ny)

par(mfrow=c(1,1))
# grafico
#opcao 1
filled.contour(leste, norte,interpol.matriz,color=terrain.colors,
               plot.axes={axis(1); axis(2);points(parana$coords)})

#opcao2
image(leste, norte, interpol.matriz,col=terrain.colors(256),
      zlim=c(min(interpol$var1.pred),max(interpol$var1.pred)))

points(parana, cex.min=1, cex.max=5, ps="quintiles", col=terrain.colors(5),add=T)

#contour(leste,norte,interpol.matriz,nlevels=56,labcex=1,col=2)


### Krigagem no geoR
#opções: SK (simples), OK (ordinária) e UK (universal)

leste = seq(100,800,10)
norte = seq(0,600,10)
grid <- expand.grid(leste,norte)

# krigagem
krigagem <- krige.conv(parana, loc=grid, 
                       krige=krige.control(type.krige="SK",trend.d="cte", 
                                           beta=mv$beta, nugget=mv$nugget, 
                                           cov.pars=mv$cov.pars))
## 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(krigagem)

# media
image(krigagem, loc=grid, coords=parana$coords, 
      col= terrain.colors(256), x.leg=c(650,700),y.leg=c(400,550), 
      vertical=TRUE)

# variancia
image(krigagem, loc=grid, coords=parana$coords, 
      values=krigagem$krige.var, col=terrain.colors(256), 
      x.leg=c(650,700), y.leg=c(400,550), vertical=TRUE)

#contour(krigagem)

# A opção "borders" precisa ser adicionada na função "krige.conv" 
# caso o seu objeto geodata não contenha o polígono.





##### Validação cruzada

valid <- xvalid(parana, model=mv)
## xvalid: number of data locations       = 143
## xvalid: number of validation locations = 143
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 
## xvalid: end of cross-validation
plot(valid, parana)

#EQM de predição
EQM.pred <- mean(valid$error^2)