Please also provide the following information in addition to your question:

- Operating System: Win10
- brms Version: 2.3.1

Hi all,

I’m fitting a multilevel Bayesian non-linear mixed-effects model using both brms and rstanarm(for comparison). However, I found brms ran 30 times slower (over 12,000s each chain) than rstanarm (400s each chain) while failed getting to convergence, the rstanarm model converged perfectly due to both trace plot and R-hat.

The nonlinear function is the logistic function with 4 parameters: A, B, xmid, scal. Within the data set, there are 3 categorical predictors, we call them par1, par2, par3. The package rstanarm is automatically re-scaling the priors, so the final prior is approximately the same between the two models.

The Bayesian model I would like to fit is:

1st level: y~N(nonlin(A, B, xmid, scal), sigma_y^2)

2nd level: A~N(A_0, sigma_{A, par1}^2 + sigma_{A, par2}^2 + sigma_{A, par3}^2), B~…, xmid~…, scal~…

3rd level: A_0~N(mu_A, sigma_A^2), sigma_{A, par1}~half-cauchy(0, scale_A),…

The brms model is:

model_brms <- brm(bf(y ~ A + (B-A)/(1+exp((xmid-x)/scal)),

cbind(A,B,xmid,scal) ~ 1 + (1|par1) +(1|par2) + (1|par3), nl =TRUE),

data=data,

prior = c(prior(normal(30,3), nlpar = A, coef = Intercept),

prior(normal(20,3), nlpar = B,coef = Intercept),

prior(normal(60,3), nlpar = xmid,coef = Intercept),

prior(normal(10,3), nlpar = scal,coef = Intercept),

prior(cauchy(0,3.5), class=sigma),

prior(cauchy(0,3.5), class = sd, nlpar = A),

prior(cauchy(0,3.5), class = sd, nlpar = B),

prior(cauchy(0,3.5), class = sd, nlpar = xmid),

prior(cauchy(0,3.5), class = sd, nlpar = scal)),

inits = c(A=30,B=20,xmid=60,scal=10),

cores = 3,

algorithm = ‘sampling’,

iter = 10000,

chains = 3,

warmup = 5000,

control = list(adapt_delta = 0.9, max_treedepth = 20)

)

The rstanarm model is:

model_arm <- stan_nlmer(y ~ SSfpl(x,A,B,xmid,scal) ~

A + B + xmid + scal +

(A|par1) + (B|par1) + (xmid|par1) + (scal|par1) +

(A|par2) + (B|par2) + (xmid|par2) + (scal|par2) +

(A|par3) + (B|par3) + (xmid|par3) + (scal|par3),

data=data,

prior = normal(location = c(30, 20, 60, 10), scale = c(1, 1, 1, 1)),

prior_aux = cauchy(0,1),

algorithm = ‘sampling’,

adapt_delta = 0.9,

iter = 10000,

chains = 3,

warmup = 5000

)

I did not fit the model with correlated random effects, because my data set only has less than 400 data points. Please correct me if there’s anything going wrong in my process of specifying the models! I know it is hard for nonlinear models to converge when there are too many predictors, but I cannot figure out why the stan_nlmer function in rstanarm can converge with the same amount of parameters and almost the same priors.

By the way, it seems like there’s no way to specify random slope in rstanarm, is that correct?

Thanks,

Zoey