Grouped nonlinear fit to breathtest data with brms

brmsgrouped,R

# See https://discourse.mc-stan.org/t/how-to-translate-a-stan-nonlinear-fit-to-brms/12764/2
suppressPackageStartupMessages(library(brms))
library(breathtestcore)
library(dplyr)
library(ggplot2)

set.seed(4711)
data("usz_13c_d", package = "breathtestcore")
data = usz_13c_d 
data = usz_13c_d[usz_13c_d$patient_id < "sub_18",]

# Plot data
ggplot2::ggplot(data, aes(x = minute, y = pdr, color = group)) +
  facet_wrap(~patient_id) +
  geom_point() + geom_smooth(se = FALSE)

form = 
  bf(pdr ~ m*beta/tempt*(1-exp(-minute / tempt))^(beta-1)*exp(-minute/tempt),
     m + tempt + beta ~ group + (1|ID1|patient_id), nl = TRUE)

prior_m_intercept = 4500
prior_tempt_intercept = 90
prior_beta_intercept = 2
lkj = 2
stan_vars = stanvar(prior_m_intercept) + 
  stanvar(prior_tempt_intercept) +
  stanvar(prior_beta_intercept) +
  stanvar(lkj)


prior = c(
  # Interepts
  prior(normal(prior_m_intercept, prior_m_intercept/10), nlpar = "m", coef = "Intercept"),
  prior(normal(prior_beta_intercept, prior_beta_intercept/5), nlpar = "beta", 
        coef = "Intercept"),
  prior(normal(prior_tempt_intercept, prior_tempt_intercept/10), nlpar = "tempt", 
        coef = Intercept),
  
  # Differences
  prior(normal(0, prior_m_intercept/10), nlpar = "m"), 
  prior(normal(0, prior_beta_intercept/10), nlpar = "beta"),
  prior(normal(0, prior_tempt_intercept/3), nlpar = "tempt"),
  
  # sd
  prior(normal(1000, 200), class = sd, nlpar = "m", group = patient_id),
  prior(normal(20, 20), class = sd, nlpar = "tempt", group = patient_id),
  prior(normal(0.3, 0.1), class = sd, nlpar = "beta", group = patient_id),
  prior(lkj(lkj), class = cor)
)

fit1 = NULL
chains = 1
nlevels = length(unique(data$group))
options(mc.cores = chains)

init = replicate( chains, list(
  # First is for intercept, following for differences
  b_m = c(rnorm(1, prior_m_intercept, prior_m_intercept/10), rep(0,nlevels - 1)), 
  b_beta = c(rnorm(1, prior_beta_intercept, prior_beta_intercept/10), rep(0, nlevels - 1)), 
  b_tempt = c(rnorm(1, prior_tempt_intercept, prior_tempt_intercept/2), rep(0 ,nlevels - 1))
), simplify = FALSE)

if (!exists("fit1") || is.null(fit1)) {
  fit1 = brm(form, prior = prior, data = data, family = gaussian(),
             init = init, 
             save_all_pars = TRUE,
             stanvars = stan_vars,
             iter = 1000, chains = chains, 
             control = list(adapt_delta = 0.95))
} else {
  fit1 = update(fit1, newdata = data, iter = 500, chains = chains, init = init)
}


pdf("fit2.pdf")
pp_check(fit1)
plot(fit1, ask = FALSE)
conds <- data.frame(patient_id = unique(data$patient_id), group = unique(data$group) )
marginal_effects(fit1, conditions = conds, re_formula = NULL, ask = FALSE)
dev.off()