R/predict_MCMC.R In bioinactivation: Mathematical Modelling of (Dynamic) Microbial Inactivation

Documented in predict_inactivation_MCMCsample_dynaFitsample_IsoFitsample_MCMCfit

#' Dynamic Prediction Intervals from a Monte Carlo Adjustment
#'
#' Given a model adjustment of a dynamic microbial inactivation process
#' performed using any of the functions in \code{bioinactivation} calculates
#' probability intervals at each time point using a Monte Carlo method.
#'
#' @param fit_object An object of classes \code{FitInactivationMCMC},
#'        \code{IsoFitInactivation} or \code{FitInactivation}.
#'
#' @param temp_profile data frame with discrete values of the temperature for
#'        each time. It must have one column named \code{time} and another named
#'        \code{temperature} providing discrete values of the temperature at
#'        time points.
#'
#' @param n_simulations a numeric indicating how many Monte Carlo simulations
#'        to perform. \code{100} by default.
#'
#' @param times numeric vector specifying the time points when results are
#'        desired. If \code{NULL}, the times in \code{MCMC_fit$best_prediction} #' are used. \code{NULL} by default. #' #' @param quantiles numeric vector indicating the quantiles to calculate in #' percentage. By default, it is set to c(2.5, 97.5) which generates a #' prediction interval with confidence 0.95. If \code{NULL}, the quantiles #' are not calculated and all the simulations are returned. #' #' @param additional_pars Additional parameters not included in the adjustment (e.g. #' the initial number of microorganism in an isothermal fit). #' #' @importFrom stats quantile median #' #' @return A data frame of class \code{PredInactivationMCMC}. On its first column, #' time at which the calculation has been made is indicated. #' If \code{quantiles = NULL}, the following columns contain the #' results of each simulation. Otherwise, the second and third columns #' provide the mean and median of the simulations at the given time #' point. Following columns contain the quantiles of the results. #' #' @export #' predict_inactivation_MCMC <- function(fit_object, temp_profile, n_simulations = 100, times = NULL, quantiles = c(2.5, 97.5), additional_pars = NULL){ if (is.FitInactivationMCMC(fit_object)){ MCMC_sample <- sample_MCMCfit(fit_object, times, n_simulations) par_sample <- MCMC_sample$par_sample
simulation_model <- MCMC_sample$simulation_model known_pars = MCMC_sample$known_pars
times <- MCMC_sample$times } else if (is.IsoFitInactivation(fit_object)){ isoFit_sample <- sample_IsoFit(fit_object, times, n_simulations) par_sample <- isoFit_sample$par_sample
simulation_model <- isoFit_sample$simulation_model known_pars = isoFit_sample$known_pars
times <- isoFit_sample$times } else if (is.FitInactivation(fit_object)) { dynaFit_sample <- sample_dynaFit(fit_object, times, n_simulations) par_sample <- dynaFit_sample$par_sample
simulation_model <- dynaFit_sample$simulation_model known_pars = dynaFit_sample$known_pars
times <- dynaFit_sample$times } else { stop("fit_object must be a of class FitInactivationMCMC") } result_matrix <- matrix(nrow = length(times), ncol = n_simulations) for (i in 1:n_simulations){ this_pars <- unlist(par_sample[i, ]) this_prediction <- predict_inactivation(simulation_model, times, c(this_pars, known_pars, additional_pars), temp_profile ) result_matrix[, i] <- this_prediction$simulation$N } ## Take means and quantiles from the results if (is.null(quantiles)) { out_matrix <- as.data.frame(cbind(times, result_matrix)) } else { quantile_matrix <- apply(result_matrix, 1, quantile, probs = quantiles/100, na.rm = TRUE) median_matrix <- apply(result_matrix, 1, median, na.rm = TRUE) mean_matrix <- rowMeans(result_matrix, na.rm = TRUE) out_matrix <- as.data.frame(cbind(times, mean = mean_matrix, median = median_matrix, t(quantile_matrix))) } class(out_matrix) <- c("PredInactivationMCMC", class(out_matrix)) return(out_matrix) } #' #' Random sample of the parameters of a FitInactivationMCMC object #' #' Function to be called by \code{\link{predict_inactivation_MCMC}}. Produces #' a random sample of the parameters calculated on the iterations of the Monte #' Carlo simulation. #' #' @param MCMC_fit An object of class \code{FitInactivationMCMC} as generated #' by \code{\link{fit_inactivation_MCMC}}. #' @param n_simulations a numeric indicating how many Monte Carlo simulations #' to perform. #' @param times numeric vector specifying the time points when results are #' desired. If \code{NULL}, the times in \code{MCMC_fit$best_prediction}
#'        are used.
#'
#' @return Returns a list with 4 components.
#'         \itemize{
#'              \item par_sample: data frame with the parameter sample.
#'              \item simulation_model: model key for the simulation
#'              \item known_pars: parameters of the model known
#'              \item times: points where the calculation is made.
#'              }
#'
#' @importFrom dplyr sample_n
#'
sample_MCMCfit <- function(MCMC_fit, times, n_simulations){

## Extract fitted parameters and identify known pars

pars <- as.data.frame(MCMC_fit$modMCMC$pars)
fit_pars <- names(pars)

best_pars <- MCMC_fit$best_prediction$model_parameters
good_indexes <- !(names(best_pars) %in% fit_pars)
known_pars <- unlist(best_pars[good_indexes])

## Extrat the rest of the model data

simulation_model <- MCMC_fit$best_prediction$model

if (is.null(times)) {
times <- MCMC_fit$best_prediction$simulation$time } # temp_values <- MCMC_fit$best_prediction$temp_approximations$temp(times)
#     temp_profile <- data.frame(time = times,
#                                temperature = temp_values)

## Sample from pars and make the simulations

par_sample <- sample_n(pars, n_simulations, replace = TRUE)

## Return the results

out <- list(par_sample = par_sample,
simulation_model = simulation_model,
known_pars = known_pars,
times = times
)
return(out)
}

#'
#' Random sample of the parameters of a IsoFitInactivation object
#'
#' Function to be called by \code{\link{predict_inactivation_MCMC}}. Produces
#' a random sample of the parameters adjusted from isothermal experiments.
#'
#' It is assumed that the parameters follow a Multivariate Normal distribution
#' with the mean and covariance matrix estimated by the adjustment. The function
#' produces a random sample using the function \code{\link{mvrnorm}}.
#'
#' @param iso_fit An object of class \code{FitInactivationMCMC} as generated
#' @param n_simulations a numeric indicating how many Monte Carlo simulations
#'        to perform.
#' @param times numeric vector specifying the time points when results are
#'        desired. If \code{NULL}, an equispaced interval between \code{0}
#'        and the maximum time of the observations with length \code{50}
#'        is used.
#'
#' @return Returns a list with 4 components.
#'         \itemize{
#'              \item par_sample: data frame with the parameter sample.
#'              \item simulation_model: model key for the simulation
#'              \item known_pars: parameters of the model known
#'              \item times: points where the calculation is made.
#'              }
#'
#' @importFrom MASS mvrnorm
#' @importFrom stats coef
#' @importFrom stats vcov
#'
#'
sample_IsoFit <- function(iso_fit, times, n_simulations){

## Take the sample of the parameters

par_sample <- as.data.frame(mvrnorm(n_simulations,
coef(iso_fit$nls), vcov(iso_fit$nls))
)

## Remove rows where any parameter is below 0

bad_index <- apply(par_sample, 1, function(x) any(x<0))

## Identify known pars

fit_pars <- names(par_sample)

best_pars <- iso_fit$parameters good_indexes <- !(names(best_pars) %in% fit_pars) known_pars <- unlist(best_pars[good_indexes]) ## Extrat the rest of the model data simulation_model <- iso_fit$model

if (is.null(times)) {
times <- seq(0, max(iso_fit$data$time), length = 50)
}

## Return the results

out <- list(par_sample = as.data.frame(par_sample),
simulation_model = simulation_model,
known_pars = known_pars,
times = times
)
return(out)

}

#'
#' Random sample of the parameters of a FitInactivation object
#'
#' Function to be called by \code{\link{predict_inactivation_MCMC}}. Produces
#' a random sample of the parameters adjusted from dynamic experiments with non
#' linear regression.
#'
#' It is assumed that the parameters follow a Multivariate Normal distribution
#' with the mean the parameters estimated by \code{\link{modFit}}. The unscaled
#' covariance matrix returned by \code{\link{modFit}} is used.
#'
#' The function produces a random sample using the function \code{\link{mvrnorm}}.
#'
#' @param dynamic_fit An object of class \code{FitInactivationMCMC} as generated
#' @param n_simulations a numeric indicating how many Monte Carlo simulations
#'        to perform.
#' @param times numeric vector specifying the time points when results are
#'        desired. If \code{NULL}, the times in \code{MCMC_fit$best_prediction} #' are used. #' #' @return Returns a list with 4 components. #' \itemize{ #' \item par_sample: data frame with the parameter sample. #' \item simulation_model: model key for the simulation #' \item known_pars: parameters of the model known #' \item times: points where the calculation is made. #' } #' #' @importFrom MASS mvrnorm #' @importFrom stats coef #' @importFrom stats vcov #' sample_dynaFit <- function(dynamic_fit, times, n_simulations){ ## Take the sample of the parameters sP <- summary(dynamic_fit) par_sample <- as.data.frame(mvrnorm(n_simulations, coef(dynamic_fit$fit_results),
sP$cov.unscaled) ) ## Remove rows where any parameter is below 0 bad_index <- apply(par_sample, 1, function(x) any(x<0)) par_sample <- par_sample[!bad_index, ] ## Identify known pars fit_pars <- names(par_sample) best_pars <- dynamic_fit$best_prediction$model_parameters good_indexes <- !(names(best_pars) %in% fit_pars) known_pars <- unlist(best_pars[good_indexes]) ## Extrat the rest of the model data simulation_model <- dynamic_fit$best_prediction$model if (is.null(times)) { times <- seq(0, max(dynamic_fit$data\$time), length = 50)
}

## Return the results

out <- list(par_sample = as.data.frame(par_sample),
simulation_model = simulation_model,
known_pars = known_pars,
times = times
)
return(out)

}


Try the bioinactivation package in your browser

Any scripts or data that you put into this service are public.

bioinactivation documentation built on Aug. 1, 2019, 5:05 p.m.