geoNimble

Here's the geoBUGS manual. Core features include:

• manipulation of adjacency matrices/data structures for CAR style models
• import of shapefiles/polygons
• CAR models:
  • car.normal (ICAR)
  • car.proper (CAR) models
• Gaussian process (GP) spatial models
• specialized models: multivariate CAR and Poisson-based model

I propose that we focus, at least initially, on handling basic CAR/ICAR and Gaussian process models with, as needed, basic handling of adjacency information.

CAR models are tricky because they involve cycles in the graph. Some things we probably want to do

• allow a user to specify a full (improper) precision and handle that (with our new pivoted (sparse) Cholesky) with multivariate MCMC samplers
  • we probably want to provide wrapping so a user can provide adjacency info and we generate the (sparse) precision
  • even without sparse matrices, if we allow a scaled precision where the matrix is fixed, we should only need to do the Cholesky once
• special vectorized scalar samplers on the elements of the spatial field
  • for the latter, we'll need to generate both conjugate and nonconjugate versions that use the adjacency info to figure out dependencies - I think this will simply involve new samplers that are designed for our new dcar.normal and dcar.proper (or whatever we name them). We also need to use getDependencies in such a way that we can get individual dependencies of individual elements of the spatial field - getDependencies in this context is subject of a Github issue.
  • another possibility is that we explicitly create a vector of scalar calculate nodeFunctions for the elements in the spatial field; not sure if this would be useful to have rather than embedding everything in the sampler functions. This would allow us to use all of our current scalar and blocked samplers, but not clear how this would be implemented.

GP models shouldn't take much work - just some additional wrapping of functionality we already have.

Questions for Andrew Thomas / MRC BUGS Team

• Density evaluation for (improper) car.normal distribution?
• Exactly how winBUGS imposes the "zero mean" constraint for car.normal() distribution?
• What features from geoBUGS were the most widely used?
• More detail about sampler assignments:
  • Truly just RWMH for continuous nodes, slice for discrete nodes?
• Implementation for conjugate samplers:
  • Using "two-point method"?
  • Efficient handling of dmnorm-dmnorm conjugacy?
• Said they have AD capability, but aren't using it extensively. What usages have they implemented?

Design decisions for dcar_normal() sampling in NIMBLE:

1. When a component of the dcar_normal process has neighbors, everything is well-defined. We apply one of:

   • a conjugate sampler (in the case of only dnorm dependents),
   • a scalar RW sampler (in the case of non-conjugate dependents), or otherwise,
   • a posterior predictive sampler (in the case of no dependents).

2. When a component has 0 neighbors, i.e., itself makes up it's own island consisting of one region (itself), then the conditional prior density is entirely uninformative, and can be thought of as either dflat(), or dnorm(mean=anything, sd=infinity). Thus, the conditional posterior is entirely determined by the dependent nodes (if any). In these cases, BUGS (arbitrarily) sets and fixes the value of these 0-neighbor island components to 0. I don't agree this this decision. Instead, we proceed as:

   • (in the conjugate case, of only dnorm dependents) sampling the component according to the correct conjugate posterior, using the conjugate dependents.
   • (in the case of non-conjugate dependents) sampling the component according to RW-MH, with posterior determined by the dependents.
   • (in the case of no dependents) the posterior is entirely un-defined. In this case, we literally do nothing (simply, return()). Do not change the value of the component, ever. Also, by virture of having no neighbors, this component does not contribute anything to the joint density of the CAR node, hence it's value (even if it's NA) does not matter, and will never be changed.

3. We also handle initial NA values differently from WinBUGS. WinBUGS mandates that all 0-neighbor components begin with initial value NA, then fixes these components at 0. We have different behaviour. We allow 0-neighbor components to have initial value NA, but do not require it. If they have initial value NA, they are handled as:

   • (in the conjugate case, of only dnorm dependents) the conditional posterior is independent of the initial value of the component, so a sample from the conjugate posterior over-writes the NA initial value on the first MCMC iteration.
   • (in the case of non-conjugate dependents) the RW sampling behaviour requires a valid initial value, from which to begin the RW. Even though the conditional posterior is independent of the current value, we still need an initial value just as a starting position for the RW sampling (e.g., first proposed jump needs to be centered around some non-NA value. If an NA is detected, it is over-written with 0 on the first MCMC iteration, and the RW sampling begins with its first proposal in the vicinity of 0, and proceeds from there.
   • (in the case of no dependents) the NA value is initially over-written with a 0, which is necessary for the zero-mean constraint to be calculated (since this calculation now includes 0-neighbor islands). This value of 0 will never be changed by the sampling of this component, but will end up being changed when the zero-mean constraint is imposed (on each MCMC iteration).

4. Regarding calculation of the "number of islands" for an adjacency structure, which is the number of constraints that must be imposed on the improper density, we follow WinBUGS. Imagine a undirected graph, with vertices being the components of the CAR distribution, and connecting edges representing the neighboring relationships between components. Count the total number of disjoint groupings of vertices. We call this c, the number of constraints. Clearly, any components with 0-neighbors constitute their own, small, island. But likewise, a group of 2 components whose only neighbors are one another, constitute another island. c is used in the conditional density calculation of the CAR node, as:

p(tau | x1,x2,...,xN) ~ tau^((N-c)/2) * exp(-tau/2 * sum_{i,j} w_ij * (xi-xj)^2),

where the sum is taken over all pairs of neighboring components xi, xj. This is the approach and calculation used by WinBUGS.

Furthermore, the value of c can instead be provided as a parameter to the dcar_normal() distribution in BUGS code, e.g. x[1:N] ~ dcar_normal(..., c = 2). In this case, the user-provided value will always be used for density evaluations, and the nimble calculation of the number of islands never takes place.

5. With regards to imposing a zero-mean constraint on the CAR process values:

   • WinBUGS imposes a zero-mean constraint on all components of the CAR process (including those with 0-neighbors, since they're fixed to equal 0). In WinBUGS, this constraint is not applied separately to each disjoint group of regions comprising its own island, but rather to all components of the CAR distribution, simultaneously. This approach is consistent with having a single dflat() intercept term used for all regions, as WinBUGS requires.
   • By default, we do not impose any such zero-mean constraint, and thus the CAR intercept (or intercepts, in the case of multiple disjoint islands) are implicit within the CAR distribution.
   • However, by specifying the distribution parameter dcar_normal(..., zero_mean = 1), then MCMC sampling will impose a zero-mean constraint, by subtracting away the mean of all CAR process components, at the conclusion of the component-wise sampling that takes place on each MCMC iteration. The calculation of this mean, and subtraction, now includes all 0-neighbor components (which themselves undergo sampling). In this case of imposing the zero-mean constraint, users are advised to include a separate intercept term to complement the (zero-mean) CAR process.

6. The dcar_normal distribution accepts vector parameters num to indicate the number of neighbors of each region (of length N), adj to provide the indices of the (necessarily symmetric) neighbors of each region (of length sum(num)), and weights to provide the unnormalised weight (also necessarily symmetric) for each neighboring relationship (same length as adj). If weights is ommitted, it assumes a default value of all 1's. This same information can also be specified in two other formats, and converted into the num, adj, weights format:

   • neighborList and weightList. neighborList is a list of length N, where each element is a numeric vector giving the indices of the neighbors of that region, and weightList is a list of length N, where each element is a numeric vector giving the weights associated with the neighboring relationships of that region.
   • weightMatrix, a symmetric N x N adjacency matrix, with diagonal elements equal to 0, and off-diagonal elements giving the weight of neighboring relationships.

These alternatives can be converted into the num, adj, weights format using as.carAdjacency(), as as.carAdjacency(neighborList, weightList) or as.carAdjacency(weightMatrix).