################## # Section 13.3.1 # ################## library(blmeco) data(redstart) mod <- lm(log1p(counts) ~ elevation + I(elevation^2), data = redstart) library(sp) library(gstat) rs <- rstandard(mod) # Create a SpatialPointsDataFrame spdata <- data.frame(resid = rs, x = redstart$x, y = redstart$y) coordinates(spdata) <- c("x", "y") bubble(spdata, main = "Residuals") jpeg(filename = "Figures/Figure13-1_bubble.jpg", width = 6000, height = 6000, units = "px", res=1200) bubble(spdata, main = "Residuals", col=c("blue", "orange")) dev.off() vario <- variogram(resid ~ 1, spdata, cutoff = 50000) plot(vario, pch = 16, cex = 1.5) jpeg(filename = "figures/Figure13-2_vario_sample.jpg", width = 6000, height = 6000, units = "px", res=1200) vario <- variogram(resid ~ 1, spdata, cutoff = 50000) plot(vario, pch = 16, cex = 1.5) dev.off() ################## # Section 13.3.2 # ################## library(lattice) jpeg(filename = "figures/Figure13-3_vario_models.jpg", width = 12000, height = 6000, units = "px", res=1200) show.vgms(models = c("Exp", "Gau", "Sph"), nugget = 0.2, sill = 0.8, strip = strip.custom(factor.levels = c("Exponential", "Gaussian", "Spherical"))) dev.off() ################## # Section 13.3.3 # ################## library(spBayes) # Standardize the covariates redstart$elevation.z <- as.numeric(scale(redstart$elevation)) redstart$elevation.z2 <- redstart$elevation.z^2 # Rescale the coordinates redstart$x <- redstart$x / 1000 redstart$y <- redstart$y / 1000 m0 <- glm(counts ~ elevation.z + elevation.z2, data = redstart, family = poisson) beta.starting <- coefficients(m0) beta.tuning <- t(chol(vcov(m0))) n.batch <- 1000 batch.length <- 100 n.samples <- n.batch*batch.length # Fit the spatial model mod <- spGLM(counts ~ elevation.z + elevation.z2, data = redstart, family = "poisson", coords = as.matrix(redstart[,c("x", "y")]), starting = list("beta" = beta.starting, "phi" = 0.5, "sigma.sq" = 1, "w" = 0), tuning = list("beta" = beta.tuning, "phi" = 0.5, "sigma.sq" = 0.1, "w" = 0.01), priors = list("beta.Flat", "phi.Unif" = c(0.01, 50), "sigma.sq.IG" = c(2, 1)), amcmc = list("n.batch" = n.batch, "batch.length" = batch.length, "accept.rate" = 0.43), cov.model = "exponential", verbose = TRUE, n.report = 10) burn.in <- 0.5*n.samples plot(window(mod$p.beta.theta.samples, start = burn.in)) summary(window(mod$p.beta.theta.samples, start = burn.in)) library(fields) # Generate the knots locations to cover the whole landscape knotlocs <- cover.design(redstart[,c("x", "y")], 30) knotlocs <- cbind(knotlocs$design[,1], knotlocs$design[,2]) plot(knotlocs, pch = 16, col = "blue", cex = 2) points(redstart[,c("x", "y")], pch = 3, col = "red") jpeg(filename = "../figures/Figure13-4_knots.jpg", width = 6000, height = 6000, units = "px", res=1200) plot(knotlocs, pch = 16, col = "blue", cex = 2, las = 1, xlab = "X-coordinates", ylab = "Y-coordinates") points(redstart[,c("x", "y")], pch = 3, col = "orange") dev.off() mod.pp <- spGLM(counts ~ elevation + elevation2, data = redstart, family = "poisson", coords = as.matrix(redstart[,c("x", "y")]), knots = as.matrix(knotlocs), starting = list("beta" = beta.starting, "phi" = 0.5, "sigma.sq" = 1, "w" = 0), tuning = list("beta" = beta.tuning, "phi" = 0.5, "sigma.sq" = 0.1, "w" = 0.01), priors = list("beta.Flat", "phi.Unif" = c(0.01, 20), "sigma.sq.IG" = c(1.5, 0.5)), amcmc = list("n.batch" = n.batch, "batch.length" = batch.length, "accept.rate" = 0.43), cov.model = "exponential", verbose = TRUE, n.report = 10) plot(window(mod.pp$p.beta.theta.samples, start = burn.in)) summary(window(mod.pp$p.beta.theta.samples, start = burn.in)) # Extract median predictions at sample points (fitted values) Xb <- cbind(1, redstart$elevation, redstart$elevation2) %*% t(mod.pp$p.beta.theta.samples[burn.in:n.samples,1:3]) yhat <- exp(Xb + mod.pp$p.w.samples[,burn.in:n.samples]) yhat.median <- apply(yhat, 1, median) # Compute the Pearson residuals resi <- (redstart$counts - yhat.median)/sqrt(yhat.median) ################## # Section 13.3.4 # ################## # MCMC settings ni <- 10000 nt <- 1 nb <- 5000 nc <- 3 # Specify model in BUGS language sink("SpatialExp.txt") cat(" model { # likelihood for (i in 1:n) { y[i] ~ dpois(lambda[i]) log(lambda[i]) <- alpha + beta1*elev[i] + beta2*elev2[i] + W[i] muW[i] <- 0 } # Spatial exponential W[1:n] ~ spatial.exp(muW[], xcoord[], ycoord[], tauSp, phi, 1) # Priors for the fixed effects alpha ~ dnorm(0, 0.01) beta1 ~ dnorm(0, 0.04) beta2 ~ dnorm(0, 0.04) # Priors for the spatial random effect tauSp <- pow(sdSp, -2) sdSp ~ dunif(0, 5) phi ~ dunif(0.01, 5) } ",fill = TRUE) sink() # Bundle data bugs.data <- list(n = nrow(redstart), y = redstart$counts, xcoord = redstart$x, ycoord = redstart$y, elev = redstart$elevation.z, elev2 = redstart$elevation.z^2) # Initial values inits <- function() {list(alpha = runif(1), beta1 = runif(1), beta2 = runif(1), phi = runif(1), W = rep(0, nrow(redstart)))} # Parameters monitored params <- c("alpha", "beta1", "beta2", "phi", "tauSp") library(dclone) # Call OpenBUGS from R (215 mins) cl <- makePSOCKcluster(3) parallel::clusterExport(cl, varlist = "redstart") samples <- bugs.parfit(cl, data = bugs.data, params = params, model = "SpatialExp.txt", inits = inits, n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, DIC = FALSE, program = "openbugs", seed = 1:3, save.history = FALSE) stopCluster(cl) plot(samples) summary(samples) library(R2OpenBUGS) mod <- bugs(bugs.data, params, model = "SpatialExp.txt", inits = inits, n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE)