Just getting back to this. Maybe you can help clear up some confusion. I’m trying to specify a multinomial model where each behavior category K varies as a (polynomial) function of age, and the correlations between the random effects also vary as a function of age. I put something together using the functions you provided and the appendix of the paper you linked. The model fits OK, but I am confused by the product of the functions that you shared. It looks like the sdcovsqrt2cov function should return the cholesky factor of a co-variance matrix, but it does not (it returns a matrix with values in the upper triangle, and the matrix is not symmetric). Here’s my code:

functions{

matrix covsqrt2corsqrt(matrix mat, int invert){ //converts from lower triangular unconstrained matrix to cor sqrt matrix

matrix[rows(mat),cols(mat)] o;

vector[rows(mat)] s;

o=mat;

for(i in 1:rows(o)){ //set upper tri to lower

for(j in min(i+1,rows(mat)):rows(mat)){

o[j,i] = inv_logit(o[j,i] )*2-1;

o[i,j] = o[j,i];

}

o[i,i]=1; // change to adjust prior for correlations

}

if(invert==1) o = inverse(o);

for(i in 1:rows(o)){

s[i] = inv_sqrt(o[i,] * o[,i]);

if(is_inf(s[i])) s[i]=0;

}

o= diag_pre_multiply(s,o);

return o;

}

matrix sdcovsqrt2cov(matrix mat, int cholesky){ //converts from lower triangular unconstrained and diag sd to cov or cholesky cov

matrix[rows(mat),rows(mat)] out;

int invert;

invert = 0; //change integer for marginal or partial prior

if(rows(mat) > 1){

out = covsqrt2corsqrt(mat,invert);

out = diag_pre_multiply(diagonal(mat), out);

}

if(rows(mat)==1) out[1,1] = mat[1,1];

if(cholesky==0) out = out * out’;

return(out);

}

}

data{

int N; // number of obs

int N_id; // number of individuals

int y[N]; // multinomial outcome

int id[N];

int K; // number of outcome behavioral categories

int ndim; // K-1

real age_z[N]; // standardized age

real age_zq[N]; // standardized age^2

real age_id[N_id]; // age for each individual

real age_idq[N_id]; // age^2 for each individual

}

parameters{

real a[K-1]; // intercepts for each behavior, minus reference category

real bA[K-1]; // parameters for age

real bQ[K-1]; // parameters for age^2

matrix[N_id,K-1] v_id; // matrix of random effects

vector<lower=0>[K-1] sigma_id; // stddev of random effects

vector[(ndimndim-ndim)/2] pcorrvecmu; // pop mean of partial correlation vectors

vector[(ndimndim-ndim)/2] pcorrvec_age; // pop mean of effect of age on correlations

vector[(ndim*ndim-ndim)/2] pcorrvec_age2; // pop mean of effect of age sq on correlations

}

transformed parameters{

vector[(ndimndim-ndim)/2] p_merged[N_id]; // combined correlation parameters

matrix[ndim,ndim] pcovmat[N_id]; // lower cov matrix plus sigmas on diag

matrix[ndim,ndim] cholcovmat[N_id]; // cholesky form

{

int counter;

for(n in 1:N_id){

counter=0;

// p_merged is a vector for the global correlation parameters, plus the age parameters

p_merged[n] = pcorrvecmu + pcorrvec_ageage_id[n] + pcorrvec_age2*age_idq[n];

// putting sd for each behavior on the diagonal

pcovmat[n] = diag_matrix(sigma_id);

// filling matrix

for (j in 1:(ndim-1))

for (i in (j+1):ndim){

counter = counter + 1;

pcovmat[n][i,j] = p_merged[n][counter];

}

cholcovmat[n] = sdcovsqrt2cov(pcovmat[n],1);

}

}

}

model{

// priors for fixed effects

a ~ normal(0,1);

bA ~ normal(0,1);

bQ ~ normal(0,1);

// priors for unconstrained corr parameters

pcorrvecmu ~ normal(0,1);

pcorrvec_age ~ normal(0,1);

pcorrvec_age2 ~ normal(0,1);

// hyper-priors

sigma_id ~ exponential(1);

for (n in 1:N_id){

v_id[n,] ~ multi_normal_cholesky(rep_vector(0,K-1),cholcovmat[n]);

}

// Likelihood function

for ( i in 1:N ) {

vector[K] p;

for ( k in 1:(K-1) )

p[k] = a[k] + v_id[id[i],k] + bA[k]*age_z[i] + bQ[k]*age_zq[i];

p[K] = 0;

y[i] ~ categorical_logit( p );

}

}

generated quantities{

vector[N] log_lik;

for ( i in 1:N ) {

vector[K] p;

for ( k in 1:(K-1) )

p[k] = a[k] + v_id[id[i],k] + bA[k]*age_z[i] + bQ[k]*age_zq[i];

p[K] = 0;

log_lik[i] = categorical_logit_lpmf( y[i] | p );

}

}

"