Rotina no geoR

A seguir iremos apresentar algumas funções dos pacotes geoR e gstat para estimação paramétrica dos modelos de semivariograma:

require(geoR)
## Loading required package: geoR
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.7-5.2 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------

Dados simulados no geoR

summary(s100)
## Number of data points: 100 
## 
## Coordinates summary
##         Coord.X    Coord.Y
## min 0.005638006 0.01091027
## max 0.983920544 0.99124979
## 
## Distance summary
##         min         max 
## 0.007640962 1.278175109 
## 
## Data summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.1680  0.2730  1.1050  0.9307  1.6100  2.8680 
## 
## Other elements in the geodata object
## [1] "cov.model" "nugget"    "cov.pars"  "kappa"     "lambda"

Gráfico descritivo dos dados

plot(s100)

Calculando a distância máxima, e utilizando sua metade no semivariograma

dmax = max(summary(s100)$distances.summary)
dmax
## [1] 1.278175
Calculando o semivariograma amostral \[\begin{eqnarray} \hat{\gamma}(\mathbf{h}) = \frac{1}{2|N(\mathbf{h})|}\displaystyle\sum_{N(\mathbf{h})}^{}\left\{Z(\mathbf{s}_{i}) - Z(\mathbf{s}_{j})\right\}^{2} \end{eqnarray}\] 
vario100 <- variog(s100, max.dist=dmax/2, estimator.type="classical")
## variog: computing omnidirectional variogram

Gráfico do semivariograma

plot(vario100, xlab="Distância (h)", ylab="Semivariograma", bty='l',
     cex=2, pch=19, cex.lab=1.5, cex.axis=1.5)
grid(col='lightseagreen', lwd=1)

Como o objetivo é comparar diferentes métodos de estimação, iremos nos concentar apenas no modelo de covariãncia exponencial: \[\begin{eqnarray} C(h) = \sigma^{2}exp\left\{-\frac{h}{\phi}\right\} \end{eqnarray}\]

Note que, o modelo exponencial tem dois parâmetros, \(\sigma^{2}\) e \(\phi\) que serão estimados. Sendo assim, precisamos atribuir valores iniciais para estes parâmetros, uma vez que, o processo de estimação é iterativo.

ini.vals <- expand.grid('phi'=seq(0,1,l=5), 'sigma2'=seq(0,1,l=5))

Ajuste utilizando Mínimos Quadrados Ordinários (OLS). Para utilizarmos o método OLS basta atribuir o argumento (wei=‘equal’)

ols <- variofit(vario100, ini=ini.vals, fix.nug=TRUE, wei='equal',
                cov.model='exponential', messages=F)

Ajuste utilizando Mínimos Quadrados Ponderados (WLS). Pesos sugeridos por Cressie (1985). Para utilizarmos o método WLS é só atribuir o argumento (wei=‘cressie’).

wls <- variofit(vario100, ini=ini.vals, fix.nug=TRUE, wei='cressie',
                cov.model='exponential', messages=F)

A seguir iremos apresentar ajustes de modelos utilizando máxima verossimilhança. Note que, este tipo de abordagem não necessita da construção de um semivariograma amostral, uma vez que, utiliza-se dos próprios dados.

Para utilizar Máxima Verossimilhança (ML) basta colocar o argumento lik.met = ‘ML’

ml <- likfit(s100, ini=ini.vals, fix.nug=TRUE, lik.met = 'ML',
                cov.model="exponential", messages=F)

Para utilizar Máxima Verossimilhança Restrita (REML) basta colocar o argumento lik.met = ‘REML’

reml <- likfit(s100, ini=ini.vals, fix.nug=TRUE, lik.met = 'REML',
               cov.model="exponential", messages=F)

Resultados apresentados em um gráfico

plot(vario100, xlab="Distância (h)", ylab="Semivariograma", 
     cex.main=1.3, main="Modelo Exponencial", bty='l',
     cex=2, pch=19, cex.lab=1.3, cex.axis=1.5)
grid(col='lightseagreen', lwd=1)

lines(ml, col=2, lwd=2)
lines(reml, col=3, lwd=2)
lines(ols, col=4, lwd=2)
lines(wls, col=5, lwd=2)

legend('topleft', c("ML", "REML", "OLS", "WLS"), col = 2:5,
       lwd = rep(2,4), bty='n', cex=1.2)

Rotina no gstat

No gstat, os gráficos descritivos são realizados utilizando rotinas do pacote sp.

require(sp)
## Loading required package: sp

Entramos com os dados da seguinte forma:

x <- s100[[1]][,1]
y <- s100[[1]][,2]
z <- s100[[2]]

dados_gstat <- data.frame('x'=x, 'y'=y, 'z'=z)

Escrevendo quem são as coordenadas

coordinates(dados_gstat) = ~ x + y
head(dados_gstat)
##               coordinates          z
## 1   (0.8071267, 0.945446)  0.9171888
## 2  (0.5499981, 0.6832649)  1.1483234
## 3  (0.3408055, 0.4585089)  1.0327563
## 4  (0.1370993, 0.4720083)  0.1219548
## 5 (0.04418569, 0.1223202)  0.6152988
## 6 (0.02781636, 0.8037459) -0.5506065

Resumo descritivo dos dados

summary(dados_gstat)
## Object of class SpatialPointsDataFrame
## Coordinates:
##           min       max
## x 0.005638006 0.9839205
## y 0.010910267 0.9912498
## Is projected: NA 
## proj4string : [NA]
## Number of points: 100
## Data attributes:
##        z          
##  Min.   :-1.1677  
##  1st Qu.: 0.2730  
##  Median : 1.1046  
##  Mean   : 0.9307  
##  3rd Qu.: 1.6102  
##  Max.   : 2.8679

Análise gráfica

bubble(dados_gstat, "z", main='')

Construção do semivariograma amostral

require(gstat)
## Loading required package: gstat
lzn.vgm = variogram(z~1, dados_gstat, cutoff=dmax/2)
plot(lzn.vgm, xlab="Distância (h)", ylab="Semivariograma", 
     cex=2, pch=19)

Modelos de semivariogramas presentes no gstat

vgm()
##    short                                      long
## 1    Nug                              Nug (nugget)
## 2    Exp                         Exp (exponential)
## 3    Sph                           Sph (spherical)
## 4    Gau                            Gau (gaussian)
## 5    Exc        Exclass (Exponential class/stable)
## 6    Mat                              Mat (Matern)
## 7    Ste Mat (Matern, M. Stein's parameterization)
## 8    Cir                            Cir (circular)
## 9    Lin                              Lin (linear)
## 10   Bes                              Bes (bessel)
## 11   Pen                      Pen (pentaspherical)
## 12   Per                            Per (periodic)
## 13   Wav                                Wav (wave)
## 14   Hol                                Hol (hole)
## 15   Log                         Log (logarithmic)
## 16   Pow                               Pow (power)
## 17   Spl                              Spl (spline)
## 18   Leg                            Leg (Legendre)
## 19   Err                   Err (Measurement error)
## 20   Int                           Int (Intercept)

Uma vantagem no gstat é que temos várias opções para a matriz de pesos para o MMQ.

Para o OLS basta fazer fit.method = 6.

olsExp <- fit.variogram(lzn.vgm, vgm(psill=0.8, "Exp", range=0.7),
                        fit.method = 6)
plot(lzn.vgm, olsExp, xlab="Distância (h)", ylab="Semivariograma", 
     cex=2, pch=19, lwd=2, main="OLS")

Já para o WLS, como sugerido por Cressie (1985), fit.method = 2.

wlsExp <- fit.variogram(lzn.vgm, vgm(psill=0.8, "Exp", range=0.7),
                        fit.method = 2)
plot(lzn.vgm, wlsExp, xlab="Distância (h)", ylab="Semivariograma", 
     cex=2, pch=19, lwd=2, main="WLS")

Para REML veja o que o gstat diz:

I guess there is much more to likelihood variogram fitting in package geoR, and probably also in nlme.

dados <- data.frame(dados_gstat)
reml <- fit.variogram.reml(z~1, ~x+y, dados, 
                   model = vgm(psill=wlsExp$psill, "Exp",
                               range=wlsExp$range))
plot(lzn.vgm, reml, xlab="Distância (h)", ylab="Semivariograma", 
     cex=2, pch=19, lwd=2, main="REML")

Mínimos quadrados generalizados no gstat

gls <- fit.variogram.gls(z~1, dados_gstat[1:40,], 
                         vgm(psill=wlsExp$psill, "Exp",       
                             range=wlsExp$range), trace = F, cutoff=dmax/2)
plot(lzn.vgm, gls, xlab="Distância (h)", ylab="Semivariograma", 
     cex=2, pch=19, lwd=2, main="GLS")