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 code 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) ## another example for Poisson data ## simulating data set.seed(371) cp <- expand.grid(seq(0,1,l=10),seq(0,1,l=10)) s <- grf(grid=cp, cov.pars=c(2, 0.2), cov.model="mat", kappa=1.5) y <- rpois(length(s$data), lambda=exp(0.5 + s$data)) dt <- as.geodata(cbind(cp,y,1), units.m.col=4) plot(dt) image(matrix(s$data, nc=10), asp=1, col=gray(seq(1,0.2,l=21))) text(cp[,1], cp[,2], as.character(dt$data), cex=1.6) ## Poisson model set.seed(371) #MCc <- mcmc.control(S.scale=0.025, phi.sc=0.1, n.iter=110000, burn.in=10000, MCc <- mcmc.control(S.scale=0.025, phi.sc=0.07, n.iter=50000, burn.in=5000, thin=50, phi.start=0.2) PGC <- prior.glm.control(phi.prior="exponential", phi=0.2, phi.discrete=seq(0,2,by=0.02),tausq.rel=0) OC <- output.glm.control(sim.pred=T) locs <- cbind(c(0.75, 0.27), c(0.25, 0.6)) points(locs, col=2, pch=19, cex=1.5) pkb <- pois.krige.bayes(dt, loc=locs, prior=PGC, mcmc=MCc, out=OC) names(pkb) names(pkb$post) dim(pkb$post$sim) names(pkb$pred) dim(pkb$pred$sim) par(mfrow=c(1,3), mar=c(3.5,3.5,0.5,0.5), mgp=c(1.8, 0.8, 0)) hist(pkb$post$beta$sample, prob=T, xlab=expression(beta), main="");lines(density(pkb$post$beta$sample)); abline(v=c(median(pkb$post$beta$sample), mean(pkb$post$beta$sample)), col=4) hist(pkb$post$sigma$sample, prob=T, xlab=expression(sigma), main="");lines(density(pkb$post$sigma$sample)); abline(v=c(median(pkb$post$sigma$sample), mean(pkb$post$sigma$sample)), col=4) hist(pkb$post$phi$sample, prob=T, xlab=expression(phi), main="");lines(density(pkb$post$phi$sample)); abline(v=c(median(pkb$post$phi$sample), mean(pkb$post$phi$sample)), col=4) plot(dt) image(matrix(s$data, nc=10), asp=1, col=gray(seq(1,0.2,l=21))) text(cp[,1], cp[,2], as.character(dt$data), cex=1.6) points(locs, col=2, pch=19, cex=1.5) hist(pkb$pred$sim[1,], prob=T, xlab="(0.75, 0.27)", main="");lines(density(pkb$pred$sim[1,])); abline(v=c(median(pkb$pred$sim[1,]), mean(pkb$pred$sim[1,])), col=4) hist(pkb$pred$sim[2,], prob=T, xlab="(0.25, 0.6)", main="");lines(density(pkb$pred$sim[2,])); abline(v=c(median(pkb$pred$sim[2,]), mean(pkb$pred$sim[2,])), col=4) par(par.ori) ## ## ## ## MCMC likelihood data(rongelap) args(glsm.mcmc) MCmle.input.fixed <- glsm.mcmc(rongelap, model=list(family="poisson", link = "boxcox", cov.pars=c(2.40, 340), beta = 6.2, nugget = 2.075*2.40, lambda=1), mcmc.input = mcmc.control(S.scale = 0.5,thin=20,burn.in=10000)) ## investigating mixing and convergence library(coda) MCmle.coda <- create.mcmc.coda(MCmle.input.fixed, mcmc.input = list(thin = 20, burn.in=10000)) plot(MCmle.coda) autocorr.plot(MCmle.coda) ## maximum likelihood args(prepare.likfit.glsm) mcmcobj <- prepare.likfit.glsm(MCmle.input.fixed) args(likfit.glsm) lik.boxcox.1.expon <- likfit.glsm(mcmcobj, ini.phi=10, fix.nugget.rel=TRUE)