Hello,
I am trying to model some data that has a hierarchical structure. I have data for 2 random variables that is indexed by 43 “country level” groups. I want to fit a multivariate normal model and in particular want to estimate hierarchical structures on my mean vectors and correlation and covariance matrices. Meaning I want 43 mu, Omega and Sigma, one for each group/country as well as an overall mu, Omega and Sigma.
My data is here (mvn_sample_data.csv (78.3 KB)). My Stan code is as follows:
data{
int N; // number of observations
int n_country; // number of countries
int n_random_vars; // number of random variables being included, 2 here
int country_id[N]; // country id vector
vector[n_random_vars] joint_data[N]; //data
}
parameters{
vector[n_random_vars] mu[n_country];
corr_matrix[n_random_vars] Omega [n_country];
vector<lower = 0>[n_random_vars] sigma [n_country];
}
transformed parameters {
cov_matrix[n_random_vars] Sigma[n_country];
for(i in 1:N){
Sigma[country_id[i]] = quad_form_diag(Omega[country_id[i]], sigma[country_id[i]]);
}
}
model{
for(i in 1:N){
joint_data[i] ~ multi_normal(mu[country_id[i]],Sigma[country_id[i]]);
mu[country_id[i]] ~ normal(0,1);
sigma[country_id[i]] ~ cauchy(0, 5);
Omega[country_id[i]] ~ lkj_corr(1);
}
}
The problems that I am facing are as follows:
- My current parametrization is slow to run, how can I take the mu, sigma and Omega out of the for loop as this might speed things up?
- I don’t know how to move past these numeric priors to prior distributions that would yield an overall mu, Omega and Sigma, potentially something like this
model{
for(i in 1:N){
joint_data[i] ~ multi_normal(mu[country_id[i]],Sigma[country_id[i]]);
}
mu ~ multi_normal(mu_prior, tau_mu_prior); //where mu_prior is a vector of zeros and ta_mu_prior is a matrix of 1's
sigma ~ multi_normal(sigma_prior, tau_sigma_prior); //where sigma_prior is a vector of zeros and ta_sigma_prior is a matrix of 1's
Omega ~ lkj_corr(Omega_prior);
}
The above block is just an attempt and does not work. I declare mu_prior, tau_mu_prior, sigma_prior, and tau_sigma_prior in the data block. I additionally have not been able to find ways to impose a hierarchical structure on Omega_prior as it has to be an integer greater than or equal to 1.
My R code to run the model is below:
# read in csv
year_1_df <- read.csv("mvn_sample_data.csv", header = T)
# prepping data for Stan
joint_data <- as.matrix(cbind(year_1_df$lbl,
year_1_df$fd))
N <- dim(joint_data)[1]
n_country <- length(unique(year_1_df$country_id))
n_random_vars <- dim(joint_data)[2]
country_id <- year_1_df$country_id
#fitting the model
fit <- stan(*model name here*,
data = list(joint_data, N, country_id, n_random_vars),
iter = 1000,
chains = 3)
My questions are as follows:
- How do I enable prior distributions on mu, Omega and Sigma in order to impose a hierarchical structure on them (i.e I want 43 mu, Omega and Sigma as well as overall mu, Omega and Sigma)?
- In particular how do I enable priors for Omega?
- Any tips to speed up my code would be useful too.
Any help would be much appreciated.