Hi all,

I will preamble by saying thank for reading this, and that any commenter should feel free to criticise any aspect of this modelling specification or workflow you deem naive or uninformed.

I am simulating two groups from this double sigmoid curve process:

Using the following code:

```
sw_dbl_log <- function(a, b, c, d, time){
((50)/(1+exp(-a*Time + b)))+((50)/(1+exp(-c*Time + d)))
}
nn = 300
Raw<- tibble(n = 1:(nn*2),
Group = rep(c("Control", "Experimental"), each = nn),
a = c(rnorm(nn, 0.1, 0.01),
rnorm(nn, 0.1, 0.01)),
b = c(rnorm(nn, 5, 0.01),
rnorm(nn, 5, 0.01)),
c = c(rnorm(nn, 0.1, 0.01),
rnorm(nn, 0.1, 0.002)),
d = c(rnorm(nn, 5, 0.01),
rnorm(nn, 25, 0.01))
) %>%
tidyr::expand(nesting(n,Group, a, b, c, d),
Time = seq(0, 300, by = 10)) %>%
mutate(y = sw_dbl_log(a, b, c, d, Time))
```

As expected, this produces gaussian distributions for each Group’s parameters (a:d), with the Control group’s time series resembling a single sigmoid with one growth phase. The Experimental group’s time series has two phases of growth, as intended. I am leveraging this specification because each parameter gives a clearly interpretable output, that is, a and c represent growth rates and b represents the time point of inflexion for each of the growth phases.

My fitting procedure is as follows :

```
formula <- bf(y ~ (50/(1+exp(-a*Time + b))) + (50/(1+exp(-c*Time + d))),
a + b + c + d ~ 0 + Group,
set_rescor(FALSE),
nl = TRUE)
prior1 <- c(
prior(normal(0.1, 0.1), nlpar = "a", lb = 0),
prior(normal(5, 2), nlpar = "b", lb = 0, ub = 10),
prior(normal(0.1, 0.1), nlpar = "c",lb = 0),
prior(normal(25, 5), nlpar = "d", coef = "GroupExperimental"),
prior(normal(5, 2), nlpar = "d", coef = "GroupControl")
)
fit <- brm(formula,
data = Raw,
prior = prior1,
control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = 10000,
warmup = 5000,
cores = 3,
chains = 3)
```

I’ve omitted prior fitting and predictive checks from the code above, but visual inspection looks good.

This issue I am currently having is that the model fit is poor, namely, the Rhat and ESS numbers are too low, and it’s principally due to bimodality in the posterior distribution of the parameters for either one of the groups, sometimes Control and sometimes Experimental, that is causing the chains not to converge.

I cannot seem to define why the model is fitting poorly and causing this bimodality. My first guess is that the priors are causing issues, but, as per my naïve understanding, the priors are specified pretty tightly around the original simulation distributions so I cannot see why they might be causing the problem. I must admit that this conclusion might expose my ignorance on the topic which is fueling my inability to fix the fit issues.

Should I change the priors, the model, stipulate init values or simply run for more iterations? Any help would be greatly appreciated.