Hierarchical model with custom function_very bad fit

Modeling
1

Hi!

I am trying to fit a growth model (length in function of age) by including the effect of the individual.
Each individual has repeated measures of both age and length, but not all individuals have the same amount of data. The growth model follows a custom function including 3 parameters: k, Linf and t0.
Fitting the function for just one individual works alright. However, when I try to fit a hierarchical model to the data, I get very bad predictions.
Does anybody knows what could be the issue here?

Thanks in advance for any input!


data {                      // Data block
int<lower=1> N;            // all data points

int J;                // total of individuals

vector[N] age;
vector[N] l;

int<lower=1,upper=J> ind[N];           //which individual
}

parameters {                // Parameters block
real<lower=0> k[J];           // k should be bigger that 0
real<lower=14> Linf[J];        //linf has to be bigger than 14 cm
real t0[J];    

real<lower=0> sigma_k;
real<lower=0> sigma_Linf;
real<lower=0> sigma_t0;

real<lower=0> sigma;      

real<lower=0> mu_k;
real<lower=14> mu_Linf;
real mu_t0;

}

model  {

//priors
mu_k ~ normal(0.1,0.05);
mu_Linf ~ normal(18,2);
mu_t0  ~ normal(0,2);
sigma ~ cauchy(0, 5); 
sigma_k ~ cauchy(0, 1); 
sigma_Linf ~ cauchy(0, 5); 
sigma_t0 ~ cauchy(0, 1); 

k ~ normal(mu_k,0.05);
Linf ~ normal(mu_Linf,2);
t0 ~ normal(mu_t0,2);

for (n in 1:N)
    l ~ normal(Linf[ind[n]] * exp(-k[ind[n]]*(age[n]-t0[ind[n]])),sigma);

}

generated quantities {
  vector[N] y_rep;
  
  for (n in 1:N)
    y_rep[n] = Linf[ind[n]] * exp(-k[ind[n]]*(age[n]-t0[ind[n]]));

}

Below the data:

data <- data.frame(agei = c(1, 2, 3, 0, 6, 7, 4, 5, 6, 5, 0, 1, 2, 
3, 4, 2, 3, 4, 1, 0, 5, 0, 1, 2, 4, 5, 3, 5, 3, 4, 0, 1, 2, 7, 
6, 0, 1, 3, 4, 5, 2, 6, 1, 2, 3, 0, 4, 5, 6, 7, 8, 9, 10, 6, 
3, 4, 5, 0, 1, 2, 7, 0, 1, 2, 4, 5, 6, 3, 8, 7, 1, 2, 0, 4, 5, 
6, 3, 8, 7, 5, 6, 4, 0, 1, 2, 3, 7, 8, 0, 1, 8, 9, 4, 5, 2, 3, 
6, 7, 0, 1, 3, 4, 5, 2, 6, 5, 6, 7, 0, 1, 2, 3, 4, 9, 10, 11, 
8, 13, 12, 5, 6, 4, 0, 1, 2, 3, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 
2, 3, 1, 0), length = c(6.7, 9.9, 11, 0.1, 14.3, 15, 12.5, 13.5, 
13.8, 12.5, 0.8, 6.7, 9, 10.3, 11.4, 9.6, 11.9, 13, 7.6, 0.8, 
13.8, 0.1, 7.6, 10.7, 14, 15, 12.4, 12.9, 11.6, 12.3, 0.1, 9.2, 
10.6, 13.9, 13.3, 0.1, 7.8, 10.7, 12.1, 13.1, 9.4, 14.2, 7.2, 
8.4, 9.1, 0.1, 10.2, 11.1, 11.7, 12.2, 12.9, 13.5, 13.9, 15.3, 
12.3, 13.1, 14.1, 0.1, 7.2, 10.5, 16.1, 0.1, 7.9, 10.2, 12.9, 
14.1, 14.8, 11.4, 15.9, 15.4, 6, 8, 0.1, 10.9, 11.9, 12.6, 9.8, 
14.2, 13.5, 13, 13.9, 12.1, 0.1, 7.4, 9.5, 11.1, 14.9, 15.8, 
0.1, 8.7, 13.8, 14.2, 11.6, 12.2, 10.3, 10.9, 12.8, 13.3, 0.1, 
7.1, 10.8, 12, 12.8, 9.5, 13.7, 9.8, 10.4, 11, 0.1, 7.4, 8.1, 
8.8, 9.3, 12, 12.5, 12.8, 11.5, 13.6, 13.2, 12.1, 12.7, 11.3, 
0.1, 7.8, 10, 10.5, 13.2, 11.9, 0.1, 7, 8.8, 9.4, 9.9, 10.5, 
11.2, 11.6, 9.7, 10.5, 7.6, 0.1), rep = c("001", "001", "001", 
"001", "001", "001", "001", "001", "002", "002", "002", "002", 
"002", "002", "002", "003", "003", "003", "003", "003", "003", 
"004", "004", "004", "004", "004", "004", "005", "005", "005", 
"005", "005", "005", "005", "005", "006", "006", "006", "006", 
"006", "006", "006", "007", "007", "007", "007", "007", "007", 
"007", "007", "007", "007", "007", "008", "008", "008", "008", 
"008", "008", "008", "008", "009", "009", "009", "009", "009", 
"009", "009", "009", "009", "010", "010", "010", "010", "010", 
"010", "010", "010", "010", "011", "011", "011", "011", "011", 
"011", "011", "011", "011", "012", "012", "012", "012", "012", 
"012", "012", "012", "012", "012", "013", "013", "013", "013", 
"013", "013", "013", "014", "014", "014", "014", "014", "014", 
"014", "014", "014", "014", "014", "014", "014", "014", "015", 
"015", "015", "015", "015", "015", "015", "015", "016", "016", 
"016", "016", "016", "016", "016", "016", "016", "017", "017", 
"017", "017")), .Names = c("agei", "length", "rep"), class = "data.frame", row.names = c(NA, 
-140L))

data_stan <- list(
 N = nrow(data),
 J = length(unique(data$rep)),
 l = data$length,
 age = as.numeric(data$agei),
 ind = as.integer(data$rep)
)
2

Hello!

Do you have any warning after fitting your model?

You do not use the hierarchical sigmas. In your code, you have fixed sd for the variability of k, linf et t0. You should replace the numbers by appropriate sigmas.

If you have trouble fitting your model, you might try to use non-centered parametrization, as explained in the stan manual. But it’s more about warnings or other sampling problems than very bad fit, I think.

Concerning the generated quantities block, you only compute the mean. Usually, we use normal_rng to get the posterior predictive distribution of every data point. The formulation is the same as in the likelihood in the model part.

Have a nice day!
Lucas

3

Hello,

Yes, I get the following warnings:
Warning messages:
1: There were 25 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
2: There were 478 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
3: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
4: Examine the pairs() plot to diagnose sampling problems

updated stan model:

data {                      // Data block
int<lower=1> N;            // all data points

int J;                // total of individuals

vector[N] age;
vector[N] l;

int<lower=1,upper=J> ind[N];           //which individual
}

parameters {                // Parameters block
real<lower=0> k[J];           // k should be bigger that 0
real<lower=14> Linf[J];        //linf has to be bigger than 14 cm
real t0[J];    

real<lower=0> sigma_k;
real<lower=0> sigma_Linf;
real<lower=0> sigma_t0;

real<lower=0> sigma;      

real<lower=0> mu_k;
real<lower=14> mu_Linf;
real mu_t0;

}

transformed parameters{

}

model  {

//priors
mu_k ~ normal(0.1,0.05);
mu_Linf ~ normal(18,2);
mu_t0  ~ normal(0,2);
sigma ~ cauchy(0, 5); 
sigma_k ~ cauchy(0, 1); 
sigma_Linf ~ cauchy(0, 5); 
sigma_t0 ~ cauchy(0, 1); 

k ~ normal(mu_k,sigma_k);
Linf ~ normal(mu_Linf,sigma_Linf);
t0 ~ normal(mu_t0,sigma_t0);

for (n in 1:N)
    l ~ normal(Linf[ind[n]] * exp(-k[ind[n]]*(age[n]-t0[ind[n]])),sigma);

}

generated quantities {
  vector[N] y_rep;
  
  for (n in 1:N)
    y_rep[n] = normal_rng(Linf[ind[n]] * exp(-k[ind[n]]*(age[n]-t0[ind[n]])),sigma);

}

If I want to use non-centered parametrization here, how exactly can I do this? I understand how to do this with simple linear regressions, but here I am not sure.

Thanks a lot in advance!

Cheers,

Nina

4

Problem solved!

Model is fitting a lot better now!

code below in case anybody is interested.

data {                      // Data block
int<lower=1> N;            // all data points

int J;                // total of individuals

vector[N] age;  //covariate vector
vector[N] l;    //response variable

vector[N] s;  //vector for structure, should be 1

int<lower=1,upper=J> ind[N];           //which individual
}

parameters {                // Parameters block
real<lower=0> k;           // k should be bigger that 0
real<lower=14> linf;        //linf has to be bigger than 14 cm
real t0;    

real<lower=0> sd_k;
real<lower=0> sd_linf;
real<lower=0> sd_t0;

real<lower=0> sigma;      

real z_k[J];
real z_linf[J];
real z_t0[J];

}

transformed parameters { 
 // group-level effects 
  vector[J] r_linf;
  // group-level effects 
  vector[J] r_k;
  // group-level effects 
  vector[J] r_t0;
for (j in 1:J){
  // group-level effects 
  r_linf[j] = sd_linf * z_linf[j];
  // group-level effects 
  r_k[j] = sd_k * z_k[j];
  // group-level effects 
  r_t0[j] = sd_t0 * z_t0[j];
}
} 

model  {
  
  vector[N] nlp_linf = linf * s;
  vector[N] nlp_k = k * s;
  vector[N] nlp_t0 = t0 * s;
  vector[N] mu;
  for (n in 1:N) { 
    nlp_linf[n] += r_linf[ind[n]];
    nlp_k[n] += r_k[ind[n]];
    nlp_t0[n] += r_t0[ind[n]];
    // compute non-linear predictor 
    mu[n] = (nlp_linf[n] * (1 - exp( - nlp_k[n] * (age[n] - nlp_t0[n]))));
  } 

//priors
  target += normal_lpdf(linf | 16, 3); 
  target += normal_lpdf(k | 0.3, 0.15); 
  target += normal_lpdf(t0 | 0, 0.01); 
  target += cauchy_lpdf(sigma | 0, 5);
  target += cauchy_lpdf(sd_linf | 0, 5);
  target += normal_lpdf(z_linf | 0, 1);
  target += cauchy_lpdf(sd_k | 0, 1);
  target += normal_lpdf(z_k | 0, 1);
  target += cauchy_lpdf(sd_t0 | 0, 0.01);
  target += normal_lpdf(z_t0 | 0, 1);
}
generated quantities {
  vector[N] y_rep;
  vector[N] nlp_linf = linf * s;
  vector[N] nlp_k = k * s;
  vector[N] nlp_t0 = t0 * s;
  real mu_linf;
  real mu_k;
  real mu_t0;
 
    for (n in 1:N) { 
      nlp_linf[n] += r_linf[ind[n]];
      nlp_k[n] += r_k[ind[n]];
      nlp_t0[n] += r_t0[ind[n]];
      y_rep[n] = (nlp_linf[n] * (1 - exp( - nlp_k[n] * (age[n] - nlp_t0[n]))));
  } 
  
 mu_linf = mean(nlp_linf);
 mu_k = mean(nlp_k);
 mu_t0 = mean(nlp_t0);
  


}
1 Like