# Dose Response Model with partial pooling on maximum value

**URL:** <https://discourse.mc-stan.org/t/dose-response-model-with-partial-pooling-on-maximum-value/17086>\
**Category:** Modeling\
**Tags:** fitting-issues\
**Created:** [July 29, 2020, 8:02pm UTC](https://discourse.mc-stan.org/t/dose-response-model-with-partial-pooling-on-maximum-value/17086 "2020-07-29T20:02:32Z")\
**Posts on this page:** 1\
**Showing post:** 14

<div class="post-metadata">

**Author:** ![maj684](https://avatars.discourse-cdn.com/v4/letter/m/b9bd4f/32.png) [@maj684](https://discourse.mc-stan.org/u/maj684)\
**Post date:** [September 1, 2020, 2:42am UTC](https://discourse.mc-stan.org/t/dose-response-model-with-partial-pooling-on-maximum-value/17086/14 "2020-09-01T02:42:36Z")

</div>

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](https://discourse.mc-stan.org/t/non-centered-parameterization-on-variance-parameter/584/3)).

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

 ![mu](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/b/bb1952649b24c50ee26956156c3dafc9c3b7b8cf.png)

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.

---

_[View the full topic](https://discourse.mc-stan.org/t/dose-response-model-with-partial-pooling-on-maximum-value/17086)._