require(geoR) ## ## Importing a data set from a data file (text format) ## options(width=70) cadf <- read.table("ca.dat", head=T) head(cadf) summary(cadf) names(cadf) cadf$area <- as.factor(cadf$area) require(geoR) ca <- as.geodata(cadf, coords.col=1:2, data.col=3, covar.col=4:5) class(ca) summary(ca) ca <- read.geodata("ca.dat", coords.col=1:2, data.col=3, covar.col=4:5, head=T) args(read.geodata) ca <- read.geodata("ca.dat", covar.col=4:5, head=T) names(ca) summary(ca) ca$covariate <- transform(ca$covariate, area= as.factor(area)) plot(ca) summary(ca) ## ## load data ca20 from geoR ## data(ca20) summary(ca20) names(ca20) plot(ca20) points(ca20, pt.div="quint") polygon(ca20$reg1) polygon(ca20$reg2) ## plot residuals from a l.m. fit on covariates plot(ca20, trend=~area+altitude) ## ## other data-sets ## data(parana) plot(parana) ## plot residuals from a lm fit on 1st degree polynomial on the coordinates plot(parana, trend="1st") ## plot transformed data data(Ksat) plot(Ksat) plot(Ksat, lambda=0) ## ## exploring the correlation structure ## v1 <- variog(parana, max.dist=400) v2 <- variog(parana, trend="1st", max.dist=400) par(mfrow=c(1,2), mar=c(3.5,3.5,1,1), mgp=c(2,1,0)) plot(v1) plot(v2) par(mfrow=c(2,2), mar=c(3.5,3.5,1,1), mgp=c(2,1,0)) v3 <- variog(ca20, max.dist=400) v4 <- variog(ca20, max.dist=400, trend=~area) v5 <- variog(ca20, max.dist=400, trend=~area+altitude) plot(v3) plot(v4) plot(v5) par(mfrow=c(1,1)) ## ## experimenting with variogram models ## plot(v4) v4e <- eyefit(v4) v4f <- variofit(v4, ini=v4e) lines(v4f, lty=2) args(variofit) ml <- likfit(ca20, trend=~area, ini=v4e) ## just checking --- but does not nedd to fit the variogram!!! lines(ml, col=2) ## ## spatial interpolation ## args(krige.conv) args(krige.control) KC <- krige.control(trend.d = ~area, trend.l = ~ind.reg, obj.model=ml) args(output.control) ## prediction locations gr <- pred_grid(ca20$borders, by=20) gr0 <- polygrid(gr, borders=ca20$border, bound=T) ind.reg <- numeric(nrow(gr0)) ind.reg[.geoR_inout(gr0, ca20$reg1)] <- 1 ind.reg[.geoR_inout(gr0, ca20$reg2)] <- 2 ind.reg[.geoR_inout(gr0, ca20$reg3)] <- 3 ind.reg <- as.factor(ind.reg) points(ca20) points(gr, pch=19, cex=.3) points(gr0, col=4, pch=19, cex=.3) kc <- krige.conv(ca20, loc=gr, krige=KC) par(mfrow=c(1,2)) image(kc, col=gray(seq(1,0,l=21))) image(kc, val=sqrt(kc$kr),col=gray(seq(1,0,l=21))) ## ## Bayesian analysis ## args(krige.bayes) args(prior.control) PC <- prior.control(phi.disc=seq(0,1000, l=21), tausq.rel.disc=seq(0,1, l=11), tausq.rel.prior="unif") args(model.control) MC <- model.control(trend.d=~area, trend.l=~ind.reg) args(output.control) OC <- output.control() ## run to obtain only the posterior #kb <- krige.bayes(ca20, prior=PC, model=MC, output=OC) par(mfrow=c(1,2)) plot(kb) ## run for predictive kb <- krige.bayes(ca20, prior=PC, model=MC, output=OC) ## ## GLGM's ## ## 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)) 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.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) ## 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.1, n.iter=5500, burn.in=500, thin=100, 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.15), c(0.25, 0.5)) pkb <- pois.krige.bayes(dt, loc=locs, prior=PGC, mcmc=MCc, out=OC) ## 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)