Articles · Get Started

Getting started with mvmetrics

This vignette walks through a complete evaluation workflow using simulated hydroclimatic data calibrated to the Valle del Cauca parameters from Arango Londoño (2026).

Installation

# From GitHub
install.packages("devtools")
devtools::install_github("darango/mvmetrics")

# Suggested packages
install.packages(c("energy", "ggplot2", "MASS"))
library(mvmetrics)

Data requirements

All functions expect three core objects:

Object	Dimensions	Description
obs	[N × K]	Observed values. Columns 1:K_c = continuous; column K = binary (0/1).
pred	[N × K]	Point predictions. Binary column = predicted probability ∈ (0,1).
draws	[N × K × B]	Posterior predictive draws. B = number of samples per observation.

Scale contract for Metric E: obs, pred, and draws must be in the original physical scale (not normalised). Metrics B and C require normalised inputs — normalise_matrix() handles this internally when called through evaluate_models().

Simulated example

The following generates data for K = 5 variables (Tmin, Tmax, HR, Rad, Pbin) with a Clausius-Clapeyron correlation structure, two competing models, and B = 50 posterior draws.

library(MASS); set.seed(2026)
N <- 500; K <- 5; K_c <- 4; B <- 50

# Empirical parameters (Chapter 7, Valle del Cauca)
mu_real <- c(Tmin=18.9, Tmax=30.6, HR=80.7, Rad=3.8)
sd_real <- c(Tmin=2.5,  Tmax=4.0,  HR=8.0,  Rad=0.8)

# Clausius-Clapeyron residual correlation (positive-definite)
R_cont <- matrix(c(
   1.000,  0.310, -0.220, -0.180,
   0.310,  1.000, -0.505, -0.535,
  -0.220, -0.505,  1.000,  0.420,
  -0.180, -0.535,  0.420,  1.000),
  nrow=4, byrow=TRUE)

# See use_simulation.R (included in package) for the full simulation

evaluate_models()

The main wrapper accepts a named list of models and returns a structured report.

report <- evaluate_models(
  obs        = obs,
  pred_list  = list(M1 = pred_m1, M3 = pred_m3),
  draws_list = list(M1 = draws_m1, M3 = draws_m3),
  mu_train   = mu_real,   # training means  (Metric E CV weights)
  sd_train   = sd_real,   # training SDs    (Metric E normalisation)
  K_c        = 4,
  alpha_D    = 0.5,       # marginal/dependence weight for Metric D
  alpha_E    = NULL,      # NULL = adaptive alpha* via dcor test
  run_dcor   = FALSE,     # TRUE for final analyses
  bootstrap  = TRUE,
  B_boot     = 500,       # 200 for quick runs, 500 for publication
  verbose    = TRUE
)

Reading the output

The returned mvmetrics_report object has five slots:

Slot	Description
$scores	data.frame — one row per model, columns for each metric (A, B, C, D, D_marg, D_dep, E, E_marg, E_dep, E_alpha_star)
$rankings	data.frame — rank per model per metric (1 = best = lowest score)
$recommendation	character — recommended metric family with rationale
$bootstrap_ci	data.frame — 95% CIs and CI overlap indicator (NULL if bootstrap=FALSE)
$marginal_detail	list — per-variable RMSE, Log-Loss, Accuracy for each model

Visualisation

# Side-by-side bars for all metrics
plot_metrics(report, type = "bars")

# Bootstrap forest plot (requires bootstrap=TRUE)
plot_metrics(report, type = "forest",
             metrics = c("A","B","D","D_dep","E","E_dep"))

# Ranking heatmap
plot_metrics(report, type = "heatmap")

# S_marg vs S_dep stacked bars (Metrics D and E)
plot_metrics(report, type = "decomp")

Individual metric functions

# RMSE per variable (named vector)
marginal_rmse(obs[,1:4], pred[,1:4])

# Metric D alpha sensitivity
for (a in c(0.2, 0.5, 0.8)) {
  res <- metric_D_mdd(obs, pred, draws, K_c=4, alpha=a)
  cat(sprintf("alpha=%.1f  S_dep=%.4f\n", a, res$S_dep))
}

# Metric E — inspect CV weights and adaptive alpha*
res_E <- metric_E_cvwmd(obs, pred, draws,
  mu_train=mu_real, sd_train=sd_real, K_c=4)
res_E$w_cont     # CV-derived weights per continuous variable
res_E$alpha_star # adaptive alpha* (0.5 if dcor not significant)
res_E$S_E_dep    # dependence score

Bootstrap confidence intervals

When bootstrap=TRUE, evaluate_models() calls bootstrap_metrics() internally. For two models, a pairwise CI overlap test is computed automatically.

Interpretation: Non-overlapping 95% CIs for S_dep (Metric D) and S_E_dep (Metric E) provide informal evidence that the dependence recovery difference is statistically distinguishable at the available sample size — without assuming distributional form.

# Access CI table directly
report$bootstrap_ci[, c("Model","Metric","CI_lo","CI_hi","CIs_overlap")]

# Run bootstrap standalone (e.g. for more resamples)
ci_tbl <- bootstrap_metrics(
  obs, pred_list, draws_list,
  mu_train=mu_real, sd_train=sd_real,
  K_c=4, B_boot=1000
)