Introduction

This tutorial implements and explains in detail the DL composite indicator algorithm in a fully self-contained way.
All core functions are explicitly defined inside the document to ensure full reproducibility, although the algorithm is also available in the public GitHub repository:

https://github.com/E-Jimenez-Fernandez/dl2

The methodological foundations of the DL2 approach were introduced in:

Jiménez-Fernández, E., Sánchez, A., & Ortega-Pérez, M. (2022). Dealing with weighting scheme in composite indicators: An unsupervised distance–machine learning proposal for quantitative data. Socio-Economic Planning Sciences, 83, 101339. https://doi.org/10.1016/j.seps.2022.101339

The theoretical perspective underpinning DL is interpreted as a norm in a m dimension vector space.
This geometric viewpoint provides a set of desirable axiomatic properties for the aggregation of indicators.

Fundamental properties of norm-based aggregation

A mapping \(\| \cdot \| : \mathbb{R}^m \to \mathbb{R}\) is a norm if it satisfies:

  1. Monotonicity
    If all indicators improve (componentwise), the composite indicator must not decrease.

  2. Positive homogeneity
    Scaling all indicators by a positive scalar \(\alpha\) scales the composite indicator by \(\alpha\).
    \[ \| \alpha x \| = \alpha \| x \| \]

  3. Convexity / Subadditivity (triangle inequality)
    The indicator must reward balanced profiles over extreme or compensatory ones.
    \[ \| x + y \| \le \| x \| + \| y \| \]

  4. Continuity
    Small changes in indicators must produce small changes in the composite index.

  5. Identity of indiscernibles
    \[ \| x \| = 0 \iff x = 0 \]

These properties provide a rigorous mathematical foundation ensuring fairness, stability, and interpretability of aggregation.

Equivalence of preferences across norms

A crucial insight of the theory is that all ℓₚ norms with \(p\geq 1\) induce the same preferences:
if one alternative dominates another in ℓ₂, it will also dominate it in ℓ₃, ℓ₄, …
This provides:

  • Stability of rankings across different norms,
  • Robustness to parametric choices,
  • Coherence in normative judgments.

As highlighted in the paper, changes in the value of \(p\) do affect the shape of the indifference curves but not the ordering of alternatives.

Loss of topological information (Brouwer)

An important theoretical message of the article is that no composite indicator can preserve all topological properties of the multidimensional space (Brouwer’s invariance theorem).
This means that every 1-dimensional indicator must necessarily:

  • Lose geometric information,
  • Introduce artificial proximities or distances,
  • Distort the original space.

The task of a good index is therefore not to avoid distortion (impossible), but to minimize norm-violations and maintain meaningful geometry.

Why DL2 fits the norm-based framework

DL follows this theoretical geometry by defining the composite indicator as a weighted ℓ₁ Minkowski norm of powered distances.

A crucial requirement of the DL framework is that all weights must be strictly positive. When the composite indicator is interpreted as a weighted Minkowski norm, the strict positivity of each \(w_j > 0\) is indispensable for preserving the fundamental axioms of a norm. If any weight were equal to zero, the corresponding dimension would become “invisible” to the index, violating the identity of indiscernibles, weakening monotonicity, and destroying the topological equivalence with classical \(\ell_p\) norms. This would produce a degenerate space in which some directions have zero length and the indicator ceases to behave as a proper distance. For this reason, the improved DL2 implementation enforces strictly positive convex weights at every iteration, ensuring that the resulting metric captures the full multidimensional structure of the data and respects the normative properties required for a valid composite indicator.

\[ CI_i = \sum_j w_j |z_{ij}| \]

with the following improvements relative to classical DP2:

  1. Strict convexity of weights \(w_j \ge 0, \sum w_j = 1\)
    ensuring that the induced distance is a proper norm.

  2. Positive lower bounds on weights,
    avoiding degenerate metrics (zero components break the identity of indiscernibles).

  3. Machine-learning variable importance (MARS) to estimate weights,
    instead of sequential regression or expert assignment.

  4. Adaptive weight smoothing,
    stabilizing convergence of rankings without modifying the indicator directly.

  5. Direct update of the indicator (no index relaxation),
    preserving the interpretation of each iteration as a genuine metric.

Under this design, DL2 can be interpreted as an iterative search for a stable norm consistent with the data’s geometry.

Structure of this document

The rest of the document proceeds as follows:

The algorithm in this document follows the main conceptual steps of the published DL2 methodology and, in addition, introduces implementation improvements at the algorithmic level, namely:

  1. Strict enforcement of convex weights (non-negative, sum to one) at every iteration, ensuring that the induced distance remains a genuine metric.
  2. An adaptive smoothing scheme for weights, governed by the magnitude of the change in the composite indicator across iterations, to stabilize convergence and rankings.
  3. A direct update of the composite indicator (without index relaxation), which preserves the interpretability of each iteration as a genuine distance-based measure.

The algorithm consists of the following main components:

  1. Polarity-based normalization of the original indicators.
  2. Powered distance transformation, typically with degree \(p = 2\).
  3. Fréchet-type Minkowski index as an initial composite indicator.
  4. Weight estimation via MARS (Multivariate Adaptive Regression Splines), approximating the functional relationship between the indicator and the single variables.
  5. Adaptive smoothing of convex weights based on the change in the indicator between iterations.
  6. Iterative refinement of the composite indicator until a convergence criterion based on squared error is satisfied.

In what follows, each function is defined with detailed comments and then an illustrative example is provided.

DL2 functions with detailed comments

# -------------------------------------------------------------------
# 1. Polarity-based Min–Max normalization
# -------------------------------------------------------------------
# This function implements the normalisation scheme described in
# Jiménez-Fernández et al. (2022), Section "Normalisation", where
# each single indicator is rescaled to [0, 1] according to its polarity.
#
# - For indicators with positive polarity (higher = better), we apply:
#     (x - min) / (max - min)
#   corresponding to eq. (3) in the paper.
#
# - For indicators with negative polarity (higher = worse), we apply:
#     (max - x) / (max - min)
#   corresponding to eq. (4) in the paper.
#
# Arguments:
#   x          : numeric data frame of raw single indicators.
#   polarity   : vector of column indices with positive polarity.
#                All other columns are treated as negative by default.
#   custom_max : optional vector of user-specified maxima per variable.
#   custom_min : optional vector of user-specified minima per variable.
#
# Returns:
#   A data frame where each column belongs to [0, 1] and the direction
#   of "better performance" has been harmonised across indicators.
#
normalization <- function(x,
                          polarity   = NULL,
                          custom_max = NULL,
                          custom_min = NULL) {
  x <- as.data.frame(x)
  names_var     <- colnames(x)
  names_regions <- rownames(x)
  m             <- ncol(x)
  columns       <- seq_len(m)

  if (is.null(polarity)) {
    # If no polarity is provided, all variables are treated as negative.
    # This is conservative and forces the user to explicitly identify
    # the positively oriented indicators when needed.
    pospol <- integer(0)
    negpol <- columns
  } else {
    pospol <- polarity
    negpol <- setdiff(columns, pospol)
  }

  normdata <- matrix(0, ncol = m, nrow = nrow(x))

  # Min–max transform: positive polarity
  norm_minmax <- function(v, min_val, max_val) {
    (v - min_val) / (max_val - min_val)
  }

  # Max–min transform: negative polarity
  norm_maxmin <- function(v, min_val, max_val) {
    (max_val - v) / (max_val - min_val)
  }

  # Use empirical or user-provided bounds
  max_vals <- if (is.null(custom_max)) apply(x, 2, max) else custom_max
  min_vals <- if (is.null(custom_min)) apply(x, 2, min) else custom_min

  if (length(max_vals) != m || length(min_vals) != m) {
    stop("Length of custom_max and custom_min must match the number of columns in x.")
  }

  # Apply positive polarity normalization (eq. (3))
  for (j in pospol) {
    normdata[, j] <- norm_minmax(x[, j], min_vals[j], max_vals[j])
  }

  # Apply negative polarity normalization (eq. (4))
  for (j in negpol) {
    normdata[, j] <- norm_maxmin(x[, j], min_vals[j], max_vals[j])
  }

  normdata <- as.data.frame(normdata)
  colnames(normdata) <- names_var
  rownames(normdata) <- names_regions
  normdata
}

# -------------------------------------------------------------------
# 2. Powered distances (Minkowski components)
# -------------------------------------------------------------------
# This function takes the normalized matrix and raises each entry
# to a given power 'degrees'. In the DL2 framework, this is linked to
# the ℓ_p norm structure underlying the distance:
#
#   d_ij = (z_ij)^p
#
# For p = 2, we are in the ℓ2 (Euclidean) case, as discussed in
# Section 3 of Jiménez-Fernández et al. (2022), where DL2 is defined
# as a weighted ℓ2 metric.
#
# Arguments:
#   x       : matrix or data frame (typically the normalized data).
#   degrees : exponent p used in the Minkowski-type transformation.
#
# Returns:
#   Matrix of powered distances (same dimension as x).
#
calculateDistance <- function(x, degrees = 2) {
  x <- as.matrix(x)
  x^degrees
}

# -------------------------------------------------------------------
# 3. Fréchet-type index (initial composite indicator)
# -------------------------------------------------------------------
# The Fréchet distance used in the original paper corresponds to a 
# Minkowski-type aggregation of the powered distances:
#
#   FD_i = ( sum_j d_ij )^(1/p)
#
# where d_ij are the powered components, typically z_ij^2.
# This index is used as an initial "seed" composite indicator in
# the iterative DL2 algorithm, before learning any weights.
#
# Arguments:
#   x       : matrix of powered distances.
#   degrees : exponent p used in the Minkowski aggregation.
#
# Returns:
#   A column matrix with the Fréchet index for each observation.
#
indexFrechet <- function(x, degrees = 2) {
  x <- as.matrix(x)
  iF <- (rowSums(x))^(1 / degrees)
  iF <- matrix(iF, ncol = 1)
  colnames(iF) <- "Frechet.Index"
  iF
}

# -------------------------------------------------------------------
# 4. DP2-style composite indicator (weighted L1 of powered distances)
# -------------------------------------------------------------------
# In the original DP2 method, and in the DL2 adaptation, the distance
# is constructed by aggregating (normalised) deviations to a reference
# through a weighted structure. In this implementation we work with
# a weighted L1 norm of the powered distance components:
#
#   CI_i = sum_j |x_ij * w_j|
#
# where x_ij are the powered distances and w_j are the learned
# weights. Conceptually, this matches the partially compensatory
# nature discussed in the paper.
#
# Arguments:
#   x         : matrix of powered distances.
#   xFactores : numeric vector of weights (one per column in x).
#   iteration : integer used for labelling the resulting indicator.
#
# Returns:
#   A list with:
#     DP1 : column matrix of the composite indicator for each unit.
#     mDP : weighted distance matrix.
#
calculateDP2 <- function(x, xFactores, iteration = 1) {
  x <- as.matrix(x)
  n <- nrow(x)
  m <- ncol(x)

  # Expand the weight vector into a matrix with n rows and m columns.
  # This allows element-wise multiplication between weights and x.
  coefs <- matrix(xFactores, n, m, byrow = TRUE)
  mDP   <- x * coefs

  # Weighted L1 norm across columns.
  DP1 <- apply(mDP, 1, function(row) sum(abs(row)))
  DP1 <- matrix(DP1, ncol = 1)
  colnames(DP1) <- paste0("p1distance.", iteration)

  list(DP1 = DP1, mDP = mDP)
}

# -------------------------------------------------------------------
# 5. Adaptive smoothing parameter for weights
# -------------------------------------------------------------------
# The original DL2 algorithm uses Kendall's tau to check the stability
# of rankings between consecutive iterations. In this implementation
# we complement that idea with a *numerical* convergence perspective:
#
#   delta = RMSE(D^(t), D^(t+1))
#
# and define an adaptive smoothing parameter λ_t that:
#
#   - is large when the indicator is already stable (small delta),
#   - is small when the indicator changes are still large (big delta).
#
# This modifies the recursion for the weights:
#
#   w^(t+1) = (1 - λ_t) * w^(t) + λ_t * w_new^(t+1)
#
# and acts as a stabilizing device for the weighting scheme, reducing
# oscillations and helping convergence of both values and rankings.
#
# Arguments:
#   delta      : current RMSE between consecutive indicators.
#   base_delta : reference RMSE from first iteration.
#   lambda_min : lower bound for λ.
#   lambda_max : upper bound for λ.
#
# Returns:
#   A scalar λ_t used to smooth the new weights.
#
adaptive_lambda <- function(delta,
                            base_delta,
                            lambda_min = 0.05,
                            lambda_max = 0.8) {
  if (is.na(delta) || is.na(base_delta) || base_delta == 0) {
    # Fallback to mid-range if baseline is degenerate
    return((lambda_min + lambda_max) / 2)
  }

  # Relative size of current variation compared to baseline
  ratio <- delta / base_delta
  ratio <- max(min(ratio, 1), 0)

  # Large delta -> lambda close to lambda_min (conservative updates)
  # Small delta -> lambda close to lambda_max (more aggressive updates)
  lambda <- lambda_max - (lambda_max - lambda_min) * ratio
  lambda
}

# -------------------------------------------------------------------
# 6. Variable weights via MARS (convex weighting scheme)
# -------------------------------------------------------------------
# This function corresponds to the step "Selecting the best set of 
# disjoint polynomials" and "Computing the weights of DL2" in
# Jiménez-Fernández et al. (2022).
#
# Conceptual steps:
#   1) Fit a MARS model:  Compind ~ z_1 + ... + z_m
#      where Compind is either the initial Frechet index or the
#      composite indicator from the previous iteration.
#   2) Use earth::evimp() to obtain RSS-based variable importance.
#   3) Transform these importances into a convex vector of weights:
#        w_j >= 0,  sum_j w_j = 1.
#
# This approximates the role of the Variable Importance (VI/FIRM)
# mechanism described in the paper, but using the RSS-based measure
# provided by the earth package as a practical and efficient proxy.
#
# Importantly, in this implementation we force:
#
#   - strictly non-negative weights,
#   - renormalisation to sum to one,
#   - a small positive weight for any variable with zero importance,
#     to preserve the metric property (all ω_j > 0).
#
# Arguments:
#   x       : normalized data frame of indicators.
#   Compind : current composite indicator (response).
#   degrees : max interaction degree for MARS.
#   nfold   : number of CV folds.
#   trace   : verbosity of earth().
#
# Returns:
#   A list with:
#     [[1]] : numeric vector of weights aligned with columns of x.
#     [[2]] : character vector of variable names in original order.
#     [[3]] : evimp object with raw importance measures.
#
weightFactors <- function(x,
                          Compind,
                          degrees,
                          nfold = 5,
                          trace = .5) {

  x      <- as.data.frame(x)
  dat_os <- x
  m      <- ncol(x)

  # Attach the composite indicator as response
  x <- dplyr::mutate(x, Cind = as.numeric(Compind))

  # MARS model capturing possibly nonlinear and interaction effects
  mymodel <- earth::earth(
    Cind ~ .,
    data   = x,
    degree = degrees,
    nfold  = nfold,
    trace  = trace
  )

  # RSS-based importance for each basis variable
  p2       <- earth::evimp(mymodel, trim = FALSE)
  Variable <- row.names(p2)

  weights  <- as.numeric(p2[, "rss"])
  Variable <- sub("-unused", "", Variable)

  # Degenerate case: no meaningful importance structure
  if (all(is.na(weights)) || sum(weights, na.rm = TRUE) == 0) {
    warning("All calculated weights are 0 or NA. Assigning equal weights.")
    weights <- rep(1 / m, m)
    res <- data.frame(t(weights))
    colnames(res) <- colnames(dat_os)
    res.ord <- as.numeric(res[colnames(dat_os)])
    return(list(res.ord, colnames(dat_os), p2))
  }

  # Replace zero importances with a small positive value to avoid
  # null weights which would break the metric property (ω_j > 0).
  mini <- min(weights[weights > 0], na.rm = TRUE)
  for (j in seq_along(weights)) {
    if (is.na(weights[j])) {
      stop("Missing values were found in the weights.")
    } else if (weights[j] == 0) {
      weights[j] <- mini / length(weights)
    }
  }

  # Enforce convexity: weights >= 0 and sum(weights) = 1
  weights <- pmax(weights, 0)
  if (sum(weights) == 0) {
    weights <- rep(1 / length(weights), length(weights))
  } else {
    weights <- weights / sum(weights)
  }

  # Align weights with the original column order of the indicators
  res <- data.frame(t(weights))
  colnames(res) <- Variable

  faltan <- setdiff(colnames(dat_os), colnames(res))
  if (length(faltan) > 0) {
    res[, faltan] <- 0
  }
  res <- res[, colnames(dat_os), drop = FALSE]

  res.ord   <- as.numeric(res[1, ])
  important <- list(res.ord, colnames(dat_os), p2)
  important
}

# -------------------------------------------------------------------
# 7. Main DL2 algorithm (improved implementation)
# -------------------------------------------------------------------
# This function orchestrates the entire DL2 procedure:
#
#   1) Optional polarity-based normalization (eqs. (3)–(4)).
#   2) Correlation analysis (diagnostic only, via corrgram).
#   3) Construction of powered distances (eq. (5) components).
#   4) Initial Fréchet-type index (eq. (6)).
#   5) Iterative loop:
#        a) MARS-based weight estimation (Section 3, MARS step).
#        b) Provisional indicator using raw new weights.
#        c) Computation of an RMSE-based delta.
#        d) Adaptive selection of λ_t (adaptive_lambda()).
#        e) Smoothing of weights ONLY (no relaxation on the index).
#        f) Final composite indicator with smoothed weights.
#        g) Convergence test based on squared error.
#
# Compared to the baseline algorithm in the paper:
#
#   - We enforce convex weights at every iteration, ensuring that
#     the induced distance is a proper metric (Section 4 of the paper).
#   - We smooth the *weights* rather than the index, which preserves
#     interpretability of each composite indicator as a genuine
#     distance-based benchmark.
#   - We use a numeric squared-error criterion instead of Kendall's τ,
#     though both aim at stabilizing the ranking / values of the units.
#
# Arguments:
#   x              : data frame of raw indicators.
#   polarity       : numeric vector of columns with positive polarity.
#   err            : error tolerance for stopping (sum of squared diffs).
#   iterations     : maximum number of iterations allowed.
#   degrees        : exponent used in calculateDistance and indexFrechet.
#   nfold          : number of CV folds in MARS.
#   trace          : verbosity of earth().
#   use_normat     : if TRUE, internal normalization is applied.
#   lambda_w_fixed : if not NULL, use a fixed smoothing parameter;
#                    otherwise, use adaptive_lambda().
#
# Returns:
#   A list with:
#     normat           : normalized data (if use_normat = TRUE).
#     Weights          : final convex weight vector.
#     derror           : error trajectory along iterations.
#     frechet          : initial Fréchet index.
#     corrgr           : corrgram object for diagnostics.
#     Comp.Indicator   : final composite indicator (distance-based).
#     Comp.Indicator_1 : same as above (for compatibility).
#     iteration        : number of iterations performed.
#     dist.in          : initial powered distance matrix.
#     dist.fin         : final weighted distance matrix.
#     base_delta       : baseline RMSE of the first iteration.
#
dl2 <- function(x,
                polarity,
                err,
                iterations,
                degrees,
                nfold = 5,
                trace = .5,
                use_normat = TRUE,
                lambda_w_fixed = NULL) {

  if (!is.data.frame(x)) {
    warning("The argument 'x' must be a data.frame object.")
    x <- as.data.frame(x)
  }

  # Step 1: optional normalization
  if (use_normat) {
    normat <- normalization(x, polarity)
  } else {
    normat <- x
  }

  # Step 2: correlation structure (descriptive)
  corrgr <- corrgram::corrgram(
    normat,
    order       = TRUE,
    lower.panel = corrgram::panel.shade,
    upper.panel = corrgram::panel.pie,
    text.panel  = corrgram::panel.txt,
    main        = "Correlation between features"
  )

  # Step 3: powered distances (Minkowski components)
  mDif <- calculateDistance(normat, degrees = degrees)

  # Step 4: initial Fréchet-type index (unweighted)
  dps       <- indexFrechet(mDif, degrees = degrees)
  dp2_aux   <- dps
  iteration <- 1
  error.d   <- numeric()
  weights   <- NULL
  calcDP2_res <- NULL
  base_delta  <- NA

  # Step 5: iterative refinement of weights and indicator
  repeat {
    cat("Iteration", iteration, "\n")

    # 5a) New weights from MARS using current indicator as response
    new_weights <- weightFactors(
      x       = normat,
      Compind = dps,
      degrees = degrees,
      nfold   = nfold,
      trace   = trace
    )[[1]]

    # 5b) Provisional indicator using raw new weights
    tmp_res <- calculateDP2(mDif, new_weights, iteration = iteration)
    d2_prov <- tmp_res$DP1

    # RMSE-based delta (used to control smoothing)
    error_prov <- sum((dps - d2_prov)^2)
    delta      <- sqrt(error_prov / length(dps))

    # Baseline delta from the first iteration
    if (iteration == 1) {
      base_delta <- ifelse(delta == 0, .Machine$double.eps, delta)
    }

    # 5c) Choose smoothing parameter for the weights
    if (is.null(lambda_w_fixed)) {
      lambda_w_cur <- adaptive_lambda(delta, base_delta)
    } else {
      lambda_w_cur <- lambda_w_fixed
    }

    # 5d) Smooth the weights (preserving convexity)
    if (is.null(weights)) {
      weights <- new_weights
    } else {
      weights <- (1 - lambda_w_cur) * weights + lambda_w_cur * new_weights
    }

    # 5e) Recompute the composite indicator with smoothed weights
    calcDP2_res <- calculateDP2(mDif, weights, iteration = iteration)
    d2_ite      <- calcDP2_res$DP1

    # 5f) Compute error and apply stopping rule
    error <- sum((dps - d2_ite)^2)
    error.d[iteration] <- error

    if ((error < err) || (iteration >= iterations)) {
      break
    }

    iteration <- iteration + 1

    # 5g) Update indicator WITHOUT relaxing it: we directly use the
    # new composite indicator as seed for the next iteration. This
    # keeps each step interpretable as a proper distance index.
    dps <- d2_ite
  }

  invisible(list(
    normat           = if (use_normat) normat else NULL,
    Weights          = weights,
    derror           = error.d,
    frechet          = dp2_aux,
    corrgr           = corrgr,
    Comp.Indicator   = d2_ite,
    Comp.Indicator_1 = dps,
    iteration        = iteration,
    dist.in          = mDif,
    dist.fin         = if (!is.null(calcDP2_res)) calcDP2_res$mDP else NULL,
    base_delta       = base_delta
  ))
}

Example dataset

We now illustrate the algorithm with a simple simulated data set comprising five indicators: three with positive polarity (Income, Education, Health) and two with negative polarity (Pollution, Crime).

set.seed(123)

X <- data.frame(
  Income     = rnorm(50, 30000, 5000),
  Education  = rnorm(50, 70, 10),
  Health     = rnorm(50, 0.8, 0.1),
  Pollution  = rnorm(50, 30, 5),
  Crime      = rnorm(50, 4, 1)
)

head(X)
##     Income Education    Health Pollution    Crime
## 1 27197.62  72.53319 0.7289593  33.93869 6.198810
## 2 28849.11  69.71453 0.8256884  33.84521 5.312413
## 3 37793.54  69.57130 0.7753308  31.66101 3.734855
## 4 30352.54  83.68602 0.7652457  24.95812 4.543194
## 5 30646.44  67.74229 0.7048381  29.40274 3.585660
## 6 38575.32  85.16471 0.7954972  28.59802 3.523753

Define polarity

polarity_vector <- c(1, 2, 3)

Run the DL2 algorithm (improved version)

result <- dl2(
  x          = X,
  polarity   = polarity_vector,
  err        = 1e-3,
  iterations = 100,
  degrees    = 2,
  nfold      = 5,
  trace      = 0,
  use_normat = TRUE
)

## Iteration 1 
## Iteration 2 
## Iteration 3 
## Iteration 4 
## Iteration 5 
## Iteration 6 
## Iteration 7 
## Iteration 8 
## Iteration 9 
## Iteration 10 
## Iteration 11 
## Iteration 12 
## Iteration 13

Normalized data

head(result$normat)
##      Income Education    Health Pollution     Crime
## 1 0.3400113 0.5698846 0.3233146 0.5389689 0.0000000
## 2 0.4198788 0.5071992 0.5562082 0.5430764 0.2490824
## 3 0.8524394 0.5040137 0.4349628 0.6390463 0.6923847
## 4 0.4925860 0.8179183 0.4106810 0.9335599 0.4652370
## 5 0.5067991 0.4633375 0.2652382 0.7382710 0.7343092
## 6 0.8902471 0.8508035 0.4835173 0.7736287 0.7517054

Final weights

result$Weights
## [1] 0.12048112 0.01922913 0.09614561 0.19650951 0.56763463

Convergence plot

df <- data.frame(
  Iteration = seq_along(result$derror),
  Error = result$derror
)

ggplot(df, aes(Iteration, Error)) +
  geom_line(size = 1.2) +
  geom_point() +
  theme_minimal() +
  ggtitle("DL2 Convergence Path (improved implementation)") +
  scale_y_log10()

Composite indicator

head(result$Comp.Indicator)
##      p1distance.13
## [1,]    0.08730744
## [2,]    0.14910581
## [3,]    0.46299542
## [4,]    0.35244026
## [5,]    0.45501776
## [6,]    0.57024263

Plot of the composite indicator

ci_df <- data.frame(
  Observation = seq_len(nrow(X)),
  CI = as.numeric(result$Comp.Indicator)
)

ggplot(ci_df, aes(x = Observation, y = CI)) +
  geom_point(size = 2) +
  geom_line() +
  theme_minimal() +
  ggtitle("Composite Indicator (DL2 Algorithm - improved version)")

Interpretation of Results

The numerical illustration included in this document provides a complete example of how the DL2 algorithm behaves when applied to a synthetic dataset composed of five indicators: Income, Education, Health, Pollution, and Crime. The first three indicators have positive polarity, meaning that higher values represent better performance, whereas Pollution and Crime have negative polarity, meaning that higher values worsen the composite assessment.

The discussion below interprets the results obtained in the example.

1. Normalized Data

The first output of the algorithm is the normalized data matrix. All variables have been transformed into the interval \([0,1]\) while respecting their polarity:

  • Income, Education, Health: higher values correspond to higher normalized scores.
  • Pollution, Crime: higher raw values correspond to lower normalized scores.

This harmonization is essential because it ensures that all indicators contribute in the same conceptual direction to the composite indicator. Normalization also removes the effect of heterogeneous units of measurement, which is critical in multidimensional assessment (e.g., euros, years, ppm, crime counts).

2. Final Weights

The vector returned in result$Weights represents the final convex weights estimated by the DL2 algorithm. These weights are:

  • Data-driven (no subjective input),
  • Strictly non-negative, and
  • Summed to one, thus forming a convex combination.

The DL2-weighting scheme assigns greater weight to indicators that provide more discriminatory power across observations. In the context of the example, one typically observes:

  • Indicators with greater variability and stronger alignment with the Fréchet-type distance structure (e.g., Income or Education) tend to receive higher weights.
  • Indicators that add redundant or weak information relative to others (e.g., if Pollution and Crime are weakly correlated with the main latent structure in the sample) tend to receive smaller weights.

This is consistent with the rationale in Jiménez-Fernández et al. (2022), where DL2 penalizes informational redundancy and rewards variables that meaningfully expand the multidimensional space.

The improved implementation provided here enhances the original methodology by:

  • Forcing convexity at every iteration,
  • Preventing zero-weights (which would break metric properties),
  • Introducing an adaptive smoothing mechanism that stabilizes convergence.

3. Convergence Behaviour

The plot titled “DL2 Convergence Path (improved implementation)” displays the squared error between consecutive iterations of the composite indicator on a logarithmic scale.

The key observations are:

  • The sequence of errors decreases monotonically or almost monotonically, indicating stable convergence.
  • The adaptive smoothing of weights ensures that early iterations change more slowly when the indicator is unstable, and more quickly once stability is achieved.
  • The decline of the error indicates that the algorithm is successfully reaching a fixed point where weights and indicator values are mutually consistent.

This behaviour aligns with the theoretical motivation of DL2, which searches for a stable configuration of distances and weights under the unsupervised learning structure of MARS.

4. Final Composite Indicator

The object result$Comp.Indicator contains the final composite scores for each of the 50 simulated observations.

The plot of the composite indicator shows:

  • A smooth ranking of units from best to worst according to the multidimensional performance encoded in the indicators.
  • Units with high Income, Education, and Health — combined with low Pollution and Crime — appear toward the top.
  • Units with poor socio-economic and environmental profiles appear toward the bottom.

Because DL2 is distance-based, the indicator values can be interpreted in relation to how far each unit is from an abstract “ideal point” in the normalized space, weighted according to the importance learned through MARS.

5. Summary Interpretation

Overall, the example demonstrates that:

  • DL2 produces rankings consistent with data structure and polarity,
  • The learned weights reflect the intrinsic importance of each indicator,
  • The convergence path is stable due to adaptive smoothing,
  • The final composite indicator forms a coherent distance-based assessment.

This small experiment illustrates the core ideas presented in Jiménez-Fernández et al. (2022), showing how the DL2 method integrates normalization theory, distance geometry, and unsupervised machine learning to construct a robust composite indicator.