[elevation0]    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[])
            }
         }

Data     (Click to open)

Initial values for exponential model
Inits for chain 1   Inits for chain 2    (Click to open)


Initial values for disc model
Inits for chain 1   Inits for chain 2    (Click to open)

Results for exponential model


         [elevation1]
         
         [elevation2]


Results for disc model

[elevation3]


[elevation4]

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 .