Curing divergences for a 3-parameter sigmoid model

Modeling
,
1

Hello all. I am trying to fit a nonlinear model specified as follows:
y = d + \frac{1-d}{1+\exp(-b(\log(x) - \log(c)))}
which is the 3-parameter log-logistic function for a curve with an upper limit of 1, a lower limit of d, an inflection point of c, and a slope of b. I have placed the following priors:
d\sim \mathcal{N}(0, 0.5) \\ \log(c)\sim \mathcal{N}(0, 5) \\ b\sim \mathcal{N}(0, 20)
where d\in(0, 1) and b\in\mathbb{R}_{\le 0}, enforcing decreasing curves.

In Stan, I have written:

data {
    int<lower=1> N;
    vector[N] log_x;
    vector<lower=0>[N] y;
}
parameters {
    real<lower=0, upper=1> d;
    real log_c;
    real<upper=0> b;
    real<lower=0> sigma;
}
model {
    d ~ normal(0, 0.5);
    log_c ~ normal(0, 5);
    b ~ normal(0, 20);
    sigma ~ normal(0, 1);

    vector[N] mu;
    for (i in 1:N) {
        mu[i] = d + (1 - d) * inv_logit(b*(log_x[i] - log_c));
    }

    y ~ normal(mu, sigma);
}

Testing this model, I have simulated data in R that are similar to one of my use cases:

# Simulated data
set.seed(1)
sigmoid_fxn <- function(log_x, d, log_c, b) d + (1 - d) / (1 + exp(-b*(log_x - log_c)))
data <- data.frame(x = rep(10 / 3^(7:0), each = 3))
data["log_x"] <- log(data$x)
data["y"] <- sigmoid_fxn(data[["log_x"]], 0.1, log(5), -10)
data["y"] <- data["y"] + (runif(3 * length(data["y"])) * 0.1 - 0.05)
# Model fit
library(rstan)
fit <- stan(file = "sigmoid.stan",
            data = list(N = dim(data)[1],
                        log_x = data[["log_x"]],
                        y = data[["y"]]),
            control = list(adapt_delta = 0.95),
            chains = 4, cores = 4, iter = 10000, seed = 1)

Which fails even with adapt_delta increased to 0.95.

Warning messages:
1: There were 225 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them. 
2: Examine the pairs() plot to diagnose sampling problems

Examining the pairs plot of the unconstrained variables suggests a problem:

library(bayesplot)
color_scheme_set("darkgray")
np <- nuts_params(fit)
mcmc_pairs(fit, pars = c("d", "log_c", "b", "sigma"), np = np,
		   transformations = list(d = \(d) log(d / (1 - d)),
		   				          b = \(b) log(-b),
    	       					  sigma = log))

2025-07-04_pairs-plot
2025-07-04_pairs-plot700×700 57.8 KB

It seems to me that the problem is centered on the narrow tip of the lower bound d. However, I am not sure how to move forward from here to solve this problem. Is there a reparameterization I can try to fix this problem, or another technique that I can use? I would appreciate any advice that you have–thank you.

1 Like
2

Does the model still yield a bunch of divergences if you simulate the data under the model you are fitting?

Ignoring the simulation for a second and focusing just on the stan model, I note that you have a normal likelihood but a lower bounded response. Is this what you intend?

1 Like
3

Hi, sigmoid are unfortunately very hard to fit unless your data cover the full dynamic rangeof the sigmoid. See e.g. Dose Response Model with partial pooling on maximum value - #3 by martinmodrak which provides links and some followup discussion and ideas. I also believe therecwas a StanConnect bio/med talk on the topic but I am not fully sure.

Best of luck with the model!

4

Thank you for pointing out the likelihood. From reviewing the literature, it seems that a gamma distribution is a reasonable choice. I have assumed Var(y) = \phi * \mathbb{E}(y) , where \phi is a noise parameter, and reparameterized the Stan gamma distribution as \alpha = \mathbb{E}(y) / \phi and \beta = 1 / \phi. In Stan:

...
model {
    d ~ normal(0, 0.5);
    log_c ~ normal(0, 5);
    b ~ normal(0, 20);
    phi ~ normal(0, 1);

    vector[N] mu;
    for (i in 1:N) {
        mu[i] = d + (1 - d) * inv_logit(b*(log_x[i] - log_c));
    }
	
	// gamma likelihood reparameterization
    y ~ gamma(mu / phi, 1 / phi);
}

This seems to resolve the sampling problem for curves with an inflection point towards the center of the tested dose range. I think that the problem for curves outside of this ideal range is due to the available data, rather than the model specification. Thank you for your help!

1 Like
5

Thank you for the link, this has some excellent discussion.

1 Like