Skip to content

Instantly share code, notes, and snippets.

@dmenne
Last active Jan 23, 2020
Embed
What would you like to do?
Grouped nonlinear fit to breathtest data with brms
# 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment