iCAR (intrinsic CAR) models via precision matrix

[pconn]

Hello, I've been trying to get an iCAR spatial regression model up and running in TMB. For Bayesian analysis, the relevant prior distribution for random effects is \Eta ~ MVN(0,(\tau*R)^-1), where R is a structure matrix as defined e.g. in Rue and Held 2005. Importantly, the precision matrix Q = \tau R is improper and noninvertible. In MCMC estimation this doesn't really matter - with data the full conditional distributions are proper. However, entering in a noninvertible Q matrix (as for intrinsic Gausian Markov random fields) appears to cause errors in TMB (I'm using GMRF(Q)(\Eta) to specify the probability of random fields). For the moment, I'm just adding a smallish value to the diagonal of R which makes Q proper and that seems to work okay, but wondering if there are workarounds here (guessing there are since INLA can handle ICAR specifications). Any suggestions or advice to get around my kludgish solution?

Thanks! Paul Conn

[skaug]

Hi,

The following code shows how it can technically be done in TMB, avoiding a determinant calculation implicit in GMRF(). My suggestion is however to turn your "smallish value" into a parameter (delta) and estimate that in TMB. If data tells that the MLE of delta=0, then I would consider using the ICAR. For delta > 0 you have a proper prior, which is easier conceptually.

Hans

#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_SPARSE_MATRIX(R);

  PARAMETER(tau);                      
  PARAMETER_VECTOR(u);      // Underlying dummy latent vecdtor
  PARAMETER_VECTOR(theta);

  Type f=0.0;

  // ICAR prior
  vector<Type> tmp = R*u;
  f -= Type(-0.5)*(u*tmp).sum();

  vector<Type> eta = tmp/tau;       // Doing the scaling as a post-step. 100% OK!!

  // The likelihood part..................

  return f;
}

[pconn]

Great - thanks a lot, Hans. Regards, Paul

On Thu, May 21, 2015 at 11:01 PM, Hans J. Skaug notifications@github.com wrote:

Hi,

The following code shows how it can technically be done in TMB, avoiding a determinant calculation implicit in GMRF(). My suggestion is however to turn your "smallish value" into a parameter (delta) and estimate that in TMB. If data tells that the MLE of delta=0, then I would consider using the ICAR. For delta > 0 you have a proper prior, which is easier conceptually.

Hans

include template

Type objective_function::operator() () { DATA_SPARSE_MATRIX(R);

PARAMETER(tau); PARAMETER_VECTOR(u); // Underlying dummy latent vecdtor PARAMETER_VECTOR(theta);

Type f=0.0;

// ICAR prior vector tmp = Ru; f -= Type(-0.5)(u*tmp).sum();

vector eta = tmp/tau; // Doing the scaling as a post-step. 100% OK!!

// The likelihood part..................

return f; }

— Reply to this email directly or view it on GitHub https://github.com/kaskr/adcomp/issues/94#issuecomment-104527583.

Paul Conn Research Statistician Polar Ecosystems Program National Marine Mammal Lab NOAA Alaska Fisheries Science Center 7600 Sand Point Way, Building 4 Seattle, WA 98115

tel 206.526.4235