Dose Response Model with partial pooling on maximum value

I am trying to get a dose response model ( four parameter log-logistic model) working with the aim of modeling the temperature that winegrapes get too cold and die in the winter, but my model is not behaving well.

The x data I have are simulated daily winter air temperatures (range of about 10 to -10 degrees C). I have added 30 to each temperature so values don’t drop below 30 but remain the same relative to each other. The y data are simulated winegrape hardiness (range about -25 to 0 degrees C) multiplied by -1 so that all values are positive and maximum hardiness (i.e. -25 degrees C) is the biggest value (25). I tried to simulate data to have the exact structure as the Stan model.

I have a model that works fine until I try adding a hierarchical element to the d parameter (maximum hardiness). Now I get ok parameter estimates in general, but lots of warning messages and an odd relationship between the hierarchical element (d_var_sigma) and the log probability density lp__.

This is the model:
\mu_{i}=f(x_{i},(b,c,d,e))=c+\frac{d_{var,i}-c}{1+e^{b(log(x_{i})-\tilde{e})}}
\tilde{y}_{i}\sim normal(\mu_{i},\sigma)

Where:
x is the concentration of the dose (amount of winter cold)
b is the response rate (slope)
d is the upper asymptote of the response (maximum hardiness)
b is the lower asymptote of the response (minimum hardiness)
e is the effective dose ED50 (winter temperature where cold hardiness is half way between min and max)
\tilde{e} is the log of the effective dose ED50

Priors (hardiness has been multipled with -1 to be positive, and 30 has been added to air temp)
b \sim gamma(7,0.75)
d \sim Normal(25, 10)
\sigma_{d,var} \sim gamma(2.5, 1.75)
d_{var} \sim Normal(d, \sigma_{d,var})
c \sim gamma(3,1)
\tilde{e} \sim normal(log(30), 0.1)
\sigma \sim normal(0,5)


data {
  int<lower=1> N;                           // Number of observations
  vector<lower=1>[N] x;                      // Dose values (air temperature)
  vector [N] y;  
  
  //Level 2 	
	int < lower = 1 > n_vars; 				// number of random effect levels (varieties) 
	int < lower = 1, upper = n_vars > variety[N]; // id of random effect (variety)// Response values (Winter cold hardiness). was an int in original code 
}

// Sampling space
parameters {
  real b;                             // slope
  real<lower=0> c;                     // lower asymptote
  real<lower=0> d;                     // upper asymptote
  real<lower=0> ehat;                           //  x where y is half way between c and d, but logged in equation
  real<lower=0> sigma_g;                        //error around the estimated hardiness 

  //level 2
  real <lower = 0> d_var_sigma; 		// a standard deviation of how much d (maximum hardiness) varies 
  real var_d[n_vars]; // a list of the effect of each variety on d (maximum hardiness)
  
}

// Calculate posterior
model {
  vector[N] mu_y;

  // Priors on parameters
  c ~ gamma(3, 1);                    // 
  d ~ normal(25, 10);               // centred around maximum hardiness, 10 sd
  ehat ~ normal(log(30), 0.15);   // centred around mean temperature, sd 20
  b ~ gamma(7, 0.75);                 // 
  sigma_g ~ normal(0, 5);           // 
  
  //level 2
  d_var_sigma ~gamma(2.5, 1.75);
  var_d ~ normal(0, d_var_sigma); //prior for the effect of variety on slope 
  
  //liklyhood function 
  for (i in 1:N) {
    mu_y[i] = c + ((d + var_d[variety[i]] - c) / (1 + exp(b * (log(x[i]) - ehat))));
  }
  y ~ normal(mu_y, sigma_g); // 
}


generated quantities {
  
  // Simulate model configuration from prior model (get mu_y)
  vector[N] mu_y;                      // Simulated mean data from likelyhood 
  vector[N] y_sim;                           //Simulated Data
  real<lower=0> e;
  e = exp(ehat);
  
  //liklyhood function 
  for (i in 1:N) {
    mu_y[i] = c + ((d+ var_d[variety[i]] - c) / (1 + exp(b * (log(x[i]) - ehat))));
  }
  // Simulate data from observational model
  for (n in 1:N) y_sim[n] = normal_rng(mu_y[n], sigma_g); 
}
****

This model gives the following warnings, but does a good job of estimating parameters:

Warning messages:
1: There were 35 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded 
2: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low 
3: Examine the pairs() plot to diagnose sampling problems
 
4: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess 
5: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess 

I think something is going on with the relationship between d_var_sigma and lp__, but I don't really understand enough to figure it out. 

![drc_dh_pairs|580x500](upload://q16Y93nlWrYHcw2XHj7kzitxFQ2.png) 


So far I have tried:
d_var_sigma prior: a half normal prior (0,5) and a gamma prior that is set up in the current version of the model.  
Running the model for 6000 warmup runs and a further 4000 iterations 
Max treedepth 15

I feel a bit lost with this problem, can anyone suggest a sensible way forward?
2 Likes

IIRC @martinmodrak has grappled with models like this one for a while.

1 Like

Yes, I’ve worked with sigmoids in models a bit. They can be challenging to fit.

But first think you might want to check out is “non-centered parametrization” of the hierarchical component, i.e. have parameters var_d_raw ~ normal(0,1); and then derive transformed parameter var_d = var_d_raw * d_var_sigma. (see e.g. https://mc-stan.org/bayesplot/articles/visual-mcmc-diagnostics.html for more details).

I’ve found the direct parametrization of the sigmoid problematic, primarily because if your data don’t cover the whole “dynamic range” of the sigmoid (i.e. data close to both lower and upper plateau and to the middle), some of the parameters can become basically unconstrained by data. My other guess here is that when you have a single d the data cover the whole dynamic range, while when you let d (and hence the whole sigmoid) vary by group (variety), this ceases to be the case for some groups. I would expect that you would see in the posterior is that as d_var_sigma gets larger, all of the sigmoid parameters are allowed to vary wildly.

I worked with @stemangiola on fitting sigmoids in a different context, you might want to check out the parametrization we used in a preprint: https://doi.org/10.1101/2020.03.16.993162, also discussed at Difficulties with logistic population growth model

I also wrote about some of the considerations at Hierarchical Gompertz model

Best of luck with your model!

2 Likes