Geoestatística

1 Simulando um conjunto de dados

1.1 Coordenadas

Crie um conjunto de dados com 500 pontos observados supondo que as coordenadas sejam a latitude e a longitude e que cada uma dessas coordenadas seja um número qualquer entre 0 e 1. Faça o gráfico dos pontos gerados.

#Simulando um dado geoestatistico

#criando uma grade regular
n   = 500
lat = runif(n,0,1) #gerando valores de "latitudes" 
lon = runif(n,0,1) #gerando valores de "longitudes"
s   = cbind(lat, lon) #vetor de coordenadas
plot(s, xlab="latitude", ylab="longitude")

1.2 Distância

Calcule a distância dos pontos gerados e faça um histograma dessas distâncias.

#calculando as distancias entre os pontos gerados
d = matrix(NA, n, n)
for(i1 in 1:n){
  for(i2 in 1:n){
    d[i1,i2] = sqrt((lat[i1] - lat[i2])^2 + (lon[i1] - lon[i2])^2)
    #d[i1,i2] = sqrt((s[i1,1] - s[i2,1])^2 + (s[i1,2] - s[i2,2])^2)
  }
}
par(mfrow=c(1,2))
hist(d, freq=F, main="Histograma", ylab="densidade", xlab="distâncias")
hist(d[row(d) < col(d)], freq=F, main="Histograma", ylab="densidade", xlab="distâncias")

1.3 Função de covariância

Crie funções para calcular a função de correlação e a função de covariância esférica, exponencial e Matérn. Para isso, use \(\phi=0,35\), \(V=1\) e \(\tau=0,1\). Depois desenhe a função de correlação, o covariograma e o semivariograma.

#criando a funcao de correlacao e covariancia esférica
V     = 1
tau   = 0.1
phi   = 0.35
CorEsf = function(d, phi){  
  saida = 1 - 1.5/phi*d+0.5*(d/phi)^3
  saida[d>phi] = 0
  return( saida )   
  }
CovEsf = function(d, phi, V, tau){  
  saida       = V * CorEsf(d,phi)
  saida[d==0] = V + tau
  return( saida )   
}
# usando uma funcao existente do R para a funcao de correlacao e 
# comparando com as funções criadas para conferir se estão corretas

#install.packages("geoR")
library(geoR)
u = seq(0, max(d), l=1000)
y <- cov.spatial(u, cov.model= "spherical",cov.pars=c(V, phi))
plot(y, CorEsf(u, phi))
lines(y,y, col="red", lwd=2)

#criando a funcao de correlacao e covariancia exponencial
CorExp = function(d, phi){  return( exp(-d/phi) )   }
CovExp = function(d, phi, V, tau){  
  saida       = V * CorExp(d,phi)
  saida[d==0] = V + tau
  return( saida )   
}

# usando uma funcao existente do R para a funcao de correlacao e 
# comparando com as funções criadas para conferir se estão corretas
y <- cov.spatial(u, cov.model= "exponential",cov.pars=c(V, phi))
plot(y, CorExp(u, phi))
lines(y,y, col="red", lwd=2)

#Desenhando as funcoes criadas
par(mfrow=c(2,3))
u = seq(0, max(d), l=1000)
plot(u, CorEsf(u,phi), type="l", ylab=expression(rho), 
     bty="n", lwd=2, xlab="distância", main="Fção de correlação esf.")
abline(v=phi, lty=3, lwd=2, col="red")
abline(h=0, lty=3, lwd=2, col="red")

plot(u, CovEsf(u, phi, V, tau), type="l", ylab="C(distância)", 
     bty="n", lwd=2, xlab="distância", main="Covariograma")
abline(v=phi, lty=3, lwd=2, col="red")
abline(h=0, lty=3, lwd=2, col="red")

plot(u, tau+V*(1-CorEsf(u,phi)), type="l", ylab="gamma(distância)", 
     bty="n", lwd=2, xlab="distância", main="Semivariograma")
abline(v=phi, lty=3, lwd=2, col="red")
abline(h=V+tau, lty=3, lwd=2, col="red")

plot(u, CorExp(u,phi), type="l", ylab=expression(rho), 
     bty="n", lwd=2, xlab="distância", main="Fção de correlação exp.")
abline(v=3*phi, lty=3, lwd=2, col="red")
abline(h=0.05, lty=3, lwd=2, col="red")

plot(u, CovExp(u, phi, V, tau), type="l", ylab="C(distância)", 
     bty="n", lwd=2, xlab="distância", main="Covariograma")
abline(v=3*phi, lty=3, lwd=2, col="red")
abline(h=0.05, lty=3, lwd=2, col="red")

plot(u, tau+V*(1-CorExp(u,phi)), type="l", ylab="gamma(distância)", 
     bty="n", lwd=2, xlab="distância", main="Semivariograma")
abline(v=3*phi, lty=3, lwd=2, col="red")
abline(h=V+tau, lty=3, lwd=2, col="red")

1.4 Comparando as funções de correlação

Compare as funções de correlação esférica, exponencial e Matérn.

#Comparando as diferentes funcoes de correlacao

par(mfrow=c(1,1))
u         = seq(0, max(d), l=1000)
phis      = c(0.02, 0.15,0.3,0.35)
#lambda2  = 0.5
lambda    = 0.8
j         = 1
plot(u, CorEsf(u, phis[j]), type="l", ylab=expression(rho), 
     bty="n", lwd=2, xlab="distância", main="Funções de correlação")
#abline(v=phis[j], lty=3, lwd=2, col=j)
for(j in 2:length(phis))
{
    lines(u, CorEsf(u, phis[j]), col=j, lwd=2)     
    #abline(v=phis[j], lty=3, lwd=2, col=j)

    lines(u, CorExp(u,phis[j]), col=j, lwd=2, lty=3)
    #abline(v=3*phis[j], lty=3, lwd=2, col=j)

    #  y <- cov.spatial(u, cov.model= "matern",cov.pars=c(V, phis[j]), kappa = lambda2)
    #  lines(u, y, col=j, lwd=2, lty=5, lwd=3)
    y <- cov.spatial(u, cov.model= "matern",cov.pars=c(V, phis[j]), kappa = lambda)
    lines(u, y, col=j, lwd=2, lty=5)
}
abline(h=0, lty=4, lwd=2)
abline(h=0.05, lty=4, lwd=2)
legend(x=0.5, y=0.9, legend=phis, lty=rep(1,length(phis)), col=1:4, lwd=rep(2,length(phis)), bty="n")
legend(x=0.9, y=0.9, legend=c("Esférica", "Exponencial", "Matern"), lty=c(1,3,5), lwd=rep(2,length(phis)), bty="n")

1.5 Gerando os dados

  1. Gere uma amostra de tamanho \(M=1000\) do vetor \(\boldsymbol{Z} = (Z_1, \ldots, Z_n)^T \sim N_n(\boldsymbol{0}, \boldsymbol{\Sigma})\), sendo \(\boldsymbol{0}\) um vetor de tamanho \(n=500\) com todos os elementos iguais a zero e \(\boldsymbol{\Sigma}\) uma matriz de covariância. Considere a função de covariância exponencial com \(\phi=0,25\), \(V=0,5\) e \(\tau=0\).

  2. Faça um histograma de \(Z_1\). Plote em cima desse histograma a curva da densidade de \(Z_1\).

  3. Faça um histograma de \(Z_{100}\). Plote em cima desse histograma a curva da densidade de \(Z_{100}\).

  4. Faça um histograma com todos os pontos amostrados de \(\boldsymbol{Z}\).

  5. Em seguida, calcule a correlação amostral de \(Z_i\) e \(Z_j\), para \(i=1, \ldots, n\) e \(j=1, \ldots, n\).

  6. Compare a correlação amostral com a correlação teórica e a variância amostral com a teórica.

#Gerando os dados
library(MASS)
M           = 1000
V           = 0.5
tau         = 0.03
Sigma       = CovExp(d, phi, V, tau)
amostra_Z   = mvrnorm(M, rep(0,n), Sigma)

par(mfrow=c(2,2))
#Lembre que Z_i ~ N(0, V+tau).
hist(amostra_Z[,1], freq=F, ylab="densidade", xlab="erros espaciais", main="Histograma de Z_1")
curve(dnorm(x, 0, sqrt(V+tau)), add=T, lwd=2, col="red")

hist(amostra_Z[,100], freq=F, ylab="densidade", xlab="erros espaciais", main="Histograma de Z_100")
curve(dnorm(x, 0, sqrt(V+tau)), add=T, lwd=2, col="red")

hist(amostra_Z, freq=F, ylab="densidade", xlab="erros espaciais", main="Histograma")

correlacao1 = correlacao2 = dif_cor = matrix(NA, n, n)
variancia = NULL
for(ind1 in 1:n){
  variancia[ind1] = var(amostra_Z[,ind1])
  for(ind2 in 1:n){
    correlacao1[ind1, ind2]  = cor(amostra_Z[,ind1], amostra_Z[,ind2])
    correlacao2[ind1, ind2] = CorExp(d[ind1, ind2], phi)
    dif_cor[ind1, ind2]     = correlacao1[ind1, ind2] - correlacao2[ind1, ind2]
  }
}

par(mfrow=c(2,2))

plot(correlacao1, correlacao2, xlab="correlação amostral", ylab="correlação teórica")
plot(as.vector(dif_cor[row(dif_cor)<col(dif_cor)]), ylab="dif. cor. amostral versus teórica")
abline(h=0, lwd=2, lty=3, col="red")
plot(variancia, ylab="variância amostral")
abline(h=V+tau, col="red", lwd=2, lty=3)
min(dif_cor)
## [1] -0.1275514
max(dif_cor)
## [1] 0.0926288

Suponha que haja uma única amostra para cada localização. Calcule o semivariograma empírico.

#http://www.leg.ufpr.br/geoR/geoRdoc/vignette/geoRintro.pdf

#Gerando os dados
library(gstat)
library(sp)

Z           = mvrnorm(1, rep(0,n), Sigma)
dados       = data.frame(lat, lon, Z)
dados       = as.geodata(dados)
summary(dados)
## Number of data points: 500 
## 
## Coordinates summary
##            lat         lon
## min 0.00369455 0.001442125
## max 0.99628446 0.999228003
## 
## Distance summary
##          min          max 
## 0.0006532249 1.3159330488 
## 
## Data summary
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.5312804 -0.1444343  0.2383864  0.2188596  0.5651669  1.8643142
plot(dados)

par(mfrow = c(2, 2))
points(dados, xlab = "Coord X", ylab = "Coord Y")
points(dados, xlab = "Coord X", ylab = "Coord Y",
pt.divide = "rank.prop")
points(dados, xlab = "Coord X", ylab = "Coord Y",
cex.max = 1.7, col = gray(seq(1, 0.1, l = 100)),
pt.divide = "equal")
points(dados, pt.divide = "quintile", xlab = "Coord X",
ylab = "Coord Y")

cloud1 <- variog(dados, option = "cloud", max.dist = max(d))
## variog: computing omnidirectional variogram
cloud2 <- variog(dados, option = "cloud", estimator.type = "modulus",
 max.dist = max(d))
## variog: computing omnidirectional variogram
bin1 <- variog(dados, uvec = seq(0, max(d), l = 50))
## variog: computing omnidirectional variogram
bin2 <- variog(dados, uvec = seq(0, max(d), l = 50),
estimator.type = "modulus")
## variog: computing omnidirectional variogram
par(mfrow = c(2, 2))
plot(cloud1, main = "classical estimator")
plot(cloud2, main = "modulus estimator")
plot(bin1, main = "classical estimator")
plot(bin2, main = "modulus estimator")

#Analisando
par(mfrow=c(1,1))
bin1 <- variog(dados, uvec = seq(0, max(d), l = 50))
## variog: computing omnidirectional variogram
plot(bin1, lwd=2, pch=19, bty="n", ylim=c(0, V+tau))
lines.variomodel(cov.model = "exp", cov.pars = c(V, phi), nugget = tau, 
                 max.dist = max(d), lwd = 2)
smooth <- variog(dados, option = "smooth", max.dist = max(d), 
                 n.points = n, kernel = "normal", band = 0.2)
## variog: computing omnidirectional variogram
lines(smooth, type = "l", lty = 2, lwd=2)
legend(0.4, 0.3, c("empirical", "exponential model",
 "smoothed"), lty = c(NA, 1, 2), pch=c(19,NA,NA), lwd = c(1,3, 1))

#Estimando os parametros
ml <- likfit(dados, ini = c(5, 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.
ml
## likfit: estimated model parameters:
##      beta     tausq   sigmasq       phi 
## "-0.0646" " 0.0361" " 0.3336" " 0.3681" 
## Practical Range with cor=0.05 for asymptotic range: 1.102587
## 
## likfit: maximised log-likelihood = -97.29
summary(ml)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##    beta 
## -0.0646 
## 
## Parameters of the spatial component:
##    correlation function: exponential
##       (estimated) variance parameter sigmasq (partial sill) =  0.3336
##       (estimated) cor. fct. parameter phi (range parameter)  =  0.3681
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (estimated) nugget =  0.0361
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 1 (no transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 1.102587
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
## "-97.29"      "4"  "202.6"  "219.4" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
## "-423.7"      "2"  "851.5"  "859.9" 
## 
## Call:
## likfit(geodata = dados, ini.cov.pars = c(5, 0.5))
ml2 <- likfit(dados, ini = c(5, 0.5), fix.nugget = T)
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optimize.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optimize.
## 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.
ml2
## likfit: estimated model parameters:
##     beta  sigmasq      phi 
## "0.1156" "0.2920" "0.0963" 
## Practical Range with cor=0.05 for asymptotic range: 0.2884802
## 
## likfit: maximised log-likelihood = -140.8
summary(ml2)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##   beta 
## 0.1156 
## 
## Parameters of the spatial component:
##    correlation function: exponential
##       (estimated) variance parameter sigmasq (partial sill) =  0.292
##       (estimated) cor. fct. parameter phi (range parameter)  =  0.0963
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (fixed) nugget = 0
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 1 (no transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 0.2884802
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
## "-140.8"      "3"  "287.7"  "300.3" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
## "-423.7"      "2"  "851.5"  "859.9" 
## 
## Call:
## likfit(geodata = dados, ini.cov.pars = c(5, 0.5), fix.nugget = T)
ajuste <- lm(Z~rep(1,n)-1)
summary(ajuste)
## 
## Call:
## lm(formula = Z ~ rep(1, n) - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.75014 -0.36329  0.01953  0.34631  1.64545 
## 
## Coefficients:
##           Estimate Std. Error t value Pr(>|t|)    
## rep(1, n)  0.21886    0.02528   8.658   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5653 on 499 degrees of freedom
## Multiple R-squared:  0.1306, Adjusted R-squared:  0.1289 
## F-statistic: 74.95 on 1 and 499 DF,  p-value: < 2.2e-16