## ## Exemplos de análises geoestatísticas ## ## - Consideramos algumas possíveis variações do "modelo básico" ## modelo básico : méda constante, dados gaussianos (normais), etc ## Y(x) = \mu + S(x) + e ## ## - Possíveis variações consideradas nos exemplos: ## * Médias não contantes \mu --> m(x) ## - médias como função das coordenadas ## - médias como função de covariáveis ## * variáveis (Y) não gaussianos/normais ## ## Exemplo 1: wolfcamp ## # lendo/importando as dados como um data-frame wolf <- read.table("wolfcamp.txt", head=T, row.names=1) head(wolf) summary(wolf) ## carregando o pacote geoR require(geoR) ## convertendo para o formato(classe) geodata da geoR class(wolf) wolf.gd <- as.geodata(wolf, coords.col=1:2, data.col=3) class(wolf.gd) ## OBS: pode-se ler e converter direto do arquivo texto para classe geodata wolf.gd <- read.geodata("wolfcamp.txt", head=T, coords.col=1:2, data.col=3, row.names=1) class(wolf.gd) ## Resumo dos dados summary(wolf.gd) wolf.gd$call ## Visualizando os dados plot(wolf.gd) # ver gráfico plot(wolf.gd, low=T) # ver gráfico ## No gráfico anterior nota-se uma "tendencia" com as coordenadas ## vamos entao fazer uma regressão as coordenadas e fazer um gráfico dos resíduos ## ou seja, passamos a ter uma média n~ao contante e considarar o modelo ## Y(x) = \mu(x) + S9x) + e plot(wolf.gd, low=T, trend="1st") ## var gráfico a comparar com anterior ## OBS: a função permite usar vários outros argumentos que pdoem ser explorados args(plot.geodata) ## Outras visualizações, agora apenas com os pontos points(wolf.gd) points(wolf.gd, trend="1st") points(wolf.gd, trend="1st", pt.div="quint") args(points.gd) ## ## A " tendencia dos dados tem grande impacto no variograma dos dados originais ## a seguir comparamos os variogramas: ## dos dados originais: Y(x) ## dos resíduos estimados: Y(x) - \mu(x) par(mfrow=c(1,2)) wolf.vg <- variog(wolf.gd, max.dist=200) plot(wolf.vg) wolf.vg.res <- variog(wolf.gd, trend="1st", max.dist=200) plot(wolf.vg.res) par(mfrow=c(1,1)) ## vamos portanto considerar dois (na verdade 4 levado em consideracao que cada modelo piode ou não ter o efeito espacial) ## modelos "candidatos" ## nos ajustes dos modelos vamos usar os variogramas empíricos para ter idéia de valores iniciais ## 1. Modelo de média contante ## 2. Modelo de média não contante dada para uma tendência linear nas coordenadas ## 3. Modelo de média não contante dada para uma tendência quadrática nas coordenadas ## wolfcamp wolf.lik <- list() wolf.lik$fit0 <- likfit(wolf.gd, ini=c(25000, 30), nug=5000) wolf.lik$fit1 <- likfit(wolf.gd, ini=c(3500, 30), nug=500, trend="1st") wolf.lik$fit2 <- likfit(wolf.gd, ini=c(3500, 30), nug=500, trend="2nd") wolf.lik$fit0 wolf.lik$fit1 wolf.lik$fit2 summary(wolf.lik$fit0) summary(wolf.lik$fit1) summary(wolf.lik$fit2) sapply(wolf.lik, logLik) ## log Lik sapply(wolf.lik, AIC) ## AIC sapply(wolf.lik, AIC, k=log(nrow(wolf.gd$coords))) ## BIC ## a análise das verossimilhanças, AIC e BIC destes modelos sugerem que os modelos com tendência se ajustam melhor. ## vamos optar pelo modelo con tendência linear ## note que o resultado confirma a relevancia do termo espacial mesmo com a média não constante summary(wolf.lik$fit1) ## Interpolação espacial ## para fazer a interpoalção vamos inicialmente definir as "bordas" da área ## o ideal seria ter esta informação nos dados mas na ausência disto vamos definir uma região arbitrária ## ao redor dos pontos regiao <- cbind(c(-245, -245, -80, -80, 190, 190, -245), c(-160, 0, 0, 150, 150, -160, -160)) points(wolf.gd) polygon(regiao) ## por conveniência e seguindo o padrão dos objetos geodata vamos colocar as bordas "dentro" do objeto de dados sob o nome "borders" ## com isto elas são "automaticamente" plotadas wolf.gd$borders <- regiao points(wolf.gd) ## agora vamos definir um grid de predição entre os menores e maiores valores das coordenadas das bordas gr <- expand.grid(seq(-245, 190, length=50), seq(-160, 150, length=50)) ## e conferir se o grid cobre a área points(wolf.gd) points(gr, col=2, cex=0.25, pch=19) ## OBS: se os valores não fossem conhecidos poderiam ser obtidos com: apply(wolf.gd$borders, 2, range) ## agora a predição espacial wolf.pred <- krige.conv(wolf.gd, loc=gr, krige=krige.control(obj.model=wolf.lik$fit1)) image(wolf.pred) image(wolf.pred, coords.data=T) ## includindo a localização dos dados ## mapa dos erros padrao de predição image(wolf.pred, val=sqrt(krige.var)) image(wolf.pred, coords.data=T, val=sqrt(krige.var)) ## includindo a localização dos dados