Introduction to orthoMTL: Multi-Task Survival Analysis

The Problem

In oncology and other therapeutic areas, genomic features associated with patient outcomes often have time-varying effects. A mutation that predicts early disease progression may be irrelevant for long-term survival — and vice versa. Standard survival models like the Cox proportional hazards model assume that each feature’s effect is constant over time. When this assumption is violated, important signals can be missed.

orthoMTL addresses this by reframing survival analysis as a multi-task learning problem. Instead of fitting a single model for time-to-event, we define binary classification tasks at multiple time thresholds: “Is this patient progression-free at 4 months? At 6 months? At 10? At 15?” Each threshold becomes a task, and the model learns a separate coefficient vector for each — while encouraging the coefficient vectors to be structured through an orthogonality penalty.

The Objective Function

The optimisation problem solved by orthoMTL() is:

\[ \min_W \; \frac{1}{2n}\|XW - Y\|^2_{\text{obs}} \;+\; \lambda\Bigl[\frac{1-\alpha}{2}\,\Omega_K(W)^2 \;+\; \alpha\,\|W\|_1\Bigr], \qquad \Omega_K(W)^2 = \sum_{s,t} K_{st}\,|W_s^\top W_t| \]

where:

\(X\) is the \(n \times p\) feature matrix
\(W\) is the \(p \times T\) coefficient matrix (\(T\) tasks)
\(Y\) is the \(n \times T\) label matrix (1 = progression-free, 0 = event, NA = censored)
\(\|\cdot\|_{\text{obs}}\) means the loss is only computed on non-censored entries
\(K\) is a constraint matrix encoding task relationships — nearby thresholds share more support, distant thresholds are pushed toward orthogonality
\(\lambda\) controls the overall penalty strength
\(\alpha \in [0, 1]\) is the elastic-net mixing: \(\alpha = 0\) gives the pure orthogonality penalty, \(\alpha = 1\) gives a pure Lasso

Simulated Data

We begin with simulated data where the ground truth is known. The simulate_mtl() function generates survival data with five types of time-varying effects: early, late, constant, increasing, and decreasing.

set.seed(42)

sim <- simulate_mtl(
  n = 300,
  p = 15,
  n_signals = 5,
  thresholds = c(4, 8, 14, 20),
  effect_strength = 1.2
)
sim
#> Simulated multi-task survival data
#> Patients: 300 | Features: 16 (incl. treatment)
#>   Continuous: 1 | Binary: 14 | Treatment: 1
#>   Signal features: 5 | Null features: 10 
#> Tasks: 4, 8, 14, 20 
#> Event rate: 48.3 %
#> Median survival time: 6.68 
#> ---
#> Signal features and effect types:
#>    mut_2 : early 
#>    mut_5 : late 
#>    mut_6 : constant 
#>    mut_11 : switch 
#>    mut_14 : decreasing 
#> ---
#> Treatment effect: -0.3 ( log-hazard , constant across tasks)

The signal features and their true temporal effect patterns:

# True coefficients for signal features
gt <- sim$ground_truth
signal_coefs <- gt$coefficients[gt$signal_features, ]
signal_coefs
#>           4          8         14   20
#> mut_2   1.2  0.8000000  0.4000000  0.0
#> mut_5   0.0  0.4000000  0.8000000  1.2
#> mut_6   1.2  1.2000000  1.2000000  1.2
#> mut_11  1.2 -0.1856406 -0.7595918 -1.2
#> mut_14 -1.2 -0.8000000 -0.4000000  0.0

# Effect types
gt$effect_types
#>        mut_2        mut_5        mut_6       mut_11       mut_14 
#>      "early"       "late"   "constant"     "switch" "decreasing"

Data Preparation

orthoMTL requires three matrices derived from the survival data.

Longitudinal binary labels

Each patient’s survival time and event indicator are converted into a binary label at each threshold. Patients who experienced an event before a threshold are labelled 0. Patients censored before a threshold have an unknown label (NA). All others are labelled 1 (progression-free).

thresholds <- sim$thresholds
Y <- create_longitudinal_labels(sim$SurvTime, sim$Event, thresholds)
head(Y, 10)
#>        4  8 14 20
#>  [1,]  0  0  0  0
#>  [2,]  1 NA NA NA
#>  [3,]  0  0  0  0
#>  [4,]  1 NA NA NA
#>  [5,]  0  0  0  0
#>  [6,]  1 NA NA NA
#>  [7,] NA NA NA NA
#>  [8,]  0  0  0  0
#>  [9,] NA NA NA NA
#> [10,]  1  1  0  0

Censoring indicator matrix

A binary matrix marking which labels are observed (1) versus censored (0). This is used by the solver to mask censored entries in the loss.

W <- create_indicator_matrix(Y)
cat("Proportion observed per threshold:\n")
#> Proportion observed per threshold:
colMeans(W)
#> [1] 0.8533333 0.7633333 0.6566667 0.5500000

Constraint matrix

The diffusion constraint matrix \(K\) encodes the prior that nearby thresholds should share coefficient support while distant thresholds can diverge.

K <- create_constraint_matrix(length(thresholds))
K
#>      [,1] [,2] [,3] [,4]
#> [1,]  0.5  0.0  1.0  2.0
#> [2,]  0.0  0.5  0.0  1.0
#> [3,]  1.0  0.0  0.5  0.0
#> [4,]  2.0  1.0  0.0  0.5

Model Fitting

We fit an initial model with a single set of hyperparameters.

fit <- orthoMTL(
  X = sim$X, Y = Y,
  lambda = 1e-3, step_size = 0.5,
  K = K, survival = TRUE, censored.mat = W,
  alpha = 0.5
)
summary(fit)
#> orthoMTL model: 16 features, 4 tasks ( survival )
#> Converged: yes ( 633 iterations, objective: 0.3112 )
#> Penalty: lambda = 0.001, alpha = 0.5 
#> Constraints: disjoint = FALSE 
#> ---
#> Hyperparameters:
#>    lambda = 0.001 
#>    alpha = 0.5 
#>    schedule = sqrt 
#>    step_size = 0.5 
#>    tol = 1e-05 
#>    stop_no_improve = 100 
#>    max_iter = 1e+06 
#>    seed = 
#>    disjoint = FALSE 
#>    logistic = FALSE 
#>    survival = TRUE 
#> ---
#> Coefficients: 16 x 4 matrix
#>   Sparsity: 0.0 %
#>   Range: [ -0.289 , 0.4445 ]
#> ---
#> Top 5 features (mean |coefficient| across tasks):
#>    mut_14 : 0.357 
#>    treatment : 0.249 
#>    mut_6 : 0.1932 
#>    mut_7 : 0.1819 
#>    mut_4 : 0.171 
#> ---
#> Task thresholds: 4, 8, 14, 20

The coefficient heatmap shows each feature’s weight across time thresholds. Blue indicates a protective effect (associated with longer progression-free survival); red indicates a risk effect.

plot_heatmap(fit)

Even with untuned hyperparameters, temporal patterns begin to emerge.

Cross-Validation

We use cv_orthoMTL() to search over a grid of hyperparameters and select the configuration with the highest cross-validated C-index.

n_tasks <- length(sim[["thresholds"]])
folds <- rep(1:5, length.out = nrow(sim[["X"]]))

cv_res <- cv_orthoMTL(
  X.train   = sim[["X"]],
  Y.train   = Y,
  W.train   = W,
  K         = K,
  lambdas   = c(1e-5, 1e-4, 1e-3),
  alphas    = c(0, 0.5, 1),
  stepsizes = c(1, 2, 5),
  diag_vals = c(0.5, n_tasks, 2 * n_tasks),
  survival  = TRUE,
  folds     = folds,
  n_cores   = 1,
  seed      = 42,
  verbose   = FALSE
)

print(cv_res)
#> Cross-validation results for orthoMTL
#> Configurations tested: 81 
#> Folds: 5 
#> ---
#> Best configuration:
#>   lambda   = 1e-04 
#>   alpha    = 1 
#>   stepsize = 2 
#>   diag_val = 0.5 
#>   CV cindex = 0.5871
#> ---
#> Top 5 configurations:
#>  lambda alpha stepsize diag_val  cv_score
#>   1e-04   1.0        2      0.5 0.5870769
#>   1e-04   1.0        2      4.0 0.5870769
#>   1e-04   1.0        2      8.0 0.5870769
#>   1e-04   0.5        2      0.5 0.5853762
#>   1e-05   0.0        2      0.5 0.5851236

Final Model

We refit on the full dataset using the best hyperparameters from cross-validation.

best <- cv_res$best
K_final <- K
diag(K_final) <- best$diag_val

fit_final <- orthoMTL(
  X = sim$X, Y = Y,
  lambda    = best$lambda,
  alpha     = best$alpha,
  step_size = best$stepsize,
  K         = K_final,
  survival  = TRUE,
  censored.mat = W
)
summary(fit_final)
#> orthoMTL model: 16 features, 4 tasks ( survival )
#> Converged: yes ( 155 iterations, objective: 0.3067 )
#> Penalty: lambda = 1e-04, alpha = 1 
#> Constraints: disjoint = FALSE 
#> ---
#> Hyperparameters:
#>    lambda = 1e-04 
#>    alpha = 1 
#>    schedule = sqrt 
#>    step_size = 2 
#>    tol = 1e-05 
#>    stop_no_improve = 100 
#>    max_iter = 1e+06 
#>    seed = 
#>    disjoint = FALSE 
#>    logistic = FALSE 
#>    survival = TRUE 
#> ---
#> Coefficients: 16 x 4 matrix
#>   Sparsity: 0.0 %
#>   Range: [ -0.2961 , 0.4425 ]
#> ---
#> Top 5 features (mean |coefficient| across tasks):
#>    mut_14 : 0.3619 
#>    treatment : 0.2459 
#>    mut_6 : 0.2034 
#>    mut_7 : 0.1843 
#>    mut_8 : 0.1728 
#> ---
#> Task thresholds: 4, 8, 14, 20

plot_heatmap(fit_final)

The task correlation map shows how similar the coefficient profiles are across thresholds. Low distance (red) between adjacent thresholds indicates smooth temporal evolution; high distance (blue) indicates divergent coefficient structures.

plot_correlation(fit_final)

Bootstrap Inference

To assess whether each feature’s coefficients are distinguishable from noise, we compare bootstrapped models (resampled data, real signal) against null models (permuted outcomes, no signal).

boot_res <- bootstrap_orthoMTL(
  X = sim$X, Y = Y,
  lambda    = best$lambda,
  alpha     = best$alpha,
  step_size = best$stepsize,
  K         = K_final,
  survival  = TRUE,
  censored.mat = W,
  n_repeats = 200,
  n_cores   = 1,
  verbose   = FALSE
)

print(boot_res)
#> Bootstrap inference for orthoMTL
#> Repeats: 200 (real) + 200 (null permutation)
#> Features: 16 | Tasks: 4 
#> ---
#> Real models -- objective:
#>   mean = 0.2865  range = [ 0.2314 , 0.3376 ]
#> Null models -- objective:
#>   mean = 0.338  range = [ 0.307 , 0.3783 ]

We select four features for detailed inspection: one true signal of each pattern, and one null feature.

signal_feats <- gt[["signal_features"]]
null_feats <- gt[["null_features"]]
effect_types <- gt[["effect_types"]]

selected <- c(
  signal_feats[effect_types == "switch"][1],
  signal_feats[effect_types == "constant"][1],
  signal_feats[effect_types == "early"][1],
  null_feats[1]
)

cat("Selected features:\n")
#> Selected features:
cat("  Switch signal:  ", selected[1], "(Cox should miss)\n")
#>   Switch signal:   mut_11 (Cox should miss)
cat("  Constant signal:", selected[2], "(Cox should find)\n")
#>   Constant signal: mut_6 (Cox should find)
cat("  Early signal:   ", selected[3], "(Cox may dilute)\n")
#>   Early signal:    mut_2 (Cox may dilute)
cat("  Null feature:   ", selected[4], "(neither should find)\n")
#>   Null feature:    cont_1 (neither should find)

suppressWarnings(
  plot_bootstrap(boot_res, features = selected)
)
#> [[1]]

For the true signal features, the real coefficients (coloured line) separate clearly from the null distribution (grey line). The temporal patterns are visible: the early feature’s effect is strongest at early thresholds and fades; the late feature’s effect emerges at later thresholds; the constant feature is stable across all thresholds.

For the null feature, real and null distributions overlap — the model correctly assigns it no meaningful effect.

Comparison with Cox Proportional Hazards

The Cox model assumes each feature’s effect is constant over time. We fit a Cox elastic-net model using glmnet and compare which features it detects.

library(survival)
library(glmnet)
#> Loading required package: Matrix
#> Loaded glmnet 4.1-10

surv_obj <- Surv(time = sim[["SurvTime"]], event = sim[["Event"]])

# Cross-validate alpha (mixing parameter)
alphas <- seq(0, 1, by = 0.1)
cv_scores <- numeric(length(alphas))

for (i in seq_along(alphas)) {
  set.seed(42)
  cvfit <- cv.glmnet(
    x = sim[["X"]],
    y = surv_obj,
    family = "cox",
    type.measure = "C",
    alpha = alphas[i]
  )
  cv_scores[i] <- cvfit[["cvm"]][cvfit[["index"]]["min", ]]
}

best_alpha <- alphas[which.max(cv_scores)]
cat("Best alpha:", best_alpha, "(CV C-index:", max(cv_scores), ")\n")
#> Best alpha: 0.7 (CV C-index: 0.6920214 )

# Refit with best alpha
set.seed(42)
cox_fit <- cv.glmnet(
  x = sim[["X"]],
  y = surv_obj,
  family = "cox",
  type.measure = "C",
  alpha = best_alpha
)

cox_coefs <- as.numeric(coef(cox_fit, s = "lambda.min"))
names(cox_coefs) <- colnames(sim[["X"]])

# Build comparison table
ortho_mean_abs <- apply(abs(coef(fit_final)), 1, mean)

comparison <- data.frame(
  feature = gt$signal_features,
  effect_type = as.character(gt$effect_types),
  orthoMTL_mean_abs = round(ortho_mean_abs[gt$signal_features], 4),
  cox_coef = round(cox_coefs[gt$signal_features], 4),
  cox_detected = cox_coefs[gt$signal_features] != 0,
  stringsAsFactors = FALSE
)

cat("Signal feature detection comparison:\n\n")
#> Signal feature detection comparison:
print(comparison, row.names = FALSE)
#>  feature effect_type orthoMTL_mean_abs cox_coef cox_detected
#>    mut_2       early            0.1496   0.9449         TRUE
#>    mut_5        late            0.1016   0.6400         TRUE
#>    mut_6    constant            0.2034   1.1669         TRUE
#>   mut_11      switch            0.0647   0.3777         TRUE
#>   mut_14  decreasing            0.3619  -0.7275         TRUE

The Cox proportional hazards model is a powerful and well-established tool for survival analysis. In terms of overall discrimination (C-index), Cox often matches or exceeds orthoMTL — it directly optimises the survival likelihood, while orthoMTL solves a regression problem on binary labels.

The value of orthoMTL is not in replacing Cox but in complementing it. Cox produces a single coefficient per feature — an average effect across the entire follow-up. orthoMTL produces a coefficient per feature per timepoint, revealing temporal dynamics that a single number cannot capture.

# Pick the switch feature — most interesting temporal pattern
switch_feat <- signal_feats[effect_types == "switch"][1]

# What Cox sees: one number
cat(switch_feat, "— Cox coefficient:", 
    round(cox_coefs[switch_feat], 3), "\n")
#> mut_11 — Cox coefficient: 0.378

# What orthoMTL sees: a trajectory
cat(switch_feat, "— orthoMTL coefficients:\n")
#> mut_11 — orthoMTL coefficients:
print(round(coef(fit_final)[switch_feat, ], 3))
#>      4      8     14     20 
#> -0.085  0.048  0.041  0.084

# What the truth is
cat(switch_feat, "— True coefficients (sign-aligned):\n")
#> mut_11 — True coefficients (sign-aligned):
print(round(-gt[["coefficients"]][switch_feat, ], 3))
#>      4      8     14     20 
#> -1.200  0.186  0.760  1.200

Cox reports a single number for this feature. orthoMTL reveals that its effect changes direction over time — information that could guide clinical interpretation of early versus late treatment response.

Ground Truth Recovery

Finally, we compare the estimated coefficients against the true data-generating coefficients for signal features.

true_coefs <- gt$coefficients[gt$signal_features, ]
est_coefs <- coef(fit_final)[gt$signal_features, ]

# Sign convention:
#   Simulation uses log-hazard scale (negative = protective, reduces hazard)
#   orthoMTL models P(progression-free) (positive = protective)
# Negate true coefficients to align
true_aligned <- -true_coefs

task_cors <- sapply(seq_len(ncol(true_aligned)), function(k) {
  cor(true_aligned[, k], est_coefs[, k])
})
names(task_cors) <- colnames(true_aligned)

cat("Correlation between true and estimated coefficients per threshold:\n")
#> Correlation between true and estimated coefficients per threshold:
print(round(task_cors, 3))
#>     4     8    14    20 
#> 0.978 0.926 0.680 0.653

The simulation generates coefficients on the log-hazard scale (negative = protective), while orthoMTL models the probability of being progression-free (positive = protective). The true coefficients are sign-flipped below so both heatmaps share the same interpretation: blue = protective, red = risk-increasing.

make_long <- function(mat, source_label) {
  data.frame(
    feature = rep(rownames(mat), ncol(mat)),
    threshold = rep(colnames(mat), each = nrow(mat)),
    weight = as.vector(mat),
    source = source_label,
    stringsAsFactors = FALSE
  )
}

# Normalize each matrix to [-1, 1] by dividing by its own max absolute value
normalize <- function(mat) mat / max(abs(mat), na.rm = TRUE)

combined <- rbind(
  make_long(normalize(true_aligned), "True (sign-aligned)"),
  make_long(normalize(est_coefs),    "Estimated")
)

combined$feature <- factor(combined$feature, levels = rev(rownames(true_aligned)))
combined$threshold <- factor(combined$threshold, levels = colnames(true_aligned))
combined$source <- factor(combined$source, levels = c("True (sign-aligned)", "Estimated"))



# Now both panels use the full color range
ggplot2::ggplot(combined, ggplot2::aes(threshold, feature, fill = weight)) +
  ggplot2::geom_tile() +
  ggplot2::scale_fill_gradient2(low = "red", mid = "white", high = "blue", midpoint = 0,
                       limits = c(-1, 1)) +
  ggplot2::facet_wrap(~ source) +
  ggplot2::labs(x = "Threshold (months)", y = NULL, fill = "Normalized\nCoefficient") +
  ggplot2::theme_minimal()

Conclusion

This vignette demonstrated the orthoMTL workflow for survival analysis.

The key difference from standard Cox modelling is not in overall predictive accuracy but in interpretability: orthoMTL reveals how each feature’s effect evolves across time thresholds. This is particularly relevant when:

Features have time-varying effects (e.g., early vs late resistance)
The clinical question is not just “does this feature matter?” but “when does it matter?”
Effect reversal over time (switch patterns) may be clinically meaningful but invisible to proportional hazards models

For the application of orthoMTL to real clinical data, see the SOLAR-1 analysis in Annals of Oncology (2026).

For details on individual functions, see the package help pages (?orthoMTL, ?cv_orthoMTL, ?bootstrap_orthoMTL, etc.).

Citation

If you use orthoMTL in your work, please cite:

Vervier, K., Mahé, P., d’Aspremont, A., Veyrieras, J.-B., & Vert, J.-P. (2014). On Learning Matrices with Orthogonal Columns or Disjoint Supports. ECML-PKDD 2014. https://hal.science/hal-00985654

For the survival extension:

Annals of Oncology (2026). DOI: 10.1016/j.annonc.2026.04.003