diff --git a/DESCRIPTION b/DESCRIPTION index e04dc85..d1f1985 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: PKPDmap Type: Package Title: MAP Bayesian estimates -Version: 1.0.1 -Date: 2024-06-10 +Version: 1.0.3 +Date: 2025-06-18 Author: Ron Keizer, Jasmine Hughes, Kara Woo Maintainer: Ron Keizer Description: Obtain MAP Bayesian estimates based on individual data and @@ -12,6 +12,6 @@ Depends: utils, stats Imports: PKPDsim, numDeriv, mvtnorm Suggests: testthat (>= 3.0.0) Remotes: InsightRX/PKPDsim -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Encoding: UTF-8 Config/testthat/edition: 3 diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..e5adbb6 --- /dev/null +++ b/NEWS.md @@ -0,0 +1,9 @@ +# Changes in master compared to last release (20240610) + +- many code chunks in main estimation function now separated off to their own function +- improved handling of floating points during warning check +- The print method for map_estimates() now uses + `sqrt(diag(PKPDsim::triangle_to_full(x$prior$omega$full)))` instead of + `sqrt(diag(PKPDsim::triangle_to_full(x$prior$omega)))` for etas +- fix for mixture models: ipreds for mixture models were found not to be correct +- censoring information is now returned with the fit object diff --git a/PKPDmap.Rproj b/PKPDmap.Rproj index 01e8c39..744633a 100755 --- a/PKPDmap.Rproj +++ b/PKPDmap.Rproj @@ -1,5 +1,5 @@ Version: 1.0 -ProjectId: 42e31d31-a19a-49a8-abcd-daac065467eb +ProjectId: 6835582e-ba67-4423-80dd-8a468b12df41 RestoreWorkspace: Default SaveWorkspace: Default diff --git a/R/calc_residuals.R b/R/calc_residuals.R new file mode 100644 index 0000000..46e3b19 --- /dev/null +++ b/R/calc_residuals.R @@ -0,0 +1,110 @@ +#' Function to calculate residuals from a fit, +#' based on a fitted parameter set, and data object. +#' +#' @inheritParams get_map_estimates +#' @param obj temporary storage object in `get_map_estimates()` +#' @param parameters_population parameters for population predictions +#' @param omega_full full omega matrix +#' @param transf transformation function +#' @param A_init_population init vector for model state (population) +#' @param A_init_individual init vector for model state (individual) +#' @param censoring_idx vector with indices for censoring +#' @param data_before_init data.frame with data before initial dose +#' +calc_residuals <- function( + obj, + data, + model, + parameters_population, + covariates, + regimen, + omega_full, + error, + weights, + transf, + A_init_population, + A_init_individual, + t_init = 0, + iov_bins = NULL, + output_include = c(), + int_step_size = 0.01, + censoring = NULL, + censoring_idx = NULL, + data_before_init = NULL, + ltbs = FALSE, + ... +) { + + ## Observation vectors + t_obs <- c(data_before_init$t, data$t) + obs_type <- c(data_before_init$obs_type, data$obs_type) + + ## Perform simulations for ipred and pred + suppressMessages({ + ## After fitting individual parameters, don't pass the mixture group to the simulator + ## (so `mixture=NULL`), otherwise `sim()` will use the population value for the + ## specified group, and not the individual fitted parameter. + sim_ipred <- PKPDsim::sim_ode( + ode = model, + parameters = obj$parameters, + mixture_group = NULL, + covariates = covariates, + n_ind = 1, + int_step_size = int_step_size, + regimen = regimen, + t_obs = t_obs, + obs_type = obs_type, + only_obs = TRUE, + checks = FALSE, + A_init = A_init_individual, + iov_bins = iov_bins, + output_include = output_include, + t_init = t_init, + ... + ) + }) + suppressMessages({ + sim_pred <- PKPDsim::sim_ode( + ode = model, + parameters = parameters_population, + mixture_group = NULL, + covariates = covariates, + n_ind = 1, + int_step_size = int_step_size, + regimen = regimen, + t_obs = c(data_before_init$t, t_obs), + obs_type = c(data_before_init$obs_type, data$obs_type), + only_obs = TRUE, + checks = FALSE, + iov_bins = iov_bins, + A_init = A_init_population, + t_init = t_init, + ... + ) + }) + + ## Parse the residuals from the predictions and the data, and add to obj + obj <- parse_residuals_from_predictions( + obj, + sim_ipred, + sim_pred, + data, + omega_full, + transf, + error, + censoring, + censoring_idx, + data_before_init, + weights + ) + + ## Add covariates and parameters to obj + if(output_include$covariates && !is.null(covariates)) { + obj$covariates_time <- sim_ipred[!duplicated(sim_ipred$t), names(covariates)] + } + if(output_include$parameters) { + obj$parameters_time <- sim_ipred[!duplicated(sim_ipred$t), names(parameters)] + } + + obj +} \ No newline at end of file diff --git a/R/check_inputs.R b/R/check_inputs.R new file mode 100644 index 0000000..5699477 --- /dev/null +++ b/R/check_inputs.R @@ -0,0 +1,17 @@ +#' Some basic input checking +#' +#' @inheritParams get_map_estimates +#' +check_inputs <- function(model, data, parameters, omega, regimen, censoring, type) { + if(tolower(type) %in% c("map", "pls")) { + if(is.null(model) || is.null(data) || is.null(parameters) || is.null(omega) || is.null(regimen)) { + stop("The 'model', 'data', 'omega', 'regimen', and 'parameters' arguments are required.") + } + } + if(!is.null(censoring) && !inherits(censoring, "character")) { + stop("Censoring argument requires label specifying column in dataset with censoring info.") + } + if(!("function" %in% class(model))) { + stop("The 'model' argument requires a function, e.g. a model defined using the new_ode_model() function from the PKPDsim package.") + } +} \ No newline at end of file diff --git a/R/get_individual_parameters_from_fit.R b/R/get_individual_parameters_from_fit.R new file mode 100644 index 0000000..25c7a8e --- /dev/null +++ b/R/get_individual_parameters_from_fit.R @@ -0,0 +1,22 @@ +#' Get the individual parameters from the fitted coefficients (etas) +#' +#' @param fit fit object +#' @param parameters list of population parameters +#' +get_individual_parameters_from_fit <- function( + fit, + parameters, + nonfixed = NULL, + as_eta = NULL +) { + par <- parameters + for(i in seq(nonfixed)) { + key <- nonfixed[i] + if(key %in% as_eta) { + par[[key]] <- as.numeric(fit$coef[i]) + } else { + par[[key]] <- as.numeric(as.numeric(par[[key]]) * exp(as.numeric(fit$coef[i]))) + } + } + par +} diff --git a/R/get_map_estimates.R b/R/get_map_estimates.R index 17dbb22..046e833 100755 --- a/R/get_map_estimates.R +++ b/R/get_map_estimates.R @@ -129,26 +129,17 @@ get_map_estimates <- function( output_include = list(covariates = FALSE, parameters = FALSE), ...) { - ## Handle weighting of priors, allow for some presets but can - ## also be set manually using `weight_prior` - if(is.null(weight_prior) || is.na(weight_prior)) { - weight_prior <- 1 - } - weight_prior <- weight_prior^2 # 3x larger IIV on SD scale + ## get prior weight for scaling of variance term + weight_prior_var <- parse_weight_prior(weight_prior, type) + + ## pick likelihood function to use, and perform some checks calc_ofv <- calc_ofv_map - if(weight_prior == 0) { + if(weight_prior_var == 0) { calc_ofv <- calc_ofv_ls } - if(tolower(type) %in% c("map", "pls")) { - if(is.null(model) || is.null(data) || is.null(parameters) || is.null(omega) || is.null(regimen)) { - stop("The 'model', 'data', 'omega', 'regimen', and 'parameters' arguments are required.") - } - } - if(tolower(type) == "pls") { - weight_prior <- 0.001 - ## RK: Empirically determined to be a good penalty. - ## In principle, PLS is just MAP with very flat priors - } + + ## Misc checks: + check_inputs(model, data, parameters, omega, regimen, censoring, type) if(tolower(type) %in% c("ls")) { if(is.null(model) || is.null(data) || is.null(regimen)) { stop("The 'model', 'data', and 'parameters' arguments are required.") @@ -156,43 +147,19 @@ get_map_estimates <- function( calc_ofv <- calc_ofv_ls error <- list(prop = 0, add = 1) } - if(is.null(obs_type_label)) { - data$obs_type <- 1 - } else { - data$obs_type <- data[[obs_type_label]] - } - if(!is.null(error)) { ## safety checks - if(is.null(error$prop)) error$prop <- 0 - if(is.null(error$add)) error$add <- 0 - if(is.null(error$exp)) error$exp <- 0 - } else { - warning("No residual error specified, using: 10% proportional + 0.1 additive error.") - list(prop = 0.1, add = 0.1, exp = 0) - } - if(!is.null(censoring) && !inherits(censoring, "character")) { - stop("Censoring argument requires label specifying column in dataset with censoring info.") - } - if(!("function" %in% class(model))) { - stop("The 'model' argument requires a function, e.g. a model defined using the new_ode_model() function from the PKPDsim package.") - } - if(!all(unlist(cols) %in% names(data))) { - stop("Expected column names were not found in data. Please use 'cols' argument to specify column names for independent and dependent variable.") - } - if("PKPDsim" %in% class(data)) { - if("comp" %in% names(data)) { - data <- data[data$comp == "obs",] - # data <- data[!duplicated(data$t),] - data$evid <- 0 - } - } - colnames(data) <- tolower(colnames(data)) - sig <- round(-log10(int_step_size)) - if("evid" %in% colnames(data)) { - data <- data[data$evid == 0,] + error <- parse_error(error) + if(is.null(attr(model, "cpp")) || !attr(model, "cpp")) { + warning("Provided model is not PKPDsim model, using generic likelihood function.") + ll_func <- ll_func_generic } + + ## Parse input data + data <- parse_input_data( + data = data, + cols = cols, + obs_type_label = obs_type_label + ) data_before_init <- data.frame() - data <- data[order(data$t, data$obs_type),] - y_orig <- data$y if(!allow_obs_before_dose) { filt_before_init <- data$t < min(regimen$dose_times) data_before_init <- data[filt_before_init,] @@ -203,82 +170,37 @@ get_map_estimates <- function( data <- data[!filt_before_init,] } } - t_obs <- data$t - if(any(duplicated(paste(t_obs, data$obs_type, sep = "_")))) { - message("Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument.") - } - if(!is.null(weights)) { - if(length(weights) != length(t_obs)) { - stop("Vector of weights of different size than observation vector!") - } - } else { - weights <- rep(1, length(data$y)) - } - if(sum(unlist(error)) == 0) { - stop("No residual error model specified, or residual error is 0.") - } + + ## Parse observation times and weights + t_obs <- parse_t_obs(data) + weights <- parse_weights(weights, t_obs) ## Log-transform both sides? if(is.null(ltbs)) { ltbs <- FALSE if(!is.null(attr(model, "ltbs")) && attr(model, "ltbs")) ltbs <- TRUE } - if(ltbs) { - transf <- function(x) log(x) - } else { - transf <- function(x) x - } - - if(is.null(attr(model, "cpp")) || !attr(model, "cpp")) { - ll_func <- ll_func_generic - } + transf <- ifelse(ltbs, function(x) log(x), function(x) x) ## check if fixed parameter actually in parameter list - if(length(intersect(fixed, names(parameters))) != length(fixed)) { - warning("Warning: not all fixed parameters were found in parameter set!\n") - } - fixed <- names(parameters)[names(parameters) %in% fixed] - if(length(fixed) == 0) { - fixed <- NULL - } - - ## fix etas - nonfixed <- names(parameters)[is.na(match(names(parameters), fixed))] - n_nonfix <- length(nonfixed) - eta <- list() - for(i in seq(nonfixed)) { - eta[[paste0("eta", sprintf("%02d", i))]] <- 0 - } - fix <- NULL - - if(inherits(omega, "matrix")) { - omega_full <- omega # dummy om matrix - } else { - omega_full <- PKPDsim::triangle_to_full(omega) # dummy om matrix - } - om_nonfixed <- PKPDsim::triangle_to_full(omega) - if(nrow(om_nonfixed) < (length(parameters) - length(fixed))) { - msg <- "Provided omega matrix is smaller than expected based on the number of model parameters. Either fix some parameters or increase the size of the omega matrix.\n" - msg <- c(msg, - paste0("Non-fixed omegas: ", paste(om_nonfixed, collapse=", "), "\n"), - paste0("Parameters: ", paste(parameters, collapse=", "), "\n"), - paste0("Fixed: ", paste(fix, collapse=","), "\n")) - stop(msg) - } + fixed <- parse_fixed(fixed, parameters) - ## check if censoring code needs to be used - censoring_idx <- NULL - if(!is.null(censoring)) { - if(any(data[[tolower(censoring)]] != 0)) { - censoring_idx <- data[[tolower(censoring)]] != 0 - if(verbose) message("One or more values in data are censored, including censoring in likelihood.") - } else { - if(verbose) message("Warning: censoring specified, but no censored values in data.") - } - } + ## Omega matrix and etas + omega <- parse_omega_matrix( + omega, + parameters, + fixed + ) + ## Check if censoring code needs to be used + censoring_idx <- parse_censoring( + censoring = censoring, + data = data, + verbose = verbose + ) + ################################################# - ## create simulation design up-front: + ## create simulation design template ################################################# mixture_group <- 1 suppressMessages({ @@ -304,7 +226,6 @@ get_map_estimates <- function( ) }) - omega_full_est <- omega_full[1:n_nonfix, 1:n_nonfix] mixture_obj <- NULL if(!is.null(attr(model, "mixture"))) { mixture <- attr(model, "mixture") @@ -317,7 +238,7 @@ get_map_estimates <- function( par_mix[[mix_par]] <- mix_par_values[i] fits[[i]] <- mle_wrapper( minuslogl = ll_func, - start = eta, + start = omega$eta, method = method, optimizer = optimizer, control = control, @@ -330,14 +251,14 @@ get_map_estimates <- function( model = model, regimen = regimen, error = error, - omega_full = omega_full_est / weight_prior, - omega_inv = solve(omega_full_est / weight_prior), - omega_eigen = sum(log(eigen(omega_full_est)$values)), - nonfixed = nonfixed, + omega_full = omega$est / weight_prior_var, + omega_inv = solve(omega$est / weight_prior_var), + omega_eigen = sum(log(eigen(omega$est / weight_prior_var)$values)), + nonfixed = omega$nonfixed, transf = transf, weights = weights, - weight_prior = weight_prior, - sig = sig, + weight_prior = weight_prior_var, + sig = round(-log10(int_step_size)), as_eta = as_eta, censoring_idx = censoring_idx, censoring_label = censoring, @@ -372,7 +293,7 @@ get_map_estimates <- function( output <- tryCatch({ fit <- mle_wrapper( ll_func, - start = eta, + start = omega$eta, method = method, optimizer = optimizer, control = control, @@ -385,14 +306,14 @@ get_map_estimates <- function( model = model, regimen = regimen, error = error, - nonfixed = nonfixed, + nonfixed = omega$nonfixed, transf = transf, - omega_full = omega_full_est / weight_prior, - omega_inv = solve(omega_full_est / weight_prior), - omega_eigen = sum(log(eigen(omega_full_est / weight_prior)$values)), + omega_full = omega$est / weight_prior_var, + omega_inv = solve(omega$est / weight_prior_var), + omega_eigen = sum(log(eigen(omega$est / weight_prior_var)$values)), weights = weights, - weight_prior = weight_prior, - sig = sig, + weight_prior = weight_prior_var, + sig = round(-log10(int_step_size)), as_eta = as_eta, censoring_idx = censoring_idx, censoring_label = censoring, @@ -410,68 +331,34 @@ get_map_estimates <- function( }) if("error" %in% class(output)) return(output) } - cf <- fit$coef - par <- parameters - for(i in seq(nonfixed)) { - key <- nonfixed[i] - if(key %in% as_eta) { - par[[key]] <- as.numeric(cf[i]) - } else { - par[[key]] <- as.numeric(as.numeric(par[[key]]) * exp(as.numeric(cf[i]))) - } - } - obj <- list(fit = fit, mixture = mixture_obj) + + ## Create output object with individual parameters + obj <- list( + fit = fit, + mixture = mixture_obj, + parameters = get_individual_parameters_from_fit( + fit, + parameters, + omega$nonfixed, + as_eta + ) + ) ################################################# - ## Non-parametric estimation + ## Non-parametric estimation (hybrid, based on + ## MAP fit as initial estimates) ################################################# if(type == "np_hybrid") { - obj$parameters_map <- par ## keep MAP estimates - np_settings <- replace_list_elements(np_settings_default, np_settings) - if(is.null(np_settings$error)) { # if no specific error magnitude is specified for NP, just use same as used for MAP - np_settings$error <- error - } - if(np_settings$grid_adaptive) { # do a first pass with a broad grid - pars_grid <- create_grid_around_parameters( - parameters = par, - span = np_settings$grid_span_adaptive, - exponential = np_settings$grid_exponential_adaptive, - grid_size = np_settings$grid_size_adaptive) - np <- get_np_estimates( - parameter_grid = pars_grid, - error = np_settings$error, - model = model, - regimen = regimen, - data = data$y, - t_obs = data$t, - covariates = covariates, - weights = weights - ) - # take the estimates with highest probability as starting point for next grid - tmp <- np$prob[order(-np$prob$like),][1,] - par <- as.list(tmp[1:length(np$parameters)]) - } - pars_grid <- create_grid_around_parameters( - parameters = par, - span = np_settings$grid_span, - exponential = np_settings$grid_exponential, - grid_size = np_settings$grid_size) - np <- get_np_estimates( - parameter_grid = pars_grid, - error = np_settings$error, + obj <- np_fit_wrapper( + obj = obj, model = model, regimen = regimen, - data = data$y, - t_obs = data$t, + data = data, covariates = covariates, - weights = weights + weights = weights, + np_settings = np_settings ) - for(i in 1:length(par)) { - par[[i]] <- np$parameters[[i]] - } - obj$np <- list(prob = np$prob) } - obj$parameters <- par ################################################# ## Add g.o.f. info @@ -483,7 +370,7 @@ get_map_estimates <- function( dose = regimen$dose_amts[1], interval = regimen$interval[1], model = model, - parameters = obj$parameters, + parameters = obj$parameters, # individual covariates = covariates, map = steady_state_analytic$map, n_transit_compartments = PKPDsim::ifelse0(steady_state_analytic$n_transit_compartments, FALSE), @@ -494,7 +381,7 @@ get_map_estimates <- function( dose = regimen$dose_amts[1], interval = regimen$interval[1], model = model, - parameters = parameters, + parameters = parameters, # population covariates = covariates, map = steady_state_analytic$map, n_transit_compartments = PKPDsim::ifelse0(steady_state_analytic$n_transit_compartments, FALSE), @@ -504,107 +391,55 @@ get_map_estimates <- function( A_init_pred <- A_init A_init_ipred <- A_init } - suppressMessages({ - ## After fitting individual parameters, don't pass the mixture group to the simulator - ## (so `mixture=NULL`), otherwise `sim()` will use the population value for the - ## specified group, and not the individual fitted parameter. - sim_ipred <- PKPDsim::sim_ode( - ode = model, - parameters = par, - mixture_group = NULL, - covariates = covariates, - n_ind = 1, - int_step_size = int_step_size, - regimen = regimen, - t_obs = c(data_before_init$t, t_obs), - obs_type = c(data_before_init$obs_type, data$obs_type), - only_obs = TRUE, - checks = FALSE, - A_init = A_init_ipred, - iov_bins = iov_bins, - output_include = output_include, - t_init = t_init, - ... - ) - }) - suppressMessages({ - sim_pred <- PKPDsim::sim_ode( - ode = model, - parameters = parameters, - mixture_group = NULL, - covariates = covariates, - n_ind = 1, - int_step_size = int_step_size, - regimen = regimen, - t_obs = c(data_before_init$t, t_obs), - obs_type = c(data_before_init$obs_type, data$obs_type), - only_obs = TRUE, - checks = FALSE, - iov_bins = iov_bins, - A_init = A_init_pred, - t_init = t_init, - ... - ) - }) - ipred <- sim_ipred$y - pred <- sim_pred$y - w_ipred <- sqrt(error$prop[data$obs_type]^2 * transf(ipred)^2 + error$add[data$obs_type]^2) - w_pred <- sqrt(error$prop[data$obs_type]^2 * transf(pred)^2 + error$add[data$obs_type]^2) - if(!isTRUE(all.equal(data$t, sim_ipred$t)) && all(data$obs_type == sim_ipred$obs_type)) { - warning("Mismatch in times and observation typese between input data and predictions. Be careful interpreting results from fit.") - } - y <- data$y - prob <- list(par = c(mvtnorm::pmvnorm(cf, mean=rep(0, length(cf)), - sigma = omega_full[1:length(cf), 1:length(cf)])), - data = stats::pnorm(transf(y) - transf(ipred), mean = 0, sd = w_ipred)) - obj$res <- (transf(y) - transf(pred)) - obj$weights <- c(rep(0, length(data_before_init$t)), weights) - obj$wres <- (obj$res / w_pred) * obj$weights - obj$cwres <- obj$res / sqrt(abs(cov(transf(pred), transf(y_orig)))) * c(rep(0, nrow(data_before_init)), obj$weights) - # Note: in NONMEM CWRES is on the population level, so can't really compare. NONMEM calls this CIWRES, it seems. - obj$ires <- (transf(y_orig) - transf(ipred)) - obj$iwres <- (obj$ires / w_ipred) - if(is.null(censoring)) { - obj$censoring <- rep(0, length(y)) - } else { - obj$censoring <- data[[censoring]] - if(any(censoring_idx)) { # turn probabilities into "IWRES"-equivalent for censored data - obj$iwres[censoring_idx] <- calc_res_from_prob(prob$data[censoring_idx]) - obj$ires[censoring_idx] <- obj$iwres[censoring_idx] * w_ipred[censoring_idx] - ## if we would calculate the likelihood for the data population parameters given the data, - ## we could also calculate the the equivalents for cwres, wres, and res. However we - ## currently don't have a need to calculate. So setting to NA to avoid wrong interpretation. - obj$cwres[censoring_idx] <- NA_real_ - obj$wres[censoring_idx] <- NA_real_ - obj$res[censoring_idx] <- NA_real_ - } - } - obj$iwres_weighted <- obj$iwres * obj$weights - obj$pred <- pred - obj$ipred <- ipred - obj$prob <- prob - obj$dv <- y_orig - obj$obs_type <- sim_ipred$obs_type - if(output_include$covariates && !is.null(covariates)) { - obj$covariates_time <- sim_ipred[!duplicated(sim_ipred$t), names(covariates)] - } - if(output_include$parameters) { - obj$parameters_time <- sim_ipred[!duplicated(sim_ipred$t), names(parameters)] - } + obj <- calc_residuals( + obj = obj, + data = data, + model = model, + parameters_population = parameters, + covariates = covariates, + int_step_size = int_step_size, + regimen = regimen, + omega_full = omega$full, + error = error, + weights = weights, + transf = transf, + A_init_population = A_init_pred, + A_init_individual = A_init_ipred, + iov_bins = iov_bins, + output_include = output_include, + t_init = t_init, + censoring = censoring, + censoring_idx = censoring_idx, + data_before_init = data_before_init, + ltbs = ltbs, + ... + ) } + + ## Add variance-covariance matrix to output object obj$vcov_full <- get_varcov_matrix( obj$fit$vcov, - fallback = omega_full + fallback = omega$full ) if(inherits(obj$vcov_full, "matrix")) { obj$vcov <- obj$vcov_full[t(!upper.tri(obj$vcov_full))] } - obj$mahalanobis <- get_mahalanobis(y, ipred, w_ipred, ltbs) + + ## Add Mahalanobis distance to output obj + obj$mahalanobis <- get_mahalanobis( + data$y, + obj$ipred, + obj$w_ipred, + ltbs + ) + + ## Add prior to output obj obj$prior <- list( parameters = parameters, omega = omega, fixed = fixed ) + class(obj) <- c(class(obj), "map_estimates") obj } diff --git a/R/np_fit_wrapper.R b/R/np_fit_wrapper.R new file mode 100644 index 0000000..e197084 --- /dev/null +++ b/R/np_fit_wrapper.R @@ -0,0 +1,67 @@ +#' Wrapper for non-parametric fit +#' +#' @param obj fit object +#' @inheritParams get_map_estimates +#' +np_fit_wrapper <- function( + obj, + model, + regimen, + data, + covariates = NULL, + weights = 1, + np_settings = list() +) { + par <- obj$parameters + obj$parameters_map <- par ## keep copy of MAP estimates + np_settings <- replace_list_elements( + np_settings_default, + np_settings + ) + if(is.null(np_settings$error)) { # if no specific error magnitude is specified for NP, just use same as used for MAP + np_settings$error <- error + } + if(np_settings$grid_adaptive) { # do a first pass with a broad grid + pars_grid <- create_grid_around_parameters( + parameters = par, + span = np_settings$grid_span_adaptive, + exponential = np_settings$grid_exponential_adaptive, + grid_size = np_settings$grid_size_adaptive + ) + np <- get_np_estimates( + parameter_grid = pars_grid, + error = np_settings$error, + model = model, + regimen = regimen, + data = data$y, + t_obs = data$t, + covariates = covariates, + weights = weights + ) + # take the estimates with highest probability as starting point for next grid + tmp <- np$prob[order(-np$prob$like),][1,] + par <- as.list(tmp[1:length(np$parameters)]) + } + pars_grid <- create_grid_around_parameters( + parameters = par, + span = np_settings$grid_span, + exponential = np_settings$grid_exponential, + grid_size = np_settings$grid_size + ) + np <- get_np_estimates( + parameter_grid = pars_grid, + error = np_settings$error, + model = model, + regimen = regimen, + data = data$y, + t_obs = data$t, + covariates = covariates, + weights = weights + ) + for(i in 1:length(par)) { + par[[i]] <- np$parameters[[i]] + } + obj$np <- list(prob = np$prob) + obj$parameters <- par + obj +} \ No newline at end of file diff --git a/R/parse_censoring.R b/R/parse_censoring.R new file mode 100644 index 0000000..0535bfc --- /dev/null +++ b/R/parse_censoring.R @@ -0,0 +1,21 @@ +#' Check if censoring functionality needs to be used, +#' and create index for censored data +#' +#' @inheritParams get_map_estimates +#' +parse_censoring <- function( + censoring, + data, + verbose = FALSE +) { + censoring_idx <- NULL + if(!is.null(censoring)) { + if(any(data[[tolower(censoring)]] != 0)) { + censoring_idx <- data[[tolower(censoring)]] != 0 + if(verbose) message("One or more values in data are censored, including censoring in likelihood.") + } else { + if(verbose) message("Warning: censoring specified, but no censored values in data.") + } + } + censoring_idx +} \ No newline at end of file diff --git a/R/parse_error.R b/R/parse_error.R new file mode 100644 index 0000000..9a2abbf --- /dev/null +++ b/R/parse_error.R @@ -0,0 +1,17 @@ +#' Parse specified error magnitudes +#' +#' @inheritParams get_map_estimates +#' +parse_error <- function(error) { + if(!is.null(error)) { ## ensure all error types are specified and not NULL + if(is.null(error$prop)) error$prop <- 0 + if(is.null(error$add)) error$add <- 0 + if(is.null(error$exp)) error$exp <- 0 + if(sum(unlist(error)) == 0) { + stop("No residual error model specified, or residual error is 0.") + } + } else { + stop("No residual error specified, cannot perform fit.") + } + error +} diff --git a/R/parse_fixed.R b/R/parse_fixed.R new file mode 100644 index 0000000..a238242 --- /dev/null +++ b/R/parse_fixed.R @@ -0,0 +1,16 @@ +#' Parse fixed parameters +#' +#' @inheritParams get_map_estimates +#' +parse_fixed <- function(fixed, parameters) { + if(is.null(fixed)) return() + fixed <- unique(fixed) + if(length(intersect(fixed, names(parameters))) != length(fixed)) { + warning("Warning: not all fixed parameters were found in parameter set!\n") + } + fixed <- names(parameters)[names(parameters) %in% fixed] + if(length(fixed) == 0) { + fixed <- NULL + } + fixed +} \ No newline at end of file diff --git a/R/parse_input_data.R b/R/parse_input_data.R new file mode 100644 index 0000000..62bd2fd --- /dev/null +++ b/R/parse_input_data.R @@ -0,0 +1,31 @@ +#' Parse input data for fit +#' +#' @inheritParams get_map_estimates +#' +parse_input_data <- function( + data, + cols = list(x = "t", y = "y"), + obs_type_label = NULL +) { + colnames(data) <- tolower(colnames(data)) + if(!all(tolower(unlist(cols)) %in% names(data))) { + stop("Expected column names were not found in data. Please use 'cols' argument to specify column names for independent and dependent variable.") + } + if(inherits(data, "PKPDsim")) { + # when fitting directly from PKPDsim-generated data: + if("comp" %in% names(data)) { + data <- data[data$comp == "obs",] + data$evid <- 0 + } + } + ## Parse data to include label for obs_type + if(is.null(obs_type_label)) { + data$obs_type <- 1 + } else { + data$obs_type <- data[[obs_type_label]] + } + if("evid" %in% colnames(data)) { + data <- data[data$evid == 0,] + } + data[order(data$t, data$obs_type),] +} \ No newline at end of file diff --git a/R/parse_omega_matrix.R b/R/parse_omega_matrix.R new file mode 100644 index 0000000..b66938b --- /dev/null +++ b/R/parse_omega_matrix.R @@ -0,0 +1,42 @@ +#' Parse omega matrix into full matrix, and +#' perform some checks +#' +#' @inheritParams get_map_estimates +#' +parse_omega_matrix <- function( + omega, + parameters, + fixed +) { + nonfixed <- names(parameters)[is.na(match(names(parameters), fixed))] + n_nonfix <- length(nonfixed) + eta <- list() + for(i in seq(nonfixed)) { + eta[[paste0("eta", sprintf("%02d", i))]] <- 0 + } + if(inherits(omega, "matrix")) { + omega_full <- omega # dummy om matrix + } else { + omega_full <- PKPDsim::triangle_to_full(omega) # dummy om matrix + } + om_nonfixed <- PKPDsim::triangle_to_full(omega) + if(nrow(om_nonfixed) < (length(parameters) - length(fixed))) { + msg <- "Provided omega matrix is smaller than expected based on the number of model parameters. Either fix some parameters or increase the size of the omega matrix.\n" + msg <- c(msg, + paste0("Non-fixed omegas: ", paste(om_nonfixed, collapse=", "), "\n"), + paste0("Parameters: ", paste(parameters, collapse=", "), "\n"), + paste0("Fixed: ", paste(fixed, collapse=","), "\n")) + stop(msg) + } + omega_full_est <- omega_full[ + 1:n_nonfix, + 1:n_nonfix + ] + list( + full = omega_full, # full omega block + est = omega_full_est, # full omega block but only for non-fixed parameters, + fixed = fixed, # vector of fixed parameters + nonfixed = nonfixed, # vector of non-fixed (estimaterd) parameters + eta = eta # list of etas + ) +} \ No newline at end of file diff --git a/R/parse_residuals_from_predictions.R b/R/parse_residuals_from_predictions.R new file mode 100644 index 0000000..38f4861 --- /dev/null +++ b/R/parse_residuals_from_predictions.R @@ -0,0 +1,69 @@ +#' Parse simulated data and extract residuals +#' +#' @inheritParams calc_residuals +#' @param sim_ipred data frame with simulated individual predictions +#' @param sim_pred data frame with simulated population predictions +#' +parse_residuals_from_predictions <- function( + obj, + sim_ipred, + sim_pred, + data, + omega_full, + transf, + error, + censoring, + censoring_idx, + data_before_init, + weights +) { + ipred <- sim_ipred$y + pred <- sim_pred$y + w_ipred <- sqrt(error$prop[data$obs_type]^2 * transf(ipred)^2 + error$add[data$obs_type]^2) + w_pred <- sqrt(error$prop[data$obs_type]^2 * transf(pred)^2 + error$add[data$obs_type]^2) + if(!isTRUE(all.equal(data$t, sim_ipred$t)) && all(data$obs_type == sim_ipred$obs_type)) { + warning("Mismatch in times and observation typese between input data and predictions. Be careful interpreting results from fit.") + } + y <- data$y + cf <- obj$fit$coef + prob <- list( + par = c( + mvtnorm::pmvnorm( + cf, + mean=rep(0, length(cf)), + sigma = omega_full[1:length(cf), 1:length(cf)] + ) + ), + data = stats::pnorm(transf(y) - transf(ipred), mean = 0, sd = w_ipred) + ) + obj$res <- (transf(y) - transf(pred)) + obj$weights <- c(rep(0, length(data_before_init$t)), weights) + obj$wres <- (obj$res / w_pred) * obj$weights + obj$cwres <- obj$res / sqrt(abs(cov(transf(pred), transf(y)))) * c(rep(0, nrow(data_before_init)), obj$weights) + # Note: in NONMEM CWRES is on the population level, so can't really compare. NONMEM calls this CIWRES, it seems. + obj$ires <- (transf(y) - transf(ipred)) + obj$iwres <- (obj$ires / w_ipred) + obj$w_ipred <- w_ipred + if(is.null(censoring)) { + obj$censoring <- rep(0, length(y)) + } else { + obj$censoring <- data[[censoring]] + if(any(censoring_idx)) { + obj$iwres[censoring_idx] <- calc_res_from_prob(prob$data[censoring_idx]) + obj$ires[censoring_idx] <- obj$iwres[censoring_idx] * w_ipred[censoring_idx] + ## if we would calculate the likelihood for the data population parameters given the data, + ## we could also calculate the the equivalents for cwres, wres, and res. However we + ## currently don't have a need to calculate. So setting to NA to avoid wrong interpretation. + obj$cwres[censoring_idx] <- NA_real_ + obj$wres[censoring_idx] <- NA_real_ + obj$res[censoring_idx] <- NA_real_ + } + } + obj$iwres_weighted <- obj$iwres * obj$weights + obj$pred <- pred + obj$ipred <- ipred + obj$prob <- prob + obj$dv <- y + obj$obs_type <- sim_ipred$obs_type + obj +} \ No newline at end of file diff --git a/R/parse_t_obs.R b/R/parse_t_obs.R new file mode 100644 index 0000000..06cfc7f --- /dev/null +++ b/R/parse_t_obs.R @@ -0,0 +1,11 @@ +#' Parse observation times +#' +#' @inheritParams get_map_estimates +#' +parse_t_obs <- function(data) { + t_obs <- data$t + if(any(duplicated(paste(t_obs, data$obs_type, sep = "_")))) { + message("Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument.") + } + t_obs +} diff --git a/R/parse_weight_prior.R b/R/parse_weight_prior.R new file mode 100644 index 0000000..5b6a435 --- /dev/null +++ b/R/parse_weight_prior.R @@ -0,0 +1,20 @@ +#' Parse the weight of the prior, make sure +#' it is on the variance scale. +#' @param weight_prior weight of the prior specified on SD scale +#' @param type type of fit, e.g. `"map"` or `"pls"` +#' +parse_weight_prior <- function( + weight_prior = 1, + type = "map" +) { + if(is.null(weight_prior) || is.na(weight_prior)) { + weight_prior <- 1 + } + weight_prior_var <- weight_prior^2 + if(tolower(type) == "pls") { + ## RK: Empirically determined to be a good penalty. + ## In principle, PLS is just MAP with very flat priors + weight_prior_var <- 0.001 + } + weight_prior_var +} \ No newline at end of file diff --git a/R/parse_weights.R b/R/parse_weights.R new file mode 100644 index 0000000..1c5dbcd --- /dev/null +++ b/R/parse_weights.R @@ -0,0 +1,14 @@ +#' Parse weights argument +#' +#' @inheritParams get_map_estimates +#' +parse_weights <- function(weights, t_obs) { + if(!is.null(weights)) { + if(length(weights) != length(t_obs)) { + stop("Vector of weights of different size than observation vector!") + } + } else { + weights <- rep(1, length(t_obs)) + } + weights +} \ No newline at end of file diff --git a/R/print.map_estimates.R b/R/print.map_estimates.R index bbb8fd8..7995798 100755 --- a/R/print.map_estimates.R +++ b/R/print.map_estimates.R @@ -7,7 +7,7 @@ #' print.map_estimates <- function(x, ...) { eta <- x$fit$coef - om <- sqrt(diag(PKPDsim::triangle_to_full(x$prior$omega))) + om <- sqrt(diag(PKPDsim::triangle_to_full(x$prior$omega$full))) res_table <- data.frame( individual = unlist(x$parameters), population = unlist(x$prior$parameters) diff --git a/man/calc_residuals.Rd b/man/calc_residuals.Rd new file mode 100644 index 0000000..002eb5a --- /dev/null +++ b/man/calc_residuals.Rd @@ -0,0 +1,88 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calc_residuals.R +\name{calc_residuals} +\alias{calc_residuals} +\title{Function to calculate residuals from a fit, +based on a fitted parameter set, and data object.} +\usage{ +calc_residuals( + obj, + data, + model, + parameters_population, + covariates, + regimen, + omega_full, + error, + weights, + transf, + A_init_population, + A_init_individual, + t_init = 0, + iov_bins = NULL, + output_include = c(), + int_step_size = 0.01, + censoring = NULL, + censoring_idx = NULL, + data_before_init = NULL, + ltbs = FALSE, + ... +) +} +\arguments{ +\item{obj}{temporary storage object in `get_map_estimates()`} + +\item{data}{data data.frame with columns `t` and `y` (and possibly evid)} + +\item{model}{model, created using `PKPDsim::new_ode_model()`} + +\item{parameters_population}{parameters for population predictions} + +\item{covariates}{list of covariates, each one created using `PKPDsim::new_coviarate()`} + +\item{regimen}{regimen} + +\item{omega_full}{full omega matrix} + +\item{error}{residual error, specified as list with arguments `add` and/or +`prop` specifying the additive and proportional parts} + +\item{weights}{vector of weights for error. Length of vector should be same +as length of observation vector. If NULL (default), all weights are equal. +Used in both MAP and NP methods. Note that `weights` argument will also +affect residuals (residuals will be scaled too).} + +\item{transf}{transformation function} + +\item{A_init_population}{init vector for model state (population)} + +\item{A_init_individual}{init vector for model state (individual)} + +\item{t_init}{initialization time before first dose, default 0.} + +\item{iov_bins}{bins for inter-occasion variability. Passed unchanged to +PKPDsim.} + +\item{output_include}{passed to PKPDsim::sim_ode(), returns covariates and +parameter values over time in return object. Only invoked if `residuals` +option is `TRUE`.} + +\item{int_step_size}{integrator step size passed to PKPDsim} + +\item{censoring}{label for column specifying censoring. If value in dataset +in this column is < 0 then censoring is assumed 0 then >ULOQ.} + +\item{censoring_idx}{vector with indices for censoring} + +\item{data_before_init}{data.frame with data before initial dose} + +\item{ltbs}{log-transform both sides? (`NULL` by default, meaning that it +will be picked up from the PKPDsim model. Can be overridden with `TRUE`). +Note: `error` should commonly only have additive part.} + +\item{...}{parameters passed on to `sim_ode()` function} +} +\description{ +Function to calculate residuals from a fit, +based on a fitted parameter set, and data object. +} diff --git a/man/check_inputs.Rd b/man/check_inputs.Rd new file mode 100644 index 0000000..5e420b4 --- /dev/null +++ b/man/check_inputs.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/check_inputs.R +\name{check_inputs} +\alias{check_inputs} +\title{Some basic input checking} +\usage{ +check_inputs(model, data, parameters, omega, regimen, censoring, type) +} +\arguments{ +\item{model}{model, created using `PKPDsim::new_ode_model()`} + +\item{data}{data data.frame with columns `t` and `y` (and possibly evid)} + +\item{parameters}{list of parameters} + +\item{omega}{between subject variability, supplied as vector specifiying the +lower triangle of the covariance matrix of random effects} + +\item{regimen}{regimen} + +\item{censoring}{label for column specifying censoring. If value in dataset +in this column is < 0 then censoring is assumed 0 then >ULOQ.} + +\item{type}{estimation type, options are `map`, `ls`, and `np_hybrid`} +} +\description{ +Some basic input checking +} diff --git a/man/get_individual_parameters_from_fit.Rd b/man/get_individual_parameters_from_fit.Rd new file mode 100644 index 0000000..7fd95bf --- /dev/null +++ b/man/get_individual_parameters_from_fit.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_individual_parameters_from_fit.R +\name{get_individual_parameters_from_fit} +\alias{get_individual_parameters_from_fit} +\title{Get the individual parameters from the fitted coefficients (etas)} +\usage{ +get_individual_parameters_from_fit( + fit, + parameters, + nonfixed = NULL, + as_eta = NULL +) +} +\arguments{ +\item{fit}{fit object} + +\item{parameters}{list of population parameters} +} +\description{ +Get the individual parameters from the fitted coefficients (etas) +} diff --git a/man/np_fit_wrapper.Rd b/man/np_fit_wrapper.Rd new file mode 100644 index 0000000..de85296 --- /dev/null +++ b/man/np_fit_wrapper.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/np_fit_wrapper.R +\name{np_fit_wrapper} +\alias{np_fit_wrapper} +\title{Wrapper for non-parametric fit} +\usage{ +np_fit_wrapper( + obj, + model, + regimen, + data, + covariates = NULL, + weights = 1, + np_settings = list() +) +} +\arguments{ +\item{obj}{fit object} + +\item{model}{model, created using `PKPDsim::new_ode_model()`} + +\item{regimen}{regimen} + +\item{data}{data data.frame with columns `t` and `y` (and possibly evid)} + +\item{covariates}{list of covariates, each one created using `PKPDsim::new_coviarate()`} + +\item{weights}{vector of weights for error. Length of vector should be same +as length of observation vector. If NULL (default), all weights are equal. +Used in both MAP and NP methods. Note that `weights` argument will also +affect residuals (residuals will be scaled too).} + +\item{np_settings}{list with settings for non-parametric estimation (if +selected), containing any of the following: `error`, `grid_span`, grid_size`, + `grid_exponential`} +} +\description{ +Wrapper for non-parametric fit +} diff --git a/man/parse_censoring.Rd b/man/parse_censoring.Rd new file mode 100644 index 0000000..660fbc6 --- /dev/null +++ b/man/parse_censoring.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_censoring.R +\name{parse_censoring} +\alias{parse_censoring} +\title{Check if censoring functionality needs to be used, +and create index for censored data} +\usage{ +parse_censoring(censoring, data, verbose = FALSE) +} +\arguments{ +\item{censoring}{label for column specifying censoring. If value in dataset +in this column is < 0 then censoring is assumed 0 then >ULOQ.} + +\item{data}{data data.frame with columns `t` and `y` (and possibly evid)} + +\item{verbose}{show more output} +} +\description{ +Check if censoring functionality needs to be used, +and create index for censored data +} diff --git a/man/parse_error.Rd b/man/parse_error.Rd new file mode 100644 index 0000000..eba79a7 --- /dev/null +++ b/man/parse_error.Rd @@ -0,0 +1,15 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_error.R +\name{parse_error} +\alias{parse_error} +\title{Parse specified error magnitudes} +\usage{ +parse_error(error) +} +\arguments{ +\item{error}{residual error, specified as list with arguments `add` and/or +`prop` specifying the additive and proportional parts} +} +\description{ +Parse specified error magnitudes +} diff --git a/man/parse_fixed.Rd b/man/parse_fixed.Rd new file mode 100644 index 0000000..a24720b --- /dev/null +++ b/man/parse_fixed.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_fixed.R +\name{parse_fixed} +\alias{parse_fixed} +\title{Parse fixed parameters} +\usage{ +parse_fixed(fixed, parameters) +} +\arguments{ +\item{fixed}{fix a specific parameters, supplied as vector of strings} + +\item{parameters}{list of parameters} +} +\description{ +Parse fixed parameters +} diff --git a/man/parse_input_data.Rd b/man/parse_input_data.Rd new file mode 100644 index 0000000..036f80a --- /dev/null +++ b/man/parse_input_data.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_input_data.R +\name{parse_input_data} +\alias{parse_input_data} +\title{Parse input data for fit} +\usage{ +parse_input_data(data, cols = list(x = "t", y = "y"), obs_type_label = NULL) +} +\arguments{ +\item{data}{data data.frame with columns `t` and `y` (and possibly evid)} + +\item{cols}{column names} + +\item{obs_type_label}{column name in `data` referring to observation type. +Can be used for specification of different residual error models for +differing observation types (e.g. venous and capillary or parent and +metabolite), Residual error should then be specified as list of vectors, +e.g. `list(prop = c(0.2, 0.1), add = c(1, 2))`.} +} +\description{ +Parse input data for fit +} diff --git a/man/parse_omega_matrix.Rd b/man/parse_omega_matrix.Rd new file mode 100644 index 0000000..7d1d37a --- /dev/null +++ b/man/parse_omega_matrix.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_omega_matrix.R +\name{parse_omega_matrix} +\alias{parse_omega_matrix} +\title{Parse omega matrix into full matrix, and +perform some checks} +\usage{ +parse_omega_matrix(omega, parameters, fixed) +} +\arguments{ +\item{omega}{between subject variability, supplied as vector specifiying the +lower triangle of the covariance matrix of random effects} + +\item{parameters}{list of parameters} + +\item{fixed}{fix a specific parameters, supplied as vector of strings} +} +\description{ +Parse omega matrix into full matrix, and +perform some checks +} diff --git a/man/parse_residuals_from_predictions.Rd b/man/parse_residuals_from_predictions.Rd new file mode 100644 index 0000000..c683d9c --- /dev/null +++ b/man/parse_residuals_from_predictions.Rd @@ -0,0 +1,51 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_residuals_from_predictions.R +\name{parse_residuals_from_predictions} +\alias{parse_residuals_from_predictions} +\title{Parse simulated data and extract residuals} +\usage{ +parse_residuals_from_predictions( + obj, + sim_ipred, + sim_pred, + data, + omega_full, + transf, + error, + censoring, + censoring_idx, + data_before_init, + weights +) +} +\arguments{ +\item{obj}{temporary storage object in `get_map_estimates()`} + +\item{sim_ipred}{data frame with simulated individual predictions} + +\item{sim_pred}{data frame with simulated population predictions} + +\item{data}{data data.frame with columns `t` and `y` (and possibly evid)} + +\item{omega_full}{full omega matrix} + +\item{transf}{transformation function} + +\item{error}{residual error, specified as list with arguments `add` and/or +`prop` specifying the additive and proportional parts} + +\item{censoring}{label for column specifying censoring. If value in dataset +in this column is < 0 then censoring is assumed 0 then >ULOQ.} + +\item{censoring_idx}{vector with indices for censoring} + +\item{data_before_init}{data.frame with data before initial dose} + +\item{weights}{vector of weights for error. Length of vector should be same +as length of observation vector. If NULL (default), all weights are equal. +Used in both MAP and NP methods. Note that `weights` argument will also +affect residuals (residuals will be scaled too).} +} +\description{ +Parse simulated data and extract residuals +} diff --git a/man/parse_t_obs.Rd b/man/parse_t_obs.Rd new file mode 100644 index 0000000..eff45cb --- /dev/null +++ b/man/parse_t_obs.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_t_obs.R +\name{parse_t_obs} +\alias{parse_t_obs} +\title{Parse observation times} +\usage{ +parse_t_obs(data) +} +\arguments{ +\item{data}{data data.frame with columns `t` and `y` (and possibly evid)} +} +\description{ +Parse observation times +} diff --git a/man/parse_weight_prior.Rd b/man/parse_weight_prior.Rd new file mode 100644 index 0000000..1b72306 --- /dev/null +++ b/man/parse_weight_prior.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_weight_prior.R +\name{parse_weight_prior} +\alias{parse_weight_prior} +\title{Parse the weight of the prior, make sure +it is on the variance scale.} +\usage{ +parse_weight_prior(weight_prior = 1, type = "map") +} +\arguments{ +\item{weight_prior}{weight of the prior specified on SD scale} + +\item{type}{type of fit, e.g. `"map"` or `"pls"`} +} +\description{ +Parse the weight of the prior, make sure +it is on the variance scale. +} diff --git a/man/parse_weights.Rd b/man/parse_weights.Rd new file mode 100644 index 0000000..4450de8 --- /dev/null +++ b/man/parse_weights.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parse_weights.R +\name{parse_weights} +\alias{parse_weights} +\title{Parse weights argument} +\usage{ +parse_weights(weights, t_obs) +} +\arguments{ +\item{weights}{vector of weights for error. Length of vector should be same +as length of observation vector. If NULL (default), all weights are equal. +Used in both MAP and NP methods. Note that `weights` argument will also +affect residuals (residuals will be scaled too).} +} +\description{ +Parse weights argument +} diff --git a/tests/testthat/output/print.map_estimates_1.txt b/tests/testthat/output/print.map_estimates_1.txt index 2d95344..e33c504 100644 --- a/tests/testthat/output/print.map_estimates_1.txt +++ b/tests/testthat/output/print.map_estimates_1.txt @@ -1,9 +1,9 @@ > PKPDmap:::print.map_estimates(tmp1) - individual population eta relative -CL 6.699894 7.67 -0.6711100 --o--|----- -V 87.000059 97.70 -0.3391079 ----o|----- -Q 3.421152 3.00 0.4154129 -----|o---- -V2 50.000000 50.00 NA + individual population eta relative +CL 6.699894 7.67 -0.671110 --o--|----- +V 87.000059 97.70 -3.668013 o-----|----- +Q 3.421152 3.00 1.313651 -----|--o-- +V2 50.000000 50.00 NA dv ipred pred res iwres 1 1119.10 1116.5710 996.0314 123.068635 0.01922780 diff --git a/tests/testthat/test_as_eta.R b/tests/testthat/test-as_eta.R similarity index 100% rename from tests/testthat/test_as_eta.R rename to tests/testthat/test-as_eta.R diff --git a/tests/testthat/test_calc_res_from_prob.R b/tests/testthat/test-calc_res_from_prob.R similarity index 100% rename from tests/testthat/test_calc_res_from_prob.R rename to tests/testthat/test-calc_res_from_prob.R diff --git a/tests/testthat/test-check_inputs.R b/tests/testthat/test-check_inputs.R new file mode 100644 index 0000000..adf5d1d --- /dev/null +++ b/tests/testthat/test-check_inputs.R @@ -0,0 +1,220 @@ +test_that("check_inputs passes with valid MAP inputs", { + # Create mock objects + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + # Should not throw any errors + expect_no_error(check_inputs(model, data, parameters, omega, regimen, NULL, "MAP")) +}) + +test_that("check_inputs passes with valid PLS inputs", { + # Create mock objects + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + # Should not throw any errors + expect_no_error(check_inputs(model, data, parameters, omega, regimen, NULL, "PLS")) +}) + +test_that("check_inputs passes with valid pls inputs (lowercase)", { + # Create mock objects + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + # Should not throw any errors + expect_no_error(check_inputs(model, data, parameters, omega, regimen, NULL, "pls")) +}) + +test_that("check_inputs passes with valid MAP inputs (lowercase)", { + # Create mock objects + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + # Should not throw any errors + expect_no_error(check_inputs(model, data, parameters, omega, regimen, NULL, "map")) +}) + +test_that("check_inputs passes with valid censoring argument", { + # Create mock objects + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + # Should not throw any errors with character censoring + expect_no_error(check_inputs(model, data, parameters, omega, regimen, "cens_col", "MAP")) +}) + +test_that("check_inputs passes with other type values", { + # Create mock objects + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + # Should not throw any errors for non-MAP/PLS types + expect_no_error(check_inputs(model, data, parameters, omega, regimen, NULL, "other_type")) +}) + +test_that("check_inputs fails when model is NULL for MAP type", { + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + expect_error( + check_inputs(NULL, data, parameters, omega, regimen, NULL, "MAP"), + "The 'model', 'data', 'omega', 'regimen', and 'parameters' arguments are required." + ) +}) + +test_that("check_inputs fails when data is NULL for MAP type", { + model <- function() {} + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + expect_error( + check_inputs(model, NULL, parameters, omega, regimen, NULL, "MAP"), + "The 'model', 'data', 'omega', 'regimen', and 'parameters' arguments are required." + ) +}) + +test_that("check_inputs fails when parameters is NULL for MAP type", { + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + expect_error( + check_inputs(model, data, NULL, omega, regimen, NULL, "MAP"), + "The 'model', 'data', 'omega', 'regimen', and 'parameters' arguments are required." + ) +}) + +test_that("check_inputs fails when omega is NULL for MAP type", { + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + regimen <- list(dose = 100, interval = 12) + + expect_error( + check_inputs(model, data, parameters, NULL, regimen, NULL, "MAP"), + "The 'model', 'data', 'omega', 'regimen', and 'parameters' arguments are required." + ) +}) + +test_that("check_inputs fails when regimen is NULL for MAP type", { + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + + expect_error( + check_inputs(model, data, parameters, omega, NULL, NULL, "MAP"), + "The 'model', 'data', 'omega', 'regimen', and 'parameters' arguments are required." + ) +}) + +test_that("check_inputs fails when model is not a function", { + model <- "not_a_function" + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + expect_error( + check_inputs(model, data, parameters, omega, regimen, NULL, "MAP"), + "The 'model' argument requires a function, e.g. a model defined using the new_ode_model\\(\\) function from the PKPDsim package." + ) +}) + +test_that("check_inputs fails when model is a list (not a function)", { + model <- list(not_a_function = TRUE) + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + expect_error( + check_inputs(model, data, parameters, omega, regimen, NULL, "MAP"), + "The 'model' argument requires a function, e.g. a model defined using the new_ode_model\\(\\) function from the PKPDsim package." + ) +}) + +test_that("check_inputs fails when censoring is not NULL and not character", { + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + # Test with numeric censoring + expect_error( + check_inputs(model, data, parameters, omega, regimen, 123, "MAP"), + "Censoring argument requires label specifying column in dataset with censoring info." + ) + + # Test with logical censoring + expect_error( + check_inputs(model, data, parameters, omega, regimen, TRUE, "MAP"), + "Censoring argument requires label specifying column in dataset with censoring info." + ) + + # Test with list censoring + expect_error( + check_inputs(model, data, parameters, omega, regimen, list(cens = "col"), "MAP"), + "Censoring argument requires label specifying column in dataset with censoring info." + ) +}) + +test_that("check_inputs allows NULL censoring", { + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + # Should not throw any errors with NULL censoring + expect_no_error(check_inputs(model, data, parameters, omega, regimen, NULL, "MAP")) +}) + +test_that("check_inputs allows character censoring", { + model <- function() {} + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + # Should not throw any errors with character censoring + expect_no_error(check_inputs(model, data, parameters, omega, regimen, "censoring_column", "MAP")) +}) + +test_that("check_inputs works with different function types", { + # Test with anonymous function + model <- function(x) x + 1 + data <- data.frame(time = 1:3, dv = c(1, 2, 3)) + parameters <- list(CL = 1, V = 10) + omega <- matrix(c(0.1, 0, 0, 0.1), nrow = 2) + regimen <- list(dose = 100, interval = 12) + + expect_no_error(check_inputs(model, data, parameters, omega, regimen, NULL, "MAP")) + + # Test with named function + test_model <- function() {} + expect_no_error(check_inputs(test_model, data, parameters, omega, regimen, NULL, "MAP")) +}) \ No newline at end of file diff --git a/tests/testthat/test_create_iov_object.R b/tests/testthat/test-create_iov_object.R similarity index 100% rename from tests/testthat/test_create_iov_object.R rename to tests/testthat/test-create_iov_object.R diff --git a/tests/testthat/test-get_individual_parameters_from_fit.R b/tests/testthat/test-get_individual_parameters_from_fit.R new file mode 100644 index 0000000..c9e8a1c --- /dev/null +++ b/tests/testthat/test-get_individual_parameters_from_fit.R @@ -0,0 +1,98 @@ +test_that("get_individual_parameters_from_fit handles basic case", { + # Create test data + fit <- list(coef = c(0.1, 0.2)) # Two coefficients + parameters <- list(CL = 10, V = 20) + nonfixed <- c("CL", "V") + as_eta <- NULL + + result <- get_individual_parameters_from_fit(fit, parameters, nonfixed, as_eta) + + # Check that parameters are correctly transformed + expect_equal(result$CL, 10 * exp(0.1)) + expect_equal(result$V, 20 * exp(0.2)) +}) + +test_that("get_individual_parameters_from_fit handles as_eta parameters", { + # Create test data with some parameters as direct etas + fit <- list(coef = c(0.1, 0.2)) + parameters <- list(CL = 10, V = 20) + nonfixed <- c("CL", "V") + as_eta <- c("CL") # CL is direct eta + + result <- get_individual_parameters_from_fit(fit, parameters, nonfixed, as_eta) + + # Check that CL is direct eta and V is transformed + expect_equal(result$CL, 0.1) + expect_equal(result$V, 20 * exp(0.2)) +}) + +test_that("get_individual_parameters_from_fit handles all as_eta parameters", { + # Create test data with all parameters as direct etas + fit <- list(coef = c(0.1, 0.2)) + parameters <- list(CL = 10, V = 20) + nonfixed <- c("CL", "V") + as_eta <- c("CL", "V") + + result <- get_individual_parameters_from_fit(fit, parameters, nonfixed, as_eta) + + # Check that all parameters are direct etas + expect_equal(result$CL, 0.1) + expect_equal(result$V, 0.2) +}) + +test_that("get_individual_parameters_from_fit handles no nonfixed parameters", { + # Create test data with no nonfixed parameters + fit <- list(coef = numeric(0)) + parameters <- list(CL = 10, V = 20) + nonfixed <- character(0) + as_eta <- NULL + + result <- get_individual_parameters_from_fit(fit, parameters, nonfixed, as_eta) + + # Check that parameters remain unchanged + expect_equal(result$CL, 10) + expect_equal(result$V, 20) +}) + +test_that("get_individual_parameters_from_fit handles negative coefficients", { + # Create test data with negative coefficients + fit <- list(coef = c(-0.1, -0.2)) + parameters <- list(CL = 10, V = 20) + nonfixed <- c("CL", "V") + as_eta <- NULL + + result <- get_individual_parameters_from_fit(fit, parameters, nonfixed, as_eta) + + # Check that parameters are correctly transformed with negative coefficients + expect_equal(result$CL, 10 * exp(-0.1)) + expect_equal(result$V, 20 * exp(-0.2)) +}) + +test_that("get_individual_parameters_from_fit handles zero coefficients", { + # Create test data with zero coefficients + fit <- list(coef = c(0, 0)) + parameters <- list(CL = 10, V = 20) + nonfixed <- c("CL", "V") + as_eta <- NULL + + result <- get_individual_parameters_from_fit(fit, parameters, nonfixed, as_eta) + + # Check that parameters remain unchanged with zero coefficients + expect_equal(result$CL, 10) + expect_equal(result$V, 20) +}) + +test_that("get_individual_parameters_from_fit handles mixed parameter types", { + # Create test data with mixed parameter types (some as_eta, some transformed) + fit <- list(coef = c(0.1, 0.2, 0.3)) + parameters <- list(CL = 10, V = 20, KA = 30) + nonfixed <- c("CL", "V", "KA") + as_eta <- c("V") # Only V is direct eta + + result <- get_individual_parameters_from_fit(fit, parameters, nonfixed, as_eta) + + # Check that parameters are correctly handled + expect_equal(result$CL, 10 * exp(0.1)) + expect_equal(result$V, 0.2) + expect_equal(result$KA, 30 * exp(0.3)) +}) \ No newline at end of file diff --git a/tests/testthat/test_get_mahalanobis.R b/tests/testthat/test-get_mahalanobis.R similarity index 100% rename from tests/testthat/test_get_mahalanobis.R rename to tests/testthat/test-get_mahalanobis.R diff --git a/tests/testthat/test_get_map_estimates.R b/tests/testthat/test-get_map_estimates.R similarity index 100% rename from tests/testthat/test_get_map_estimates.R rename to tests/testthat/test-get_map_estimates.R diff --git a/tests/testthat/test_get_map_estimates.ss.R b/tests/testthat/test-get_map_estimates.ss.R similarity index 100% rename from tests/testthat/test_get_map_estimates.ss.R rename to tests/testthat/test-get_map_estimates.ss.R diff --git a/tests/testthat/test_get_varcov_matrix.R b/tests/testthat/test-get_varcov_matrix.R similarity index 100% rename from tests/testthat/test_get_varcov_matrix.R rename to tests/testthat/test-get_varcov_matrix.R diff --git a/tests/testthat/test-parse_censoring.R b/tests/testthat/test-parse_censoring.R new file mode 100644 index 0000000..48b4b41 --- /dev/null +++ b/tests/testthat/test-parse_censoring.R @@ -0,0 +1,36 @@ +test_that("parse_censoring handles NULL censoring", { + data <- data.frame(dv = c(1, 2, 3)) + result <- parse_censoring(NULL, data) + expect_null(result) +}) + +test_that("parse_censoring handles no censored values", { + data <- data.frame(dv = c(1, 2, 3), cens = c(0, 0, 0)) + result <- parse_censoring("cens", data) + expect_null(result) + + # Test with verbose messaging + expect_message( + parse_censoring("cens", data, verbose = TRUE), + "Warning: censoring specified, but no censored values in data." + ) +}) + +test_that("parse_censoring handles censored values", { + data <- data.frame(dv = c(1, 2, 3), cens = c(0, 1, 0)) + result <- parse_censoring("cens", data) + expect_equal(result, c(FALSE, TRUE, FALSE)) + + # Test with verbose messaging + expect_message( + parse_censoring("cens", data, verbose = TRUE), + "One or more values in data are censored, including censoring in likelihood." + ) +}) + +test_that("parse_censoring is case insensitive", { + # Test with lowercase column name + data <- data.frame(dv = c(1, 2, 3), cens = c(0, 1, 0)) + result <- parse_censoring("CENS", data) + expect_equal(result, c(FALSE, TRUE, FALSE)) +}) diff --git a/tests/testthat/test-parse_error.R b/tests/testthat/test-parse_error.R new file mode 100644 index 0000000..fa74ac3 --- /dev/null +++ b/tests/testthat/test-parse_error.R @@ -0,0 +1,53 @@ +test_that("parse_error handles NULL input", { + expect_error( + parse_error(NULL), + "No residual error specified, cannot perform fit." + ) +}) + +test_that("parse_error handles missing error components", { + # Test with only prop specified + error <- list(prop = 0.1) + result <- parse_error(error) + expect_equal(result$prop, 0.1) + expect_equal(result$add, 0) + expect_equal(result$exp, 0) + + # Test with only add specified + error <- list(add = 0.2) + result <- parse_error(error) + expect_equal(result$prop, 0) + expect_equal(result$add, 0.2) + expect_equal(result$exp, 0) + + # Test with only exp specified + error <- list(exp = 0.3) + result <- parse_error(error) + expect_equal(result$prop, 0) + expect_equal(result$add, 0) + expect_equal(result$exp, 0.3) +}) + +test_that("parse_error handles all zeros", { + error <- list(prop = 0, add = 0, exp = 0) + expect_error( + parse_error(error), + "No residual error model specified, or residual error is 0." + ) +}) + +test_that("parse_error handles valid error specifications", { + # Test with all components specified + error <- list(prop = 0.1, add = 0.2, exp = 0.3) + result <- parse_error(error) + expect_equal(result$prop, 0.1) + expect_equal(result$add, 0.2) + expect_equal(result$exp, 0.3) + + # Test with mixed zero and non-zero components + error <- list(prop = 0.1, add = 0, exp = 0.3) + result <- parse_error(error) + expect_equal(result$prop, 0.1) + expect_equal(result$add, 0) + expect_equal(result$exp, 0.3) +}) \ No newline at end of file diff --git a/tests/testthat/test-parse_fixed.R b/tests/testthat/test-parse_fixed.R new file mode 100644 index 0000000..9804732 --- /dev/null +++ b/tests/testthat/test-parse_fixed.R @@ -0,0 +1,147 @@ +test_that("parse_fixed returns correct fixed parameters when all exist", { + parameters <- list(CL = 1, V = 10, KA = 2) + fixed <- c("CL", "V") + + result <- parse_fixed(fixed, parameters) + expect_equal(result, c("CL", "V")) +}) + +test_that("parse_fixed returns correct fixed parameters when some exist", { + parameters <- list(CL = 1, V = 10, KA = 2) + fixed <- c("CL", "V", "NONEXISTENT") + + # Should warn about missing parameters + expect_warning( + result <- parse_fixed(fixed, parameters), + "Warning: not all fixed parameters were found in parameter set!" + ) + + # Should return only existing parameters + expect_equal(result, c("CL", "V")) +}) + +test_that("parse_fixed returns NULL when no fixed parameters exist", { + parameters <- list(CL = 1, V = 10, KA = 2) + fixed <- c("NONEXISTENT1", "NONEXISTENT2") + + # Should warn about missing parameters + expect_warning( + result <- parse_fixed(fixed, parameters), + "Warning: not all fixed parameters were found in parameter set!" + ) + + # Should return NULL when no parameters match + expect_null(result) +}) + +test_that("parse_fixed returns NULL when fixed is empty", { + parameters <- list(CL = 1, V = 10, KA = 2) + fixed <- character(0) + + result <- parse_fixed(fixed, parameters) + expect_null(result) +}) + +test_that("parse_fixed handles case sensitivity correctly", { + parameters <- list(CL = 1, V = 10, KA = 2) + fixed <- c("cl", "V") # Note: "cl" is lowercase + + # Should warn about missing parameters (case sensitive) + expect_warning( + result <- parse_fixed(fixed, parameters), + "Warning: not all fixed parameters were found in parameter set!" + ) + + # Should return only exact matches + expect_equal(result, "V") +}) + +test_that("parse_fixed handles single parameter correctly", { + parameters <- list(CL = 1, V = 10, KA = 2) + fixed <- "CL" + + result <- parse_fixed(fixed, parameters) + expect_equal(result, "CL") +}) + +test_that("parse_fixed handles single non-existent parameter", { + parameters <- list(CL = 1, V = 10, KA = 2) + fixed <- "NONEXISTENT" + + # Should warn about missing parameter + expect_warning( + result <- parse_fixed(fixed, parameters), + "Warning: not all fixed parameters were found in parameter set!" + ) + + # Should return NULL + expect_null(result) +}) + +test_that("parse_fixed handles parameters with special characters", { + parameters <- list("CL_1" = 1, "V_central" = 10, "KA_abs" = 2) + fixed <- c("CL_1", "V_central") + + result <- parse_fixed(fixed, parameters) + expect_equal(result, c("CL_1", "V_central")) +}) + +test_that("parse_fixed handles mixed existing and non-existing parameters", { + parameters <- list(CL = 1, V = 10, KA = 2) + fixed <- c("CL", "NONEXISTENT", "V") + + # Should warn about missing parameters + expect_warning( + result <- parse_fixed(fixed, parameters), + "Warning: not all fixed parameters were found in parameter set!" + ) + + # Should return existing parameters in order they appear in fixed + expect_equal(result, c("CL", "V")) +}) + +test_that("parse_fixed handles duplicate parameters in fixed", { + parameters <- list(CL = 1, V = 10, KA = 2) + fixed <- c("CL", "V", "CL") # CL appears twice + + result <- parse_fixed(fixed, parameters) + expect_equal(result, c("CL", "V")) +}) + +test_that("parse_fixed handles empty parameters list", { + parameters <- list() + fixed <- c("CL", "V") + + # Should warn about missing parameters + expect_warning( + result <- parse_fixed(fixed, parameters), + "Warning: not all fixed parameters were found in parameter set!" + ) + + # Should return NULL + expect_null(result) +}) + +test_that("parse_fixed handles NULL fixed parameter", { + parameters <- list(CL = 1, V = 10, KA = 2) + + # NULL should be treated as empty character vector + expect_null(parse_fixed(NULL, parameters)) +}) + +test_that("parse_fixed handles numeric parameters", { + parameters <- list(CL = 1, V = 10, KA = 2) + fixed <- c("CL", "V") + + result <- parse_fixed(fixed, parameters) + expect_equal(result, c("CL", "V")) +}) + +test_that("parse_fixed handles character parameters", { + parameters <- list(CL = "value1", V = "value2", KA = "value3") + fixed <- c("CL", "V") + + result <- parse_fixed(fixed, parameters) + expect_equal(result, c("CL", "V")) +}) + diff --git a/tests/testthat/test-parse_input_data.R b/tests/testthat/test-parse_input_data.R new file mode 100644 index 0000000..2457a2f --- /dev/null +++ b/tests/testthat/test-parse_input_data.R @@ -0,0 +1,135 @@ +test_that("parse_input_data handles regular data frame input", { + # Create test data + data <- data.frame( + T = c(2, 1, 3), + OBS_TYPE = c(2, 1, 1), + Y = c(10, 20, 30) + ) + + result <- parse_input_data(data) + + # Check column names are lowercase + expect_equal(names(result), c("t", "obs_type", "y")) + + # Check sorting + expect_equal(result$t, c(1, 2, 3)) + expect_equal(result$obs_type, c(1, 1, 1)) +}) + +test_that("parse_input_data handles PKPDsim object", { + # Create mock PKPDsim object + data <- structure( + data.frame( + t = c(2, 1, 3), + comp = c("obs", "dose", "obs"), + dv = c(10, 20, 30) + ), + class = c("data.frame", "PKPDsim") + ) + + result <- parse_input_data(data, cols=list(x = "t", y = "dv")) + + # Check that only obs rows are kept + expect_equal(nrow(result), 2) + expect_true(all(result$comp == "obs")) + + # Check that evid column is added and set to 0 + expect_true("evid" %in% names(result)) + expect_true(all(result$evid == 0)) +}) + +test_that("parse_input_data handles EVID filtering", { + # Create test data with EVID column + data <- data.frame( + t = c(2, 1, 3), + obs_type = c(2, 1, 1), + dv = c(10, 20, 30), + EVID = c(0, 1, 0) + ) + + result <- parse_input_data(data, cols=list(x = "t", y = "dv")) + + # Check that only EVID=0 rows are kept + expect_equal(nrow(result), 2) + expect_true(all(result$evid == 0)) +}) + +test_that("parse_input_data handles obs_type_label argument", { + # Test with NULL obs_type_label (default behavior) + data <- data.frame( + t = c(2, 1, 3), + dv = c(10, 20, 30) + ) + result <- parse_input_data(data, cols=list(x = "t", y = "dv")) + expect_equal(result$obs_type, c(1, 1, 1)) + + # Test with custom obs_type_label + data <- data.frame( + t = c(2, 1, 3), + dv = c(10, 20, 30), + measurement_type = c(2, 1, 3) + ) + result <- parse_input_data( + data, + cols=list(x = "t", y = "dv"), + obs_type_label = "measurement_type" + ) + expect_equal(result$obs_type, c(1, 2, 3)) +}) + +test_that("parse_input_data handles PKPDsim object with obs_type_label", { + # Create mock PKPDsim object with obs_type_label + data <- structure( + data.frame( + t = c(2, 1, 3), + comp = c("obs", 1, "obs"), + dv = c(10, 20, 30), + measurement_type = c(2, 1, 3) + ), + class = c("data.frame", "PKPDsim") + ) + + result <- parse_input_data( + data, + cols=list(x = "t", y = "dv"), + obs_type_label = "measurement_type" + ) + + # Check that only obs rows are kept and obs_type is set correctly + expect_equal(nrow(result), 2) + expect_true(all(result$comp == "obs")) + expect_equal(result$obs_type, c(2, 3)) +}) + +test_that("parse_input_data handles multiple EVID values", { + # Create test data with multiple EVID values + data <- data.frame( + t = c(2, 1, 3, 4), + obs_type = c(2, 1, 1, 1), + dv = c(10, 20, 30, 40), + EVID = c(0, 1, 2, 0) + ) + + result <- parse_input_data(data, cols=list(x = "t", y = "dv")) + + # Check that only EVID=0 rows are kept + expect_equal(nrow(result), 2) + expect_true(all(result$evid == 0)) + expect_equal(result$t, c(2, 4)) +}) + +test_that("parse_input_data preserves data values after transformations", { + # Create test data + data <- data.frame( + t = c(2, 1, 3), + obs_type = c(2, 1, 1), + DV = c(10, 20, 30), + EVID = c(0, 1, 0) + ) + + result <- parse_input_data(data, cols=list(x = "t", y = "dv")) + + # Check that data values are preserved after transformations + expect_equal(result$dv, c(10, 30)) + expect_equal(result$t, c(2, 3)) +}) diff --git a/tests/testthat/test-parse_omega_matrix.R b/tests/testthat/test-parse_omega_matrix.R new file mode 100644 index 0000000..51f761f --- /dev/null +++ b/tests/testthat/test-parse_omega_matrix.R @@ -0,0 +1,91 @@ +test_that("parse_omega_matrix handles matrix input", { + # Create test matrix + omega <- matrix(c(0.1, 0.2, 0.2, 0.3), nrow = 2, ncol = 2) + parameters <- c("CL", "V") + fixed <- character(0) + + result <- parse_omega_matrix( + omega, + parameters, + fixed + ) + + # Check that matrix is returned unchanged + expect_equal(result$full, omega) + expect_equal(dim(result$full), c(2, 2)) +}) + +test_that("parse_omega_matrix handles triangle matrix input", { + # Create test triangle matrix (lower triangle) + omega <- c(0.1, 0.2, 0.3) # represents matrix c(0.1, 0.2, 0.2, 0.3) + parameters <- c("CL", "V") + fixed <- character(0) + + result <- parse_omega_matrix(omega, parameters, fixed) + + # Check that triangle matrix is converted to full matrix + expect_equal(dim(result$full), c(2, 2)) + expect_equal(result$full[1,1], 0.1) + expect_equal(result$full[2,1], 0.2) + expect_equal(result$full[1,2], 0.2) + expect_equal(result$full[2,2], 0.3) +}) + +test_that("parse_omega_matrix validates matrix size correctly", { + # Test with matching sizes + omega <- c(0.1, 0.2, 0.3) # 2x2 matrix + parameters <- c("CL", "V") + fixed <- character(0) + + result <- parse_omega_matrix(omega, parameters, fixed) + expect_equal(dim(result$full), c(2, 2)) + + # Test with fixed parameters + omega <- c(0.1, 0.2, 0.3) # 2x2 matrix + parameters <- c("CL", "V", "KA") + fixed <- "KA" + + result <- parse_omega_matrix(omega, parameters, fixed) + expect_equal(dim(result$full), c(2, 2)) +}) + +test_that("parse_omega_matrix handles error cases", { + # Test with too small omega matrix + omega <- c(0.1) # 1x1 matrix + parameters <- c("CL", "V") + fixed <- character(0) + + expect_error( + parse_omega_matrix(omega, parameters, fixed), + "Provided omega matrix is smaller than expected based on the number of model parameters" + ) + + # Test with too small omega matrix and fixed parameters + omega <- c(0.1) # 1x1 matrix + parameters <- c("CL", "V", "KA") + fixed <- "KA" + + expect_error( + parse_omega_matrix(omega, parameters, fixed), + "Provided omega matrix is smaller than expected based on the number of model parameters" + ) +}) + +test_that("parse_omega_matrix handles edge cases", { + # Test with single parameter + omega <- c(0.1) # 1x1 matrix + parameters <- c("CL") + fixed <- character(0) + + result <- parse_omega_matrix(omega, parameters, fixed) + expect_equal(dim(result$full), c(1, 1)) + expect_equal(result$full[1,1], 0.1) + + # Test with all parameters fixed + omega <- c(0.1) # 1x1 matrix + parameters <- c("CL", "V") + fixed <- c("CL", "V") + + result <- parse_omega_matrix(omega, parameters, fixed) + expect_equal(dim(result$full), c(1, 1)) +}) diff --git a/tests/testthat/test-parse_t_obs.R b/tests/testthat/test-parse_t_obs.R new file mode 100644 index 0000000..c566a3b --- /dev/null +++ b/tests/testthat/test-parse_t_obs.R @@ -0,0 +1,165 @@ +test_that("parse_t_obs returns observation times correctly", { + data <- data.frame(t = c(0, 1, 2, 3), obs_type = c("PK", "PK", "PD", "PD")) + + result <- parse_t_obs(data) + expect_equal(result, c(0, 1, 2, 3)) +}) + +test_that("parse_t_obs handles single observation", { + data <- data.frame(t = 0, obs_type = "PK") + + result <- parse_t_obs(data) + expect_equal(result, 0) +}) + +test_that("parse_t_obs handles empty data frame", { + data <- data.frame(t = numeric(0), obs_type = character(0)) + + result <- parse_t_obs(data) + expect_equal(result, numeric(0)) +}) + +test_that("parse_t_obs handles duplicate times with different obs_type", { + data <- data.frame(t = c(0, 1, 1, 2), obs_type = c("PK", "PK", "PD", "PD")) + + # Should not message about duplicates since obs_type differs + expect_silent(result <- parse_t_obs(data)) + expect_equal(result, c(0, 1, 1, 2)) +}) + +test_that("parse_t_obs detects and messages about duplicates with same obs_type", { + data <- data.frame(t = c(0, 1, 1, 2), obs_type = c("PK", "PK", "PK", "PD")) + + # Should message about duplicates + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(0, 1, 1, 2)) +}) + +test_that("parse_t_obs detects duplicates with same time and obs_type", { + data <- data.frame(t = c(0, 1, 1, 1, 2), obs_type = c("PK", "PK", "PK", "PD", "PD")) + + # Should message about duplicates + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(0, 1, 1, 1, 2)) +}) + +test_that("parse_t_obs handles multiple duplicate groups", { + data <- data.frame( + t = c(0, 1, 1, 2, 2, 2, 3), + obs_type = c("PK", "PK", "PK", "PD", "PD", "PK", "PK") + ) + + # Should message about duplicates + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(0, 1, 1, 2, 2, 2, 3)) +}) + +test_that("parse_t_obs handles all identical times and obs_types", { + data <- data.frame(t = c(1, 1, 1), obs_type = c("PK", "PK", "PK")) + + # Should message about duplicates + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(1, 1, 1)) +}) + +test_that("parse_t_obs handles numeric obs_type", { + data <- data.frame(t = c(0, 1, 1, 2), obs_type = c(1, 1, 1, 2)) + + # Should message about duplicates + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(0, 1, 1, 2)) +}) + +test_that("parse_t_obs handles character obs_type with spaces", { + data <- data.frame(t = c(0, 1, 1, 2), obs_type = c("PK", "PK", "PK", "PD")) + + # Should message about duplicates + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(0, 1, 1, 2)) +}) + +test_that("parse_t_obs handles NA values in obs_type", { + data <- data.frame(t = c(0, 1, 1, 2), obs_type = c("PK", NA, NA, "PD")) + + # Should message about duplicates (NA values are treated as identical) + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(0, 1, 1, 2)) +}) + +test_that("parse_t_obs handles NA values in time", { + data <- data.frame(t = c(0, NA, NA, 2), obs_type = c("PK", "PK", "PK", "PD")) + + # Should message about duplicates (NA values are treated as identical) + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(0, NA, NA, 2)) +}) + +test_that("parse_t_obs handles mixed data types in obs_type", { + data <- data.frame(t = c(0, 1, 1, 2), obs_type = c("PK", 1, 1, "PD")) + + # Should message about duplicates + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(0, 1, 1, 2)) +}) + +test_that("parse_t_obs handles decimal times", { + data <- data.frame(t = c(0.5, 1.0, 1.0, 2.5), obs_type = c("PK", "PK", "PK", "PD")) + + # Should message about duplicates + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(0.5, 1.0, 1.0, 2.5)) +}) + +test_that("parse_t_obs handles negative times", { + data <- data.frame(t = c(-1, 0, 0, 1), obs_type = c("PK", "PK", "PK", "PD")) + + # Should message about duplicates + expect_message( + result <- parse_t_obs(data), + "Duplicate times were detected in data. Estimation will proceed but please check that data is correct. For putting more weight on certain measurements, please use the `weights` argument." + ) + expect_equal(result, c(-1, 0, 0, 1)) +}) + +test_that("parse_t_obs handles large datasets without duplicates", { + # Create a large dataset with no duplicates + n <- 1000 + data <- data.frame( + t = 1:n, + obs_type = rep(c("PK", "PD"), length.out = n) + ) + + # Should not message about duplicates + expect_silent(result <- parse_t_obs(data)) + expect_equal(result, 1:n) +}) diff --git a/tests/testthat/test-parse_weight_prior.R b/tests/testthat/test-parse_weight_prior.R new file mode 100644 index 0000000..6c4f7e9 --- /dev/null +++ b/tests/testthat/test-parse_weight_prior.R @@ -0,0 +1,69 @@ +test_that("parse_weight_prior handles default values", { + # Test default weight_prior + result <- parse_weight_prior() + expect_equal(result, 1) # 1^2 = 1 + + # Test default weight_prior with explicit type + result <- parse_weight_prior(type = "map") + expect_equal(result, 1) +}) + +test_that("parse_weight_prior handles NULL and NA inputs", { + # Test NULL weight_prior + result <- parse_weight_prior(NULL) + expect_equal(result, 1) + + # Test NA weight_prior + result <- parse_weight_prior(NA) + expect_equal(result, 1) + + # Test NULL weight_prior with explicit type + result <- parse_weight_prior(NULL, type = "map") + expect_equal(result, 1) +}) + +test_that("parse_weight_prior handles different weight values", { + # Test with weight = 2 + result <- parse_weight_prior(2) + expect_equal(result, 4) # 2^2 = 4 + + # Test with weight = 0.5 + result <- parse_weight_prior(0.5) + expect_equal(result, 0.25) # 0.5^2 = 0.25 + + # Test with weight = 0 + result <- parse_weight_prior(0) + expect_equal(result, 0) # 0^2 = 0 +}) + +test_that("parse_weight_prior handles PLS type", { + # Test with lowercase "pls" + result <- parse_weight_prior(type = "pls") + expect_equal(result, 0.001) + + # Test with uppercase "PLS" + result <- parse_weight_prior(type = "PLS") + expect_equal(result, 0.001) + + # Test with mixed case "Pls" + result <- parse_weight_prior(type = "Pls") + expect_equal(result, 0.001) + + # Test with weight_prior and PLS type + result <- parse_weight_prior(2, type = "pls") + expect_equal(result, 0.001) # Should ignore weight_prior for PLS +}) + +test_that("parse_weight_prior handles MAP type with different cases", { + # Test with lowercase "map" + result <- parse_weight_prior(2, type = "map") + expect_equal(result, 4) + + # Test with uppercase "MAP" + result <- parse_weight_prior(2, type = "MAP") + expect_equal(result, 4) + + # Test with mixed case "Map" + result <- parse_weight_prior(2, type = "Map") + expect_equal(result, 4) +}) \ No newline at end of file diff --git a/tests/testthat/test-parse_weights.R b/tests/testthat/test-parse_weights.R new file mode 100644 index 0000000..90b80f6 --- /dev/null +++ b/tests/testthat/test-parse_weights.R @@ -0,0 +1,204 @@ +test_that("parse_weights returns default weights when weights is NULL", { + t_obs <- c(0, 1, 2, 3) + + result <- parse_weights(NULL, t_obs) + expect_equal(result, c(1, 1, 1, 1)) +}) + +test_that("parse_weights returns default weights for single observation", { + t_obs <- 0 + + result <- parse_weights(NULL, t_obs) + expect_equal(result, 1) +}) + +test_that("parse_weights returns default weights for empty observation vector", { + t_obs <- numeric(0) + + result <- parse_weights(NULL, t_obs) + expect_equal(result, numeric(0)) +}) + +test_that("parse_weights returns provided weights when valid", { + t_obs <- c(0, 1, 2, 3) + weights <- c(1.5, 2.0, 0.5, 1.0) + + result <- parse_weights(weights, t_obs) + expect_equal(result, weights) +}) + +test_that("parse_weights returns provided weights for single observation", { + t_obs <- 0 + weights <- 2.5 + + result <- parse_weights(weights, t_obs) + expect_equal(result, 2.5) +}) + +test_that("parse_weights returns provided weights for empty observation vector", { + t_obs <- numeric(0) + weights <- numeric(0) + + result <- parse_weights(weights, t_obs) + expect_equal(result, numeric(0)) +}) + +test_that("parse_weights handles integer weights", { + t_obs <- c(0, 1, 2) + weights <- c(1L, 2L, 3L) + + result <- parse_weights(weights, t_obs) + expect_equal(result, c(1, 2, 3)) +}) + +test_that("parse_weights handles decimal weights", { + t_obs <- c(0, 1, 2, 3) + weights <- c(0.1, 0.5, 1.5, 2.0) + + result <- parse_weights(weights, t_obs) + expect_equal(result, weights) +}) + +test_that("parse_weights handles zero weights", { + t_obs <- c(0, 1, 2) + weights <- c(0, 1, 0) + + result <- parse_weights(weights, t_obs) + expect_equal(result, c(0, 1, 0)) +}) + +test_that("parse_weights handles negative weights", { + t_obs <- c(0, 1, 2) + weights <- c(-1, 0, 1) + + result <- parse_weights(weights, t_obs) + expect_equal(result, c(-1, 0, 1)) +}) + +test_that("parse_weights handles large weights", { + t_obs <- c(0, 1, 2) + weights <- c(1000, 2000, 3000) + + result <- parse_weights(weights, t_obs) + expect_equal(result, c(1000, 2000, 3000)) +}) + +test_that("parse_weights handles very small weights", { + t_obs <- c(0, 1, 2) + weights <- c(1e-10, 1e-5, 1e-3) + + result <- parse_weights(weights, t_obs) + expect_equal(result, c(1e-10, 1e-5, 1e-3)) +}) + +test_that("parse_weights handles all equal weights", { + t_obs <- c(0, 1, 2, 3, 4) + weights <- rep(2.5, 5) + + result <- parse_weights(weights, t_obs) + expect_equal(result, rep(2.5, 5)) +}) + +test_that("parse_weights handles mixed positive and negative weights", { + t_obs <- c(0, 1, 2, 3) + weights <- c(1, -1, 0, 0.5) + + result <- parse_weights(weights, t_obs) + expect_equal(result, c(1, -1, 0, 0.5)) +}) + +test_that("parse_weights handles large observation vectors", { + n <- 1000 + t_obs <- 1:n + weights <- rep(1.5, n) + + result <- parse_weights(weights, t_obs) + expect_equal(result, rep(1.5, n)) +}) + +test_that("parse_weights handles large observation vectors with NULL weights", { + n <- 1000 + t_obs <- 1:n + + result <- parse_weights(NULL, t_obs) + expect_equal(result, rep(1, n)) +}) + +test_that("parse_weights fails when weights vector is shorter than t_obs", { + t_obs <- c(0, 1, 2, 3) + weights <- c(1, 2) # Only 2 weights for 4 observations + + expect_error( + parse_weights(weights, t_obs), + "Vector of weights of different size than observation vector!" + ) +}) + +test_that("parse_weights fails when weights vector is longer than t_obs", { + t_obs <- c(0, 1) + weights <- c(1, 2, 3, 4) # 4 weights for 2 observations + + expect_error( + parse_weights(weights, t_obs), + "Vector of weights of different size than observation vector!" + ) +}) + +test_that("parse_weights fails when weights is empty but t_obs is not", { + t_obs <- c(0, 1, 2) + weights <- numeric(0) + + expect_error( + parse_weights(weights, t_obs), + "Vector of weights of different size than observation vector!" + ) +}) + +test_that("parse_weights fails when weights is not empty but t_obs is", { + t_obs <- numeric(0) + weights <- c(1, 2, 3) + + expect_error( + parse_weights(weights, t_obs), + "Vector of weights of different size than observation vector!" + ) +}) + +test_that("parse_weights handles NA values in weights", { + t_obs <- c(0, 1, 2) + weights <- c(1, NA, 3) + + result <- parse_weights(weights, t_obs) + expect_equal(result, c(1, NA, 3)) +}) + +test_that("parse_weights handles Inf values in weights", { + t_obs <- c(0, 1, 2) + weights <- c(1, Inf, 3) + + result <- parse_weights(weights, t_obs) + expect_equal(result, c(1, Inf, 3)) +}) + +test_that("parse_weights handles NaN values in weights", { + t_obs <- c(0, 1, 2) + weights <- c(1, NaN, 3) + + result <- parse_weights(weights, t_obs) + expect_equal(result, c(1, NaN, 3)) +}) + +test_that("parse_weights preserves weights data type", { + t_obs <- c(0, 1, 2) + + # Test with numeric weights + weights_numeric <- c(1.0, 2.0, 3.0) + result_numeric <- parse_weights(weights_numeric, t_obs) + expect_equal(result_numeric, weights_numeric) + + # Test with integer weights + weights_integer <- c(1L, 2L, 3L) + result_integer <- parse_weights(weights_integer, t_obs) + expect_equal(result_integer, c(1, 2, 3)) # Note: integers become numeric +}) + diff --git a/tests/testthat/test_print.PKPDmap.R b/tests/testthat/test-print.PKPDmap.R similarity index 100% rename from tests/testthat/test_print.PKPDmap.R rename to tests/testthat/test-print.PKPDmap.R