Hierarchical prior (for partial pooling) on correlation matrices?

General
27

Thank you for all your input - Seriously, a life saver.

Also for future me (or anyone who wants to do this simply), I wrote two functions that helps compute the per-group correlation matrices (or choleskies).
It takes a global matrix, a local diff matrix, and alpha.
If the local matrix is smaller than the global matrix, then it uses the inner-most subset from the global matrix.
E.g., if you have a 9x9 global matrix, but a 6x6 local matrix, then it will use global[1:6,1:6]; this means the function requires that the local matrix is literally a top-right-corner subset of the global matrix.

	matrix convex_combine_cholesky(matrix global_chol_cor, matrix local_chol_cor, real alpha){
		int dim = rows(local_chol_cor);
		int global_dim = rows(global_chol_cor);
		matrix[global_dim,global_dim] global_cor = multiply_lower_tri_self_transpose(global_chol_cor);
		matrix[dim,dim] local_cor = multiply_lower_tri_self_transpose(local_chol_cor);
		matrix[dim,dim] global_cor_use;
		matrix[dim,dim] combined_chol_cor;
		
		if(dim < rows(global_cor)){
			global_cor_use = global_cor[1:dim,1:dim];
		} else {
			global_cor_use = global_cor;
		}
		combined_chol_cor = cholesky_decompose((1 - alpha)*global_cor_use + (alpha)*local_cor);
		return(combined_chol_cor);
	}
	
	matrix convex_combine(matrix global_cor, matrix local_cor, real alpha){
		int dim = rows(local_cor);
		matrix[dim,dim] global_cor_use;
		matrix[dim,dim] combined_cor;
		if(dim < rows(global_cor)){
			global_cor_use = global_cor[1:dim,1:dim];
		} else {
			global_cor_use = global_cor;
		}
		combined_cor = (1 - alpha)*global_cor_use + alpha*local_cor;
		return(combined_cor);
	}

To use it, you need to declare a global matrix, which should be the same dimension as the largest group matrix; group-specific ‘diff’ matrices; and alpha.
It looks like this:

parameters {
    cholesky_factor_corr[9] R_global_L;
    cholesky_factor_corr[9] R_diff_L[4];
    real<lower=0,upper=1> alpha;
}
transformed parameters {
    cholesky_factor_corr[9] R_group_L[4];
    for(g in 1:4){
        R_group_L[g] = convex_combine_cholesky(R_global_L,R_diff_L[g],alpha);
    }
}
model {
    alpha ~ beta(2,2) 
    R_global_L ~ lkj_corr_cholesky(3);
    for(g in 1:4){
        R_diff_L[g] ~ lkj_corr_cholesky(1/alpha);
        data[g] ~ multi_normal_cholesky(muVector, R_group_L[g]);
    }

}
5 Likes