Bayesian kriging and spatial prediction:Surface elevation

The data for this example are taken from Diggle and Riberio (2000) (and came originally from Davis (1973)). The data file contains the variables height, x and y, giving surface elevation at each of 52 locations (x, y) within a 310-foot square. The unit of distance is 50 feet, whereas one unit of height represents 10 feet of elevation. A Gaussian kriging model can be fitted to these data in OpenBUGS using either the spatial.exp or spatial.disc distributions. The data file also contains a set of locations x.pred and y.pred representing a 15 x 15 grid of points at which we wish to predict surface elevation. Predictions can be obtained using either the spatial.pred or spatial.unipred predictive distributions in OpenBUGS

Model
model {

# Spatially structured multivariate normal likelihood
      # exponential correlation function
   height[1:N] ~ spatial.exp(mu[], x[], y[], tau, phi, kappa)
      # disc correlation function
#   height[1:N] ~ spatial.disc(mu[], x[], y[], tau, alpha)

for(i in 1:N) {
   mu[i] <- beta
}

# Priors
beta ~ dflat()
tau ~ dgamma(0.001, 0.001)
sigma2 <- 1/tau

# priors for spatial.exp parameters
      # prior range for correlation at min distance (0.2 x 50 ft) is 0.02 to 0.99
   phi ~ dunif(0.05, 20)
# prior range for correlation at max distance (8.3 x 50 ft) is 0 to 0.66
   kappa ~ dunif(0.05,1.95)

# priors for spatial.disc parameter
      # prior range for correlation at min distance (0.2 x 50 ft) is 0.07 to 0.96
#   alpha ~ dunif(0.25, 48)
# prior range for correlation at max distance (8.3 x 50 ft) is 0 to 0.63

# Spatial prediction

# Single site prediction
for(j in 1:M) {
   height.pred[j] ~ spatial.unipred(beta, x.pred[j], y.pred[j], height[])
}

# Only use joint prediction for small subset of points, due to length of time it takes to run
for(j in 1:10) { mu.pred[j] <- beta }
   height.pred.multi[1:10] ~ spatial.pred(mu.pred[], x.pred[1:10], y.pred[1:10], height[])
            }
         }


Initial values for exponential model

   

Initial values for disc model

   

The GeoBUGS map tool can be used to produce an approximate map of the posterior mean predicted surface elevation. It is not possible to map point locations using GeoBUGS. However, a map file called elevation is already loaded in the GeoBUGS Map Tool --- this is a 15 x 15 grid with the centroids of the grid cell corresponding to the locations of each of the prediction points. You can use this to produce a map of the posterior mean (or other posterior summaries) of the vector of predicted values height.pred.

Results for exponential model

         
         
         

Results for disc model