Trabalho 4

Estimação em um modelo geoestatístico gaussiano

Aluno:Fernando G. Moro

Disciplina: CE089-Estatística Computacional II

Descrição e Objetivos com o Trabalho

O trabalho consiste em cinco etapas:

  1. Apresentar a teoria básica do modelo proposto

  2. Simular dados do modelo

  3. Construir a função de log-verossimilhança e implementar o algorítimo de estimação (Newton-Rhapson)

  4. Ver o comportamento da função de verossimilhança

  5. Aplicar a um conjunto de dados reais

O trabalho tem por objetivo exemplificar a implementação de um processo de estimação de um modelo estatístico.

Estrutura do modelo geoestatístico gaussiano

Seja \(\textbf{Z}\) o vetor de observações da variável resposta, \(\textbf{X}\) a matrix das covariáveis e \(\gamma\) o efeito aleatório das distâncias das coordenadas \((x,y)\), então temos:

\[ \begin{align} \textbf{Z} = \textbf{X}\beta + \gamma \end{align} \]

Onde a distribuição de \(\textbf{Z}\) é dada por:

\[ \begin{align} \textbf{Z} \sim NM(\textbf{X} \beta,\tau^2 \textbf{I}) \end{align} \]

Note que a distribuição do vetor de respostas \(\textbf{Z}\) considera-se iid.

Já o vetor \(\gamma\) contém os efeitos aleatórios associados a cada observação da minha amostra, ou seja:

\[ \begin{align} \gamma \sim NM(\textbf{0},\sigma^2 \Sigma_{\gamma}) \end{align} \]

Sendo \(\sigma^2\) a variância associada a minha variável latente e \(\Sigma_{\gamma}\) a matriz de correlações associada às distâncias entre as coordenadas de cada ponto de observação no espaço. No caso, a distância utilizada é a euclidiana, representada por \(\textbf{h}=||x_i - y_j||\). Porém, para podermos definir uma correlação temos várias estruturas definidas na literatura, para este trabalho iremos considerar a estrutura de correlação exponencial, dada por:

\[ \begin{align} \rho (\textbf{h},\phi)= exp(\textbf{h}/\phi) \end{align} \]

Onde o parâmetro \(\phi\) controla a intensidade de decaimento da correlação de acordo com a medida de distância \(\textbf{h}\).

Podemos resumir o modelo em:

\[ \begin{align} \textbf{Z} \sim NM(\textbf{X}\beta,\Sigma) \end{align} \]

Onde \(\Sigma\) fica dada por:

\[\Sigma = \tau^2 \textbf{I} + \sigma^2 \Sigma_{\gamma}\]

Simulando dados do modelo proposto

Implementação da matriz \(\Sigma\):

Matr<-function(sigma, t, phi, Dist){
    MSig<-as.matrix(sigma^2*exp(-Dist/phi))
    diag(MSig)<- sigma^2 + t^2
    return(MSig)
}

Implementando Função para simular dados do modelo:

random.geo<-function(bheta,sigma,t,phi,n){
    cords<-data.frame(cx=runif(n),cy=runif(n))
    Cdist<-dist(cords,diag=T,upper=T)
    MSig<-Matr(sigma=sigma,t=t,phi=phi,Dist=Cdist)
    Z<-rnorm(n)
    Y<-bheta+crossprod(chol(MSig),Z)
    return(cbind(cords,Y=Y))
}

Simulando dados do modelo:

sim<-random.geo(bheta=0,sigma=6,t=2,phi=4,n=125)
head(sim)
##         cx     cy       Y
## 1 0.005668 0.7681  0.3886
## 2 0.039049 0.2816 -7.8218
## 3 0.875802 0.1629  0.6647
## 4 0.135953 0.5762 -3.0131
## 5 0.803781 0.6524 -1.9263
## 6 0.923122 0.1453 -3.2844

Implementação e estimação de parâmetros

Utilizando regressão linear como exemplo “didático”

reglin<- lm(Y ~ cx + cy,data=sim)
summary(reglin)
## 
## Call:
## lm(formula = Y ~ cx + cy, data = sim)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.364 -1.544 -0.307  1.338  5.956 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -7.862      0.519  -15.15  < 2e-16 ***
## cx             3.463      0.640    5.41  3.1e-07 ***
## cy             3.774      0.728    5.18  8.7e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.24 on 122 degrees of freedom
## Multiple R-squared:  0.339,  Adjusted R-squared:  0.328 
## F-statistic: 31.3 on 2 and 122 DF,  p-value: 1.05e-11

Questão: Quais são os possíveis problemas de ajustar um modelo linear?

  1. Não considerar a dependência ou correlação entre as observações no espaço.

  2. As coordenadas observadas representam uma superfície de possíveis pontos no espaço, de forma que não podem ser considerados como efeitos fixos.

  3. Em consequência das demais teremos inferências errôneas.

Escrevendo a função de Log-verossimilhança:

\[ \begin{align} ll(\sigma,\tau,\phi)=\frac{-n}{2log(2 \pi)}-\frac{1}{2}log|\Sigma|-\frac{1}{2}(\textbf{Z}-\textbf{X}\beta)^T \Sigma^{-1} (\textbf{Z}-\textbf{X}\beta) \end{align} \]

Implementando a Log-verossimilhança

ll.geo <- function(theta, dados){
s<-exp(theta[1])
t<-exp(theta[2])
phi<-exp(theta[3])
U <- dist(dados[,1:2], diag=TRUE, upper=TRUE)
D<-as.vector(rep(1,times=nrow(dados)))
Sigma<-Matr(sigma=s,t=t,phi=phi,Dist=U)
bheta<- (t(D)%*%solve(Sigma)%*%D)%*%(t(D)%*%solve(Sigma)%*%dados$Y)
ll = -length(dados$Y)/2*log(2*pi)-1/2*log(det(Sigma))-1/2*t(dados$Y-D%*%bheta)%*%solve(Sigma)%*%(dados$Y-D%*%bheta)
return(-ll)
}

Equação do algorítimo Newton-Rhapson:

\[ \begin{align} \textbf{X} = \textbf{X} - \textbf{H}^{-1} \textbf{G} \end{align} \]

Onde \(\textbf{H}\) representa a matriz de derivadas segundas ou matriz hessiana e \(\textbf{G}\) representa o vetor gradiente.

Implementando o algorítimo temos:

#Chamando pacote para as derivadas numéricas
require(rootSolve)
## Loading required package: rootSolve
require(numDeriv)
## Loading required package: numDeriv
## 
## Attaching package: 'numDeriv'
## 
## The following object is masked from 'package:rootSolve':
## 
##     hessian
maxiter <- 500; i <- 1         
tol <- 1e-4; error <- 100*tol  

theta <- matrix(NA, nrow=1, ncol=3)
theta[1,] <- c(1,1,1)


while(i <= maxiter & error>tol){
  theta <- rbind(theta, rep(NA,3))
  theta[i+1,] <- theta[i,]-
    solve(hessian(ll.geo,x=as.numeric(theta[i,]), dados=sim))%*%t(gradient(ll.geo,x=as.numeric(theta[i,]),dados=sim))
  error <- sum(abs((theta[i+1,]-theta[i,])/theta[i+1,]))
  i <- i+1
  print(cbind(sigma=exp(theta[i,1]),tau=exp(theta[i,2]),phi=exp(theta[i,3]),error, iter=i))
}
##      sigma   tau  phi error iter
## [1,] 3.795 1.871 2.27 1.067    2
##      sigma   tau   phi error iter
## [1,]  4.09 1.794 1.275 2.504    3
##      sigma   tau   phi  error iter
## [1,] 4.027 1.816 1.459 0.3892    4
##      sigma  tau   phi   error iter
## [1,] 4.017 1.82 1.479 0.04073    5
##      sigma  tau   phi     error iter
## [1,] 4.017 1.82 1.479 0.0002611    6
##      sigma  tau   phi   error iter
## [1,] 4.017 1.82 1.479 1.6e-06    7

Calculando a estimativa de \(\beta\): Estimador de máxima verossimilhança:

\[ \begin{align} \beta= (\textbf{X}^{T}\Sigma^{-1}\textbf{X})(\textbf{X}^{T}\Sigma^{-1}\textbf{Z}) \end{align} \]

Implementando temos:

U <- dist(sim[,1:2], diag=TRUE, upper=TRUE)
D<-as.vector(rep(1,times=nrow(sim)))
Sigma<-Matr(sigma=exp(theta[max(i),1]),t=exp(theta[max(i),2]),phi=exp(theta[max(i),3]),Dist=U)
bheta<- (t(D)%*%solve(Sigma)%*%D)%*%(t(D)%*%solve(Sigma)%*%sim$Y)
bheta
##          [,1]
## [1,] -0.02764

Gráficos da superfície de verossimilhança

  1. Superfície de verossimilhança entre \(\sigma\) e \(\tau\)
ll.V<-Vectorize(ll.geo2 <- function(s,t,phi, dados){
  U <- dist(dados[,1:2], diag=TRUE, upper=TRUE)
  D<-as.vector(rep(1,times=nrow(dados)))
  Sigma<-Matr(sigma=s,t=t,phi=phi,Dist=U)
  bheta<- (t(D)%*%solve(Sigma)%*%D)%*%(t(D)%*%solve(Sigma)%*%dados$Y)
  ll = -length(dados$Y)/2*log(2*pi)-1/2*log(det(Sigma))-1/2*t(dados$Y-D%*%bheta)%*%solve(Sigma)%*%(dados$Y-D%*%bheta)
  return(ll)
},c("s","t"))

Considerando sequência de valores para construir a superfície

s<-seq(exp(theta[max(i),1])-0.5*exp(theta[max(i),1]),exp(theta[max(i),1])+0.5*exp(theta[max(i),1]),len=70)
t<-seq(exp(theta[max(i),2])-0.5*exp(theta[max(i),2]),exp(theta[max(i),2])+0.5*exp(theta[max(i),2]),len=70)
phi<-exp(theta[max(i),3])

Construindo a superfície:

llST<-outer(s,t,ll.V,dados=sim,phi=phi)
m<-max(llST)
contour(s, t, llST, levels=seq(m-200, m-1, by=1),xlab=expression(sigma),ylab=expression(tau),lty=2)

plot of chunk unnamed-chunk-9

  1. Superfície de verossimilhança entre \(\sigma\) e \(\phi\)
ll.V2<-Vectorize(ll.geo2 <- function(s,t,phi, dados){
  U <- dist(dados[,1:2], diag=TRUE, upper=TRUE)
  D<-as.vector(rep(1,times=nrow(dados)))
  Sigma<-Matr(sigma=s,t=t,phi=phi,Dist=U)
  bheta<- (t(D)%*%solve(Sigma)%*%D)%*%(t(D)%*%solve(Sigma)%*%dados$Y)
  ll = -length(dados$Y)/2*log(2*pi)-1/2*log(det(Sigma))-1/2*t(dados$Y-D%*%bheta)%*%solve(Sigma)%*%(dados$Y-D%*%bheta)
  return(ll)
},c("s","phi"))

Considerando sequência de valores para construir a superfície

s<-seq(exp(theta[max(i),1])-0.5*exp(theta[max(i),1]),exp(theta[max(i),1])+0.5*exp(theta[max(i),1]),len=70)
phi<-seq(exp(theta[max(i),3])-0.5*exp(theta[max(i),3]),exp(theta[max(i),3])+0.5*exp(theta[max(i),3]),len=70)
t<-exp(theta[max(i),2])

Construindo a superfície:

llST2<-outer(s,phi,ll.V2,dados=sim,t=t)
m<-max(llST2)
contour(s, phi, llST2, levels=seq(m-50, m-0.1, by=0.1),xlab=expression(sigma),ylab=expression(phi),col=2,lty=6)

plot of chunk unnamed-chunk-12

  1. Superfície de verossimilhança entre \(\tau\) e \(\phi\)
ll.V3<-Vectorize(ll.geo2 <- function(s,t,phi, dados){
  U <- dist(dados[,1:2], diag=TRUE, upper=TRUE)
  D<-as.vector(rep(1,times=nrow(dados)))
  Sigma<-Matr(sigma=s,t=t,phi=phi,Dist=U)
  bheta<- (t(D)%*%solve(Sigma)%*%D)%*%(t(D)%*%solve(Sigma)%*%dados$Y)
  ll = -length(dados$Y)/2*log(2*pi)-1/2*log(det(Sigma))-1/2*t(dados$Y-D%*%bheta)%*%solve(Sigma)%*%(dados$Y-D%*%bheta)
  return(ll)
},c("t","phi"))

Considerando sequência de valores para construir a superfície

t<-seq(exp(theta[max(i),2])-0.5*exp(theta[max(i),2]),exp(theta[max(i),2])+0.5*exp(theta[max(i),2]),len=70)
phi<-seq(exp(theta[max(i),3])-0.5*exp(theta[max(i),3]),exp(theta[max(i),3])+0.5*exp(theta[max(i),3]),len=70)
s<-exp(theta[max(i),1])

Construindo a superfície:

llST3<-outer(t,phi,ll.V3,dados=sim,s=s)
m<-max(llST3)
contour(t, phi, llST3, levels=seq(m-20, m-1, by=0.5),xlab=expression(tau),ylab=expression(phi),col="blue",lty=2)

plot of chunk unnamed-chunk-15

Aplicando a um conjunto de dados reais

Importando conjunto de dados:

require(geoR)
## Loading required package: geoR
## Loading required package: sp
## Loading required package: MASS
## --------------------------------------------------------------
##  Analysis of geostatistical data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.7-4.1 (built on 2012-06-29) is now loaded
## --------------------------------------------------------------
data(parana)


hist(parana$data,xlab="Acúmulo",ylab="Densidade",prob=T,main="Média de chuva",col=4)

plot of chunk unnamed-chunk-16

Selecionando as coordenadas:

dados<-data.frame(cx=parana$coords[,1], cy=parana$coords[,2],Y=((parana$data-mean(parana$data))/sd(parana$data)))
head(dados)
##      cx    cy       Y
## 1 403.0 164.5  0.5521
## 2 501.7 428.8 -1.2815
## 3 556.3 445.3 -1.8707
## 4 573.4 447.0 -1.9457
## 5 702.4 272.3 -1.9317
## 6 668.5 261.7 -1.6696

Estimação dos parâmetros:

require(rootSolve)
require(numDeriv)
maxiterd <- 500; id <- 1         
told<- 1e-4; errord <- 100*told  
thetad <- matrix(NA, nrow=1, ncol=3)
thetad[1,] <- c(0.1,-1,6)

while(id <= maxiterd & errord>told){
  thetad <- rbind(thetad, rep(NA,3))
  thetad[id+1,] <- thetad[id,]-
    solve(hessian(ll.geo,x=as.numeric(thetad[id,]), dados=dados))%*%t(gradient(ll.geo,x=as.numeric(thetad[id,]),dados=dados))
  errord <- sum(abs((thetad[id+1,]-thetad[id,])/thetad[id+1,]))
  id <- id+1
  print(cbind(sigmapad=exp(thetad[id,1]),taupad=exp(thetad[id,2]),phiest=exp(thetad[id,3]),errord, iterd=id))
}
##      sigmapad taupad phiest errord iterd
## [1,]    1.073 0.3152    625 0.6188     2
##      sigmapad taupad phiest errord iterd
## [1,]    1.145 0.3293  726.9 0.5405     3
##      sigmapad taupad phiest errord iterd
## [1,]    1.238 0.3294  854.8 0.3914     4
##      sigmapad taupad phiest errord iterd
## [1,]    1.376 0.3296   1066  0.363     5
##      sigmapad taupad phiest errord iterd
## [1,]    1.592 0.3299   1439 0.3552     6
##      sigmapad taupad phiest errord iterd
## [1,]    1.703   0.33   1649  0.145     7
##      sigmapad taupad phiest   errord iterd
## [1,]    1.709   0.33   1661 0.007694     8
##      sigmapad taupad phiest    errord iterd
## [1,]    1.709   0.33   1661 5.321e-05     9

Estimando \(\beta\):

U <- dist(dados[,1:2], diag=TRUE, upper=TRUE)
D<-as.vector(rep(1,times=nrow(dados)))
Sigma<-Matr(sigma=exp(thetad[max(id),1]),t=exp(thetad[max(id),2]),phi=exp(thetad[max(id),3]),Dist=U)
bheta<- (t(D)%*%solve(Sigma)%*%D)%*%(t(D)%*%solve(Sigma)%*%dados$Y)
bheta
##          [,1]
## [1,] -0.09653

Gráficos de superfície de Verossimilhança

  1. Superfície de verossimilhança entre \(\sigma\) e \(\tau\)

Considerando sequência de valores para construir a superfície

s<-seq(exp(thetad[max(id),1])-0.5*exp(thetad[max(id),1]),exp(thetad[max(id),1])+0.5*exp(thetad[max(id),1]),len=70)
t<-seq(exp(thetad[max(id),2])-0.5*exp(thetad[max(id),2]),exp(thetad[max(id),2])+0.5*exp(thetad[max(id),2]),len=70)
phi<-exp(thetad[max(id),3])

Construindo a superfície:

llST<-outer(s,t,ll.V,dados=dados,phi=phi)
m<-max(llST)
contour(s, t, llST, levels=seq(m-80, m-0.5, by=0.5,xlab=expression(sigma),ylab=expression(tau),lty=2))
## Warning: extra arguments 'xlab', 'ylab', 'lty' will be disregarded

plot of chunk unnamed-chunk-21

  1. Superfície de verossimilhança entre \(\sigma\) e \(\phi\)

Considerando sequência de valores para construir a superfície

s<-seq(exp(thetad[max(id),1])-0.5*exp(thetad[max(id),1]),exp(thetad[max(id),1])+0.5*exp(thetad[max(id),1]),len=150)
phi<-seq(exp(thetad[max(id),3])-0.5*exp(thetad[max(id),3]),exp(thetad[max(id),3])+0.5*exp(thetad[max(id),3]),len=150)
t<-exp(thetad[max(id),2])

Construindo a superfície:

llST2<-outer(s,phi,ll.V2,dados=dados,t=t)
m<-max(llST2)
contour(s, phi, llST2, levels=seq(m-80, m, by=0.2),xlab=expression(sigma),ylab=expression(phi),col=2,lty=2)

plot of chunk unnamed-chunk-23

  1. Superfície de verossimilhança entre \(\tau\) e \(\phi\)

Considerando sequência de valores para construir a superfície

t<-seq(exp(thetad[max(id),2])-0.5*exp(thetad[max(id),2]),exp(thetad[max(id),2])+0.5*exp(thetad[max(id),2]),len=150)
phi<-seq(exp(thetad[max(id),3])-0.5*exp(thetad[max(id),3]),exp(thetad[max(id),3])+0.5*exp(thetad[max(id),3]),len=150)
s<-exp(thetad[max(id),1])

Construindo a superfície:

llST3<-outer(t,phi,ll.V3,dados=dados,s=s)
m<-max(llST3)
contour(t, phi, llST3, levels=seq(m-80, m, by=0.2),xlab=expression(tau),ylab=expression(phi),col="blue",lty=2)

plot of chunk unnamed-chunk-25

O pacote geoR possui já funções implementadas para este tipo de modelagem. Iremos utilizar o mesmo a fim de comparação com o modelo utilizado.

Modelando com a função likfit:

datageo<-as.geodata(dados)
ml <- likfit(datageo, trend="cte", ini=c(1,1),lambda=1)
## ---------------------------------------------------------------
## 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.

Saída do modelo:

summary(ml)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood 
## 
## Parameters of the mean component (trend):
##  beta 
## -0.54 
## 
## Parameters of the spatial component:
##    correlation function: exponential
##       (estimated) variance parameter sigmasq (partial sill) =  2.91
##       (estimated) cor. fct. parameter phi (range parameter)  =  1659
##    anisotropy parameters:
##       (fixed) anisotropy angle = 0  ( 0 degrees )
##       (fixed) anisotropy ratio = 1
## 
## Parameter of the error component:
##       (estimated) nugget =  0.109
## 
## Transformation parameter:
##       (fixed) Box-Cox parameter = 1 (no transformation)
## 
## Practical Range with cor=0.05 for asymptotic range: 4969
## 
## Maximised Likelihood:
##    log.L n.params      AIC      BIC 
##  "-92.5"      "4"    "193"    "205" 
## 
## non spatial model:
##    log.L n.params      AIC      BIC 
##   "-202"      "2"    "409"    "415" 
## 
## Call:
## likfit(geodata = datageo, trend = "cte", ini.cov.pars = c(1, 
##     1), lambda = 1)

Parâmetros estimados:

cbind(tau=sqrt(as.numeric(ml[2])),sigma=sqrt(as.numeric(ml$cov.pars[1])),phi=as.numeric(ml$cov.pars[2]),beta=as.numeric(ml$beta))
##       tau sigma  phi    beta
## [1,] 0.33 1.707 1659 -0.5403

Valores preditos e gráficos de superfície:

gx <- seq(50,800, by=10)
gy <- seq(-100,650, by=10)
loc0 <- expand.grid(gx,gy)

KC <- krige.control(trend.d="1st", trend.l="1st", obj=ml)
kc <- krige.conv(datageo, loc=loc0, krige= KC, bor=parana$borders)
## 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

Gráfico de superfície:

image(kc, col=terrain.colors(12))
contour(kc,add=T,lty=2,col=4)

plot of chunk unnamed-chunk-30