# Need help with first hierarchical model

Hi,
I am trying to fit a bunch of curves (J) to a parametric form of the form f(x) = A*((x - B)/(100 - B))^(-C)
Below is the stan code for fitting these curves individually and I get good fits.

data {
int<lower=4> N;
int<lower=0> J;
vector[N] y;
vector[N] x;
int curveId[N];
vector[J] minX;
vector[J] minY;
}

parameters {
vector<lower = 0>[J] C;
vector<lower=0>[J] A;
vector<lower=0>[J] B;
real<lower=0> eps;
}
transformed parameters {

}
model {

vector[N] y_hat;
for (i in 1:N){
y_hat[i] = (A[cureveId[i]]*((x[i]-B[cureveId[i]])/(100. - B[cureveId[i]]))^(-C[cureveId[i]]));
}

for (j in 1:J){
C[j] ~ normal(1,5);
B[j] ~ skew_normal(minX[j],5,-2);
A[j] ~ normal(minY[j],5);
}
y ~ normal(y_hat,eps);
}
â€śâ€ť"

Now I am trying to add a group level predictor z for each curve J. I would like to extend it to a hierarchical model where I would also like to model A,B,C which depend on z. I am very new to this stuff so I was just trying to start with only modeling B first. Here is my attempt at it but it diverges. Any help or pointers are appreciated. I have been trying to follow the Radon example

pystan_model_code_heir = â€śâ€ť"
data {
int<lower=4> N;
int<lower=0> J;
vector[N] y;
vector[N] x;
int curveId[N];
vector[J] minX;
vector[J] minY;
vector[J] z;
}

parameters {
vector<lower = 0.25,upper=5>[J] C;
vector<lower=0>[J] A;
real<lower=0> eps;

// B model
real<upper=0> beta_B;
real alpha_B;
}
transformed parameters {
vector<lower=0,upper=100>[J] B;
vector[N] y_hat;

for (j in 1:J){
B[j] = exp(alpha_B + beta_B*z[j]);
}

for (i in 1:N){
y_hat[i] = (A[curveId[i]]*((x[i]-B[curveId[i]])/(100. - B[curveId[i]]))^(-C[curveId[i]]));
}
}
model {

for (j in 1:J){
C[j] ~ normal(1,5)T[0.25,5];
A[j] ~ normal(minY[j],5);
}
y ~ normal(y_hat,eps);
}
â€śâ€ť"

Thanks

Parameters defined as transformed
parameters are not sampled.

So you should have your J individual Bs defined in the parameter block.

Then itâ€™s the conditional distribution of the individual level parameters given the hierarchical parameter that you have, nit the value itseld.

So youâ€™ll want to replace

B[j] = exp(alpha_B + beta_B*z[j]);

With something like

B[j] ~ log_normal(alpha_B + beta_B*z[j], some sorta standard deviation);

1 Like

Thank you very much for the clarification!! I made the changes and the results look good.

1 Like