Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
| — |
projetos:camada1 [2007/01/19 11:53] (atual) abtmartins criada |
||
|---|---|---|---|
| Linha 1: | Linha 1: | ||
| + | %<<echo=F>>= | ||
| + | %require(geoR) | ||
| + | %rm(list=ls(all=TRUE)) | ||
| + | %@ | ||
| + | |||
| + | \subsection{Análise da variável alumínio medida na camada~1} | ||
| + | |||
| + | % ******************************************************************** | ||
| + | % | ||
| + | \subsubsection{Análise exploratória dos dados} | ||
| + | % | ||
| + | % ******************************************************************** | ||
| + | |||
| + | \noindent Criando o objeto geodata para Alumínio. Estamos indicando que as coordenadas $(x,y)$ dos pontos amostrados encontram-se na primeira e segunda coluna do banco de dados (arquivo Quimicos.txt), a variável Alumínio encontra-se na terceira coluna e a covariável região encontra-se na coluna de número~12. | ||
| + | |||
| + | <<>>= | ||
| + | Alum<-as.geodata(Quimicos,coords.col=1:2,data.col=3,covar.col=12) | ||
| + | @ | ||
| + | |||
| + | \noindent O comando a seguir atribui à variável {\sf Alum\$borders} as coordenadas que se encontram no arquivo {\sf borda.txt}. | ||
| + | |||
| + | <<>>= | ||
| + | Alum$borders <- borda | ||
| + | @ | ||
| + | %$ | ||
| + | |||
| + | \noindent Descrevendo os dados rapidamente | ||
| + | <<>>= | ||
| + | summary(Alum) | ||
| + | @ | ||
| + | |||
| + | \noindent Visualizando graficamente os dados originais, ou seja, sem qualquer transformação ou arranjo. | ||
| + | <<fig=T>>= | ||
| + | plot(Alum) | ||
| + | @ | ||
| + | |||
| + | \noindent Pelo histograma podemos observar que os dados não seguem uma distribuição normal, sendo necessária uma transformação nos dados. Para efetuarmos essa transformação, precisaremos identificar o valor de $\lambda$ e aplicar o método de Box-Cox. A figura \ref{fig:boxcoxAl} dá uma indicação gráfica de quais valores de $\lambda$ pertencem ao intervalo de 95\% de confiança. No comando seguinte, por conveniência visual, o gráfico foi construído em um intervalo de -1 a 0.5. O resultado foi armazenado na variável {\sf BC1} (sem o efeito da região) e {\sf BC2} (com o efeito da região). | ||
| + | |||
| + | <<fig=F>>= | ||
| + | par(mfrow=c(1,2)) | ||
| + | boxcox(Alum,lam = seq(-1,0.5, l=100))-> BC1 | ||
| + | boxcox(Alum, trend=~regiao, lam = seq(-1,0.5, l=100))-> BC2 | ||
| + | par(mfrow=c(1,1)) | ||
| + | @ | ||
| + | |||
| + | \begin{figure}[!ht] | ||
| + | \SweaveOpts{width=7,height=4} | ||
| + | \centering | ||
| + | <<fig=T, echo=F>>= | ||
| + | par(mfrow=c(1,2)) | ||
| + | boxcox(Alum,lam = seq(-1,0.5, l=100))-> BC1 | ||
| + | boxcox(Alum, trend=~regiao, lam = seq(-1,0.5, l=100))-> BC2 | ||
| + | par(mfrow=c(1,1)) | ||
| + | @ | ||
| + | \vspace{-0.5cm} | ||
| + | \caption{Log-verossimilhança do parâmetro $\lambda$ para a transformação de Box-Cox (esquerda) sem efeito da região e (direita) com efeito da região} | ||
| + | \label{fig:boxcoxAl} | ||
| + | \end{figure} | ||
| + | |||
| + | \noindent Considerando um nível de 99\% de confiança decidimos adotar $\lambda=0$ devido a proximidade dos valores originais a zero, usando a log-normal e evitando assim problemas no truncamento. | ||
| + | |||
| + | \noindent Revendo, graficamente, os dados incluindo um efeito da região temos: | ||
| + | |||
| + | \begin{figure}[!ht] | ||
| + | \SweaveOpts{width=5.5,height=5.5} | ||
| + | <<fig=T>>= | ||
| + | plot(Alum, trend=~regiao) | ||
| + | @ | ||
| + | \end{figure} | ||
| + | |||
| + | % | ||
| + | % ******************************************************************** | ||
| + | % | ||
| + | \subsubsection{Análise Espacial} | ||
| + | % | ||
| + | % ******************************************************************** | ||
| + | \noindent Inspecionando os dados espacialmente (Figura \ref{fig:postplotAl} ). Neste ponto, consideramos o valor de $\lambda=0$ e estudaremos o efeito da região no modelo. | ||
| + | |||
| + | <<fig=F>>= | ||
| + | par(mar=c(3.5,3,0,0.5),mgp=c(2,0.8,0)) | ||
| + | points(Alum, lambda=0,cex.min=.3,cex.max=1, axes=F, ann=F,col=gray(seq(0.9,0,l=length(Alum$data)))) | ||
| + | polygon(A1) | ||
| + | polygon(A2) | ||
| + | par(mar=c(3.5,3.5,0.5,0.5),mgp=c(2,0.8,0)) | ||
| + | @ | ||
| + | %$ | ||
| + | |||
| + | <<fig=F>>= | ||
| + | par(mar=c(3.5,3,0,0.5),mgp=c(2,0.8,0)) | ||
| + | points(Alum, trend=~regiao, lambda=0,cex.min=.3,cex.max=1, axes=F, ann=F,col=gray(seq(0.9,0,l=length(Alum$data)))) | ||
| + | polygon(A1) | ||
| + | polygon(A2) | ||
| + | par(mar=c(3.5,3.5,0.5,0.5),mgp=c(2,0.8,0)) | ||
| + | @ | ||
| + | %$ | ||
| + | |||
| + | \begin{figure}[!ht] | ||
| + | \centering | ||
| + | \SweaveOpts{width=6,height=3.5} | ||
| + | \setkeys{Gin}{width=0.90\textwidth} | ||
| + | <<fig=T,echo=FALSE>>= | ||
| + | par(mfrow=c(1,2), mar=c(1,2,1,1),mgp=c(2,0.8,0)) | ||
| + | points(Alum, lambda=0,cex.min=.3,cex.max=1, axes=F, ann=F,col=gray(seq(0.9,0,l=length(Alum$data)))) | ||
| + | polygon(A1) | ||
| + | polygon(A2) | ||
| + | |||
| + | points(Alum, trend=~regiao, lambda=0,cex.min=.3,cex.max=1, axes=F, ann=F,col=gray(seq(0.9,0,l=length(Alum$data)))) | ||
| + | polygon(A1) | ||
| + | polygon(A2) | ||
| + | par(mfrow=c(1,2), mar=c(2,2,1,1),mgp=c(2,0.8,0)) | ||
| + | @ | ||
| + | %$ | ||
| + | \vspace{-1.5cm} | ||
| + | \caption{Locação dos dados na área de estudo e suas respectivas regiões, sem e com a covariável} | ||
| + | \label{fig:postplotAl} | ||
| + | \end{figure} | ||
| + | |||
| + | \noindent A próxima ação é a construção do semivariograma experimental a partir dos dados transformados, visando uma análise exploratória do comportamento da correlação espacial. Foi considerado também o efeito das diferentes regiões. A figura \ref{fig:variogramaAl} sugere uma dependência espacial, com valores iniciais de efeito pepita de 0,6, variância total de 1,4 e um alcance de 150 metros. | ||
| + | |||
| + | <<fig=F, results=hide>>= | ||
| + | Alum.v<-variog(Alum,max.dist=350,lambda=0, trend=~regiao) | ||
| + | plot(Alum.v, pch=19,cex=0.8) | ||
| + | @ | ||
| + | \begin{figure}[!ht] | ||
| + | \centering | ||
| + | \SweaveOpts{width=6,height=6} | ||
| + | \setkeys{Gin}{width=0.60\textwidth} | ||
| + | <<fig=T,echo=FALSE, results=hide>>= | ||
| + | Alum.v<-variog(Alum,max.dist=350,lambda=0, trend=~regiao) | ||
| + | plot(Alum.v, pch=19,cex=0.8) | ||
| + | @ | ||
| + | \vspace{-.5cm} | ||
| + | \caption{Semivariograma experimental com agrupamento de semivariâncias em intervalos de distâncias} | ||
| + | \label{fig:variogramaAl} | ||
| + | \end{figure} | ||
| + | |||
| + | % ************************************************************************** | ||
| + | % | ||
| + | \subsubsection{Ajustando um modelo teórico ao semivariograma experimental} | ||
| + | % | ||
| + | % ************************************************************************** | ||
| + | |||
| + | \noindent A idéia aqui é avaliar a aderência de uma função de correlação de Matérn ao semivariograma experimental produzido pelos dados transformados, segundo 5 valores diferentes do parâmetro de suavidade {\it kappa}. Os parâmetros aqui são estimados pelo método da máxima-verossimilhança: | ||
| + | |||
| + | \begin{itemize} | ||
| + | \item sem o efeito da covariável {\sf região}. | ||
| + | <<results=hide>>= | ||
| + | Al.ml <- list() | ||
| + | Al.ml$Al.1<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=0.5, lambda=0) | ||
| + | Al.ml$Al.2<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=1.5, lambda=0) | ||
| + | Al.ml$Al.3<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=2.5, lambda=0) | ||
| + | Al.ml$Al.4<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=3.5, lambda=0) | ||
| + | Al.ml$Al.5<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=4.5, lambda=0) | ||
| + | @ | ||
| + | %$ | ||
| + | |||
| + | \item com o efeito da covariável {\sf região}. | ||
| + | <<results=hide>>= | ||
| + | Al.ml1 <- list() | ||
| + | Al.ml1$Al.1<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=0.5, lambda=0, trend=~regiao) | ||
| + | Al.ml1$Al.2<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=1.5, lambda=0, trend=~regiao) | ||
| + | Al.ml1$Al.3<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=2.5, lambda=0, trend=~regiao) | ||
| + | Al.ml1$Al.4<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=3.5, lambda=0, trend=~regiao) | ||
| + | Al.ml1$Al.5<-likfit(Alum, ini=c(0.6,150),cov.model="matern", kappa=4.5, lambda=0, trend=~regiao) | ||
| + | @ | ||
| + | %$ | ||
| + | \end{itemize} | ||
| + | |||
| + | \noindent Assim, podemos apresentar os resultados obtidos da otimização da função log-verossimilhanca e escolher o modelo que apresentar o maior valor para a log-verossimilhança. | ||
| + | |||
| + | \begin{itemize} | ||
| + | \item[a)] Verossimilhanças sem o efeito da covariável | ||
| + | <<>>= | ||
| + | sapply(Al.ml, logLik) | ||
| + | @ | ||
| + | \item[b)] Verossimilhanças com o efeito da covariável | ||
| + | <<>>= | ||
| + | sapply(Al.ml1, logLik) | ||
| + | @ | ||
| + | \end{itemize} | ||
| + | |||
| + | \noindent O maior valor obtido foi -45.16316 que corresponde ao modelo com a presença da covariável {\sf região} e para um valor de {\it kappa} de 3,5. | ||
| + | |||
| + | \subsubsection{Construção dos mapas temáticos} | ||
| + | |||
| + | Definindo o {\it grid} para a interpolação de valores por krigagem | ||
| + | <<results=hide>>= | ||
| + | gr <- pred_grid(Alum$borders, by=10) | ||
| + | gr0<- polygrid(gr, borders=Alum$borders, bound=T) | ||
| + | ind.reg<-numeric(nrow(gr0)) | ||
| + | ind.reg[.geoR_inout(gr0,A1)]<-1 | ||
| + | ind.reg[.geoR_inout(gr0,A2)]<-2 | ||
| + | ind.reg[.geoR_inout(gr0,A3)]<-3 | ||
| + | ind.reg<-as.factor(ind.reg) | ||
| + | @ | ||
| + | |||
| + | \noindent Efetuando a krigagem | ||
| + | <<results=hide>>= | ||
| + | set.seed(512) | ||
| + | KC1<-krige.control(trend.d=~regiao, trend.l=~ind.reg, obj.model=Al.ml1$Al.4) | ||
| + | OC=output.control(n.pred=3) # define 3 simulacoes | ||
| + | |||
| + | pred<-krige.conv(Alum, loc=gr, krige=KC1,output=OC) | ||
| + | @ | ||
| + | %$ | ||
| + | |||
| + | \noindent A figura \ref{fig:mapas} é subdividida em quatro mapas. O primeiro trata-se de uma estimativa dos pontos com base na influência média de seus vizinhos, processo esse executado por uma krigagem convencional. Nos demais mapas os pontos são simulados com o modelo definido no ajuste de máxima-verossimilhança com seus respectivos parâmetros. | ||
| + | |||
| + | <<fig=F, results=hide>>= | ||
| + | par(mfrow=c(2,2), mar=c(1,2,1,1)) | ||
| + | image(pred, col=gray(seq(0.95,0.1,l=51))) | ||
| + | image(pred, val=pred$simul[,1],col=gray(seq(0.95,0.1,l=51))) | ||
| + | image(pred, val=pred$simul[,2],col=gray(seq(0.95,0.1,l=51))) | ||
| + | image(pred, val=pred$simul[,3],col=gray(seq(0.95,0.1,l=51))) | ||
| + | par(mfrow=c(1,1)) | ||
| + | @ | ||
| + | %$ | ||
| + | \begin{figure}[!ht] | ||
| + | \centering | ||
| + | \SweaveOpts{width=6,height=6} | ||
| + | \setkeys{Gin}{width=0.99\textwidth} | ||
| + | <<fig=T, echo=F, results=hide>>= | ||
| + | par(mfrow=c(2,2), mar=c(1,2,1,1)) | ||
| + | image(pred, col=gray(seq(0.95,0.1,l=51))) | ||
| + | image(pred, val=pred$simul[,1],col=gray(seq(0.95,0.1,l=51))) | ||
| + | image(pred, val=pred$simul[,2],col=gray(seq(0.95,0.1,l=51))) | ||
| + | image(pred, val=pred$simul[,3],col=gray(seq(0.95,0.1,l=51))) | ||
| + | par(mfrow=c(1,1)) | ||
| + | @ | ||
| + | %$ | ||
| + | \vspace{-.5cm} | ||
| + | \caption{Mapa de distribuição espacial da média (superior esquerdo) e de 3 simulações} | ||
| + | \label{fig:mapas} | ||
| + | \end{figure} | ||
| + | |||
| + | |||
| + | |||
| + | |||
| + | |||
