Coding a new matrix/vector whose elements are a mean of elements from other matrices/vectors

Dear STAN experts,

I’m trying to run a hierarchical model of a multinormal regression. Briefly, this model assesses the influence of external feedback on pre- to post- feedback changes in particpants’ ratings of a self-performance. Moreover, the model differentiates between cases wherein the feedback was either negative/positive/equal to participants’ own ratings via conditioning. The conditioning for the “equal” cases involves a covariance matrix, whose elements are supposed to be the mean of parameters of two other covariance matrices (i.e. for negative or positive feedback). I’m not sure that this averaged covariance matrix is produced correctly, as this version of the model yields many divergent transitions. However, when this conditioning is eliminated, and the equal cases are merged either into the negative or positive feedback conditioning, the model works fine. Hence, I suspect that I’m not averaging vectors and matrices correctly, and I couldn’t really find in the different manuals how this is supposed to be done. Currently the conditioning is coded as follows (the full model is attached separately):

int <lower=1> Nsubj ; // number of participants
int <lower=1> Ntotal ; // number of questions
int <lower=1> nYlevels ; // number of responses on scale (0-10 scale, so 11 responses)
int <lower=1> s[Ntotal] ; // index of subject at a given question
int <lower=1> k ; // num. of predictors
int yPre[Ntotal] ; // subjects' self-rating before feedback
int yJud[Ntotal] ; // feedback from external judges
int yPost[Ntotal] ; // subjects' self-rating after feedback

model {

mu_tau_pos ~ cauchy(0, 2.5) ;
omega_corr_pos ~ gamma(2,1) ;
mu_tau_neg ~ cauchy(0, 2.5) ;
omega_corr_neg ~ gamma(2,1) ;

for (i in 1:Nsubj) {

tau_pos[i,] ~ gamma((mu_tau_pos/sigma_tau_pos)^2, mu_tau_pos/(sigma_tau_pos^2));
Omega_pos[i,,] ~ lkj_corr(omega_corr_pos);

tau_neg[i,] ~ gamma((mu_tau_neg/sigma_tau_neg)^2, mu_tau_neg/(sigma_tau_neg^2));
Omega_neg[i,,] ~ lkj_corr(omega_corr_neg);


for (i in 1:Ntotal) {

if (yJud[i]-yPre [i]>0) {
[yPre[i], yJud[i], yPost[i]] ~ multi_normal([muPreSub[s[i]], muJudPosSub[s[i]], muPostSub[s[i]]],
 quad_form_diag(Omega_pos[s[i],,], tau_pos[s[i],]));

if (yJud[i]-yPreFb[i]<0) {
[yPre[i], yJud[i], yPost[i]] ~ multi_normal([muPreSub[s[i]], muJudNegSub[s[i]], muPostSub[s[i]]],
 quad_form_diag(Omega_neg[s[i],,], tau_neg[s[i],]));

if (yJud[i]-yPre [i]==0) {
real muJudZeroSub = (muJudPosSub[s[i]]+muJudNegSub[s[i]])/2 ;
vector[k] tau_zero = (tau_pos[s[i],]+tau_neg[s[i],])/2 ;
matrix[k,k] Omega_zero = (Omega_pos[s[i],,]+Omega_neg[s[i],,])/2;
[yPre[i], yJud[i], yPost[i]] ~ multi_normal([muPreSub[s[i]], muJudZeroSub[s[i]], muPostSub[s[i]]],
quad_form_diag(Omega_Zero[s[i],,], tau_Zero[s[i],]));

And if I merge the equal cases into the positive/negative, e.g. as follows:

for (i in 1:Ntotal) {


if (yJud[i]-yPre [i] **<=** 0) {

[yPre[i], yJud[i], yPost[i]] ~ multi_normal([muPreSub[s[i]], muJudPosSub[s[i]], muPostSub[s[i]]],

 quad_form_diag(Omega_pos[s[i],,], tau_pos[s[i],]));


The model works fine.

Any advice will be greatly appreciated,

Many thanks,

OfiryPreJudPost_multinormal.stan (2.9 KB)

Hi, sorry we took so long the respond.
I am not 100% sure what you aim for and I don’t think I understand the model fully, but if tau is supposed to represent the standard deviation, than you should IMHO average the varainces, e.g.:

vector[k] tau_zero = sqrt( (tau_pos[s[i],] ^ 2 + tau_neg[s[i],] ^ 2)/2 ) ; 

But I am not sure this is so important. Otherwise I don’t see any obvious flaws in the model.

I also find it a bit weird that you branch on the sign of the difference and ignore its magnitude. Why wouldn’t having a single covariance matrix and make the mean effect depend on yJud - yPre. Also, if the rating is on a scale, you might be better of using ordinal models (as brms does -

More broadly, couldn’t your model be expressed in brms, I guess it should be possible (or almost possible) and would save you a lot of headaches with coding… But that’s obviously your call.

Also note that you can use triple backticks (```) to mark chunks of source code in your post.

Hope that helps