Modeling group parameter

Modeling
1

Hi stan users!
I’ve been struggled on modeling the hierarchical model with grouped parameter these days, so I want some advice here.
The goal of my model is to figure out the group parameter of some parameters. Below is simple example of my model.

functions{
 functions{
  vector example_model(real t, vector x, real r1f, real r1r, real r2f, real r3f, real r4f, real r6f ){
    vector[8] dydt;
    
    real r5f;
    
    r5f = 0.083;
   
    dydt[1] = -r1f*x[1]*x[2] +r1r*x[3] + r6f*x[7];
    dydt[2] = -r1f*x[1]*x[2] + r1r*x[3] + r6f*x[7];
    dydt[3] = r1f*x[1]*x[2] -r1r*x[3] - r2f*x[3]*x[4];
    dydt[4] = -r3f*x[5]*x[4]-r2f*x[5]*x[4];
    dydt[5] = r2f*x[3]*x[4]-r4f*x[5];
    dydt[6] = r4f*x[5] -r5f*x[6];
    dydt[7] = r5f*x[6]-r6f*x[7];
    dydt[8] = r5f*x[7];
    
    return dydt;
  }
}

data{
  int<lower=1> T; // number of row
  int<lower=1> N; // number of column (components number) 
  array[N] vector[7] y0; // inital data 
  array[20,7] real z; // observed data
  real t0;
  array[T] real ts; // time seqeunce
}

transformed data{
// to make array vector to vector for using ode solver
  vector[8] y1 = to_vector(y0[,1]);
  vector[8] y2 = to_vector(y0[,2]);
  vector[8] y3 = to_vector(y0[,3]);
  vector[8] y4 = to_vector(y0[,4]);
  vector[8] y5 = to_vector(y0[,5]);
  vector[8] y6 = to_vector(y0[,6]);
  vector[8] y7 = to_vector(y0[,7]);
}

parameters{
  
  // hyperparameters for r1f ~ r6f (group)

  vector<lower = 0, upper = 1>[6] mu;
  vector<lower = 0, upper = 3>[6] sigma;
  vector<lower = 0, upper = 3>[6] sigma_2;

  //subject-level raw parameters for fitting

  vector<lower = 0>[7] r1f_rp;
  vector<lower = 0>[7] r3f_rp;
  vector<lower = 0>[7] r1r_rp;
  vector<lower = 0>[7] r2f_rp;
  vector<lower = 0>[7] r4f_rp;
  vector<lower = 0>[7] r6f_rp;
  
  }
  
transformed parameters{
  
  // subject-level pars
  real<lower = 0> r1f[7];
  real<lower = 0> r3f[7];
  real<lower = 0> r1r[7];
  real<lower = 0> r2f[7];
  real<lower = 0> r4f[7];
  real<lower = 0> r6f[7];
  
  for(i in 1:7){
    
    r1f[i] = mu[1] + sigma[1]*r1f_rp[i];
    r1r[i] = mu[2] + sigma[2]*r1r_rp[i];
    r2f[i] = mu[3] + sigma[3]*r2f_rp[i];
    r3f[i] = mu[4] + sigma[4]*r3f_rp[i];
    r4f[i] = mu[5] + sigma[5]*r4f_rp[i];
    r6f[i] = mu[6] + sigma[6]*r6f_rp[i];

  }

}

model{

  mu ~ normal(0,1);
  sigma ~ cauchy(3,1);
  
  sigma_2 ~ cauchy(3,1);
  
  //prior distribution
  
  r1f_rp ~ normal(1,1);
  r1r_rp ~ normal(1,1);
  r2f_rp ~ normal(1,1);
  r3f_rp ~ normal(1,1);
  r4f_rp ~ normal(1,1);
  r6f_rp ~ normal(1,1);
  

    r1f ~ normal(0,sigma_2[1]);
    r1r ~ normal(0,sigma_2[2]);
    r2f ~ normal(0,sigma_2[3]);
    r3f ~ normal(0,sigma_2[4]);
    r4f ~ normal(0,sigma_2[5]);
    r6f ~ normal(0,sigma_2[6]);

  // ode solve 
  array[T] vector[8] z_hat1 = ode_bdf(example_model, y1, t0, ts, r1f[1], r1r[1], r2f[1], r3f[1], r4f[1], r6f[1]);
  array[T] vector[8] z_hat2 = ode_bdf(example_model, y2, t0, ts, r1f[2], r1r[2], r2f[2], r3f[2], r4f[2], r6f[2]);
  array[T] vector[8] z_hat3 = ode_bdf(example_model, y3, t0, ts, r1f[3], r1r[3], r2f[3], r3f[3], r4f[3], r6f[3]);
  array[T] vector[8] z_hat4 = ode_bdf(example_model, y4, t0, ts, r1f[4], r1r[4], r2f[4], r3f[4], r4f[4], r6f[4]);
  array[T] vector[8] z_hat5 = ode_bdf(example_model, y5, t0, ts, r1f[5], r1r[5], r2f[5], r3f[5], r4f[5], r6f[5]);
  array[T] vector[8] z_hat6 = ode_bdf(example_model, y6, t0, ts, r1f[6], r1r[6], r2f[6], r3f[6], r4f[6], r6f[6]);
  array[T] vector[8] z_hat7 = ode_bdf(example_model, y7, t0, ts, r1f[7], r1r[7], r2f[7], r3f[7], r4f[7], r6f[7]);

  array[20,8] real y;
    
  // likelihood
  for(i in 1:20){
      y[i,1] = z_hat1[i*100,8];
      y[i,2] = z_hat2[i*100,8];
      y[i,3] = z_hat3[i*100,8];
      y[i,4] = z_hat4[i*100,8];
      y[i,5] = z_hat5[i*100,8];
      y[i,6] = z_hat6[i*100,8];
      y[i,7] = z_hat7[i*100,8];

    }
  
  for(i in 1:7){
    
    z[,i] ~ lognormal(log(y[,i]), 5);
    
  }
  
}

}

I have 7 data each are different in initial condition. And with these data, I implement ode solver for each data to predict the final component(x[8]). Lastly, I use these as likelihood with real data z.

What I got problem is : when I input below code, I got very low parameters(they were close to 0).


  for(i in 1:7){
    
    r1f[i] = mu[1] + sigma[1]*r1f_rp[i];
    r1r[i] = mu[2] + sigma[2]*r1r_rp[i];
    r2f[i] = mu[3] + sigma[3]*r2f_rp[i];
    r3f[i] = mu[4] + sigma[4]*r3f_rp[i];
    r4f[i] = mu[5] + sigma[5]*r4f_rp[i];
    r6f[i] = mu[6] + sigma[6]*r6f_rp[i];

  }

The purpose of above lines is each of 7 data are belong in to one normal distribution. But I think I’m doing wrong thing.

Second, the result of sampling does not give information of grouped parameter at all. I can only get raw parameters result.
As I know, transformed parameters are constraint for raw_parameters and this applies to raw_parameters when mcmc sampling is going on. Then, the result of raw_parameter refelcts the grouped parameter now?(I’m not sure that I understood correctly.)

This is my first time to modeling hierarchical, so some advice will real help for me.
Thank you for help!

2

I can’t answer your actual question without seeing the rest of the model definition. I do have a few suggestions from what I could see, but this isn’t addressing your main question.

In the ODE system, it’s better to avoid duplicated arithmetic and define intermediate values that can be reused. Something like this:

vector example_model(real t, vector x, real r1f, real r1r, real r2f, real r3f, real r4f, real r6f ){
  real r5f = 0.083;
  vector[8] dydt;
  real neg_rf1x1x2 = -rf1 * x[1] * x[2];
  dydt[1] = neg_rf1x1x2 + r1r*x[3] + r6f*x[7];
  dydt[2] = neg_rf1x1x2  + r1r*x[3] + r6f*x[7];
  dydt[3] = negrf1x1x2 * x[1]*x[2] -r1r*x[3] - r2f*x[3]*x[4];
  dydt[4] = -r3f*x[5]*x[4]-r2f*x[5]*x[4];
  dydt[5] = r2f*x[3]*x[4]-r4f*x[5];
  dydt[6] = r4f*x[5] -r5f*x[6];
  dydt[7] = r5f*x[6]-r6f*x[7];
  dydt[8] = r5f*x[7];
  return dydt;
}

and so on for the other variables that are re-used.

Then if the idea of this

r1f[i] = mu[1] + sigma[1] * r1f_rp[i];

is to manually do non-centered parameterizations, we added an affine transform that makes this much easier. This can just be

parameters {
  vector<offset=mu[1], multiplier=sigma[1]>[N] r1f;
...
model {
  r1f ~ normal(mu[1], sigma[1]);

Then if you drop the offset and multiplier you go back to the centered parameterization.

3

Then this means that I consequently tried centered parameterization because I implemented additional prior declaration. I should try without this term.

4

Unless you have a lot of data per parameter being estimated, the non-centered parameterization with offset and multiplier specified is usually better.