Example: Implementing a custom cost function
Overview
This example demonstrates how to implement a custom penalty-based cost function for mlxModelFinder, inspired by the approach described in the following paper:
Richardson, S., Irurzun Arana, I., Nowojewski, A. et al. A machine learning approach to population pharmacokinetic modelling automation. Commun Med 5, 327 (2025). https://doi.org/10.1038/s43856-025-01054-8
This cost function combines the model’s statistical fit (AIC) with penalties on various quality metrics such as relative standard errors (RSEs), omega values, shrinkage, and the condition number. It allows model selection to account for both goodness-of-fit and model stability.
Penalty Structure
Each model metric is evaluated against a set of thresholds, with penalties increasing as the metric deviates from acceptable ranges.
The function applies piecewise linear penalties defined as follows:
Parameter | Estimated value | Threshold | Penalty |
|---|---|---|---|
Fixed effects | Relative error | 10% | 0 |
30% | 20 | ||
75% | 100 | ||
1000% | 100000 | ||
Random effects | Relative error | 20% | 0 |
40% | 30 | ||
65% | 100 | ||
1000% | 100000 | ||
Random effects | Value (CV%) | 0 (0 CV%) | 1000000 |
0.001 (3 CV%) | 1000 | ||
0.01 (10 CV%) | 0 | ||
1 (131 CV%) | 0 | ||
1.65 (205 CV%) | 100 | ||
2.3 (300 CV%) | 500 | ||
3 (437 CV%) | 10000 | ||
3.5 (567 CV%) | 1000000 | ||
Random effects | Shrinkage | 40% | 0 |
100% | 200 | ||
Condition number | Condition number | 500 | 0 |
1000 | 100 | ||
10000 | 1000 | ||
1000000 | 10000000 |
Implementation in R
Below is a sample implementation of the penalized cost function.
It reads Monolix output files, extracts relevant metrics, computes penalties via linear interpolation, and returns the penalized AIC value.
calculate_penalized_aic <- function() {
# --- 1. Setup and Helper Functions ---
# Helper for piecewise-linear interpolation
# It calculates a y-value for a given x based on a set of (xp, yp) points.
piecewise_linear_interp <- function(x, xp, yp) {
stopifnot(length(xp) == length(yp), length(xp) >= 2)
if (is.na(x)) return(0)
if (x <= xp[1]) return(yp[1])
if (x >= tail(xp, 1)) return(tail(yp, 1))
i <- findInterval(x, xp)
x0 <- xp[i]; x1 <- xp[i + 1]
y0 <- yp[i]; y1 <- yp[i + 1]
y0 + (y1 - y0) * (x - x0) / (x1 - x0)
}
# --- 2. Penalty Calculation Definitions ---
# These functions define the penalty curves for different metrics.
pen_theta_rse <- function(rse_pct) {
piecewise_linear_interp(rse_pct,
xp = c(10, 30, 75, 1000),
yp = c(0, 20, 100, 100000))
}
pen_omega_rse <- function(rse_pct) {
piecewise_linear_interp(rse_pct,
xp = c(20, 40, 65, 1000),
yp = c(0, 30, 100, 100000))
}
pen_omega_value <- function(omega_sd_val) {
omega_var_val <- omega_sd_val^2
xp <- c(0.0000, 0.0010, 0.0100, 1.0000, 1.6500, 2.3000, 3.0000, 3.5000)
yp <- c(1e6, 1e3, 0, 0, 100, 500, 1e4, 1e6)
piecewise_linear_interp(omega_var_val, xp, yp)
}
pen_omega_shrink <- function(shr_pct) {
piecewise_linear_interp(shr_pct, xp = c(40, 100), yp = c(0, 200))
}
pen_condition_number <- function(cond_num) {
piecewise_linear_interp(cond_num,
xp = c(500, 1000, 10000, 1e6),
yp = c(0, 100, 1000, 1e7))
}
# --- 3. Extract Metrics from Monolix Project ---
proj_dir <- getProjectSettings()$directory
AIC <- getEstimatedLogLikelihood()$importanceSampling[["AIC"]]
if (is.null(AIC) || is.na(AIC)) {
warning("AIC from Importance Sampling not found. Using AIC from Linearization.")
AIC <- getEstimatedLogLikelihood()$linearization[["AIC"]]
if (is.null(AIC)) stop("Could not retrieve AIC value.")
}
pop_params <- getEstimatedPopulationParameters()
std_errs <- getEstimatedStandardErrors()$stochasticApproximation
is_theta <- grepl("(_pop|beta_)", names(pop_params)) & !grepl("omega", names(pop_params))
is_omega <- grepl("^omega_", names(pop_params)) & !grepl("corr|cov", names(pop_params))
theta_rse <- std_errs[is_theta, "rse"]
omega_rse <- std_errs[is_omega, "rse"]
omega_values <- pop_params[is_omega]
shrinkage_file <- file.path(proj_dir, "individualParameters", "shrinkage.txt")
shrinkage <- if (file.exists(shrinkage_file)) {
read.csv(shrinkage_file)$shrinkage_mode
} else {
warning("Shrinkage file not found. Skipping shrinkage penalty.")
numeric(0)
}
fim_file <- file.path(proj_dir, "FisherInformation", "covarianceEstimatesSA.txt")
condition_number <- if (file.exists(fim_file)) {
fim <- as.matrix(read.csv(fim_file, header = FALSE, row.names = 1))
eigen_values <- eigen(fim, only.values = TRUE)$values
max(eigen_values) / min(eigen_values)
} else {
warning("FIM file not found. Skipping condition number penalty.")
NA
}
# --- 4. Compute Penalties ---
p_theta_rse <- sum(sapply(theta_rse, pen_theta_rse))
p_omega_rse <- sum(sapply(omega_rse, pen_omega_rse))
p_omega_val <- sum(sapply(omega_values, pen_omega_value))
p_shrinkage <- sum(sapply(shrinkage, pen_omega_shrink))
p_cond_num <- pen_condition_number(condition_number)
# --- 5. Aggregate and Return Results ---
penalty_breakdown <- c(
THETA_RSE = p_theta_rse,
OMEGA_RSE = p_omega_rse,
OMEGA_value = p_omega_val,
OMEGA_shrinkage = p_shrinkage,
Condition_number = p_cond_num
)
total_penalty <- sum(penalty_breakdown, na.rm = TRUE)
penalized_AIC <- AIC + total_penalty
penalized_AIC
}Usage
Provide this function as an argument to the findModel function:
res <- findModel(
project = "theophylline_project.mlxtran",
library = "pk",
requiredFilters = list(administration = "oral", parametrization = "clearance"),
# optional: filters / algorithm / settings ...
metric_func = calculate_penalized_aic
)