What do you guys think of the idea whereby no multivariate normal is modelled, but the target is incremented by the sampling distribution expectation of each bivariate correlation separately? That way, we can model each bivariate correlation hierarchically separately in a relatively straight-forward manner. Here’s a rough-but-compiles implementation of what I mean (the one thing I didn’t have time to work out is how to add my toggle trick for centered/non-centered parameterization, so it’s currently center-parameterized for all hierarchical quantities):
data{
// num_subj: number of subj
int<lower=1> num_subj ;
// num_var: number of variables modelled as intercorrelated
int<lower=2> num_var ;
// num_obs: number of observations per variable and subject
int<lower=1> num_obs ;
// subj_obs: array of matrices, one per subject
matrix[num_obs,num_var] subj_obs[num_subj] ;
}
transformed data{
// num_cors: number of correlations implied by the number of variables
int num_cors = ((num_var-1)*num_var)/2;
// corz_se: standard error for correlations on Fisher's Z scale implied by the number of observations
real corz_se = 1/sqrt(num_obs-3) ;
// sqrt_num_obs: pre-computing this for use in likelihood
real sqrt_num_obs = sqrt(num_obs) ;
// sd_gamma_shape: pre-computing this for use in likelihood
real sd_gamma_shape = (num_obs - 1) / 2 ;
// compute the empirical means, SDs, variances, and correlations on Fisher's Z scale
vector[num_var] obs_means[num_subj] ;
vector[num_var] obs_sds[num_subj] ;
vector[num_var] obs_variances[num_subj] ;
vector[num_cors] obs_z[num_subj] ;
for(i_subj in 1:num_subj){
matrix[num_obs,num_var] this_subj_obs = subj_obs[i_subj] ; //temp variable
for(i_var in 1:num_var){
obs_means[i_subj][i_var] = mean(this_subj_obs[,i_var]);
obs_sds[i_subj][i_var] = sd(this_subj_obs[,i_var]);
obs_variances[i_subj][i_var] = sd(this_subj_obs[,i_var])^2;
}
int k = 0 ; //temp variable
vector[num_cors] obs_r ; //temp variable
// looping over lower-tri to compute correlation between each column in this_subj_obs
for(this_var in 1:(num_var-1)){
for(this_other_var in (this_var+1):num_var){
k += 1 ;
obs_r[k] = (
sum(
( this_subj_obs[,this_var] - obs_means[i_subj][this_var] )
.* ( this_subj_obs[,this_other_var] - obs_means[i_subj][this_other_var] )
)
) / (
(num_obs-1) * obs_sds[i_subj][this_var] * obs_sds[i_subj][this_other_var]
) ;
}
}
obs_z[i_subj] = atanh(obs_r) ;
}
}
parameters{
// mean_mean: mean (across subjects) mean
vector[num_cors] mean_mean ;
// mean_sd: sd (across subjects) mean
vector<lower=0>[num_cors] mean_sd ;
// mean_subj: by-subject mean
vector[num_cors] mean_subj[num_subj] ;
// logsd_mean: mean (across subjects) log-sd
vector[num_cors] logsd_mean ;
// logsd_sd: sd (across subjects) log-sd
vector<lower=0>[num_cors] logsd_sd ;
// sd_subj: by-subject sd
vector<lower=0>[num_cors] sd_subj[num_subj] ;
// corz_mean: mean (across subjects) correlations on Fisher's Z scale
vector[num_cors] corz_mean ;
// corz_sd: sd (across subjects) correlations on Fisher's Z scale
vector<lower=0>[num_cors] corz_sd ;
// corz_subj: by-subject correlations on Fisher's Z scale
vector[num_cors] corz_subj[num_subj] ;
}
model{
////
// Priors
////
// Peaked-at-zero normal(0,1) prior on mean of means
mean_mean ~ std_normal() ;
// Peaked-at-zero normal(0,1) prior on sd of means
mean_sd ~ std_normal() ;
// Peaked-at-zero normal(0,1) prior on mean of logsds
logsd_mean ~ std_normal() ;
// Peaked-at-zero normal(0,1) prior on sd of logsds
logsd_sd ~ std_normal() ;
// normal(0,.74) priors on the mean correlations
// yields flat prior on abs(cor)<.7, and diminishing to zero by abs(cor)==1
// view in R via: hist(tanh(abs(rnorm(1e7)*.74)),br=1e4,xlim=c(0,1))
corz_mean ~ normal(0,.74) ;
// Peaked-at-zero normal(0,1) prior on variability of correlations
// (would want to do PPC for this!)
corz_sd ~ std_normal() ;
// hierarchical structure:
for(i_subj in 1:num_subj){
mean_subj[i_subj] ~ normal(mean_mean,mean_sd) ;
sd_subj[i_subj] ~ lognormal(logsd_mean,logsd_sd) ;
corz_subj[i_subj] ~ normal(corz_mean,corz_sd) ;
}
////
// Likelihood
////
for(i_subj in 1:num_subj){
obs_means[i_subj] ~ normal( mean_subj[i_subj] , sd_subj[i_subj]/sqrt_num_obs ) ;
obs_variances[i_subj] ~ gamma( sd_gamma_shape , sd_gamma_shape ./ pow(sd_subj[i_subj],2) ) ;
obs_z[i_subj] ~ normal(corz_subj[i_subj],corz_se) ;
}
}
generated quantities{
// cor_mean: the upper-tri of the group-level-mean correlation matrix, flattened to a vector for efficient storage
vector[num_cors] cor_mean = tanh(corz_mean) ;
}
Now, that’s encoding a model with independent hierarchies for everything, but it’d also be possible be possible to combine things with multivariate normal hierarchical structures, from just within-quantitiy-type (ex. currently the subject-level mean on each var is modelled as independent from their mean on the other vars, but the set of subject-level means could be treated as mvn) to doing everything as one big mvn (capturing things like subjects with a high mean on var1 also tend to have a high correlation between var2 and var3).