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):

```
data{
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)