sal <- read.table("sal.dat", head=T, row.names=1) str(sal) require(geoR) salg <- as.geodata(sal, data.col=10) plot(salg) points(salg) bor <- locator(15) salg$borders <- matrix(unlist(bor), nc=2) points(salg) summary(salg) require(MASS) boxcox(salg) plot(salg, lam=0.5, low=T) salg <- as.geodata(sal, data.col=10, covar.col=1) salg$borders <- matrix(unlist(bor), nc=2) summary(salg) boxcox(salg, trend=~COORX) plot(salg, lam=0.5, low=T, trend=~COORX) plot(variog(salg, max.dist=20000)) ##practicalRange(10, cov.model="mat", kapp=2.5)/10 args(krige.bayes) args(model.control) args(prior.control) ## fazendo primeiro sem predição MC <- model.control(trend.d = ~salg$covar$COORX, lam=0.5) PC <- prior.control(phi="reciprocal", phi.discrete=seq(0, 15000, length=16)) kb1 <- krige.bayes(salg, model=MC, prior=PC) ## problema numérico, cvamos redefinir as coordenadas min.x <- min(salg$borders[,1]) min.y <- min(salg$borders[,2]) salg$coords[,1] <- (salg$coords[,1] - min.x)/1000 salg$coords[,2] <- (salg$coords[,2] - min.y)/1000 salg$covariate$COORX <- (salg$covariate$COORX - min.x)/1000 salg$borders[,1] <- (salg$borders[,1] - min.x)/1000 salg$borders[,2] <- (salg$borders[,2] - min.y)/1000 summary(salg) MC <- model.control(trend.d = ~COORX, lam=0.5) #PC <- prior.control(phi.prior="reciprocal", phi.discrete=seq(0, 15, length=16)) #PC <- prior.control(phi.prior="reciprocal", phi.discrete=seq(0, 25, length=26)) PC <- prior.control(phi.prior="square", phi.discrete=seq(0, 30, length=31)) OC <- output.control(n.post=10000) kb1 <- krige.bayes(salg, model=MC, prior=PC, output=OC) plot(kb1) #PC <- prior.control(phi.prior="square", phi.discrete=seq(0, 30, length=31), # tausq.rel.disc=seq(0, 2, l=21), tausq.rel.prior="unif") PC <- prior.control(phi.prior="square", phi.discrete=seq(0, 30, length=16), tausq.rel.disc=seq(0, 1, l=11), tausq.rel.prior="unif") kb1 <- krige.bayes(salg, model=MC, prior=PC, output=OC) par(mfrow=c(1,2)) plot(kb1) names(kb1) names(kb1$post) names(kb1$post$sample) par(mfrow=c(2,3), mar=c(3,3,.5, .5)) lapply(kb1$post$sample, hist, main="") par(mfrow=c(2,3), mar=c(3,3,.5, .5)) lapply(kb1$post$sample, function(x) plot(density(x))) par(mfrow=c(2,3), mar=c(3,3,.5, .5)) lapply(kb1$post$sample, function(x){hist(x, prob=T, main=""); lines(density(x))}) ## verificando se a covariavel é relevante pelo intervalo de credibilidade ## (95% baseado em quantis) quantile(kb1$post$sample$beta1, prob=c(0.025, 0.975)) (kb1$post$joint) (kb1$post$phi) (kb1$post$tau) (kb1$post$sigmasq) (kb1$post$beta) head(kb1$post$sample) dim(kb1$post$sample) ## prediçao gr <- pred_grid(salg$borders, by = 0.5) par(mfrow=c(1,1)) points(salg) points(gr, pch=21, cex=0.3, col=2) ind <- .geoR_inout(gr, salg$borders) gr0 <- gr[ind,] points(gr0, pch=21, cex=0.3, col=4) dim(gr) dim(gr0) MC <- model.control(trend.d = ~COORX, trend.l = ~gr0[,1], lam=0.5) OC <- output.control(n.post=10000, n.pred=1000) ## o comando a seguir pode ser MUUUUUUUUUIIIIIITOOOO demorado para rodar! kb.p <- krige.bayes(salg, loc=gr, borders=salg$borders, model=MC, prior=PC, output=OC)