I am attempting to deal with a short time series with four different data sets of surveys. Some data sets have repeated measures on the same respondents. I am using a normal regression model but when there are repeated measurements the likelihood is multivariate normal. I have reason to believe from data exploration that the variance is dependent on demographic characteristics. I want to put LKJ priors on the covariance matrices. When I run the stan code I get the following error:
“Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: multi_normal_lpdf: Covariance matrix is not symmetric. Covariance matrix[1,2] = 0, but Covariance matrix[2,1] = 2.9525 (in ‘modelc141bfc6bac_stantestpolicyv2’ at line 101)”
I suspect I have not specified the prior on the covariance correctly but I’m not sure what I’m doing wrong.
data {
int<lower=0> N; // total sample size
int<lower=1> D; // dimension of non-demographic covariates
int<lower=1> K; // number of demographic covariates
int<lower=1> T; //number of time points
//group 1
int<lower=1> N_g1; // sample size group 1
vector[D] x_g1[N_g1]; // covariates group 1
vector[K] cvec_g1[N_g1]; // demographic covariates group 1
real y_g1[N_g1]; // response group 1
//group 2
int<lower=1> N_g2; // sample size group 2
matrix[D,2] x_g2[N_g2]; // covariates group 2
vector[K] cvec_g2[N_g2]; // demographic covariates group 2
vector[2] y_g2[N_g2]; // response group 2
//group 3
int<lower=1> N_g3; // sample size group 3
matrix[D,6] x_g3[N_g3]; // covariates group 3
vector[K] cvec_g3[N_g3]; // demographic covariates group 3
vector[6] y_g3[N_g3]; // response group 3
// group 4
int<lower=1> N_g4; // sample size group 4
vector[D] x_g4[N_g4]; // covariates group 4
vector[K] cvec_g4[N_g4]; // demographic covariates group 4
real y_g4[N_g4]; // response group 4
// vector of demographic group
int<lower=1> C_g1[N_g1];
int<lower=1> C_g2[N_g2];
int<lower=1> C_g3[N_g3];
int<lower=1> C_g4[N_g4];
}
// The parameters accepted by the model. Our model
parameters {
vector[K] alpha; // incenterepts
matrix [D,T] beta; // regresion coefficients
vector<lower=0>[K] sigma_alpha; //variance of demographic coefficients
vector<lower=0>[D] sigma_beta; //variance of non-demographic coefficients
vector<lower=0>[32] sigma_error_s1; // variance of normal for s=1 for each covariance
vector<lower=0>[32] sigma_error_s4;
cholesky_factor_corr[2] group2cor[32]; // covariance of matrix for people observed 2
cholesky_factor_corr[6] group3cor[32]; // covariance of matrix for people observed 6
vector<lower=0>[2] sigmag2[32];
vector<lower=0>[6] sigmag3[32];
real<lower=0> group1var[32];
real<lower=0> group4var[32];
//real mu;
//real<lower=0> sigma;
}
transformed parameters{
//int my_subset[2] = {1,7};
cholesky_factor_cov[2] covg2[32];
cholesky_factor_cov[6] covg3[32];
for(k in 1:32) covg2[k] = diag_pre_multiply(sigmag2[k],group2cor[k]);
for(k in 1:32) covg3[k] = diag_pre_multiply(sigmag3[k],group3cor[k]);
}
// The model to be estimated. We model the output
// 'y' to be normally distributed with mean 'mu'
// and standard deviation 'sigma'.
model {
//y ~ normal(mu, sigma);
for(i in 1:K){
sigma_alpha[i] ~ uniform(0,1);
}
for(i in 1:D){
sigma_beta[i] ~ uniform(0,1);
}
for(i in 1:K){
alpha[i] ~ normal(0, sigma_alpha[i]);
}
for(i in 1:D){
beta[i,:] ~ normal(0, sigma_beta[i]);
}
for(i in 1:32){
group2cor[i] ~ lkj_corr_cholesky(1);
group3cor[i] ~ lkj_corr_cholesky(1);
sigmag2[i] ~ uniform(0,1);
sigmag3[i] ~ uniform(0,1);
group1var[i] ~ uniform(0,1);
group4var[i] ~ uniform(0,1);
}
for(i in 1:N_g1){
int c = C_g1[i];
y_g1[i] ~ normal(dot_product(cvec_g1[i], alpha)+dot_product(x_g1[i],beta[:,1]), group1var[c]);
}
for(i in 1:N_g2){
int c = C_g2[i];
vector[2] mu2;
mu2[1] = dot_product(cvec_g2[i], alpha)+x_g2[i][:,1]'*beta[:, 1];
mu2[2] = dot_product(cvec_g2[i], alpha)+x_g2[i][:,2]'*beta[:, 7];
y_g2[i] ~ multi_normal(mu2, covg2[c]);
}
for(i in 1:N_g3){
int c = C_g3[i];
vector[6] mu3;
mu3[1] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,2]'*beta[:, 2];
mu3[2] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,3]'*beta[:, 3];
mu3[3] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,4]'*beta[:, 4];
mu3[4] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,5]'*beta[:, 5];
mu3[5] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,6]'*beta[:,6];
mu3[6] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,7]'*beta[:,7];
y_g3[i] ~ multi_normal(mu3, covg3[c]);
}
for(i in 1:N_g4){
int c = C_g4[i];
y_g4[i] ~ normal(dot_product(cvec_g4[i], alpha)+dot_product(x_g4[i],beta[:,1]), group1var[c]);
}
}