diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/D24_2014.R | 23 | ||||
| -rw-r--r-- | R/endpoints.R | 1 | ||||
| -rw-r--r-- | R/mkinds.R | 7 | ||||
| -rw-r--r-- | R/plot.mixed.mmkin.R | 2 | ||||
| -rw-r--r-- | R/saem.R | 50 |
5 files changed, 74 insertions, 9 deletions
diff --git a/R/D24_2014.R b/R/D24_2014.R index aaed029c..8a077ff1 100644 --- a/R/D24_2014.R +++ b/R/D24_2014.R @@ -7,6 +7,9 @@ #' context of pesticide registrations, as the use of the data may be #' constrained by data protection regulations. #' +#' The acronyms used for the 2,4-D transformation products are DCP for +#' 2,4-dichlorophenol and DCA for 2,4-dichloroanisole. +#' #' Data for the first dataset are from p. 685. Data for the other four #' datasets were used in the preprocessed versions given in the kinetics #' section (p. 761ff.), with the exception of residues smaller than 1 for DCP @@ -33,5 +36,25 @@ #' DCP = mkinsub("SFO", to = "DCA"), #' DCA = mkinsub("SFO")) #' print(m_D24_2) +#' D24_2014_data <- lapply(D24_2014$ds, function(x) x$data) +#' names(D24_2014_data) <- sapply(D24_2014$ds, function(x) x$title) +#' f_D24_2014_const <- mmkin( +#' models = list( +#' "SFO-SFO-SFO" = m_D24, +#' "DFOP-SFO-SFO" = m_D24_2), +#' data = D24_2014_data, +#' quiet = TRUE) +#' print(f_D24_2014_const) +#' plot(f_D24_2014_const[, 3]) +#' f_D24_2014_tc <- update(f_D24_2014_const, error_model = "tc") +#' print(f_D24_2014_tc) +#' plot(f_D24_2014_tc[, 3]) +#' # For dataset 3, the best fit is obtained using constant +#' # variance and the DFOP-SFO-SFO model +#' AIC( +#' f_D24_2014_const[["SFO-SFO-SFO", 3]], +#' f_D24_2014_const[["DFOP-SFO-SFO", 3]], +#' f_D24_2014_tc[["SFO-SFO-SFO", 3]], +#' f_D24_2014_tc[["DFOP-SFO-SFO", 3]]) #' } "D24_2014" diff --git a/R/endpoints.R b/R/endpoints.R index c3e0a0d3..70a9eef3 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -45,6 +45,7 @@ endpoints <- function(fit, covariates = NULL, covariate_quantile = 0.5) { if (!is.null(fit$covariate_models)) { if (is.null(covariates)) { + # Use covariate quantiles if no explicit covariates are given covariates = as.data.frame( apply(fit$covariates, 2, quantile, covariate_quantile, simplify = FALSE)) @@ -3,7 +3,7 @@ #' @description #' At the moment this dataset class is hardly used in mkin. For example, #' mkinfit does not take mkinds datasets as argument, but works with dataframes -#' such as the on contained in the data field of mkinds objects. Some datasets +#' such as the one contained in the data field of mkinds objects. Some datasets #' provided by this package come as mkinds objects nevertheless. #' #' @importFrom R6 R6Class @@ -87,6 +87,11 @@ print.mkinds <- function(x, data = FALSE, ...) { #' Time normalisation factors are initialised with a value of 1 for each #' dataset if no data are supplied. #' +#' @note +#' Currently, no functions making use of the defined class structure +#' are available in this package. Refer to [D24_2014] for an example +#' dataset in this structure, with some example evaluations. +#' #' @examples #' #' mdsg <- mkindsg$new("Experimental X", experimental_data_for_UBA_2019[6:10]) diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index f05f1110..058e01a8 100644 --- a/R/plot.mixed.mmkin.R +++ b/R/plot.mixed.mmkin.R @@ -157,7 +157,7 @@ plot.mixed.mmkin <- function(x, ifelse(length(x$covariate_models) == 1, "Covariate", "Covariates"), rownames(covariates)) - } + } degparms_pop <- parms(x, covariates = covariates) pop_curves <- TRUE } @@ -50,6 +50,9 @@ utils::globalVariables(c("predicted", "std")) #' initialisation of [saemix::saemixModel] is used for omega.init. #' @param covariates A data frame with covariate data for use in #' 'covariate_models', with dataset names as row names. +#' @param center_covariates Either NA, for no centering, or your +#' preferred function for calculating the center, currently either +#' median or mean. #' @param covariate_models A list containing linear model formulas with one explanatory #' variable, i.e. of the type 'parameter ~ covariate'. Covariates must be available #' in the 'covariates' data frame. @@ -119,9 +122,11 @@ utils::globalVariables(c("predicted", "std")) #' plot(f_saem_dfop_sfo) #' summary(f_saem_dfop_sfo, data = TRUE) #' -#' # The following takes about 6 minutes -#' f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve", -#' nbiter.saemix = c(200, 80)) +#' # The following took about 6 minutes and currently fails with +#' # "Error in deSolve::lsoda(y = odeini, times = outtimes_deSolve, func = lsoda_func, : +#' # illegal input detected before taking any integration steps - see written message" +#' #f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve", +#' # nbiter.saemix = c(200, 80)) #' #' #anova( #' # f_saem_dfop_sfo, @@ -148,6 +153,7 @@ saem.mmkin <- function(object, omega.init = "auto", covariates = NULL, covariate_models = NULL, + center_covariates = c(NA, "median", "mean"), no_random_effect = NULL, error.init = c(1, 1), nbiter.saemix = c(300, 100), @@ -158,6 +164,18 @@ saem.mmkin <- function(object, { call <- match.call() transformations <- match.arg(transformations) + center_covariates <- match.arg(center_covariates) + if (is.na(center_covariates)) { + covariate_centers <- NA + centered_covariates <- NA + } else { + center_covariate_func <- get(center_covariates) + centered_covariates <- data.frame(lapply(covariates, + function(x) x - center_covariate_func(x))) + rownames(centered_covariates) <- rownames(covariates) + covariate_centers <- sapply(covariates, center_covariate_func) + } + m_saemix <- saemix_model(object, verbose = verbose, error_model = error_model, degparms_start = degparms_start, @@ -166,12 +184,14 @@ saem.mmkin <- function(object, transformations = transformations, covariance.model = covariance.model, omega.init = omega.init, - covariates = covariates, + covariates = if (is.na(center_covariates)) covariates else centered_covariates, covariate_models = covariate_models, error.init = error.init, no_random_effect = no_random_effect, ...) - d_saemix <- saemix_data(object, covariates = covariates, verbose = verbose) + d_saemix <- saemix_data(object, + covariates = if (is.na(center_covariates)) covariates else centered_covariates, + verbose = verbose) fit_failed <- FALSE FIM_failed <- NULL @@ -227,6 +247,8 @@ saem.mmkin <- function(object, transform_rates = object[[1]]$transform_rates, transform_fractions = object[[1]]$transform_fractions, covariates = covariates, + centered_covariates = centered_covariates, + covariate_centers = covariate_centers, covariate_models = covariate_models, sm = m_saemix, so = f_saemix, @@ -817,7 +839,8 @@ update.saem.mmkin <- function(object, ..., evaluate = TRUE) { #' each column corresponding to a row of the data frame holding the covariates #' @param covariates A data frame holding covariate values for which to #' return parameter values. Only has an effect if 'ci' is FALSE. -parms.saem.mmkin <- function(object, ci = FALSE, covariates = NULL, ...) { +parms.saem.mmkin <- function(object, ci = FALSE, covariates = NULL, ...) +{ cov.mod <- object$sm@covariance.model n_cov_mod_parms <- sum(cov.mod[upper.tri(cov.mod, diag = TRUE)]) n_parms <- length(object$sm@name.modpar) + @@ -839,12 +862,25 @@ parms.saem.mmkin <- function(object, ci = FALSE, covariates = NULL, ...) { names(estimate) <- rownames(conf.int) + # If confidence intervals are requested, we do not care about covariates if (ci) { return(conf.int) } else { if (is.null(covariates)) { return(estimate) } else { + # If there are covariates in the model, calculate parameter + # estimates for the corresponding covariate values passed + # in 'covariates'. + covariate_centers <- object$covariate_centers + parm_covariates <- covariates + if (!is.na(covariate_centers[1])) { + # If covariates were centered by the saem function, we need to center + # user specified covariates as well + for (i in names(covariates)) { + parm_covariates[i] <- covariates[i] - object$covariate_centers[i] + } + } est_for_cov <- matrix(NA, nrow = length(object$sm@name.modpar), ncol = nrow(covariates), dimnames = (list(object$sm@name.modpar, rownames(covariates)))) @@ -856,7 +892,7 @@ parms.saem.mmkin <- function(object, ci = FALSE, covariates = NULL, ...) { beta_degparm_name <- paste0("beta_", covariate, "(", deg_parm_name, ")") est_for_cov[deg_parm_name, ] <- estimate[deg_parm_name] + - covariates[[covariate]] * estimate[beta_degparm_name] + parm_covariates[[covariate]] * estimate[beta_degparm_name] } else { est_for_cov[deg_parm_name, ] <- estimate[deg_parm_name] } |