It would be useful to see your data in order to run the model and check the diagnostics first hand. However, here is an alternative formulation, albeit using a sigmoidal emax rather than log-logistic, but I believe they are similar in what they offer in terms of their parametric form and I have used the sigmoid emax previously so am more more familiar with this form. I haven’t used a random effect on the parameters for these models previously, but have introduced one here on \beta_2 using a non-centered parameterisation (you might want to check the manual on this, or Non centered parameterization on variance parameter).
data {
int<lower=0> N;
int<lower=1> Nc;
real y[N];
vector[N] x;
int c[N];
vector[4] mu0;
vector[5] sig0;
}
parameters {
// upper and lower asymptote
real b1;
real<lower=0> b3; // the ed50
real b4; // governs the shape of the curve
real<lower=0> xi;
vector[Nc] b2_z;
real<lower=0> nu;
}
transformed parameters{
real mu[N];
vector[Nc] b2;
b2 = b2_z * nu;
for(i in 1:N){
mu[i] = b1 + b2[c[i]] * pow(x[i], b4) * inv(pow(b3, b4) + pow(x[i], b4));
}
}
model {
target += normal_lpdf(b1 | mu0[1], sig0[1]);
target += normal_lpdf(b2_z | 0, 1);
target += normal_lpdf(nu | 0, 5);
target += normal_lpdf(b3 | mu0[3], sig0[3]) - normal_lccdf(0 | mu0[3], sig0[3]);
target += normal_lpdf(b4 | mu0[4], sig0[4]);
target += normal_lpdf(xi | 0, sig0[5]);
target += normal_lpdf(y | mu, xi);
}
generated quantities {
}
Based on your description, I simulated some data as follows:
sigemax <- function(x, a, b, c, d){
a + b * (x^d) / (c^d + x^d)
}
# data simulation
N <- 50
Nc <- 5
n <- N / Nc
set.seed(6)
x <- seq(20, 40, len = n)
x <- rep(x, len = N)
a <- 0
b <- rnorm(Nc, 0, 2)
c <- 30
d <- 20
k <- rep(1:Nc, each = n)
y <- sigemax(x, a, 20 + b[k], c, d)
y <- y + rnorm(N, 0, 1)
plot(x, y, xlim = c(20, 40), col = k)
and then compiled and ran the sampler as usual
mod1 <- rstan::stan_model("emax1.stan", verbose = F)
ld <- list(N = N,
Nc = Nc,
c = k,
y = y,
x = x,
mu0 = c(0, 25, 30, 20),
sig0 = c(2, 4, 3, 10, 4)
)
library(rstan)
library(scales)
# mc (many cores please)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
fit1 <- rstan::sampling(object = mod1,
data = ld,
chains = 8,
thin = 1,
iter = 5000,
refresh = 1000,
control=list(adapt_delta=0.9, max_treedepth=12))
you can check on the results with
params <- rstan::extract(fit1)
bayesplot::mcmc_dens(fit1, pars = c("b1", "b3", "b4"), facet_args = list(nrow = 3))
bayesplot::mcmc_dens(fit1, pars = c("xi"), facet_args = list(nrow = 1))
bayesplot::mcmc_dens(fit1, pars = paste0("b2[", 1:5,"]"), facet_args = list(nrow = 5))
20 + b
and
xnew <- seq(20, 40, by = 1)
mu <- lapply(1:Nc, function(i){
sapply(seq_along(xnew), function(j){
sigemax(xnew[j], params$b1, params$b2[, i], params$b3, params$b4)
})
})
lwr <- sapply(mu, function(z){
apply(z, 2, function(w) quantile(w, probs = c(0.1)))
})
upr <- sapply(mu, function(z){
apply(z, 2, function(w) quantile(w, probs = c(0.9)))
})
mumu <- sapply(mu, function(w){colMeans(w)})
plot(1, xlab = "x", ylab = "mu", xlim = range(xnew), ylim = range(mumu))
for(i in 1:Nc){
lines(xnew, mumu[, i], col = i) # scales::alpha("black", 0.7))
lines(xnew, lwr[, i], col = i, lty = 2)
lines(xnew, upr[, i], col = i, lty = 2)
}
which seem at least somewhat reasonable and I didn’t have any issues raised by stan. However, I will say that I am aware that these models are quite sensitive to the priors and so I have made them relatively informative here based on the description of your study.
I’d also add that I have found dynamic linear models (DLMs) to be particularly useful for dose-response modelling as they essentially offer smoothing without enforcing a parametric form. Scott Berry outlines some basic DLMs in his 2011 text on Bayesian adaptive clinical trials.