
#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$:

```{r,echo=TRUE}
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:

```{r,echo=TRUE}
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:

```{r,echo=TRUE}
sim<-random.geo(bheta=0,sigma=6,t=2,phi=4,n=125)
head(sim)
```

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

Utilizando regressão linear como exemplo "didático"

```{r,echo=TRUE}
reglin<- lm(Y ~ cx + cy,data=sim)
summary(reglin)
```

**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.

1. 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.

1. 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**

```{r,echo=TRUE}
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:

```{r,est1,echo=TRUE}
#Chamando pacote para as derivadas numéricas
require(rootSolve)
require(numDeriv)
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))
}
```

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:

```{r,echo=TRUE}
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
```


###Gráficos da superfície de verossimilhança###


a. Superfície de verossimilhança entre $\sigma$ e $\tau$

```{r,echo=TRUE}
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

```{r,echo=TRUE}
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:

```{r,echo=TRUE}
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)
```

b. Superfície de verossimilhança entre $\sigma$ e $\phi$

```{r,echo=TRUE}
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

```{r,echo=TRUE}
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:

```{r,echo=TRUE}
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)
```


c. Superfície de verossimilhança entre $\tau$ e $\phi$

```{r,echo=TRUE}
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

```{r,echo=TRUE}
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:

```{r,echo=TRUE}
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)
```

###Aplicando a um conjunto de dados reais###


Importando conjunto de dados:

```{r,echo=TRUE}
require(geoR)
data(parana)


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

Selecionando as coordenadas:

```{r,echo=TRUE}
dados<-data.frame(cx=parana$coords[,1], cy=parana$coords[,2],Y=((parana$data-mean(parana$data))/sd(parana$data)))
head(dados)
```


Estimação dos parâmetros:

```{r,echo=TRUE}
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))
}
```

Estimando $\beta$:


```{r,echo=TRUE}
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
```


Gráficos de superfície de **Verossimilhança**


a. Superfície de verossimilhança entre $\sigma$ e $\tau$

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

```{r,echo=TRUE}
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:

```{r,echo=TRUE}
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))
```

b. Superfície de verossimilhança entre $\sigma$ e $\phi$

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

```{r,echo=TRUE}
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:

```{r,echo=TRUE}
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)
```


c. Superfície de verossimilhança entre $\tau$ e $\phi$


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

```{r,echo=TRUE}
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:

```{r,echo=TRUE}
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)
```

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**:

```{r,echo=TRUE}
datageo<-as.geodata(dados)
ml <- likfit(datageo, trend="cte", ini=c(1,1),lambda=1)
```

Saída do modelo:

```{r,echo=TRUE}
summary(ml)
```

Parâmetros estimados:

```{r,echo=TRUE}
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))
```


Valores preditos e gráficos de superfície:

```{r,echo=TRUE,warning=FALSE}
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)
```

Gráfico de superfície:

```{r,echo=TRUE,warning=FALSE}
image(kc, col=terrain.colors(12))
contour(kc,add=T,lty=2,col=4)
```








