Nothing
#' Fitting of dynamic inactivation with MCMC
#'
#' Fits the parameters of an inactivation model to experimental using
#' the Markov Chain Monte Carlo fitting algorithm implemented in
#' the function \code{\link{modMCMC}} of the package \code{\link{FME}}.
#'
#' @param experiment_data data frame with the experimental data to be adjusted.
#' It must have a column named \dQuote{time} and another one named
#' \dQuote{N}.
#' @param simulation_model character identifying the model to be used.
#' @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 starting_points starting values of the parameters for the adjustment.
#' @param upper_bounds named numerical vector defining the upper bounds of the
#' parameters for the adjustment.
#' @param lower_bounds named numerical vector defining the lower bounds of the
#' parameters for the adjustment.
#' @param known_params named numerical vector with the fixed (i.e., not
#' adjustable) model parameters.
#' @param minimize_log logical. If \code{TRUE}, the adjustment is based on the
#' minimization of the error of the logarithmic count.
#' @param ... other arguments for \code{\link{modMCMC}}.
#' @param tol0 numeric. Observations at time 0 make Weibull-based models singular.
#' The time for observatins taken at time 0 are changed for this value.
#'
#' @importFrom dplyr mutate_
#' @importFrom dplyr %>%
#' @importFrom dplyr select_
#' @importFrom FME modMCMC
#' @importFrom lazyeval interp
#'
#' @return A list of class \code{FitInactivationMCMC} with the following items:
#' \itemize{
#' \item \code{modMCMC}: a list of class \code{modMCMC} with the
#' information of the adjustment process.
#' \item best_prediction: a list of class \code{SimulInactivation}
#' with the prediction generated by the best predictor.
#' \item data: a data frame with the data used for the fitting.
#' }
#'
#' @export
#'
#' @examples
#' ## EXAMPLE 1 ------
#' data(dynamic_inactivation) # The example data set is used.
#'
#' get_model_data() # Retrieve the valid model keys.
#'
#' simulation_model <- "Peleg" # Peleg's model will be used
#'
#' model_data <- get_model_data(simulation_model)
#' model_data$parameters # Set the model parameters
#'
#' dummy_temp <- data.frame(time = c(0, 1.25, 2.25, 4.6),
#' temperature = c(70, 105, 105, 70)) # Dummy temp. profile
#'
#' ## Set known parameters and initial points/bounds for unknown ones
#'
#' known_params = c(temp_crit = 100)
#'
#' starting_points <- c(n = 1, k_b = 0.25, N0 = 1e+05)
#' upper_bounds <- c(n = 2, k_b = 1, N0 = 1e6)
#' lower_bounds <- c(n = 0.5, k_b = 0.1, N0 = 1e4)
#'
#' MCMC_fit <- fit_inactivation_MCMC(dynamic_inactivation, simulation_model,
#' dummy_temp, starting_points,
#' upper_bounds, lower_bounds,
#' known_params,
#' niter = 100)
#' # It is recommended to increase niter
#'
#' plot(MCMC_fit)
#' goodness_of_fit(MCMC_fit)
#'
#' ## END EXAMPLE 1 -----
#'
fit_inactivation_MCMC <- function(experiment_data, simulation_model, temp_profile,
starting_points, upper_bounds, lower_bounds,
known_params, ...,
minimize_log = TRUE, tol0 = 1e-5) {
#- Check of the model parameters
# check_model_params(simulation_model, known_params, starting_points, TRUE)
#- Gather the information
model_data <- get_model_data(simulation_model)
if (minimize_log) {
data_for_fit <- mutate_(experiment_data,
logN = interp(~log10(N), N=as.name("N")))
data_for_fit <- select_(data_for_fit, as.name("time"), as.name("logN"))
} else {
data_for_fit <- select_(experiment_data, as.name("time"), as.name("N"))
}
#- Add a small tolerance to data at time 0 to avoid singularities
data_0 <- data_for_fit$time < tol0
data_for_fit[data_0, "time"] <- tol0
#- Call the fitting function
Fit <- modMCMC(f = get_prediction_cost, p = unlist(starting_points),
temp_profile = temp_profile,
data_for_fit = data_for_fit,
lower = unlist(lower_bounds), upper = unlist(upper_bounds),
known_params = known_params,
simulation_model = simulation_model,
...)
#- Prepare the output
prediction_time <- seq(tol0, max(data_for_fit$time), length=100)
best_prediction <- predict_inactivation(simulation_model,
prediction_time,
as.list(c(Fit$bestpar, known_params)),
temp_profile
)
out <- list(modMCMC = Fit,
best_prediction = best_prediction,
data = data_for_fit)
class(out) <- c("FitInactivationMCMC", class(out))
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.