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))