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.
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?
