 # 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,