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/mkinfit.R | 6 | ||||
| -rw-r--r-- | R/mkinpredict.R | 20 | ||||
| -rw-r--r-- | R/multistart.R | 4 | ||||
| -rw-r--r-- | R/parplot.R | 19 | ||||
| -rw-r--r-- | R/plot.mixed.mmkin.R | 4 | ||||
| -rw-r--r-- | R/saem.R | 50 |
9 files changed, 112 insertions, 22 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/mkinfit.R b/R/mkinfit.R index c851fddb..52053685 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -21,7 +21,8 @@ utils::globalVariables(c("name", "time", "value")) #' "FOMC", "DFOP", "HS", "SFORB", "IORE"). If a shorthand name is given, a #' parent only degradation model is generated for the variable with the #' highest value in \code{observed}. -#' @param observed A dataframe with the observed data. The first column called +#' @param observed A dataframe or an object coercible to a dataframe +#' (e.g. a \code{tibble}) with the observed data. The first column called #' "name" must contain the name of the observed variable for each data point. #' The second column must contain the times of observation, named "time". #' The third column must be named "value" and contain the observed values. @@ -292,6 +293,9 @@ mkinfit <- function(mkinmod, observed, # Get the names of observed variables obs_vars <- names(mkinmod$spec) + # Coerce observed data to a dataframe + observed <- as.data.frame(observed) + # Subset observed data with names of observed data in the model and remove NA values observed <- subset(observed, name %in% obs_vars) observed <- subset(observed, !is.na(value)) diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 4c6d7862..2c865cb7 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -172,6 +172,15 @@ mkinpredict.mkinmod <- function(x, } if (solution_type == "deSolve") { + + # We need to make sure that 0 is included in outtimes, + # even if we do not have observed data for time zero, + # because otherwise the initial values are not applied + # to time zero but to the time of the first observation + # in the data + if (min(outtimes) > 0) outtimes_deSolve <- c(0, outtimes) + else outtimes_deSolve <- outtimes + if (!is.null(x$cf) & use_compiled[1] != FALSE) { if (!is.null(x$symbols) & use_symbols) { @@ -182,7 +191,7 @@ mkinpredict.mkinmod <- function(x, out <- deSolve::lsoda( y = odeini, - times = outtimes, + times = outtimes_deSolve, func = lsoda_func, initfunc = "initpar", dllname = x$dll_info[["name"]], @@ -209,7 +218,7 @@ mkinpredict.mkinmod <- function(x, } out <- deSolve::ode( y = odeini, - times = outtimes, + times = outtimes_deSolve, func = mkindiff, parms = odeparms, method = method.ode, @@ -219,6 +228,13 @@ mkinpredict.mkinmod <- function(x, ... ) } + + # Now we need to remove time zero, in case we did not + # have observations for it + if (min(outtimes) > 0) { + out <- out[-1, ] + } + n_out_na <- sum(is.na(out)) if (n_out_na > 0 & na_stop) { cat("odeini:\n") diff --git a/R/multistart.R b/R/multistart.R index aeea2d81..4ba82a43 100644 --- a/R/multistart.R +++ b/R/multistart.R @@ -40,7 +40,7 @@ #' f_mmkin <- mmkin("DFOP", dmta_ds, error_model = "tc", cores = 7, quiet = TRUE) #' f_saem_full <- saem(f_mmkin) #' f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16) -#' parplot(f_saem_full_multi, lpos = "topleft") +#' parplot(f_saem_full_multi, lpos = "topleft", las = 2) #' illparms(f_saem_full) #' #' f_saem_reduced <- update(f_saem_full, no_random_effect = "log_k2") @@ -50,7 +50,7 @@ #' library(parallel) #' cl <- makePSOCKcluster(12) #' f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cluster = cl) -#' parplot(f_saem_reduced_multi, lpos = "topright", ylim = c(0.5, 2)) +#' parplot(f_saem_reduced_multi, lpos = "topright", ylim = c(0.5, 2), las = 2) #' stopCluster(cl) #' } multistart <- function(object, n = 50, diff --git a/R/parplot.R b/R/parplot.R index 3da4b51a..a33112a5 100644 --- a/R/parplot.R +++ b/R/parplot.R @@ -35,9 +35,6 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA, scale = c("best", "median"), lpos = "bottomleft", main = "", ...) { - oldpar <- par(no.readonly = TRUE) - on.exit(par(oldpar, no.readonly = TRUE)) - orig <- attr(object, "orig") orig_parms <- parms(orig) start_degparms <- orig$mean_dp_start @@ -59,11 +56,10 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA, selected <- which(ll > llmin) selected_parms <- all_parms[selected, ] - par(las = 1) if (orig$transformations == "mkin") { degparm_names_transformed <- names(start_degparms) degparm_index <- which(names(orig_parms) %in% degparm_names_transformed) - orig_parms[degparm_names_transformed] <- backtransform_odeparms( + orig_degparms <- backtransform_odeparms( orig_parms[degparm_names_transformed], orig$mmkin[[1]]$mkinmod, transform_rates = orig$mmkin[[1]]$transform_rates, @@ -74,14 +70,17 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA, transform_fractions = orig$mmkin[[1]]$transform_fractions) degparm_names <- names(start_degparms) - names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index])) + orig_parms_back <- orig_parms + orig_parms_back[degparm_index] <- orig_degparms + names(orig_parms_back)[degparm_index] <- degparm_names + orig_parms <- orig_parms_back selected_parms[, degparm_names_transformed] <- t(apply(selected_parms[, degparm_names_transformed], 1, backtransform_odeparms, orig$mmkin[[1]]$mkinmod, transform_rates = orig$mmkin[[1]]$transform_rates, transform_fractions = orig$mmkin[[1]]$transform_fractions)) - colnames(selected_parms)[1:length(degparm_names)] <- degparm_names + colnames(selected_parms)[degparm_index] <- degparm_names } start_errparms <- orig$so@model@error.init @@ -99,6 +98,12 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA, # Boxplots of all scaled parameters selected_scaled_parms <- t(apply(selected_parms, 1, function(x) x / parm_scale)) + i_negative <- selected_scaled_parms <= 0 + parms_with_negative_scaled_values <- paste(names(which(apply(i_negative, 2, any))), collapse = ", ") + if (any(i_negative)) { + warning("There are negative values for ", parms_with_negative_scaled_values, " which are set to NA for plotting") + } + selected_scaled_parms[i_negative] <- NA boxplot(selected_scaled_parms, log = "y", main = main, , ylab = "Normalised parameters", ...) diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index d6c3d0de..058e01a8 100644 --- a/R/plot.mixed.mmkin.R +++ b/R/plot.mixed.mmkin.R @@ -93,7 +93,7 @@ plot.mixed.mmkin <- function(x, nrow.legend = ceiling((length(i) + 1) / ncol.legend), rel.height.legend = 0.02 + 0.07 * nrow.legend, rel.height.bottom = 1.1, - pch_ds = 1:length(i), + pch_ds = c(1:25, 33, 35:38, 40:41, 47:57, 60:90)[1:length(i)], col_ds = pch_ds + 1, lty_ds = col_ds, frame = TRUE, ... @@ -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] } |