Dear Stan users/developers.
I am trying to write a code for a multivariate, multilevel regression. To start, I am simulating data from only two dependent (y1, y2) and two independent variables (x1, x2) for 40 participants, with 20 measurements each I attach herewith the Stan code and the code for simulation of data. I’m having divergent transitions, although I’ve used multinormal_cholesky for the multivariate outcome variable. Note also that to simplify things I assume that the coefficients for each outcome are independent, such that multivariate normal distributions are not used for level-1 parameters. Any help with eliminating the divergent transitions would be appreciated!
Thanks,
Isaac
data{
int<lower=0> Nvars; //2 independent and 2 dependent variables in our simulation
int<lower=0> N; //Total number of data points
int<lower=0> Nsubj; //Number of subjects
int subj[N];
matrix[N,Nvars] x;
vector[Nvars] y [N];
}
parameters{
cholesky_factor_corr[Nvars] L_Omega; //For the correlations between outcome variables
vector<lower=0>[Nvars] L_sigma; //Residual SD of outcome variables
matrix[Nvars,Nsubj] beta0; //Intercepts at subject level, separately for each outcome variable
matrix[Nvars,Nsubj] beta [Nvars]; //Slopes at the subject level, for each outcome and predictor variables (i.e. same variable, with the former being at time t and the latter at time t-1)
vector[Nvars] Beta0; //Level-two intercepts
vector[Nvars] Beta [Nvars]; //Level-two slopes
vector<lower=0>[Nvars] Tau0; //Random effect of intercepts
matrix<lower=0> [Nvars,Nvars] Tau; //Random effect of slopes
}
transformed parameters{
vector[Nvars] mu[N];
matrix [Nvars,Nvars] Sigma;
Sigma = diag_pre_multiply(L_sigma, L_Omega);
for (i in 1:N){
for (v in 1:Nvars){
mu[i,v]=beta0[v,subj[i]]+dot_product(to_vector(beta[1:Nvars,v,subj[i]]),x[i,1:Nvars]);
}
}
}
model{
//Priors
L_Omega~lkj_corr_cholesky(1);
L_sigma~cauchy(0,5);
Tau0~cauchy(0,5);
to_vector(Tau)~cauchy(0,5);
to_vector(beta0)~normal(0,5);
for (v in 1:Nvars){to_vector(beta[v,,])~normal(0,5);}
for (i in 1:N){
y[i,1:Nvars]~multi_normal_cholesky(mu[i,1:Nvars],Sigma);
}
for (s in 1:Nsubj){
for (v in 1:Nvars){
beta0[v,s]~normal(Beta0[v],Tau0[v]);
beta[1:Nvars,v,s]~normal(Beta[1:Nvars,v],Tau[1:Nvars,v]);}
}
}
'
And here is the code for the simulations:
sigma=matrix(c(1,0.5,0.3,0.4,
0.5,1,0.6,0.7,
0.3,0.6,1,0,
0.4,0.7,0,1),nrow=4)
colnames(sigma)<-c("y1","y2","x1","x2")
rownames(sigma)<-c("y1","y2","x1","x2")
library(MASS)
Data=mvrnorm(n = 20,mu = c(10,20,0,0),Sigma = sigma,empirical = F)
for (i in 2:40){
Data<-rbind.data.frame(Data,mvrnorm(n = 20,mu = c(10,20,0,0),Sigma = sigma,empirical = F))
}
Data$subj<-rep(1:40,each=20)
#Data<-Data[sample(rownames(Data),900),]
y=Data[,1:2]
x=Data[,3:4]
subj=Data$subj
datalist=list(y=y,
x=x,
N=nrow(x),
Nvars=2,
Nsubj=max(subj),
subj=subj)
