## par.ori <- par(no.readonly=T) ## importando dados camg <- read.table("regioes.dat", head=TRUE) head(camg) dim(camg) str(camg) require(geoR) #options(width=120) ## escolhendo uma variável para análise e duas como covariáveis ctc40 <- as.geodata(camg, coords.col=1:2, data.col=9, covar.col=c(3, 11)) summary(ctc40) ## verificando possível relação com covariáveis with(camg, plot(ctc40 ~ altitude)); with(camg, lines(lowess(ctc40 ~ altitude))) with(camg, plot(ctc40 ~ regiao)) ## gráficos com a inclusão de diferentes covariáveis plot(ctc40, low=T) plot(ctc40, low=T, trend=~regiao) plot(ctc40, low=T, trend=~altitude) plot(ctc40, low=T, trend=~altitude+regiao) ## verificando necessidade de transformação da resposta require(MASS) par(mfrow=c(2,2), mar=c(3,3,0, 0.5)) boxcox(ctc40) boxcox(ctc40, trend=~altitude) boxcox(ctc40, trend=~regiao) boxcox(ctc40, trend=~regiao+altitude) ## variogramas para diferentes modelos de média par(mfrow=c(2,2), mar=c(3,3,0.3, 0.3), mgp=c(1.8, .7, 0)) plot(variog(ctc40, max.dist=900)) plot(variog(ctc40, max.dist=900, trend=~regiao)) plot(variog(ctc40, max.dist=900, trend=~altitude)) plot(variog(ctc40, max.dist=900, trend=~regiao+altitude)) ## identificando possíveis dados atípicos/influentes par(par.ori) with(camg, boxplot(ctc40)) points(ctc40) points(camg[with(camg, which(ctc40 > 200)), 1:2], col=2, pch=19) ## gerando novo objeto excluindo atipico(s) ctc40.1 <- subset(ctc40, data < 200) class(ctc40.1) summary(ctc40.1) ## refazendo variogramas sem dado(s) atípico(s) par(mfrow=c(2,2), mar=c(3,3,0.3, 0.3), mgp=c(1.8, 0.7, 0)) plot(variog(ctc40.1, max.dist=900)) plot(variog(ctc40.1, max.dist=900, trend=~regiao)) plot(variog(ctc40.1, max.dist=900, trend=~altitude)) plot(variog(ctc40.1, max.dist=900, trend=~regiao+altitude)) ## Modelos: (diferentes funções de correlação e modelos para média (covariáveis) ## com todos os dados ctc40.ml0 <- list() ctc40.ml0$kappa0.5 <- likfit(ctc40, ini=c(600, 200), nug=100) ctc40.ml0$kappa1.5 <- likfit(ctc40, ini=c(600, 200), nug=100, kappa=1.5) ctc40.ml0$kappa2.5 <- likfit(ctc40, ini=c(600, 200), nug=100, kappa=2.5) sapply(ctc40.ml0, logLik) ctc40.ml.r <- list() ctc40.ml.r$kappa0.5 <- likfit(ctc40, ini=c(600, 200), nug=100, trend=~regiao) ctc40.ml.r$kappa1.5 <- likfit(ctc40, ini=c(600, 200), nug=100, kappa=1.5, trend=~regiao) ctc40.ml.r$kappa2.5 <- likfit(ctc40, ini=c(600, 200), nug=100, kappa=2.5, trend=~regiao) sapply(ctc40.ml.r, logLik) ctc40.ml.ra <- list() ctc40.ml.ra$kappa0.5 <- likfit(ctc40, ini=c(600, 200), nug=100, trend=~regiao+altitude) ctc40.ml.ra$kappa1.5 <- likfit(ctc40, ini=c(600, 200), nug=100, kappa=1.5, trend=~regiao+altitude) ctc40.ml.ra$kappa2.5 <- likfit(ctc40, ini=c(600, 200), nug=100, kappa=2.5, trend=~regiao+altitude) sapply(ctc40.ml.ra, logLik) ctc40.ml0 ctc40.ml.r ctc40.ml.ra ## comparando ajustes sapply(ctc40.ml0, logLik) sapply(ctc40.ml.r, logLik) sapply(ctc40.ml.ra, logLik) ## valores de referencia para o dobro da diferença em verossimilhanças qchisq(0.95, df=2) qchisq(0.95, df=1) ## outras medidas (uteis para modelos não encaixados sapply(ctc40.ml0, AIC) sapply(ctc40.ml.r, AIC) sapply(ctc40.ml.ra, AIC) sapply(ctc40.ml0, function(x) x$BIC) sapply(ctc40.ml.r, function(x) x$BIC) sapply(ctc40.ml.ra, function(x) x$BIC) ## repetindo acima para os dados de ctc40.1: ctc40.1.ml0 <- list() ctc40.1.ml0$kappa0.5 <- likfit(ctc40.1, ini=c(600, 200), nug=100) ctc40.1.ml0$kappa1.5 <- likfit(ctc40.1, ini=c(600, 200), nug=100, kappa=1.5) ctc40.1.ml0$kappa2.5 <- likfit(ctc40.1, ini=c(600, 200), nug=100, kappa=2.5) sapply(ctc40.1.ml0, logLik) ctc40.1.ml.r <- list() ctc40.1.ml.r$kappa0.5 <- likfit(ctc40.1, ini=c(600, 200), nug=100, trend=~regiao) ctc40.1.ml.r$kappa1.5 <- likfit(ctc40.1, ini=c(600, 200), nug=100, kappa=1.5, trend=~regiao) ctc40.1.ml.r$kappa2.5 <- likfit(ctc40.1, ini=c(600, 200), nug=100, kappa=2.5, trend=~regiao) sapply(ctc40.1.ml.r, logLik) ctc40.1.ml.ra <- list() ctc40.1.ml.ra$kappa0.5 <- likfit(ctc40.1, ini=c(600, 200), nug=100, trend=~regiao+altitude) ctc40.1.ml.ra$kappa1.5 <- likfit(ctc40.1, ini=c(600, 200), nug=100, kappa=1.5, trend=~regiao+altitude) ctc40.1.ml.ra$kappa2.5 <- likfit(ctc40.1, ini=c(600, 200), nug=100, kappa=2.5, trend=~regiao+altitude) sapply(ctc40.1.ml.ra, logLik) ## comparando ajustes sapply(ctc40.1.ml0, logLik) sapply(ctc40.1.ml.r, logLik) sapply(ctc40.1.ml.ra, logLik) ## modelo adotando: ctc40.1.ml.ra$kappa0.5 ## verificando se ajuste não está sendo influenciado pela escolha de valores iniciais options(geoR.messages=F) likfit(ctc40.1, ini=c(600, 50), nug=100, trend=~regiao+altitude) likfit(ctc40.1, ini=c(600, 150), nug=100, trend=~regiao+altitude) likfit(ctc40.1, ini=c(600, 250), nug=100, trend=~regiao+altitude) likfit(ctc40.1, ini=c(600, 300), nug=100, trend=~regiao+altitude) likfit(ctc40.1, ini=c(600, 350), nug=100, trend=~regiao+altitude) likfit(ctc40.1, ini=c(600, 350), nug=100, trend=~regiao+altitude) ## PROBLEMAS DE CONVERGENCIA: ## mudar: criterios de convergencia e/ou fc de minimização e/ou metodo likfit(ctc40.1, ini=c(600, 50), nug=100, trend=~regiao+altitude, lik.m="RML") likfit(ctc40.1, ini=c(600, 150), nug=100, trend=~regiao+altitude, lik.m="RML") likfit(ctc40.1, ini=c(600, 200), nug=100, trend=~regiao+altitude, lik.m="RML") likfit(ctc40.1, ini=c(600, 250), nug=100, trend=~regiao+altitude, lik.m="RML") ## grid de valores iniciais: ini.grid <- cbind(seq(400, 800, by=25), seq(0, 400, by=25)) ini.grid options(geoR.messages=T) ctc40.ml.ra <- likfit(ctc40.1, ini=ini.grid, nug=seq(400, 0, by=-25), trend=~regiao+altitude) ctc40.ml.ra ## alguns diagnósticos: ## envelope de variograma (empírico e ajustado) ## validação cruzada xvalid() ## Elegendo (com desconforto!!!!!) um modelo para prosseguir ## predição espacial ## pegando as bordas da regiao como um todo de outro objeto: data(ca20) ctc40$borders <- ca20$borders plot(ctc40) ## malha de pontos de predição args(pred_grid) gr <- pred_grid(ctc40$borders, by=20) par(par.ori) points(ctc40) points(gr, col=2, pch=19, cex=0.3) args(krige.conv) args(krige.control) ## construindo as covariáveis nas localizações de predição ## regiões ctc40$reg1 <- ca20$reg1 ctc40$reg1 ctc40$reg2 <- ca20$reg2 ctc40$reg2 ctc40$reg3 <- ca20$reg3 ctc40$reg3 points(ctc40) points(gr, col=2, pch=19, cex=0.3) polygon(ctc40$reg1) polygon(ctc40$reg2) polygon(ctc40$reg3) gr0 <- gr[.geoR_inout(gr, ctc40$borders),] points(gr0, col=4, pch=19, cex=0.3) dim(gr) dim(gr0) cov.regiao <- numeric(nrow(gr0)) cov.regiao[.geoR_inout(gr0, ctc40$reg1)] <- "R1" cov.regiao[.geoR_inout(gr0, ctc40$reg2)] <- "R2" cov.regiao[.geoR_inout(gr0, ctc40$reg3)] <- "R3" cov.regiao <- as.factor(cov.regiao) table(cov.regiao) ## interpolando altitude utilizando splines (pacote akima) install.packages("akima") require(akima) head(camg) cov.altitude <- interpp(camg[,1], camg[,2], camg[,3], gr0[,1], gr0[,2], linear=F, extrap=T)$z cov.altitude KC <- krige.control(type="OK", obj.model=ctc40.ml.ra, trend.d=~regiao+altitude, trend.l=~cov.regiao+cov.altitude) length(cov.regiao) length(cov.altitude) dim(gr0) ctc40.kc <- krige.conv(ctc40, loc=gr0, krige=KC) image(ctc40.kc, gr, ctc40$borders, col=terrain.colors(21), x.leg=c(4950, 5450), y.leg=c(4770, 4830)) polygon(ctc40$reg1) polygon(ctc40$reg2) ## explorando mais os resultados: image(ctc40.kc, loc=gr, bor=ctc40$borders, col=terrain.colors(21), val=sqrt(ctc40.kc$krige.var), x.leg=c(4950, 5450), y.leg=c(4770, 4830)) polygon(ctc40$reg1) polygon(ctc40$reg2) args(output.control) OC <- output.control(thres=150, quantile=c(0.1, 0.9)) ctc40.kc <- krige.conv(ctc40, loc=gr0, krige=KC, out=OC) names(ctc40.kc) dim(ctc40.kc$simula) ## mapa de probabilidades: P(S(x)|y > 150) image(ctc40.kc, loc=gr, bor=ctc40$borders, col=terrain.colors(21), val=1-ctc40.kc$prob,x.leg=c(4950, 5450), y.leg=c(4770, 4830)) polygon(ctc40$reg1) polygon(ctc40$reg2) ## quantis names(ctc40.kc$quant) ## mapa do 90o percentil -- "baixos confiáveis": P(S(x)|y < a) = 0.90 image(ctc40.kc, loc=gr, bor=ctc40$borders, col=terrain.colors(21), val=ctc40.kc$quantile$q90,x.leg=c(4950, 5450), y.leg=c(4770, 4830)) polygon(ctc40$reg1) polygon(ctc40$reg2) ## mapa do 10o percentil -- "altos confiáveis": P(S(x)|y < a) = 0.10 image(ctc40.kc, loc=gr, bor=ctc40$borders, col=terrain.colors(21), val=ctc40.kc$quantile$q10,x.leg=c(4950, 5450), y.leg=c(4770, 4830)) polygon(ctc40$reg1) polygon(ctc40$reg2) ## outros resultados (outros funcionais) names(ctc40.kc) dim(ctc40.kc$simul) length(ctc40.kc$mean.sim) length(ctc40.kc$sim.means) ## histograma de média geral hist(ctc40.kc$sim.means) ## histograma da proporcao da área superior a certo valor (150) p150 <- apply(ctc40.kc$simul, 2, function(x) mean(x > 150)) hist(p150, prob=T); lines(density(p150)) quantile(p150, probs=c(0.025, 0.975)) ## avaliação do erro de Monte Carlo (simulacoes) names(ctc40.kc) with(ctc40.kc, summary(mean.simulations - predict)) ## poderia avaliar outras quantidades!! ## se necessário modifique argumento n.predictive em output.control() ## default é 1000