[forest0]    Poisson-gamma spatial moving average
                  (convolution) model: Distribution of
                  hickory trees in Duke forest


Ickstadt and Wolpert (1998) and Wolpert and Ickstadt (1998) analyse the spatial distribution of hickory trees (a subdominant species) in a 140 metre x 140 metre research plot in Duke Forest, North Carolina, USA, with the aim of investigating the regularity of the distribution of trees: spatial clustering would reveal that the subdominant species is receding from, or encroaching into, a stand recovering from a disturbance, while regularity would suggest a more stable and undisturbed recent history.

The data comprise counts of trees, Y
i , in each of i=1,...,16 30m x 30m quadrats (each having area A i = 900 m 2 ) on a 4x4 grid covering the research plot (excluding a 10m boundary region to minimise edge effects), together with the x and y co-ordinates of the centroid of each plot (relative to the south-west corner of the research plot). Ickstadt and Wolpert (1998) fit a Poisson-gamma spatial moving average model as follows:

          Y i ~ Poisson( m i ) i = 1,...,16
          m i = A i * l i
          l i = S j k ij g j

where g j can be thought of as latent unobserved risk factors associated with locations or sub-regions of the study region indexed by j, and k ij is a kernel matrix of 'weights' which depend on the distance or proximity between latent location j and quadrat i of the study region (see Appendix 2 for further details of this convolution model). Ickstadt and Wolpert (1998) take the locations of the latent g j to be the same as the partition of the study region used for the observed data i.e. j = 1,....16 with latent sub-region j being the same as quadrat i of the study region. Note that it is not necessary for the latent partition to be the same as the observed partition - see Hudersfield . The g j are assigned independent gamma prior distributions

      
   g j ~ Gamma( a , b ) j = 1,....,16. 

Ickstadt and Wolpert (1998) consider two alternative kernel functions: an adjacency-based kernel and a distance-based kernel. The adjacency based kernel is defined as:

         k
ij = 1 if i = j
         k
ij = exp( q 2 )/n i if j is a nearest neighbour of i (only north-south and east-west
         neighbours considered) and n
i is the number of neighbours of area i
         k
ij = 0 otherwise

The distance based kernel is defined as:

k
ij = exp( - [d ij / exp( q 2 )] 2 ) where d ij is the distance between the centroid of areas i and j

For both kernels, the parameter exp(
q 2 ) reflects the degree of spatial dependence or clustering in the data, with large values of q 2 suggesting spatial correlation and irregular distribution of hickory trees.

Ickstadt and Wolpert (1998) elicit expert prior opinion concerning the unknown hyperparameters
q 0 = log( a ), q 1 = - log( b ) and q 2 and translate this information into Normal prior distributions for q 0 , q 1 and q 2 .

This model may be implemented in WinBUGS 1.4 using the
pois.conv distribution. The code is given below.

Model

          model {

          # likelihood
              for (i in 1:N) {
                Y[i] ~ dpois.conv(mu[i,])
                for (j in 1:J) {
                   mu[i, j] <- A[i] * lambda[i, j]
                   lambda[i, j] <- k[i, j] * gamma[j]
                }
             }

          # priors (see Ickstadt and Wolpert (1998) for details of prior elicitation)
             for (j in 1:J) {
                gamma[j] ~ dgamma(alpha, beta)
             }
             alpha <- exp(theta0)
             beta <- exp(-theta1)
             
             theta0 ~ dnorm(-0.383, 1)
             theta1 ~ dnorm(-5.190, 1)
             theta2 ~ dnorm(-1.775, 1)            # prior on theta2 for adjacency-based kernel
             #   theta2 ~ dnorm(1.280, 1)         # prior on theta2 for distance-based kernel


          # compute adjacency-based kernel
          # Note N = J in this example (necessary for adjacency-based kernel)
              for (i in 1:N) {
                k[i, i] <- 1
                for (j in 1:J) {
          # distance between areas i and j
                   d[i, j] <- sqrt((x[i] - x[j])*(x[i] - x[j]) + (y[i] - y[j])*(y[i] - y[j]))
          # (only needed to help compute nearest neighbour weights;
          # alternatively, read matrix w from file)
                   w[i, j] <- step(30.1 - d[i, j])         # nearest neighbour weights:
          # areas are 30 sq m, so any pair of areas with centroids > 30m apart are not
         # nearest neighbours (0.1 added to account for numeric imprecision in d)
                }
                for (j in (i+1):J) {
                   k[i, j] <- w[i, j] * exp(theta2) / (sum(w[i,]) - 1)
                   k[j, i] <- w[j, i] * exp(theta2) / (sum(w[j,]) - 1)
                }
             }

          # alternatively, compute distance-based kernel
          #   for (i in 1:N) {
          #      k[i, i] <- 1
          #      for (j in 1:J) {
          # distance between areas i and j
          #         d[i, j] <- sqrt((x[i] - x[j])*(x[i] - x[j]) + (y[i] - y[j])*(y[i] - y[j]))
          #      }
          #      for (j in (i+1):J) {
          #         k[i, j] <- exp(-pow(d[i, j]/exp(theta2), 2))
          #         k[j, i] <- exp(-pow(d[j, i]/exp(theta2), 2))
          #      }
          #   }

          # summary quantities for posterior inference
             for (i in 1:N) {
          # smoothed density of hickory trees (per sq metre) in area i
                density[i] <- sum(lambda[i, ])
          # observed density of hickory trees (per sq metre) in area i
                obs.density[i] <- Y[i]/A[i]
             }
          # large values indicate strong spatial dependence;
          # spatial.effect -> 0 indicates spatial independence
             spatial.effect <- exp(theta2)
          }

Data click here to open data

Inits click here to open initial values

Results

Assuming adjacency-based kernel (equivalent to Ickstadt and Wolpert's model M S ):    

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   spatial.effect   0.2856   0.1793   0.009916   0.04003   0.2588   0.7088   10001   10000
   theta0   -0.6121   0.4562   0.02313   -1.519   -0.6084   0.228   10001   10000
   theta1   -5.111   0.5088   0.0179   -6.046   -5.133   -4.05   10001   10000
   theta2   -1.487   0.7497   0.04143   -3.218   -1.352   -0.3441   10001   10000

Ickstadt and Wolpert report a posterior mean (sd) of 0.281 (0.156) for the spatial effect, exp(
q 2 ), from their analysis using a 4x4 partition of the study region (Table 1.3).


Assuming distance-based kernel (equivalent to Ickstadt and Wolpert's model M D ):
   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   spatial.effect   8.042   7.474   0.5676   0.5371   4.684   23.9   20001   20000
   theta0   -0.297   0.5198   0.02996   -1.533   -0.2463   0.5766   20001   20000
   theta1   -5.276   0.522   0.01989   -6.197   -5.308   -4.157   20001   20000
   theta2   1.563   1.108   0.06682   -0.6216   1.544   3.174   20001   20000
   
   
Ickstadt and Wolpert report a posterior mean (sd) of 7.449 (6.764) for the spatial effect, exp(
q 2 ), from their analysis using a 4x4 partition of the study region (Table 1.3), and posterior means of -0.02, -5.28 and 1.54 respectively for theta0, theta1 and theta2. Ickstadt and Wolpert noted that the posterior for the spatial effect was multi-modal, as seen in the density plot below:

   
   
MISSING FIGURE
   

Note however that convergence of this distance-based model is not very good, and the results reported above are not particularly stable (the adjacency-based model appears to converge much better). This may be partly due to the rather strong prior information about theta2 for the distance-based model which, as noted by Ickstadt and Wolpert, appears to conflict with the data.