par.ori <- par(no.readonly=TRUE) ## ## GLGM's ## require(geoR) ## binary data example ## simulating data set.seed(34) locs <- seq(0,1,l=401) s <- grf(grid=cbind(locs,1), cov.pars=c(5, 0.1), cov.model="matern", kappa=1.5) p <- exp(s$data)/(1+exp(s$data)) plot(locs, p, ylim=c(0,1), ty="l") ind <- seq(1,401,by=8) y <- rbinom(p[ind], size=1, prob=p[ind]) plot(locs[ind],y, xlab="locations", ylab="data") lines(locs,p) d.bin <- as.geodata(cbind(s$coords[ind,], y)) d.bin$units.m <- rep(1, length(d.bin$data)) ## geoRglm código de análise require(geoRglm) args(binom.krige) args(krige.glm.control) model <- krige.glm.control(cov.pars=c(5, 0.1), cov.model="matern", kappa=1.5, beta =0) args(mcmc.control) #mcmc.c <- mcmc.control(S.scale=0.1, thin=100, burn.in = 30000) mcmc.c <- mcmc.control(S.scale=0.25, thin=100, burn.in = 30000) args(output.glm.control) out <- output.glm.control(sim.predict=F) set.seed(217) pred <- binom.krige(d.bin, loc=cbind(locs,1), krige=model, mcmc.in=mcmc.c, output=out) lines(locs, pred$pred, lty=2) lines(locs, pred$pred - 1.96*sqrt(pred$kr), lty=3) lines(locs, pred$pred + 1.96*sqrt(pred$kr), lty=3)