Scaling iCAR

[John Clare]

Hi,

I had a question regarding scaling iCAR (car_normal) effects in nimble. This is something I'd like to implement in a couple of different settings: a BYM model where it would be helpful to parameterize the the marginal variance of spatial and non-spatial RE to be roughly equal/equivalent to one divided by car_normal's tau argument (or at least have consistent priors for each of these), and a latent factor model where it would be preferred for the marginal variance of the car_normal effects to be roughly 1. 

Inla approaches this by rescaling the precision matrix--there are function (at bottom) that work to scale things by the geometric mean of the matrix diagonal following a sum to zero constraint so that the iCAR values are approximately ~N(mu [or 0], 1/tau). My understanding is that it would require custom functions to use a scaled precision matrix like this directly in nimble, and that nimble's iCAR sum to zero is a little more deterministic (subtracting the mean after  the fact). What I''ve seen done in other software applications (https://doi.org/10.1016/j.sste.2019.100301) is to derive the scaling factor in inla, pass it as a value or constant to the software, and rescale centered iCAR effects after the fact. I'm curious whether something like the stuff at the bottom here would work? (I think, at best, it might work only if there are not multiple disconnected regions).

I recognize this is maybe a little outside of the scope of the list, so particular thanks for any thoughts on this.

Thanks,

John

###Let Q be a precision matrix, N be number of areas, 

##adj.matrix be full adjacency matrix, or sparseMatrix

###in R 

###precision matrix for iCAR

Q=  INLA::Diagonal(N, rowSums(adj.matrix)) - adj.matrix

###function to scale precision matrix 
inla.scale.model.new = function(Q, constr = NULL, eps = sqrt(.Machine$double.eps))
{
  ## returns scaled Q,  the scaled marginal variances

  ## and the scaling factor--only modification is returning 'fac'
 
  marg.var = rep(0, nrow(Q))
  Q = inla.as.sparse(Q)
  g = inla.read.graph(Q)
  for(k in seq_len(g$cc$n)) {
    i = g$cc$nodes[[k]]
    n = length(i)
    QQ = Q[i, i, drop=FALSE]
    if (n == 1) {
      QQ[1, 1] = 1
      marg.var[i] = 1
    } else {
      cconstr = constr
      if (!is.null(constr)) {
        ## the GMRFLib will automatically drop duplicated constraints; how convenient...
        cconstr$A = constr$A[, i, drop = FALSE]
        eeps = eps
      } else {
        eeps = 0
      }
      res = inla.qinv(QQ + Diagonal(n) * max(diag(QQ)) * eeps, constr = cconstr) 
      fac = exp(mean(log(diag(res)))) ###this is the scaling factor--geometric mean...
      QQ = fac * QQ
      marg.var[i] = diag(res)/fac
    }
    Q[i, i] = QQ
  }
  return (list(Q=Q, var = marg.var, fac=fac))
}

scale_factor<-inla.scale.model.new(Q, constr=list(A=matrix(1, 1, N), e=0))$fac

###in nimbleCode

rho[1:X]~dcar_normal(....tau=1...., zero_mean=1)

rhostar[1:X] <- rho[1:X]/sqrt(scale_factor)

## in this case, want marginal variance of rhostar to be about 1

[Chris Paciorek]

Hi John,

A few thoughts here, though I'm not particularly familiar with the scaling approach from a statistical perspective.

1) It sounds like scale_factor is a constant in terms of the model, so there's no problem with calculating it separately however you want and then using it in the model as your last few lines of code show. Is that "all" that you are asking or am I missing something?

2) If using a. proper CAR, you could in principle set up the precision matrix manually in nimble and run MCMC but it will be inefficient because it wouldn't take advantage of the sparsity of the precision matrix. You couldn't do that with an improper CAR because our linear algebra calculations are not set up to deal with the fact that the matrix is singular. Not sure this is relevant for you but wanted to mention it.

3) On the off-chance it's useful, note that nimbleRcall() will allow you to call out to R. So, for example if you needed to repeatedly calculate scale_factor as part of the model calculations, that's a way you could do it. Similarly, if one needed sparse matrix manipulations you could do that within R via nimbleRcall. 

Not sure if I'm fully following what you're asking. Feel free to iterate here.

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/125edaca-94b9-4476-ac59-b170306e5b8dn%40googlegroups.com.

[John Clare]

Hi Chris,

Thanks a bunch. Mostly was interested in getting a better handle on whether I was missing potential options with nimble, and trying to understand if the hacky solution I proposed made any sense with respect to how nimble works through iCAR stuff.

Yeah, scale_factor is just a constant that is not updated. I think INLA works with the scaled precision matrix directly and employs a sum-to-zero constraint I haven't looked closely at. Applications in Stan I've seen have used a pairwise difference formulation of the joint distribution and a soft sum-to-zero constraint (e.g., mean(rho)~N(0, sd=some small number) but rescale rho (the iCAR vector) after the fact by the scaling factor. The scaling does not make sense if the vector (or each connected component) is not roughly centered on zero, so partially just thinking out loud about whether this makes sense given nimble's order or operations--I think sampling the conditional distribution for rho[i], making density calculations, and then subtracting the mean post-hoc? I think further subsequent scaling by a constant would work ok provided the graph consists of one connected component. Am I reading correctly that if there were two connected components (e.g., two islands or something), then zero_mean = 1 makes little sense (would have to center each component separately after the density statement)?

Thanks,

John

---

From: christophe...@gmail.com <christophe...@gmail.com>
Sent: Friday, February 25, 2022 6:59 PM
To: John Clare <jcl...@wisc.edu>
Cc: nimble-users <nimble...@googlegroups.com>
Subject: Re: Scaling iCAR

 

[Chris Paciorek]

hi John,

Yes, I think that's right that the zero_mean=1 makes little sense if one has more than one connected component. As far as order of operations, if zero_mean = 1, the centering happens after having updated each of the individual components and before going on to sample other parts of the model (i.e., at the end of running sampler_CAR_normal).

-chris