Objective 2 — Cox-LASSO Prognostic Signature
Judy Ngobeni
2026-09-22
1 Overview
This file runs the full Objective 2 pipeline on the current analysis set: 23 patients with usable survival follow-up (13 deaths, 10 censored), out of the 64 patients in the Objective 1 harmonized multi-omics dataset (1,221 proteins, 29 metabolites, 1,250 features total).
Every step below was first run and checked interactively; this file
consolidates those checks into one reproducible script, with a fixed
random seed and the tie-handling method
(cox.ties = "efron", matching survival::coxph)
set explicitly throughout. Numbers here may differ very slightly from
earlier exploratory console output, because the console runs used
glmnet’s default tie method before this was fixed.
If your supervisor’s response on the 41 excluded patients changes the
analysis set, only the “Build the survival dataset” section below needs
updating; every later step re-runs on whatever setB and
X/y it produces.
Caution on knit time: the permutation test (500
shuffles) and the two stability-selection loops (200 resamples each) are
the slow steps, together taking roughly 15–20 minutes on a laptop.
cache = TRUE is set above so that later edits to text or
downstream chunks do not force a full re-run.
2 Load data
harm <- readRDS("harmonized_multiomics_dataset.rds")
clin <- readRDS("clinical_clean.rds")
standardize_id <- function(id) {
id <- str_trim(id)
prefix <- str_extract(id, "^[A-Za-z]+")
number <- str_extract(id, "\\d+")
paste0(toupper(prefix), sprintf("%03d", as.integer(number)))
}3 Build the survival dataset
surv_full <- clin %>%
mutate(id_std = standardize_id(linked_id)) %>%
filter(id_std %in% harm$id_std) %>%
mutate(event = case_when(
dead_alive == 0 ~ 1L, # dead_alive is inverted in the source file: 0 = dead
dead_alive == 1 ~ 0L, # 1 = alive
TRUE ~ NA_integer_
))
setB <- surv_full %>%
mutate(followup_use = case_when(
!is.na(followup) ~ followup,
# P034: Excel #VALUE! error in source file; recovered from test/last-visit dates
id_std == "P034" ~ as.numeric(date_of_last_visit_parsed - date_of_test_parsed),
TRUE ~ NA_real_
)) %>%
filter(!is.na(followup_use), !is.na(event)) %>%
mutate(stage_grp = if_else(pathology == "RPC", "Resectable", "Advanced"))
cat("Analysis set:", nrow(setB), "patients,", sum(setB$event), "events\n")## Analysis set: 23 patients, 13 eventstable(setB$pathology, setB$event, dnn = c("pathology", "event"))## event
## pathology 0 1
## LAPC 2 3
## MPC 0 2
## RPC 8 8# Descriptive follow-up
survfit(Surv(followup_use, 1 - event) ~ 1, data = setB) # reverse KM: median follow-up## Call: survfit(formula = Surv(followup_use, 1 - event) ~ 1, data = setB)
##
## n events median 0.95LCL 0.95UCL
## [1,] 23 10 358 253 NAsurvdiff(Surv(followup_use, event) ~ stage_grp, data = setB)## Call:
## survdiff(formula = Surv(followup_use, event) ~ stage_grp, data = setB)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## stage_grp=Advanced 7 5 5.02 7.62e-05 0.000145
## stage_grp=Resectable 16 8 7.98 4.80e-05 0.000145
##
## Chisq= 0 on 1 degrees of freedom, p= 1The 41 patients excluded from setB were deliberately
left out: of them, 13 have a recorded death with no usable follow-up
start/end, 4 are recorded alive with no last-visit date, and 24 have no
recorded vital status at all. Adding the recorded deaths without being
able to censor the survivors on the same basis would bias the analysis
towards poor survival, so they are excluded here.
4 Assemble the feature matrix
dat <- setB %>%
select(id_std, followup_use, event) %>%
inner_join(harm, by = "id_std")
X <- dat %>% select(-id_std, -followup_use, -event) %>% as.matrix()
rownames(X) <- dat$id_std
y <- cbind(time = dat$followup_use, status = dat$event)
stopifnot(nrow(X) == 23, sum(!is.finite(X)) == 0, sum(apply(X, 2, var) == 0) == 0)
dim(X)## [1] 23 12505 Helper functions
# Stratified k-fold split: keeps events spread evenly across folds
strat_folds <- function(status, k = 5) {
fid <- integer(length(status))
for (s in unique(status)) {
idx <- which(status == s)
fid[idx] <- sample(rep(1:k, length.out = length(idx)))
}
fid
}
# Cross-validated deviance gain of the best model over a no-feature model,
# averaged over nrep independent fold splits
cv_gain <- function(X, y, nrep = 10, k = 5) {
lam <- glmnet(X, y, family = "cox", cox.ties = "efron")$lambda
res <- lapply(1:nrep, function(r) {
fid <- strat_folds(y[, "status"], k)
cv.glmnet(X, y, family = "cox", lambda = lam, foldid = fid, cox.ties = "efron")$cvm
})
L <- min(lengths(res))
m <- rowMeans(sapply(res, function(v) v[1:L]))
m[1] - min(m)
}
# Honest cross-validated C-index: feature selection is repeated inside every
# training fold, so this is not inflated by in-sample selection
cv_cindex <- function(X, y, nrep = 20, k = 5, seed_offset = 500) {
sapply(1:nrep, function(r) {
set.seed(seed_offset + r)
fid <- strat_folds(y[, "status"], k)
lp <- numeric(nrow(X))
for (f in 1:k) {
tr <- fid != f
cvf <- cv.glmnet(X[tr, ], y[tr, ], family = "cox", nfolds = 4, cox.ties = "efron")
lp[!tr] <- as.numeric(predict(cvf, X[!tr, ], s = "lambda.min"))
}
concordance(Surv(y[, "time"], y[, "status"]) ~ lp, reverse = TRUE)$concordance
})
}6 Feature selection: repeated cross-validation
The penalty strength (lambda) is chosen by 5-fold cross-validation, repeated 100 times with independent fold splits, because a single split is too noisy with only 13 events.
lam <- glmnet(X, y, family = "cox", cox.ties = "efron")$lambda
res <- lapply(1:100, function(r) {
set.seed(1000 + r)
fid <- strat_folds(y[, "status"], 5)
cv.glmnet(X, y, family = "cox", lambda = lam, foldid = fid, cox.ties = "efron")$cvm
})
L <- min(lengths(res))
mean_cvm <- rowMeans(sapply(res, function(v) v[1:L]))
lam_best <- lam[which.min(mean_cvm)]
c(null_deviance = mean_cvm[1], best_deviance = min(mean_cvm),
gain = mean_cvm[1] - min(mean_cvm))## null_deviance best_deviance gain
## 3.5958439 3.4186484 0.1771955fit_all <- glmnet(X, y, family = "cox", lambda = lam, cox.ties = "efron")
b <- as.matrix(coef(fit_all, s = lam_best))
b <- b[b[, 1] != 0, , drop = FALSE]
b## 1
## dag1 -0.4147170
## gstk1 0.1806352
## sptan1 0.45635367 Stability selection
The model is refitted 200 times, each on a random 70% of patients (stratified by event status), and we count how often each feature is selected. Features chosen only occasionally are not credible signature members.
set.seed(11)
nsub <- 200
counts <- setNames(numeric(ncol(X)), colnames(X))
ev <- which(y[, "status"] == 1); ce <- which(y[, "status"] == 0)
for (i in 1:nsub) {
idx <- c(sample(ev, round(0.7 * length(ev))), sample(ce, round(0.7 * length(ce))))
f <- glmnet(X[idx, ], y[idx, ], family = "cox", lambda = lam_best, cox.ties = "efron")
counts <- counts + (as.numeric(coef(f)) != 0)
}
freq <- sort(counts / nsub, decreasing = TRUE)
head(freq, 10)## sptan1 dag1 gstk1 man2a1 mcam wasf2 anxa2 diaph1 g3bp1 pcyox1
## 0.705 0.355 0.335 0.240 0.230 0.185 0.170 0.170 0.120 0.110top <- head(freq, 10)
barplot(top, horiz = TRUE, las = 1, xlim = c(0, 1),
col = ifelse(top >= 0.6, "#1F9E93", "#B0B0B0"),
xlab = "Selected in resamples (proportion)",
main = "Stability selection (200 resamples)")
abline(v = 0.6, lty = 2, col = "grey40")
8 Permutation test
The survival outcomes are shuffled among patients 500 times, and the entire cross-validated fitting procedure is repeated on each shuffle. This shows how good a result can look by chance alone, given the same feature matrix and sample size.
set.seed(1)
obs_gain <- cv_gain(X, y, nrep = 10)
obs_gain## [1] 0.1734671set.seed(3)
n_perm <- 500
perm_gain <- replicate(n_perm, cv_gain(X, y[sample(nrow(X)), ], nrep = 10))
perm_p <- (1 + sum(perm_gain >= obs_gain)) / (1 + length(perm_gain))
perm_p## [1] 0.05788423quantile(perm_gain, c(0.5, 0.9, 0.95, 0.99))## 50% 90% 95% 99%
## 0.00000000 0.08419071 0.17896396 0.36435471hist(perm_gain, breaks = 30, col = "#D9CBEA", border = "white",
main = "Permutation test: shuffled deviance gain",
xlab = "Cross-validated deviance gain")
abline(v = obs_gain, col = "#5B3A8C", lwd = 2)
text(obs_gain, par("usr")[4] * 0.9, "observed", pos = 4, col = "#5B3A8C")
9 Honest cross-validated performance
Feature selection is repeated inside every training fold below, so this C-index is not inflated by having chosen features on the same patients it is tested on (unlike a plain in-sample Cox summary).
cidx <- cv_cindex(X, y, nrep = 20)
summary(cidx)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4968 0.5815 0.5918 0.6016 0.6313 0.702510 Sensitivity analysis: excluding P009
P009 is the retained metabolomics outlier (flagged in Objective 1) and is one of the 13 deaths in this set. This checks whether the SPTAN1 finding depends on this one patient.
keep <- rownames(X) != "P009"
X2 <- X[keep, ]; y2 <- y[keep, ]
table(y2[, "status"])##
## 0 1
## 10 12lam2 <- glmnet(X2, y2, family = "cox", cox.ties = "efron")$lambda
res2 <- lapply(1:100, function(r) {
set.seed(1000 + r)
fid <- strat_folds(y2[, "status"], 5)
cv.glmnet(X2, y2, family = "cox", lambda = lam2, foldid = fid, cox.ties = "efron")$cvm
})
L2 <- min(lengths(res2))
m2 <- rowMeans(sapply(res2, function(v) v[1:L2]))
lam2_best <- lam2[which.min(m2)]
c(null = m2[1], best = min(m2), gain = m2[1] - min(m2))## null best gain
## 3.4411011 3.1791917 0.2619093f2 <- glmnet(X2, y2, family = "cox", lambda = lam2, cox.ties = "efron")
b2 <- as.matrix(coef(f2, s = lam2_best))
b2[b2[, 1] != 0, , drop = FALSE]## 1
## dag1 -0.315854368
## gstk1 0.117034207
## mylk 0.007592185
## sptan1 0.619883244set.seed(21)
ev2 <- which(y2[, "status"] == 1); ce2 <- which(y2[, "status"] == 0)
counts2 <- setNames(numeric(ncol(X2)), colnames(X2))
for (i in 1:200) {
idx <- c(sample(ev2, round(0.7 * length(ev2))), sample(ce2, round(0.7 * length(ce2))))
f <- glmnet(X2[idx, ], y2[idx, ], family = "cox", lambda = lam2_best, cox.ties = "efron")
counts2 <- counts2 + (as.numeric(coef(f)) != 0)
}
freq_noP009 <- sort(counts2 / 200, decreasing = TRUE)
head(freq_noP009, 8)## sptan1 anxa2 dag1 g3bp1 gstk1 mcam wasf2 man2a1
## 0.765 0.305 0.275 0.250 0.210 0.210 0.160 0.13511 Do the metabolites alone carry any signal?
The full model above searches 1,221 proteins against only 29 metabolites, and the LASSO penalty tends to favour the larger block. This tests the 29 metabolites on their own, so the proteins cannot crowd them out.
met_names <- setdiff(names(readRDS("harmonized_metabolomics.rds")), "id_std")
Xm <- X[, intersect(colnames(X), met_names)]
dim(Xm)## [1] 23 29set.seed(1)
obs_m <- cv_gain(Xm, y, nrep = 10)
set.seed(2)
perm_m <- replicate(200, cv_gain(Xm, y[sample(nrow(Xm)), ], nrep = 10))
obs_m## [1] 0(1 + sum(perm_m >= obs_m)) / (1 + length(perm_m))## [1] 112 Proportional hazards check
top_feature <- names(freq)[1] # SPTAN1, based on the stability results above
d1 <- data.frame(
time = y[, "time"], status = y[, "status"],
feature = scale(X[, top_feature])[, 1]
)
fit_cox <- coxph(Surv(time, status) ~ feature, data = d1, ties = "efron")
summary(fit_cox)## Call:
## coxph(formula = Surv(time, status) ~ feature, data = d1, ties = "efron")
##
## n= 23, number of events= 13
##
## coef exp(coef) se(coef) z Pr(>|z|)
## feature 1.0880 2.9684 0.3382 3.217 0.0013 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## feature 2.968 0.3369 1.53 5.76
##
## Concordance= 0.81 (se = 0.071 )
## Likelihood ratio test= 12.6 on 1 df, p=4e-04
## Wald test = 10.35 on 1 df, p=0.001
## Score (logrank) test = 13.09 on 1 df, p=3e-04cox.zph(fit_cox)## chisq df p
## feature 0.0446 1 0.83
## GLOBAL 0.0446 1 0.83Caution: the hazard ratio and p-value above are
optimistic, because top_feature was chosen from ~1,250
candidates using the same 23 patients it is now tested on. Report the
permutation p-value and the cross-validated C-index above as the primary
evidence; report this Cox summary only as an inflated, in-sample
figure.
13 Results summary
results <- data.frame(
Measure = c("Analysis set", "Events",
"Features selected (best lambda)",
"Top stable feature", "Stability (all 23)", "Stability (excl. P009)",
"Permutation p-value", "Cross-validated C-index (median)",
"Cross-validated C-index (range)",
"Metabolite-only permutation p-value"),
Value = c(
paste(nrow(X), "patients"),
paste(sum(y[, "status"]), "deaths"),
nrow(b),
top_feature,
sprintf("%.1f%%", freq[[top_feature]] * 100),
sprintf("%.1f%%", freq_noP009[[top_feature]] * 100),
sprintf("%.3f", perm_p),
sprintf("%.2f", median(cidx)),
sprintf("%.2f\u2013%.2f", min(cidx), max(cidx)),
sprintf("%.3f", (1 + sum(perm_m >= obs_m)) / (1 + length(perm_m)))
)
)
knitr::kable(results)| Measure | Value |
|---|---|
| Analysis set | 23 patients |
| Events | 13 deaths |
| Features selected (best lambda) | 3 |
| Top stable feature | sptan1 |
| Stability (all 23) | 70.5% |
| Stability (excl. P009) | 76.5% |
| Permutation p-value | 0.058 |
| Cross-validated C-index (median) | 0.59 |
| Cross-validated C-index (range) | 0.50–0.70 |
| Metabolite-only permutation p-value | 1.000 |
14 Save results
saveRDS(
list(
setB = setB, X = X, y = y,
lam_best = lam_best, coef_best = b,
freq = freq, freq_noP009 = freq_noP009,
obs_gain = obs_gain, perm_gain = perm_gain, perm_p = perm_p,
cidx = cidx,
obs_m = obs_m, perm_m = perm_m,
fit_cox_top_feature = fit_cox
),
"objective2_final_results.rds"
)15 Session info
sessionInfo()## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Africa/Johannesburg
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] survival_3.8-6 glmnet_5.0 Matrix_1.7-5 stringr_1.6.0 dplyr_1.2.1
##
## loaded via a namespace (and not attached):
## [1] jsonlite_2.0.0 compiler_4.6.1 tidyselect_1.2.1 Rcpp_1.1.2
## [5] jquerylib_0.1.4 splines_4.6.1 yaml_2.3.12 fastmap_1.2.0
## [9] lattice_0.22-9 R6_2.6.1 generics_0.1.4 shape_1.4.6.1
## [13] knitr_1.51 iterators_1.0.14 tibble_3.3.1 bslib_0.11.0
## [17] pillar_1.11.1 rlang_1.2.0 cachem_1.1.0 stringi_1.8.7
## [21] xfun_0.59 sass_0.4.10 otel_0.2.0 cli_3.6.6
## [25] withr_3.0.3 magrittr_2.0.5 digest_0.6.39 foreach_1.5.2
## [29] grid_4.6.1 rstudioapi_0.19.0 lifecycle_1.0.5 vctrs_0.7.3
## [33] evaluate_1.0.5 glue_1.8.1 codetools_0.2-20 rmarkdown_2.31
## [37] tools_4.6.1 pkgconfig_2.0.3 htmltools_0.5.9