par.ori <- par(no.readonly=T) options(width=70) ## lendo dados q2 <- read.table("quimicosB.txt", head=T, dec=",") ## lendo bordas da área e sub-áreas bor <- read.table("borda.txt") a1 <- read.table("area1.txt") a2 <- read.table("area2.txt") a3 <- read.table("area3.txt") head(q2) names(q2) q2$Regiao <- as.factor(q2$Regiao) summary(q2) apply(q2[,1:2], 2 , range) ## ajustando valores das coordenadas (opcional) q2$X <- (q2$X - 688000) q2$Y <- (q2$Y - 7100000) ## require(geoR) ## Escolhendo a variável cálcio e "montando" um geodata Ca <- as.geodata(q2, coords.col=1:2, data.col=4, covar.col=13:14) points(Ca) polygon(bor) polygon(a1, col=2, density=5) polygon(a2, col=3, density=5) polygon(a3, col=4, density=5) Ca$borders <- bor ## note problema com valor negativo de teor de cálcio summary(Ca) ## aqui vamos (arbitrariamente) remover os pontos com dados não positivos Ca <- subset(Ca, data > 0) summary(Ca) ## iniciando algumas análises exploratórias plot(Ca) plot(Ca, low=T) ## investigando relação com covariáveis plot(Ca$covar$Prof, Ca$data) lines(lowess(Ca$data~Ca$covar$Prof)) boxplot(Ca$data~Ca$covar$Regi, notch=T, varwidth=T) plot(Ca, trend=~Regiao) ## notas: potenciais outliers ## localizacao dos "outliers" ## faz sentido agrupar 1 e 3 (vs 2) ?? ## investigando transformação require(MASS) boxcox(Ca) boxcox(Ca, trend=~Regiao) par(mfrow=c(1,2)) boxcox(Ca, lam=seq(0,1,l=201)) boxcox(Ca, trend=~Regiao, lam=seq(0,1,l=201)) ## portanto neste caso a transformação prescrita independe do modelo adotado para média ## Modelos candidatos: ## Pendências a ser discutidas no contexto dos dados: ## - valores negativos ## - covariável região e possibilidade de agrupar valores ## - outliers globais ## - outliers locais par(par.ori) #points(Ca) #points(Ca, cex.min=0.3, cex.max=3) points(Ca, lam=0.5) par(mfrow=c(1,2)) points(Ca, cex.min=0.3, cex.max=3, lam=0.5) polygon(a1) polygon(a2) points(Ca, trend=~Regiao, cex.min=0.3, cex.max=3, lam=0.5) polygon(a1) polygon(a2) ## explorando variogramas summary(Ca) v0 <- variog(Ca, max.dist=150, lam=0.5) v1 <- variog(Ca, max.dist=150, lam=0.5, trend=~Regiao) plot(v0) plot(v1) ## experimentar com diferentes opções de variogramas empíricos!! ## seguem algumas entre muitas possíveis opções v0 <- variog(Ca, max.dist=100, lam=0.5) v1 <- variog(Ca, max.dist=100, lam=0.5, trend=~Regiao) ## vamos tb garantir a mesma escala em ambos gráficos vmax <- max(c(v0$v, v1$v)) plot(v0, ylim=c(0, vmax)) plot(v1, ylim=c(0, vmax)) v0 <- variog(Ca, max.dist=100, lam=0.5, uvec=seq(0,120, l=10) ) v1 <- variog(Ca, max.dist=100, lam=0.5, trend=~Regiao, uvec=seq(0,120, l=10)) plot(v0, ylim=c(0, vmax)) plot(v1, ylim=c(0, vmax)) v0 <- variog(Ca, lam=0.5, uvec=seq(0,250, l=8) ) v1 <- variog(Ca, lam=0.5, trend=~Regiao, uvec=seq(0,250, l=8)) plot(v0, ylim=c(0, vmax)) plot(v1, ylim=c(0, vmax)) v0 <- variog(Ca, lam=0.5, uvec=seq(0,80, l=10) ) v1 <- variog(Ca, lam=0.5, trend=~Regiao, uvec=seq(0,80, l=10)) plot(v0, ylim=c(0, vmax)) plot(v1, ylim=c(0, vmax)) ## ver tb outras opções na função variog() ## Ajustando, avaliando e escolhendo modelo(s) ## Nota: pode-se usar eyefit() sobre variograms para indicar os valores iniciais ## para o algorítmo de ajuste de modelos mods <- list() mods$m0.exp <- likfit(Ca, ini=c(0.35, 10), nug=0.10, lam=0.5) mods$m1.exp <- likfit(Ca, ini=c(0.35, 10), nug=0.10, trend=~Regiao, lam=0.5) mods$m0.exp mods$m1.exp ## um ajuste mais cuidadoso considerando múltiplos valores iniciais inis <- cbind(cbind(seq(0.2, 0.6, l=6), seq(0, 50, l=6))) mods$m0.exp <- likfit(Ca, ini=inis, nug=seq(0, 0.4, l=5), lam=0.5) mods$m1.exp <- likfit(Ca, ini=inis, nug=seq(0, 0.4, l=5), trend=~Regiao, lam=0.5) mods$m0.exp mods$m1.exp ## visualizando ajuste ("tarja preta" aqui!!!) plot(v0, ylim=c(0, vmax)) lines(m0.exp) plot(v1, ylim=c(0, vmax)) lines(m1.exp) ## considerando mais uma potencial covariável mods$m2.exp <- likfit(Ca, ini=inis, nug=seq(0, 0.4, l=5), trend=~Regiao+Prof, lam=0.5) mods$m2.exp mods$m2.exp <- NULL ## considerando agora diferentes opções para fç de covariância Matérn mods$m0.15 <- likfit(Ca, ini=inis, nug=seq(0, 0.4, l=5), lam=0.5, kappa=1.5) mods$m1.15 <- likfit(Ca, ini=inis, nug=seq(0, 0.4, l=5), trend=~Regiao, lam=0.5, kappa=1.5) mods$m0.25 <- likfit(Ca, ini=inis, nug=seq(0, 0.4, l=5), lam=0.5, kappa=2.5) mods$m1.25 <- likfit(Ca, ini=inis, nug=seq(0, 0.4, l=5), trend=~Regiao, lam=0.5, kappa=2.5) mods$m0.exp mods$m0.15 mods$m0.25 summary(mods$m0.exp) mods$m1.exp mods$m1.15 mods$m1.25 summary(mods$m1.exp) names(mods) sapply(mods, function(x) x) sapply(mods, function(x) x$loglik) sapply(mods, function(x) with(x,c(loglik,AIC,BIC))) ## comentários: ## - poderia-se verificar a REML restrita ## - verificar impacto de outlier global nos resultados ## - idem para locais ## Note outras opções a serem investigadas! ## modelo escolhido será m1.exp cafit <- likfit(Ca, ini=c(0.35, 10), nug=0.10, trend=~Regiao, lam=0.5, comp=T) cafit$model.comp names(cafit) args(likfit) ## interpolação espacial par(par.ori) points(Ca) ## definindo o grid de predição gr <- pred_grid(bor, by=10) dim(gr) points(gr, pch=19, cex=0.25) gr0 <- polygrid(gr, border=bor) dim(gr0) points(gr0, pch=19, cex=0.25, col=2) ## construindo a covariável do modelo area.gr <- numeric(nrow(gr0)) area.gr[.geoR_inout(gr0, a1)] <- 1 area.gr[.geoR_inout(gr0, a2)] <- 2 area.gr[.geoR_inout(gr0, a3)] <- 3 table(area.gr) area.gr <- as.factor(area.gr) ## interpolação args(krige.conv) args(krige.control) KC <- krige.control(trend.d=~Regiao, trend.l=~area.gr, obj.mod=cafit) ## interpolação ca.kc <- krige.conv(Ca, loc=gr0, krige=KC) names(ca.kc) ## passando argumentos para garantir que a fc image() consiga ## montar a estrutura retangular da dados image(ca.kc, loc=gr, border=bor) ## ou..... ca.kc <- krige.conv(Ca, loc=gr, border=bor, krige=KC) image(ca.kc) image(ca.kc, x.leg=c(600, 1000), y.leg=c(89550,89600)) image(ca.kc, col=gray(seq(1,0.15, l=51)), x.leg=c(600, 1000), y.leg=c(89550,89600)) polygon(a1) polygon(a2) image(ca.kc, val=sqrt(ca.kc$krige.var),col=gray(seq(1,0.15, l=51)), x.leg=c(600, 1000), y.leg=c(89550,89600)) polygon(a1) polygon(a2) ## opções extras de saídas args(output.control) OC <- output.control(n.pred=500, simulation=T, mean.var=T, quantile=c(0.1, 0.5, 0.75), thres=1.8, sim.m=T, sim.vars=T) ca.kc <- krige.conv(Ca, loc=gr0, krige=KC, output=OC) names(ca.kc) attach(ca.kc) dim(simulations) length(mean.simulations) dim(quantiles.simulations) all.equal(quantiles.simulations[,2],apply(simulations, 1, median)) ## mapa da medianas image(ca.kc, value=ca.kc$quant[,2],col=gray(seq(1,0.15, l=51))) polygon(a1) polygon(a2) ## verificar outras saidas ## outros resumos de interesse podem ser obtidos da simulacao. ## exemplo: distr. da prop. da area com teores abaixo de 1.5 ca.prop1.5 <- apply(ca.kc$simulations, 2, function(x) mean(x < 1.5)) mean(ca.prop1.5) hist(ca.prop1.5) names(ca.kc) ## legend.kriging() ## agora volte e defina um grid mais fino ## Análise Bayesiana args(krige.bayes) args(model.control) MC <- model.control(trend.d=~Regiao, trend.l=~area.gr, lam=0.5) args(prior.control) PC <- prior.control(phi.prior="reciprocal", phi.discrete=seq(0,100,l=11), tausq.rel.pr="rec", tausq.rel.disc=seq(0,1, l=11)) ca.kb <- krige.bayes(Ca, model=MC, prior=PC) PC <- prior.control(phi.prior="reciprocal", phi.discrete=seq(0,100,l=11), tausq.rel.pr="unif", tausq.rel.disc=seq(0,1, l=11)) ca.kb <- krige.bayes(Ca, model=MC, prior=PC) PC <- prior.control(phi.prior="reciprocal", phi.discrete=seq(0,40,l=51)) ca.kb <- krige.bayes(Ca, model=MC, prior=PC) par(par.ori) plot(ca.kb) hist(ca.kb) ## experimentar MUITO com outras prioris etc etc ## PS: ## "caldo de galinha!" contruindo bordas com criatividade :) ##via locator points(Ca, bor=NULL) bor <- locator(type="l") bor <- matrix(unlist(bor), nc=2) ## via poligono convexo require(MASS) bor <- ap$coords[chull(ap$coords),] polygon(bor) ## e nao esqueça do mapa na tela!!!