How to translate a Stan nonlinear fit to brms

Interfaces brms
1

I am trying to demonstrate how a Stan nonlinear fit can be “simplified” by brms, and failed miserably. Data set is from https://github.com/dmenne/breathtestcore, Stan in https://github.com/dmenne/breathteststan

library(ggplot)
library(breathteststan) # on cran
library(breathtestcore)
library(brms)

data(usz_13c)
# Use a subset
solids = str_detect(usz_13c$patient_id, "00") &  str_detect(usz_13c$group, "solid")
bt = usz_13c[solids,] 
# cleanup_data standardizes columns
# Stan does not converge when samples are too dense, so do subsampling
bt = cleanup_data(breathtestcore::subsample_data(bt,15))

# Group parameter (patient/normal) is not used here, but that's the subject 
# for another day

#plot(null_fit(bt))

# https://github.com/dmenne/breathteststan/blob/master/inst/stan/breath_test_1.stan
# My simple stan model. 1000 iterations, about 30 seeconds on my computer
fit = stan_fit(bt)
plot(fit)
cf = coef(fit)
cf[cf$parameter %in% c("m","k","beta"),c(2,3,5)]

# Simplified model, incorporated constant "dose" parameter into m
fit_model = brm( bf(pdr ~ m*k*beta*(1-exp(-k*minute))^(beta-1)*exp(-k*minute),
                    m  ~ 1 + (1|patient_id),
                    k  ~ 1 + (1|patient_id),
                    beta ~ 1 + (1|patient_id), 
                    nl = TRUE), 
                 data = bt, 
                 chains = 1,
                 iter = 1,
                 family = gaussian(),
                 control = list(adapt_delta = 0.9, max_treedepth = 20),
                 inits = list(list(m = 2000, k = 0.02, beta = 2 )),
                 prior = c(
                   prior(normal(2000,400), nlpar = "m"),
                   prior(lognormal(-5, 1), nlpar = "k"),
                   prior(normal(2, 0.6), nlpar = beta)
                 ))
update(fit_model, iter = 1000)
# No convergence.....
1 Like
2

Sorry, Paul, with your help some years ago we had already solved this.

See https://gist.github.com/dmenne/db4bacb7128f2922aaa75bcd9c1f8766

3

So no more action required?

4

No. I still don’t understand why the simple problem does not work, but I am checking.

Closed and thanks.

5

Krrr… I forgot the dummy ID … Done this before.

m  ~ 1 + (1|ID|patient_id),
k  ~ 1 + (1|ID|patient_id),
beta ~ 1 + (1|ID|patient_id),