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 );
}
}
"