rm(list=ls(all=TRUE)) setwd("C:/Users/Gabriel/Desktop/Geoestatistica/Simula") # Gerando dados com distribuição de Poisson require(geoR) require(MASS) set.seed(126) data(parana) s <- grf(grid=parana$coords, cov.pars=c(1.2, 70), mean=1, cov.model="exp") # gerando dados com correlaçao espacial c(sigma2,phi) y <- rpois(length(s$data), lambda=exp(s$data)) # gerando dados poison dt <- as.geodata(cbind(parana$coords,y,1), units.m.col=4) # geodadta dt$borders <- parana$borders # x11() # plot(dt) # Modelo normal com transformação log dt2 <- dt dt2$data <- dt2$data+.5 boxcox(dt2) # logo tranf log ef <- eyefit(variog(dt2, lam=0, max.dist=400)) lf <- likfit(dt2, lam=0, ini=ef) grd <- expand.grid(seq(min(parana$borders[,1]),max(parana$borders[,1]),l=100), seq(min(parana$borders[,2]),max(parana$borders[,2]),l=100)) # points(parana); points(grd, pch='.') kr <- krige.conv(dt2, loc=grd, krige=krige.control(obj.model=lf)) # Modelo Poisson require(geoRglm) model1 <- list(cov.pars=c(1,50), beta=1, family="poisson") # c(sigma2=1,phi=50) mcmc1 <- mcmc.control(S.scale=0.38, n.iter=10000, burn.in=1000) # passo para test1 <- glsm.mcmc(dt, model=model1, mcmc.input=mcmc1) # mudar passao até ter 60, 38 é bom premodel <- prepare.likfit.glsm(test1) model2 <- likfit.glsm(premodel, ini.phi=40, cov.model='exp') mcmc2 <- mcmc.control(S.scale=0.40, n.iter=10000, burn.in=1000) test2 <- glsm.mcmc(dt, model=model2, mcmc.input=mcmc2) kr.pois <- glsm.krige(test2, loc=grd) # image(kr.pois) kr.pois2 <- kr kr.pois2$predict <- kr.pois$predict kr.pois2$krige.var <- kr.pois$krige.var zl <- range(c(range(kr$predict),range(kr.pois$predict))) par(mfrow=c(1,2)) image(kr, col=heat.colors(100), zlim=zl); contour(kr, add=T, zlim=zl, breaks=10) image(kr.pois2, col=heat.colors(100), zlim=zl); contour(kr.pois2, add=T, zlim=zl, breaks=10) par(mfrow=c(1,1)) fm <- lm(kr.pois2$predict ~ I(kr$predict) + I(kr$predict^2)) anova(fm); summary(fm) require(car); Anova(fm) f <- function(a,b,c,x) a+b*x+c*x^2 plot(kr$predict, kr.pois2$predict, cex=.6, xlim=c(0,20), ylim=c(0,20)) abline(0, 1, col=2, lwd=2, lty=2) curve(f(coef(fm)[1], coef(fm)[2], coef(fm)[3], x), add=T, col=3, lwd=2) lf$parameters model2$parameters