Binomial model with hierarchical randoms

Interfaces brms
1

I try to fit a model being a product of two logit functions (then each from 0 to 1) multiplied by a scale parameter (prior being uniform from 0 to 1) and applying this model to observations with a binomial link.
When I don’t have random factors or only one random factor, it is ok. When I have two random factors, I get this error:
Exception: binomial_lpmf: Probability parameter[1] is 9.06077, but must be in the interval [0, 1]
[the value 9.06077 change].
I don’t understand how it is possible because the product of two logits * a value from 0 to 1 should be between 0 to 1.

Here is a “minimal” reproducible example.

library(brms)

# fake data
set.seed(1)
x <- sample(x = 0:20, size=100, replace = TRUE)
Total <- x + sample(x = 10:20, size=100, replace = TRUE)
y <- runif(100, min=20, max=35)
Ran1 <- c(rep(x = "K", 50), rep(x = "L", 50))
Ran2 <- c(rep(x = c("A", "B"), 25), rep(x = c("C", "D"), 25))

dta <- data.frame(x=x,
                  Total=Total,
                  y=y,
                  Ran1=Ran1, 
                  Ran2=Ran2)

prior <- set_prior(prior = "uniform(20, 35)", class = "b", coef="", 
                   nlpar = "P", lb = 20, ub = 35, check = TRUE) + 
  set_prior(prior = "uniform(0, 12)", class = "b", coef="", 
            nlpar = "DeltaP", lb = 0, ub = 12, check = TRUE) + 
  set_prior(prior = "uniform(0.0001, 4)", class = "b", coef="", 
            nlpar = "SLow", lb = 0.0001, ub = 4, check = TRUE) + 
  set_prior(prior = "uniform(0.0001, 4)", class = "b", coef="", 
            nlpar = "SHigh", lb = 0.0001, ub = 4, check = TRUE) + 
  set_prior(prior = "uniform(0.01, 0.99)", class = "b", coef="", 
            nlpar = "MaxHS", lb = 0.01, ub = 0.99, check = TRUE)

options(mc.cores = 1)

control <- list(
  # adapt_engaged = TRUE,
  adapt_delta = 0.99, #increased from default of 0.8
  # stepsize = 0.05, # 0.05 default
  max_treedepth = 20
)


logit_Ran12 <-brm(bf(x | trials(Total) ~ 
                                          (1/(1+exp((4*abs(SLow))*(P-y)))) * 
                                          (1/(1+exp((4*-abs(SHigh))*(P + abs(DeltaP) -y)))) * MaxHS, 
                                        P ~ (1 + 1 | Ran1 / Ran2), 
                                        SLow ~ (1 + 1 | Ran1 / Ran2), 
                                        SHigh ~ (1 + 1 |Ran1 / Ran2), 
                                        DeltaP ~ (1 + 1 | Ran1 / Ran2), 
                                        MaxHS ~ (1 + 1 | Ran1 / Ran2), 
                                        nl = TRUE),
                                     family = binomial(link="identity"), 
                                     data = dta, 
                                     save_pars = save_pars(all = TRUE), 
                                     control = control, 
                                     sample_prior = TRUE, 
                                     prior = prior,
                                     warmup = 50, 
                                     thin = 1, 
                                     iter = 100, 
                                     chains = 1, 
                                     cores=1,
                                     seed = 123, 
                                     init=rep(list(list(b_P_Intercept=28.93, 
                                                        b_SLow_Intercept=0.05, 
                                                        b_SHigh_Intercept=-0.05, 
                                                        b_DeltaP_Intercept=2, 
                                                        b_MaxHS_Intercept=0.3,  
                                                        sd_Ran1__P_Intercept=0.1, 
                                                        sd_Ran1__SLow_Intercept=0.1, 
                                                        sd_Ran1__DeltaP_Intercept=0.1, 
                                                        sd_Ran1__SHigh_Intercept=0.1,
                                                        sd_Ran1__MaxHS_Intercept=0.1,
                                                        "sd_Ran1:Ran2__P_Intercept"=0.1, 
                                                        "sd_Ran1:Ran2__SLow_Intercept"=0.1, 
                                                        "sd_Ran1:Ran2__DeltaP_Intercept"=0.1, 
                                                        "sd_Ran1:Ran2__SHigh_Intercept"=0.1, 
                                                        "sd_Ran1:Ran2__MaxHS_Intercept"=0.1)),
                                              1))

The results begin by:
SAMPLING FOR MODEL ‘anon_model’ NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: binomial_lpmf: Probability parameter[1] is -0.0489592, but must be in the interval [0, 1] (in ‘anon_model’, line 208, column 4 to column 44)

And many others after

Thanks a lot

brms 2.21.0
R 4.4.0

2

I continue to explore this model trying to understand what I do wrong.
The problem comes from :
MaxHS ~ (1 + 1 | Ran1 / Ran2)
Based on my understanding, the prior of MaxHS is applied before the random effects are applied. Then MaxHS is no more constraint to be between 0 and 1 after application of random effect.
I solve it using:
set_prior(prior = “uniform(-10, 10)”, class = “b”, coef=“”,
nlpar = “MaxHS”, lb = -10, ub = 10, check = TRUE)
and
(1/(1+exp((4abs(SLow)) (P-y)))) *
(1/(1+exp((4*-abs(SHigh))*(P + abs(DeltaP) -y)))) * (1/(1+exp(MaxHS)))
Then this function generates output only in the range [0-1]