Multivariate growth outcomes with multivariate Gaussian

Modeling
,
1

Hi all,

I am attempting to specify a model that approximates the underlying data-generating process of several traits related to human growth. Let x be a scalar independent variable and y be a vector of observed continuous variables indexed by k = 1,2,...,K. I assume the overall response vector is distributed as y_n ~ N (h(x), \sigma^2(x)) where h(x) is a function for the mean and \sigma^2(x) a function for the standard deviation. The mean is a K length response vector where h_k(x) = a_k*x^r_k. For the sd I assume linearity where \sigma(x) = s_k[1+\kappa_k*x]. Essentially each continuous response has a parameter vector [a, r, b, s, \kappa]. The ultimate goal is to translate this up where K is quite large and is a mix of continuous and ordinal traits, but currently I have simulated some data and am trying to test the model with 2 continuous responses.

I attach a toy model that does run, but throws MANY divergence errors, so I am curious as to the best way to code out what I’ve described above. Note, I have not included the SD function here because when I do, the diag_pre_multiply() function throws an error. I’m sure I could supply better priors and such, but I am most concerned about the most efficient way to model the mean and sd function and then the decomposition of \Sigma. From research, I know noise increases with x (demonstrated with the sd function) and I know there is dependence between response traits.

My simulated input data looks something like this:

[1] 2
[1] 1000
[1] 0.75311542 1.15415400 0.06741298 2.67019476 3.96470039 3.32218096
          [,1]      [,2]
[1,]  49.12179  50.41544
[2,]  66.51528  66.02949
[3,] -15.19807  37.58481
[4,] 130.97413  96.35293
[5,] 152.44352 104.65998
[6,] 138.78272 106.90922

data{
  //multivariate outcomes
  int N; // # of individuals
  int K; // # of responses
  vector[K] y[N]; // vector of responses per individual
  //predictors
  real X[N]; //age
}
parameters{
  vector[K] a;
  vector[K] r;
  vector[K] b;
  vector[K] s_scale;
  vector[K] kappa;
  //cholesky factors for rho
  cholesky_factor_corr[K] L_Omega;
}
transformed parameters{
  vector[K] mu[N];
  matrix[K,K] Sigma;
  for(i in 1:N){
    for(k in 1:K){
      mu[i,k] = a[k]*X[i]^(1-r[k])+b[k];
      }
  }
  Sigma = diag_pre_multiply(s_scale,L_Omega);
}
model{
  for(k in 1:K){
    a[K] ~ gamma(0.001,0.001);
	  r[K] ~ gamma(0.001,0.001);
	  b[K] ~ normal(0,1000);
	  kappa[K] ~ gamma(0.001,0.001);
	  s_scale[k] ~ gamma(4,1);
  }
  L_Omega ~ lkj_corr_cholesky(1);
	y ~ multi_normal_cholesky(mu, Sigma);
}