Multivariate CAR modelling based on paper of 'Hierarchical multivariate CAR models for spatio-temporally correlated survival data'

Hey there,

I am sort of new to stan. I am using stan to fit a Multivariate CAR model (area spatial model for two diseases ), the MCAR can be express a multivariate Gaussian distribution


where n is the number of subareas, p=2 represents that there are 2 diseases. And

R_1 and R_2 are Cholesky factorization of diag(\omega_{i+}) - \alpha_1 W and diag(\omega_{i+}) - \alpha_2 W respectively. And \alpha_1 and \alpha_2 are two different propriety parameters for the two diseases. W is the ajacency matrix and diag(\omega_{i+}) denotes a diagonal matrix with entries equal to the number of neighours for each subarea.

I write a stan model to mimic the process,

data {
   int<lower=0> N;// number of subarea
   int<lower=0> p;// number of disease, in this case p =2 
   int<lower=0> y1[N];  //response of disease 1
   int<lower=0> y2[N]; // response of disease 2
   
   matrix[N, N] D_W; // diagonal matrix with neighbour numbers as entries
   matrix[N, N] W;   // Adjacency matrix. 
   matrix[p,p] LAMBDA0; // hyperparameter for LAMBDA, fixed, fed as data.
}

// The parameters accepted by the model. 
parameters {
    matrix[p,p] LAMBDA;   // inverse covariance matrix among the p=2 disease

    vector[p*N] phi; // Multivariate CAR (MCAR) effect vector for p=2 diseases.
    real<lower=0,upper=1> a1; // propriety parameter for disease 1
    real<lower=0,upper=1> a2; // propriety parameter for disease 2
}

transformed parameters {
  //
  vector[N] zeros;
  cholesky_factor_cov[N] R1; // cholesky factorization for disease 1
  cholesky_factor_cov[N] R2; // cholesky factorization for disease 2
  matrix[N*p, N*p] B_P;   // precision matrix of MCAR
  
  zeros = rep_vector(0, N);
  
  R1 = cholesky_decompose(D_W-a1*W);
  R2 = cholesky_decompose(D_W-a2*W);
  B_P[1:N,1:N] = LAMBDA[1,1]*R1*R1';
  B_P[1:N,(1+N):(N+N)] = LAMBDA[1,2]*R1*R2';
  B_P[(1+N):(N+N),1:N] = LAMBDA[2,1]*R2*R1';
  B_P[(1+N):(N+N),(1+N):(N+N)] = LAMBDA[2,2]*R2*R2';

}


// The model to be estimated. We model the output
// 'y' to be normally distributed with mean 'mu'
// and standard deviation 'sigma'.
model {
  for ( i in 1:N)
  {
    y1[i] ~ poisson(1000*exp(phi[i]));
    y2[i] ~ poisson(1000*exp(phi[i+N]));
  }
  // MCAR prior for spatial random effects accross two disease.
  phi ~ multi_normal_prec(zeros, B_P);
  // prior for LAMBDA
  LAMBDA ~ wishart(2, LAMBDA0 );
  // prior for propriety parameters (disease specific, should close to 1.
  // a's=0 means independent random effect across subareas.)
  a1 ~ beta(18,2);
  a2 ~ beta(18,2);
}

in the transformed parameter part, compute B_{P} using
cholesky_decompose function. However, I encounter an error saying " Exception: multi_normal_prec_lpdf: the Precision matrix is not symmetric. Precision matrix[1,5] = -3.0333, but Precision matrix[5,1] = 0.444777 ". I am doubting that the issue might arise in the computing of B_{P}. I am sure It’s correct (symmetric precision matrix) when I generated data in R as follows

library(rstan)  
library(nimble)
library(R2jags)
library(runjags)
library(mcmcplots)
library(ngspatial)
options(mc.cores = parallel::detectCores())  
rstan_options(auto_write = TRUE)
#source("Functions.R")

dat <- list()

###############################################generate the spatial structure and simulated data
###############################################
###############################################
###############################################
sim<-list()
#sim$N<-100 # Number of observations.
sim$N <- 4
# Set-up adjacency structure for ICAR model.
sim$A<-adjacency.matrix(sqrt(sim$N))
sim$n_adj<-rowSums(sim$A)
sim$D<-diag(sim$n_adj)
sim$adj<-as.carAdjacency(sim$A)$adj



a1 = rbeta(1,18,2)
a2 = rbeta(1,18,2)
a1;a2

QR1 = chol(sim$D-a1*sim$A) 
QR2 = chol(sim$D-a2*sim$A) 

R = Matrix::bdiag(QR1,QR2) %>% as.matrix()

Lam = rWishart(n=1, df=2, Sigma=diag(10,nrow=2))
Lam = Lam[,,1]

######precision matrix for MCAR 
B_P = t(R)%*% kronecker(Lam,diag(1,nrow=sim$N))%*% R


phi0 = MASS::mvrnorm(1,mu=rep(0,2*sim$N),Sigma = solve(B_P))

sim$phi = phi0[1:sim$N]
sim$eta = phi0[(sim$N+1):(sim$N+sim$N)]


B_P_2 = matrix(NA,nrow = 8, ncol = 8)
B_P_2[1:4,1:4] = Lam[1,1]*t(QR1)%*%(QR1)
B_P_2[1:4,5:8] = Lam[1,2]*t(QR1)%*%(QR2)
B_P_2[5:8,1:4] = Lam[2,1]*t(QR2)%*%(QR1)
B_P_2[5:8,5:8] = Lam[2,2]*t(QR2)%*%(QR2)

B_P_2

#############simulate data
txtstring <- '
data{
# Likelihood:
for (i in 1:N){
    y1[i]~dpois(1000*mu1[i])
    
    y2[i]~dpois(1000*mu2[i])
    
    log(mu1[i])<-  phi[i] 
    
    log(mu2[i])<-  eta[i] 
  }
}
model{
fake <- 0
}
'

# parameters are treated as data for the simulation step
data<-list(N=sim$N, phi=sim$phi, eta=sim$eta)

# run jags
out <- run.jags(txtstring, data = data, monitor=c("y1","y2"),sample=1, n.chains=1, summarise=FALSE)

Simulated <- coda::as.mcmc(out)
Simulated


dim(Simulated)

dat = as.vector(Simulated)
dat

##################fit the simulated data
###compile the model
model = stan_model(file='MCAR_Simplified.stan')
### setting initials 
inits_chain1 <- list(phi=rep(0.001,2*sim$N), a1=0.8329573, a2 = 0.9409176, LAMBDA = Lam)
inits_chain2 <- list(phi=rep(-0.001,2*sim$N), a1=0.8329573, a2 = 0.9409176, LAMBDA = Lam)

####fit the full model use stan
MCAR_Simplified_sample = sampling(model, 
                               data=list(N=sim$N,
                                         p=2,
                                         y1 = dat[1:N],
                                         y2=dat[(N+1):(2*N)],
                                         D_W = sim$D,
                                         W = sim$A,
                                         LAMBDA0 = diag(x=10,nrow=2)), 
                               #warmup=1000,
                               iter=1000, 
                               chains = 2,
                               control = list(adapt_delta = 0.8, 
                                              max_treedepth = 10))

Can anyone give me some suggestion? Thanks very much!

Writing fast so let me know if this doesn’t make sense

1 Like

Hi Daniel,
I followed your advice, but the same problem existed.