bmn {BayesDissolution}R Documentation

Bayesian Multivariate Normal Model for Dissolution Profile Modeling


This function implements the Bayesian multivariate normal model described in Pourmohamad et al (2022).


bmn(dis_data, B = 10000)



A data frame containing the dissolution data. The first column of the data frame should denote the group labels identifying whether a given dissolution belongs to the "reference" or "test" formulation group. For a given dissolution run, the remaining columns of the data frame contains the individual run's dissolution measurements sorted in time.


A positive integer specifying the number of posterior samples to draw. By default B is set to 10000.


The function returns a list of B posterior samples for the following parameters:


You should always check MCMC diagnostics on the posterior samples before drawing conclusions.


Novick, S., Shen, Y., Yang, H., Peterson, J., LeBlond, D., and Altan, S. (2015). Dissolution Curve Comparisons Through the F2 Parameter, a Bayesian Extension of the f2 Statistic. Journal of Biopharmaceutical Statistics, 25(2):351-371.

Pourmohamad, T., Oliva Aviles, C.M., and Richardson, R. Gaussian Process Modeling for Dissolution Curve Comparisons. Journal of the Royal Statistical Society, Series C, 71(2):331-351.


### dis_data comes loaded with the package
### We set B = 1000 to obtain 1000 posterior samples
B <- 1000
post <- bmn(dis_data, B = B)

### We can check how well the posterior samples of the means are mixing by
### plotting the individual chains by time point
burnin <- B * 0.1     # A 10% burn-in
post$mu_R <- post$muR[,-(1:burnin)]
post$mu_T <- post$muT[,-(1:burnin)]

N <- B - burnin      # Number of posterior samples after burn-in
chains <- data.frame(samples = rep(c(1:N, 1:N), each = ncol(dis_data) - 1),
                     group = rep(c("muR", "muT"), each = (ncol(dis_data) - 1) * N),
                     timepoint = paste("Time Point", rep(1:(ncol(dis_data) - 1), 2 * N)),
                     values = c(c(post$mu_R), c(post$mu_T)))

g <- ggplot2::ggplot(chains, ggplot2::aes(samples, values)) +
                     ggplot2::geom_line() +
                     ggplot2::labs(x = "Iterations", y = "Posterior Sample Values") +
                     ggplot2::facet_wrap(group ~ timepoint) +
                     ggplot2::theme(text = ggplot2::element_text(size = 16))

### If we want to calculate the Pr(f2 > 50)
post$f2<- post$f2[-(1:burnin)]
prob <- sum(post$f2 > 50) / (B - burnin)

### Or if we want calculate a 95% credible interval for f2
alpha <- 0.05
f2_cred <- c(quantile(post$f2, alpha / 2),quantile(post$f2, 1 - alpha / 2))

[Package BayesDissolution version 0.1.0 Index]