Fits exponential beta curves to 13C breath test series data using Bayesian Stan methods. See https://menne-biomed.de/blog/breath-test-stan/ for a comparison between single curve, mixed-model population and Bayesian methods.
stan_fit(
data,
dose = 100,
sample_minutes = 15,
student_t_df = 10,
chains = 2,
iter = 1000,
seed = 4711
)Data frame or tibble as created by cleanup_data,
with mandatory columns patient_id, group, minute and pdr.
It is recommended to run all data through cleanup_data which
will insert dummy columns for patient_id and minute if the
data are distinct, and report an error if not. It is possible to fit single curves,
but this does not make fully use of the stabilization of the estimated t50.
Instead of fitting single curves, use a set of 20+ of the typical data in your
clinic, and append the data set from the single patient as a mix-in.
Dose of acetate or octanoate. Currently, only one common dose for all records is supported.
If mean sampling interval is < sampleMinutes, data are subsampled using a spline algorithm
When student_t_df < 10, the student distribution is used to model the residuals. Recommended values to model typical outliers are from 3 to 6. When student_t_df >= 10, the normal distribution is used.
Number of chains for Stan
Number of iterations for each Stan chain
Optional seed for rstan
Name of model; use names(stanmodels) for other models.
A list of classes "breathteststanfit" and "breathtestfit" with elements
coef Estimated parameters as data frame in a key-value format with
columns patient_id, group, parameter, method and value.
Has an attribute AIC.
data The effectively analyzed data. If density of points
is too high, e.g. with BreathId devices, data are subsampled before fitting.
stan_fit The Stan fit for use with shinystan::launch_shiny
or extraction of chains.
Base methods coef, plot, print; methods from package
broom: tidy, augment.
# This needs some time !!!
library(breathtestcore)
cmdstanr::check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)
#> Warning: The 'fix' argument is deprecated as of CmdStanR 1.0.0 and will be removed in a future release.
library(dplyr, quietly = TRUE, warn.conflicts = FALSE)
d = breathtestcore::simulate_breathtest_data(n_records = 3) # default 3 records
data = breathtestcore::cleanup_data(d$data)
# Use more than 80 iterations and 4 chains for serious fits
fit = stan_fit(data, chains = 1, iter = 80)
#> Error in registerNames(names, package, ".__global__", add): The namespace for package "breathteststan" is locked; no changes in the global variables list may be made.
plot(fit) # calls plot.breathtestfit
#> Error: object 'fit' not found
# Extract coefficients and compare these with those
# used to generate the data
options(digits = 2)
cf = coef(fit)
#> Error: object 'fit' not found
cf %>%
filter(grepl("m|k|beta", parameter)) %>%
select(-method, -group) %>%
tidyr::spread(parameter, value) %>%
inner_join(d$record, by = "patient_id") %>%
select(
patient_id,
m_in = m.y,
m_out = m.x,
beta_in = beta.y,
beta_out = beta.x,
k_in = k.y,
k_out = k.x
)
#> Error: object 'cf' not found
# For a detailed analysis of the fit, use the shinystan library
# shinystan::launch_shinystan(fit$stan_fit)
# \donttest{
# The following plots are somewhat degenerate because
# of the few iterations in stan_fit
library(bayesplot)
#> This is bayesplot version 1.15.0
#> - Online documentation and vignettes at mc-stan.org/bayesplot
#> - bayesplot theme set to bayesplot::theme_default()
#> * Does _not_ affect other ggplot2 plots
#> * See ?bayesplot_theme_set for details on theme setting
color_scheme_set("viridisD")
options(device.height = 2, device.width = 4)
drws = fit$stan_fit$draws(variables = c("beta[1]", "beta[2]", "beta[3]"))
#> Error: object 'fit' not found
mcmc_dens(drws)
#> Error: object 'drws' not found
mcmc_hist(fit$stan_fit$draws(variables = c("beta[1]", "beta[2]", "beta[3]")))
#> Error: object 'fit' not found
mcmc_intervals(fit$stan_fit$draws(variables = c("beta[1]", "beta[2]", "beta[3]")))
#> Error: object 'fit' not found
mcmc_trace(fit$stan_fit$draws(variables = c("beta[1]", "beta[2]", "beta[3]")))
#> Error: object 'fit' not found
mcmc_areas(fit$stan_fit$draws(variables = c("beta[1]", "beta[2]", "beta[3]")))
#> Error: object 'fit' not found
# }