From e53b8edd69d96dcc455e3bbae9085cd5d81287f4 Mon Sep 17 00:00:00 2001 From: zhiz Date: Thu, 14 Mar 2024 18:11:41 +0100 Subject: [PATCH] fix crossrefs --- _quarto.yml | 7 +++++++ docs/search.json | 23 ++++++++++++++++++++++- docs/survomics.html | 14 +++++++------- survomics.qmd | 8 ++++---- 4 files changed, 40 insertions(+), 12 deletions(-) diff --git a/_quarto.yml b/_quarto.yml index cbec697..d312dea 100644 --- a/_quarto.yml +++ b/_quarto.yml @@ -58,4 +58,11 @@ format: grid: sidebar-width: 350px margin-width: 350px + + # pdf: + # geometry: + # - bottom=17mm + # - left=15mm + # - right=15mm + # fontsize: 16pt \ No newline at end of file diff --git a/docs/search.json b/docs/search.json index b3056f4..cd24e76 100644 --- a/docs/search.json +++ b/docs/search.json @@ -25,7 +25,7 @@ "href": "survomics.html#tcga-omics-data", "title": "Supplemental information for ‘Tutorial on survival modeling with applications to omics data’", "section": "TCGA omics data", - "text": "TCGA omics data\nWe use function GDCquery() to query and use GDCdownload() and GDCprepare() to download TCGA omics data from one cancer type (breast cancer). The argument data.category in function GDCquery() specifies the type of omics data, such as \"Copy Number Variation\", \"DNA Methylation\", \"Transcriptome Profiling\", \"Simple Nucleotide Variation\". Note that the downloaded omics data are accompanied by metadata including survival outcomes, clinical and demographic variables. The accompanied metadata are almost the same as the clinical data downloaded via GDCquery_clinic() in the previous section but here only corresponding to one cancer type.\n\n# download TCGA breast cancer (BRCA) mRNA-Seq data using GDC api method\nquery <- TCGAbiolinks::GDCquery(\n project = \"TCGA-BRCA\",\n data.category = \"Transcriptome Profiling\",\n data.type = \"Gene Expression Quantification\",\n workflow.type = \"STAR - Counts\",\n experimental.strategy = \"RNA-Seq\",\n sample.type = c(\"Primary Tumor\")\n)\nTCGAbiolinks::GDCdownload(query = query, method = \"api\")\ndat <- TCGAbiolinks::GDCprepare(query = query)\n\nSummarizedExperiment::assays(dat)$unstranded[1:5, 1:2]\n\n TCGA-A7-A26E-01B-06R-A277-07 TCGA-A2-A0CU-01A-12R-A034-07\nENSG00000000003.15 691 1429\nENSG00000000005.6 20 73\nENSG00000000419.13 335 1674\nENSG00000000457.14 1292 1018\nENSG00000000460.17 536 450\nIt is recommended to use DESeq2 or TMM normalization method for RNA-seq data before further statistical analysis (Y. Zhao et al. 2021). Here we demonstrate how to use the R/Bioconductor package DESeq2 (Love, Huber, and Anders 2014) to normalize the RNA count data.\n\nmeta <- colData(dat)[, c(\"project_id\", \"submitter_id\", \"age_at_diagnosis\", \"ethnicity\", \"gender\", \"days_to_death\", \"days_to_last_follow_up\", \"vital_status\", \"paper_BRCA_Subtype_PAM50\", \"treatments\")]\nmeta$treatments <- unlist(lapply(meta$treatments, function(xx) {\n any(xx$treatment_or_therapy == \"yes\")\n}))\ndds <- DESeq2::DESeqDataSetFromMatrix(assays(dat)$unstranded, colData = meta, design = ~1)\ndds2 <- DESeq2::estimateSizeFactors(dds)\nRNA_count <- DESeq2::counts(dds2, normalized = TRUE)\nRNA_count[1:5, 1:2]\n\n TCGA-A7-A26E-01B-06R-A277-07 TCGA-A2-A0CU-01A-12R-A034-07\nENSG00000000003.15 1899.76848 1419.51789\nENSG00000000005.6 54.98606 72.51561\nENSG00000000419.13 921.01656 1662.89219\nENSG00000000457.14 3552.09968 1011.24507\nENSG00000000460.17 1473.62649 447.01403\nTo perform survival analysis with both clinical/demographic variables and omics data, in the following code we extract female breast cancer patients with their corresponding survival outcomes, clinical/demographic variables and RNA-seq features.\n\nmeta$time <- apply(meta[, c(\"days_to_death\", \"days_to_last_follow_up\")], 1, max, na.rm = TRUE) / 365.25\nmeta$status <- meta$vital_status\nmeta$age <- meta$age_at_diagnosis / 365.25\nclin <- subset(meta, gender == \"female\" & !duplicated(submitter_id) & time > 0 & !is.na(age))\nclin <- clin[order(clin$submitter_id), ]\nRNA_count <- RNA_count[, rownames(clin)]\n\n\n\n\n\n\n\nNote\n\n\n\n\nBioconductor might provide an outdated version of TCGAbiolinks. Here, we use the GitHub version TCGAbiolinks_2.29.6. If you encounter some issues when using this tutorial, please check your installed TCGAbiolinks version. If necessary, you can re-install the package from its GitHub repository. Otherwise, download the data from and load the dat object with: load(\"TCGA_data.rda\").\nThe package TCGAbiolinks cannot retrieve any proteomics or metabolomics data. It is always useful to look at your data first, in particular the data type and dimensions (i.e. numbers of rows and columns for a data frame or matrix)." + "text": "TCGA omics data\nWe use function GDCquery() to query and use GDCdownload() and GDCprepare() to download TCGA omics data from one cancer type (breast cancer). The argument data.category in function GDCquery() specifies the type of omics data, such as \"Copy Number Variation\", \"DNA Methylation\", \"Transcriptome Profiling\", \"Simple Nucleotide Variation\". Note that the downloaded omics data are accompanied by metadata including survival outcomes, clinical and demographic variables. The accompanied metadata are almost the same as the clinical data downloaded via GDCquery_clinic() in the previous section but here only corresponding to one cancer type.\n\n# download TCGA breast cancer (BRCA) mRNA-Seq data using GDC api method\nquery <- TCGAbiolinks::GDCquery(\n project = \"TCGA-BRCA\",\n data.category = \"Transcriptome Profiling\",\n data.type = \"Gene Expression Quantification\",\n workflow.type = \"STAR - Counts\",\n experimental.strategy = \"RNA-Seq\",\n sample.type = c(\"Primary Tumor\")\n)\nTCGAbiolinks::GDCdownload(query = query, method = \"api\")\ndat <- TCGAbiolinks::GDCprepare(query = query)\n\nSummarizedExperiment::assays(dat)$unstranded[1:5, 1:2]\n\n TCGA-A7-A26E-01B-06R-A277-07 TCGA-A2-A0CU-01A-12R-A034-07\nENSG00000000003.15 691 1429\nENSG00000000005.6 20 73\nENSG00000000419.13 335 1674\nENSG00000000457.14 1292 1018\nENSG00000000460.17 536 450\nIt is recommended to use DESeq2 or TMM normalization method for RNA-seq data before further statistical analysis (Y. Zhao et al. 2021). Here we demonstrate how to use the R/Bioconductor package DESeq2 (Love, Huber, and Anders 2014) to normalize the RNA count data.\n\nmeta <- colData(dat)[, c(\"project_id\", \"submitter_id\", \"age_at_diagnosis\", \"ethnicity\", \"gender\", \"days_to_death\", \"days_to_last_follow_up\", \"vital_status\", \"paper_BRCA_Subtype_PAM50\", \"treatments\")]\nmeta$treatments <- unlist(lapply(meta$treatments, function(xx) {\n any(xx$treatment_or_therapy == \"yes\")\n}))\ndds <- DESeq2::DESeqDataSetFromMatrix(assays(dat)$unstranded, colData = meta, design = ~1)\ndds2 <- DESeq2::estimateSizeFactors(dds)\nRNA_count <- DESeq2::counts(dds2, normalized = TRUE)\nRNA_count[1:5, 1:2]\n\n TCGA-A7-A26E-01B-06R-A277-07 TCGA-A2-A0CU-01A-12R-A034-07\nENSG00000000003.15 1899.76848 1419.51789\nENSG00000000005.6 54.98606 72.51561\nENSG00000000419.13 921.01656 1662.89219\nENSG00000000457.14 3552.09968 1011.24507\nENSG00000000460.17 1473.62649 447.01403\nTo perform survival analysis with both clinical/demographic variables and omics data, in the following code we extract female breast cancer patients with their corresponding survival outcomes, clinical/demographic variables and RNA-seq features.\n\nmeta$time <- apply(meta[, c(\"days_to_death\", \"days_to_last_follow_up\")], 1, max, na.rm = TRUE) / 365.25\nmeta$status <- meta$vital_status\nmeta$age <- meta$age_at_diagnosis / 365.25\nclin <- subset(meta, gender == \"female\" & !duplicated(submitter_id) & time > 0 & !is.na(age))\nclin <- clin[order(clin$submitter_id), ]\nRNA_count <- RNA_count[, rownames(clin)]\n\n\n\n\n\n\n\nNote\n\n\n\n\nBioconductor might provide an outdated version of TCGAbiolinks. Here, we use the GitHub version TCGAbiolinks_2.29.6. If you encounter some issues when using this tutorial, please check your installed TCGAbiolinks version. If necessary, you can re-install the package from its GitHub repository. Otherwise, download the data from here and load the dat object with: load(\"TCGA_data.rda\").\nThe package TCGAbiolinks cannot retrieve any proteomics or metabolomics data. It is always useful to look at your data first, in particular the data type and dimensions (i.e. numbers of rows and columns for a data frame or matrix)." }, { "objectID": "survomics.html#survival-analysis-with-low-dimensional-input-data", @@ -47,5 +47,26 @@ "title": "Supplemental information for ‘Tutorial on survival modeling with applications to omics data’", "section": "Survival model validation", "text": "Survival model validation\nThe ideal evaluation of a prognostic model is based on completely independent validation data, since high-dimensional survival models built on the training data can be overfitted. If there are no independent validation data, it is recommended to use resampling-based methods for estimating the uncertainty of the model’s prediction performance. This can be done for example by repeatedly splitting the dataset to training/validation sets and evaluating a model’s performance on the different validation sets using various evaluation metrics.\n\n\n\n\n\n\nModel validation\n\n\n\nTo validate a prediction model systematically, the predictive performance of the model is commonly addressed by\n\nDiscrimination: the ability of the model to distinguish between low and high risk patients\nCalibration: the agreement between the observed and predicted survival probabilities\nOverall performance: the distance between the observed and predicted survival probabilities\n\n\n\nThe performance metrics can be time-dependent or time-independent, with the time-dependent metrics being more informative in general compared to integrated measures (i.e. evaluated across many time points). For survival data, we can assess the discriminatory power of a model (i.e. how well does it ranks patients) or how well a model is calibrated (i.e. how closely the predicted survival probabilities agree numerically with the actual survival outcomes). For example, measures such as the receiver operating characteristic (ROC) curve, the (integrated) area under time-specific ROC curves (AUC, Heagerty and Zheng (2005)) and the concordance index (C-index, Harrell et al. (1982)) are measures of discrimination, while the right-censored logarithmic loss (RCLL, Avati et al. (2020)) and the well-known Brier score (Graf et al. 1999) are used to evaluate both discrimination and calibration performance.\n\nModel evaluation (classic)\n\n\n\n\n\n\nNote\n\n\n\n‘Classic’ here refers to the use of manual R code in combination with many separate R packages which have been routinely used in academia the latest 10+ years for evaluating survival models.\n\n\nTo evaluate the performance of a statistical model, we first split the data into training and validation data sets. For example, we can randomly split the 1047 BRCA patients from TCGA into \\(80\\%\\) as training set and \\(20\\%\\) as validation set.\n\nset.seed(123)\nn <- nrow(x)\nidx <- sample(1:n, n * 0.8, replace = FALSE)\nx_train <- x[idx, ]\ny_train <- y[idx, ]\nx_validate <- x[-idx, ]\ny_validate <- y[-idx, ]\n\n\n\n\n\n\n\nImportant\n\n\n\nThe \\(20\\%\\) split of a dataset is often not considered an independent dataset and resampling-based methods should be used in such cases to provide an unbiased estimate of the predictive accuracy of a prognostic model.\n\n\n\nDiscrimination metrics\n Goodness-of-fit \nThe simplest way to demonstrate the prognostic power of a survival model is to dichotomize the prognostic scores (i.e., linear predictor \\(lp\\) in the Cox model) by median value, and then to use a log-rank test to compare the survival curves of the patients in the two groups. We use the built model to predict the prognostic scores based on the \\(20\\%\\) validation data. The following code shows the goodness-of-fit of a Lasso Cox model with the BRCA patients survival and PAM50 mRNA-Seq data from TCGA.\n\n# train a Lasso Cox model, similarly for other Cox-type models\nset.seed(123)\ncvfit <- cv.glmnet(x_train, y_train, family = \"cox\", nfolds = 5, penalty.factor = pf)\npred_lp <- predict(cvfit, newx = x_validate, s = cvfit$lambda.min, type = \"link\")\n\n# dichotomize by prognostic scores (linear predictor) by median to divide the validation patients into two groups\ngroup_dichotomize <- as.numeric(pred_lp > median(pred_lp))\n\n# draw two survival curves based on KM estimation and compare them by a log-rank test\ndat_tmp <- data.frame(time = y_validate[, 1], status = y_validate[, 2], group = group_dichotomize)\nsfit <- survfit(Surv(time, status) ~ group, data = dat_tmp)\n\nggsurv <- ggsurvplot(sfit,\n conf.int = TRUE, risk.table = TRUE,\n xlab = \"Time since diagnosis (year)\", legend = c(.2, .3),\n legend.labs = c(\"Low risk\", \"High risk\"), legend.title = \"Dichotomized groups\",\n risk.table.y.text.col = TRUE, risk.table.y.text = FALSE\n)\nggsurv$plot <- ggsurv$plot +\n annotate(\"text\", x = 2.6, y = .03, label = paste0(\"Log-rank test:\\n\", surv_pvalue(sfit)$pval.txt))\nggsurv$table <- ggsurv$table + labs(y = \"Dichotomized\\n groups\")\nggsurv\n\n\n\n\nKaplan-Meier curves of the BRCA patients data dichotomized by the median of prognostic scores (calculated from the Lasso Cox model with patients’ survival and mRNA-Seq data) into two groups. The log-rank test is to compare the two survival distributions corresponding to the two groups of patients.\n\n\n\nThe prognostic scores can also be divided into three or more groups based on quantiles and the log-rank test can be used to compare the difference of multiple survival curves.\n\ngroup <- pred_lp\ngroup[pred_lp >= quantile(pred_lp, 2 / 3)] <- 3\ngroup[pred_lp >= quantile(pred_lp, 1 / 3) & pred_lp < quantile(pred_lp, 2 / 3)] <- 2\ngroup[pred_lp < quantile(pred_lp, 1 / 3)] <- 1\n\n# draw two survival curves based on KM estimation and compare them by a log-rank test\ndat_tmp <- data.frame(time = y_validate[, 1], status = y_validate[, 2], group = group)\nsfit <- survfit(Surv(time, status) ~ group, data = dat_tmp)\n\nggsurv <- ggsurvplot(sfit,\n conf.int = TRUE, risk.table = TRUE,\n xlab = \"Time since diagnosis (year)\", legend = c(.2, .3),\n legend.labs = c(\"Low risk\", \"Middle risk\", \"High risk\"), legend.title = \"Groups\",\n risk.table.y.text.col = TRUE, risk.table.y.text = FALSE\n)\nggsurv$plot <- ggsurv$plot +\n annotate(\"text\", x = 3.5, y = .05, label = paste0(\"Log-rank test:\\n\", surv_pvalue(sfit)$pval.txt))\nggsurv\n\n\n\n\nKaplan-Meier curves of the BRCA patients data divided by 33% and 67% quantiles of prognostic scores (calculated from the Lasso Cox model with patients’ survival and mRNA-Seq data) into three groups. The log-rank test is to compare the two survival distributions corresponding to the three groups of patients.\n\n\n\n ROC curve \nThe R package risksetROC (Heagerty and Zheng 2005) can estimate a ROC curve at an evaluation time point. The following code draws a ROC curve at 5-years survival evaluation time point for the 20% TCGA validation data and based on a Lasso Cox model learned from the 80% training data.\n\nROC <- risksetROC(\n Stime = y_validate[, 1], status = y_validate[, 2],\n marker = pred_lp, predict.time = 5, method = \"Cox\",\n main = \"ROC Curve\", col = \"seagreen3\", type = \"s\",\n lwd = 2, xlab = \"1 - Specificity\", ylab = \"Sensitivity\"\n)\ntext(0.7, 0.2, paste(\"AUC =\", round(ROC$AUC, 3)))\n\n\n\n\nROC curve estimated at 5-years survival evaluation time point for the 20% TCGA validation data and based on a Lasso Cox model learned from the 80% training data. The AUC value is the area under the ROC curve. The diagonal line represents the performance of a random prediction of the outcome event with AUC = 0.5.\n\n\n\n Time-dependent AUC \nBoth time-dependent and integrated AUCs can be estimated by the R package risksetROC. We demonstrate the calculation based on both training and validation data.\n\n\n\n\n\n\nSurvival prediction in Cox models\n\n\n\nA Cox proportional hazards model (and Lasso Cox as a consequence) is a semi-parametric model, which means that it does not produce survival distribution predictions by default. However, using the function risksetROC::CoxWeights() you can transform the cv.glmnet’s output linear predictors (lp) to survival distribution predictions. This transformation internally uses the Breslow estimator for the cumulative baseline hazard.\n\n\n\n# unique event times for patients in the training and validation data sets\nutimes_train <- sort(unique(y_train[y_train[, 2] == 1, 1]))\nutimes_validate <- sort(unique(y_validate[y_validate[, 2] == 1, 1]))\n\n# markers from the estimated linear predictors of a Lasso Cox model\npred_lp_train <- predict(cvfit, newx = x_train, s = cvfit$lambda.min, type = \"link\")\npred_lp_validate <- predict(cvfit, newx = x_validate, s = cvfit$lambda.min, type = \"link\")\n\n## compute time-dependent AUC\nAUC_train <- rep(NA, length(utimes_train))\nAUC_validate <- rep(NA, length(utimes_validate))\nfor (j in 1:length(utimes_train)) {\n out <- risksetROC::CoxWeights(\n marker = pred_lp_train, Stime = y_train[, 1],\n status = y_train[, 2], predict.time = utimes_train[j]\n )\n AUC_train[j] <- out$AUC\n}\nfor (j in 1:length(utimes_validate)) {\n out <- risksetROC::CoxWeights(\n marker = pred_lp_validate, Stime = y_validate[, 1],\n status = y_validate[, 2], predict.time = utimes_validate[j]\n )\n AUC_validate[j] <- out$AUC\n}\n\n# draw the time-dependent AUC from the training and validation data sets\ndat_AUC <- data.frame(\n tAUC = c(AUC_train, AUC_validate),\n times = c(utimes_train, utimes_validate),\n group = c(rep(\"AUC_train\", length(AUC_train)), rep(\"AUC_validate\", length(AUC_validate)))\n)\nggplot(dat_AUC, aes(times, tAUC, group = group, color = group)) +\n xlab(\"Evaluation time points (year)\") +\n ylab(\"AUC\") +\n ylim(0.5, 1) +\n geom_step(direction = \"vh\") +\n theme(legend.position = c(0.7, 0.8), legend.title = element_blank())\n\n\n\n\nTime-dependent AUC based on a Lasso Cox model applied to the BRCA patients data from TCGA. The red line shows the Time-dependent AUC calculated from the 80% training data, and the green line shows the Time-dependent AUC calculated from the 20% validation data.\n\n\n\n Integrated AUC \nThe R package risksetROC (Heagerty and Zheng 2005) provides function IntegrateAUC() to estimate integrated AUC.\n\n# Before computing integrated AUC, first estimate survival probabilities at unique survival times\nsurv_prob_train <- unique(survfit(Surv(y_train[, 1], y_train[, 2]) ~ 1)$surv)\nsurv_prob_validate <- unique(survfit(Surv(y_validate[, 1], y_validate[, 2]) ~ 1)$surv)\n\n## integrated AUC (e.g. over tmax=10 years) to get concordance measure based on training data\n(iAUC_train <- risksetROC::IntegrateAUC(AUC_train, utimes_train, surv_prob_train, tmax = 10))\n\n[1] 0.6279646\n\n## integrated AUC (e.g. over tmax=10 years) to get concordance measure based on validation data\n(iAUC_validate <- risksetROC::IntegrateAUC(AUC_validate, utimes_validate, surv_prob_validate, tmax = 10))\n\n[1] 0.6318253\n Time-dependent C-index \nThe C-index is not proper for \\(t\\)-year predictions, see Blanche, Kattan, and Gerds (2019). Consider using time-dependent AUC or time-dependent Brier score instead. For a time-dependent discrimination index for survival data, see Antolini et al. (2005).\n C-index \nThe R package glmnet provides the function glmnet::Cindex() to estimate Harrell’s C-index from a “coxnet” object. The R package survAUC provides the function survAUC::UnoC() to estimated Uno’s C-index. See an example calculation for both C-indexes using a Lasso Cox model below.\n\nset.seed(123)\ncvfit <- cv.glmnet(x_train, y_train, family = \"cox\", nfolds = 5, penalty.factor = pf)\npred <- predict(cvfit, newx = x_validate, type = \"link\", s = cvfit$lambda.min)\n# Harrell's C-index\n(Cindex_Harrell <- Cindex(pred = pred[, 1], y = y_validate))\n\n[1] 0.7246466\n\n# Uno's C-index\n(Cindex_Uno <- survAUC::UnoC(y_train, y_validate, pred))\n\n[1] 0.5772041\n\n\n\nCalibration metrics\nSee a calibration plot in the following section Graphical computation.\n\n\n\nOverall metrics\n Time-dependent Brier score \nThe R package riskRegression can assess the prediction error curves of survival models based on the time-dependent Brier score. Similar to the time-dependent AUC, one needs to first calculate the linear predictors (\\(lp\\)) from a frequentist or Bayesian Cox model, and then use survival::coxph() to regress the survival outcomes on the linear predictor, which is prepared as input of riskRegression::Score() to estimate the (time-dependent) Brier score. If the survival model was fitted via the R package psbcSpeedUp, the curve of the time-dependent Brier score can be drawn by using the function psbcSpeedUp::plotBrier().\n\n## time-dependent Brier score\n\n# use the (x_train, y_train) 80% samples for training\n# and the (x_validate, y_validate) 20% samples for testing\n\ny_train_surv <- Surv(y_train[, \"time\"], y_train[, \"status\"])\ny_validate_surv <- Surv(y_validate[, \"time\"], y_validate[, \"status\"])\nset.seed(123)\ncvfit <- cv.glmnet(x_train, y_train_surv, family = \"cox\", nfolds = 5, penalty.factor = pf)\nlp_train <- predict(cvfit, newx = x_train, s = cvfit$lambda.min, type = \"link\")\nlp_validate <- predict(cvfit, newx = x_validate, s = cvfit$lambda.min, type = \"link\")\n\n# prepare data format suited for function Score() from the riskRegression package\ndata_train <- data.frame(time = y_train[, \"time\"], status = y_train[, \"status\"], lp = as.vector(lp_train))\ndata_validate <- data.frame(time = y_validate[, \"time\"], status = y_validate[, \"status\"], lp = as.vector(lp_validate))\nlasso_train <- coxph(Surv(time, status) ~ lp, data = data_train, y = TRUE, x = TRUE)\nlasso_validate <- coxph(Surv(time, status) ~ lp, data = data_validate, y = TRUE, x = TRUE)\n\n# calculate Brier scores based on both training and validation data\nBrier_train <- riskRegression::Score(list(\"Brier_train\" = lasso_train), formula = Surv(time, status) ~ 1, data = data_train, conf.int = FALSE, metrics = \"brier\", summary = \"ibs\", times = sort(unique(data_train$time)))$Brier$score\nBrier_validate <- riskRegression::Score(list(\"Brier_validate\" = lasso_validate), formula = Surv(time, status) ~ 1, data = data_validate, conf.int = FALSE, metrics = \"brier\", summary = \"ibs\", times = sort(unique(data_validate$time)))$Brier$score\nBrier_score <- rbind(Brier_train, Brier_validate)\nBrier_score <- Brier_score[Brier_score$model != \"Null model\", ]\n\nggplot(Brier_score, aes(times, Brier, group = model, color = model)) +\n xlab(\"Evaluation time points (year)\") +\n ylab(\"Brier score\") +\n geom_step(direction = \"vh\") +\n theme(legend.position = c(0.15, 0.88), legend.title = element_blank())\n\n\n\n\nTime-dependent Brier score based on a Lasso Cox model applied to the BRCA patients data from TCGA. The red line shows the Time-dependent Brier score calculated from the 80% training data, and the green line shows the Time-dependent Brier score calculated from the 20% validation data.\n\n\n\n Integrated Brier score (IBS) \nThe function riskRegression::Score() also summarizes IBS when specifying argument summary = \"ibs\". We can extract the IBS corresponding to the largest evaluation time point.\n\nBrier_validate_ibs <- Brier_validate[Brier_validate$model == \"Brier_validate\", ]\nBrier_validate_ibs$IBS[which.max(Brier_validate_ibs$times)]\n\n[1] 0.1721158\n\n\n\nUncertainty Quantification\n\n\n\n\n\n\nImportant\n\n\n\nIt is recommended to use resampling-based methods for estimating the uncertainty of the model’s performance, if there are no independent validation data for model evaluation. This can be done for example by repeatedly splitting the dataset to training/validation sets and evaluating a model’s performance on the different validation sets using various discrimination or calibration metrics.\n\n\nWe demonstrate how to randomly split the data, e.g. \\(100\\) times, train a Lasso Cox model and estimate the integrated AUC based on the validation data in each replication. For other Cox-type models, we can just replace the model fitting part cv.glmnet() (and predict()) in the for loop below. However, most of the Bayesian Cox models introduced previously are computationally time-consuming when randomly splitting the data many times.\n\n# split the data 100 times\nset.seed(123)\nk <- 100\niAUC <- rep(NA, k)\nfor (i in 1:k) {\n idx <- sample(1:n, n * 0.8, replace = FALSE)\n x_train <- x[idx, ]\n y_train <- y[idx, ]\n x_validate <- x[-idx, ]\n y_validate <- y[-idx, ]\n cvfit <- cv.glmnet(x_train, y_train, family = \"cox\", nfolds = 5, penalty.factor = pf)\n pred_lp <- predict(cvfit, newx = x_validate, s = cvfit$lambda.min, type = \"link\")\n utimes <- sort(unique(y_validate[y_validate[, 2] == 1, 1]))\n AUC <- rep(NA, length(utimes))\n for (j in 1:length(utimes)) {\n out <- CoxWeights(marker = pred_lp, Stime = y_validate[, 1], status = y_validate[, 2], predict.time = utimes[j])\n AUC[j] <- out$AUC\n }\n surv_prob <- unique(survfit(Surv(y_validate[, 1], y_validate[, 2]) ~ 1)$surv)\n iAUC[i] <- IntegrateAUC(AUC, utimes, surv_prob, tmax = 10)\n}\ndat_tmp <- data.frame(x = \"Lasso Cox\", y = iAUC)\n\nset.seed(123)\nggplot(dat_tmp, aes(x, y)) +\n geom_boxplot() +\n ylim(0.5, 1) +\n xlab(\"\") +\n ylab(\"Integrated AUC\") +\n geom_jitter(color = \"blue\", size = 0.5, alpha = 0.5)\n\n\n\n\nIntegrated AUC based on randomly split validation data 100 times. The blue dots are the 100 values of integrated AUC.\n\n\n\nSimilar to obtaining uncertainty of the integrated AUC, we can also estimate the uncertainty of the C-index for evaluating the global performance of our model’s discrimination.\n\n# split the data 100 times\nset.seed(123)\nk <- 100\nCindex_all <- data.frame(Harrell = rep(NA, k), Uno = rep(NA, k))\nfor (i in 1:k) {\n idx <- sample(1:n, n * 0.8, replace = FALSE)\n x_train <- x[idx, ]\n y_train <- y[idx, ]\n x_validate <- x[-idx, ]\n y_validate <- y[-idx, ]\n cvfit <- cv.glmnet(x_train, y_train, family = \"cox\", nfolds = 5, penalty.factor = pf)\n pred <- predict(cvfit, newx = x_validate, type = \"response\", s = cvfit$lambda.min)\n Cindex_all$Harrell[i] <- mean(apply(pred, 2, Cindex, y = y_validate))\n Cindex_all$Uno[i] <- UnoC(y_train, y_validate, pred)\n}\ndat_tmp <- data.frame(x = rep(c(\"Harrell\", \"Uno\"), each = k), y = unlist(Cindex_all))\n\nset.seed(123)\nggplot(dat_tmp, aes(x, y, col = x)) +\n geom_boxplot() +\n geom_jitter(size = 0.5, alpha = 0.5) +\n ylim(0, 1) +\n xlab(\"\") +\n ylab(\"C-index\") +\n theme(legend.position = \"none\")\n\n\n\n\nC-index (Harrell’s and Uno’s) based on randomly split validation data 100 times.\n\n\n\nThe R package c060 (Sill et al. 2014) includes wrapper functions for the glmnet algorithm and implements resampling-based methods (e.g. cross-validation and bootstrap - with and without replacement) based on the peperr package to calculate the time-dependent Brier score. c060 extends peperr package to allow mandatory features without penalization. Binder and Schumacher (2008) recommends to draw bootstrap samples without replacement (i.e. subsampling), because bootstrap samples with replacement often result in too complex models in high-dimensional settings. To use resampling by CV properly for survival data, see R. M. Simon et al. (2011). Note that resampling-based methods here are similar to splitting \\(80\\%/20\\%\\) the data many times which allows us to quantify the uncertainty of the time-dependent Brier score.\n\n## time-dependent Brier score by subsampling from the whole data\nset.seed(123)\npeperr_object <- peperr::peperr(\n response = y_surv, x = x, fit.fun = fit.glmnet,\n args.fit = list(family = \"cox\", penalty.factor = pf),\n complexity = complexity.glmnet,\n args.complexity = list(family = \"cox\", nfolds = 5, penalty.factor = pf),\n indices = resample.indices(n = n, method = \"sub632\", sample.n = 100)\n)\nc060::Plot.peperr.curves(peperr_object)\n\n\n\n\nResampling-based prediction error curves (time-dependent Brier score) a the Lasso Cox model applied to the BRCA data set from TCGA. The gray area indicates the pointwise 2.5% and 97.5% quantiles of the 100 out-of-bag bootstrap samples. The other lines show the prediction error curves of the null model (estimated by the Kaplan-Meier estimator without covariate information), the full apparent error estimates (i.e., the errors as estimated when applying the model to the entire training data set), and the .632+ bootstrap error estimates.\n\n\n\n\n\nFeature stability analysis\nTo identify stable omics features, a straightforward way is to find the overlapped omics features with nonzero coefficients among different data subsets (e.g. CV folds or resamples). The following code summarizes the Lasso Cox selected omics features which were identified at least \\(2\\) or \\(5\\) out of \\(10\\) resamples. Similarly, this approach can be applied to other Lasso-type or Bayesian Cox models that perform feature selection for identifying stable selected features.\n\n# specify the number of resamples k\nk <- 10\nbeta_all <- matrix(nrow = ncol(x), ncol = k)\nset.seed(123)\nfor (j in 1:k) {\n resample_id <- sample(1:nrow(y), nrow(y), replace = TRUE)\n resample_x <- x[resample_id, ]\n resample_y <- y[resample_id, ]\n cvfit <- cv.glmnet(resample_x, resample_y, family = \"cox\", nfolds = 5, penalty.factor = pf)\n beta_all[, j] <- as.vector(coef(cvfit, s = cvfit$lambda.min))\n}\n\n(stable_features <- colnames(x)[rowSums(beta_all != 0) >= 2])\n\n [1] \"age\" \"ethnicity\" \"ANLN\" \"UBE2T\" \"NDC80\" \"PGR\" \"ORC6\" \n [8] \"ESR1\" \"PHGDH\" \"MMP11\" \"SFRP1\" \"CCNE1\" \"BLVRA\" \"BAG1\" \n[15] \"MLPH\" \"CENPF\" \"KRT17\" \"FOXA1\" \"ACTR3B\" \"CCNB1\" \"MDM2\" \n[22] \"MYC\" \"CEP55\" \"SLC39A6\" \"GRB7\" \"NUF2\" \"EGFR\" \"MKI67\" \n[29] \"TMEM45B\" \"FGFR4\" \"MELK\" \"NAT1\" \"CXXC5\" \"BCL2\" \"GPR160\" \n[36] \"TYMS\" \"KRT5\" \"MAPT\" \"MIA\"\n\n(stable_features <- colnames(x)[rowSums(beta_all != 0) >= 5])\n\n [1] \"age\" \"ethnicity\" \"ANLN\" \"ORC6\" \"MMP11\" \"BLVRA\" \"BAG1\" \n [8] \"CCNB1\" \"EGFR\" \"TMEM45B\" \"BCL2\" \"TYMS\" \"KRT5\" \"MIA\"\nAlternatively for a Bayesian Cox model, its median probability model (MPM) can be obtained based on the coefficient estimates over MCMC iterations. The following code shows how to obtain the MPM’s coefficients of the penalized semiparametric Bayesian Cox model with Elastic Net prior run previously.\n\ngammas <- colMeans(matrix(as.numeric(EN_beta_p != 0), ncol = ncol(EN_beta_p)))\nbeta_MPM <- (gammas >= 0.5) * colMeans(EN_beta_p) / gammas\nbeta_MPM[is.na(beta_MPM)] <- 0\nbeta_MPM\n\n age ethnicity ANLN FOXC1 CDH3 UBE2T \n 1.305162e-02 5.348458e-03 -1.299443e-03 -1.857811e-02 -6.123574e-03 -5.467111e-03 \n NDC80 PGR BIRC5 ORC6 ESR1 PHGDH \n-6.652927e-03 -2.101243e-06 -1.640386e-02 -1.237153e-02 -1.077863e-02 2.483990e-02 \n CDC6 MMP11 MYBL2 SFRP1 CCNE1 BLVRA \n-9.079708e-03 -1.588726e-02 5.225344e-03 -1.383981e-02 -3.181265e-03 -2.632373e-02 \n BAG1 MLPH CDC20 CENPF KRT17 FOXA1 \n-3.913529e-02 -1.435805e-02 -2.027232e-02 -2.476495e-02 -2.871143e-02 -3.017213e-03 \n ACTR3B CCNB1 MDM2 MYC CEP55 SLC39A6 \n-2.504869e-03 -1.346817e-03 -2.156041e-02 1.431062e-02 1.421036e-02 -1.150196e-02 \n ERBB2 GRB7 KIF2C NUF2 EGFR MKI67 \n-6.347367e-03 -1.008689e-02 6.033792e-03 -2.405689e-03 -1.964927e-02 1.956661e-02 \n TMEM45B FGFR4 PTTG1 MELK NAT1 CXXC5 \n 2.736216e-02 1.842323e-03 -5.651905e-03 2.894074e-02 -2.126163e-02 2.571472e-02 \n BCL2 RRM2 GPR160 EXO1 UBE2C TYMS \n-5.140894e-03 2.881004e-02 -3.927705e-02 -1.710419e-02 -1.343832e-02 -1.884342e-02 \n KRT5 KRT14 MAPT MIA \n-2.180294e-02 -1.386489e-03 -2.587557e-02 -1.033317e-02\n\n\n\nGraphical representation\nAfter identifying stable omics features predictive of survival outcomes, we can draw a nomogram to allows the graphical calculation of survival probabilities and report a calibration plot for practitioners.\n Nomogram \nWe demonstrate a nomogram using the stable selected features from TCGA breast cancer data preprocessed previously. The R package regplot draws an enhanced regression nomogram based on the rms package.\n\n# remove patients without reporting ethnicity\nyy <- y[x[, 2] != 3, ]\nxx <- x[x[, 2] != 3, ]\n# specify the number of resamples k\nk <- 10\nbeta_all <- matrix(nrow = ncol(xx), ncol = k)\nset.seed(123)\nfor (j in 1:k) {\n resample_id <- sample(1:nrow(yy), nrow(yy), replace = TRUE)\n resample_x <- xx[resample_id, ]\n resample_y <- yy[resample_id, ]\n cvfit <- cv.glmnet(resample_x, resample_y, family = \"cox\", nfolds = 5, penalty.factor = pf)\n beta_all[, j] <- as.vector(coef(cvfit, s = cvfit$lambda.min))\n}\n\n# identify features at least 80% frequently selected\nx_stable <- data.frame(xx[, rowSums(beta_all != 0) >= k * 0.8])\nx_stable$ethnicity <- factor(x_stable$ethnicity)\nlevels(x_stable$ethnicity) <- c(\"Hispanic/latino\", \"Not hispanic/latino\")\n\ndata_tmp <- data.frame(times = yy[, \"time\"], status = yy[, \"status\"], x_stable)\nf <- cph(\n formula = Surv(times, status) ~ age + ethnicity + ANLN + BLVRA + EGFR,\n data = data_tmp, x = TRUE, y = TRUE, surv = TRUE\n)\nddist <- datadist(data_tmp)\noldoption <- options(datadist = \"ddist\")\nsurv <- Survival(f)\nnom <- nomogram(f,\n fun = list(function(x) surv(1, x), function(x) surv(3, x), function(x) surv(5, x)),\n funlabel = c(\n \"1-Year Survival Probability\",\n \"3-Year Survival Probability\",\n \"5-Year Survival Probability\"\n ),\n lp = FALSE\n)\nregplot::regplot(f,\n observation = data_tmp[1, ], failtime = c(1, 3, 5), title = \"\",\n prfail = FALSE, points = TRUE, showP = FALSE, subticks = TRUE\n)\n\n\n\n\nNomogram developed to estimate the overall survival probability for TCGA’s BRAC patients based on demographic and Lasso Cox selected mRNA features. The red coloured symbols represent one patient’s information and predicted probabilities of 1-year, 3-year and 5-year survival.\n\n\n\n Calibration plot \nA calibration plot is a straightforward visualization to show the prediction ability of the nomogram, i.e., the agreement between predicted survival probabilities from the final model and the KM estimated survival probabilities in different percentiles of the predicted values at a time point of interest. We demonstrate below calibration plots based on training and validation data sets, respectively.\n\n# Calibration at 5-year time-point\n\n# prepare suitable data format for calibration plot\nset.seed(123)\ntrain_id <- sample(1:nrow(yy), nrow(yy) * 0.8, replace = FALSE)\ndata_train <- data_tmp[train_id, ]\ndata_validate <- data_tmp[-train_id, ]\n\nddist <- datadist(data_train)\noptions(datadist = \"ddist\")\nf_train <- cph(\n formula = Surv(times, status) ~ age + ethnicity + ANLN + BLVRA + EGFR,\n data = data_train, x = TRUE, y = TRUE, surv = TRUE, time.inc = 5\n)\nf_validate <- update(f_train, data = data_validate)\ncal_train <- calibrate(f_train, u = 5, cmethod = \"KM\", m = nrow(data_train) / 4, B = 200)\ncal_validate <- calibrate(f_validate, u = 5, cmethod = \"KM\", m = nrow(data_validate) / 4, B = 200)\n\nlayout(matrix(1:2, nrow = 1))\nplot(cal_train,\n lwd = 2, lty = 1, errbar.col = \"seagreen3\",\n xlab = \"Predicted survival probability\", ylab = \"Actual survival probability\",\n xlim = c(0, 1), ylim = c(0, 1), col = \"seagreen3\", subtitles = FALSE\n)\ntitle(main = \"Calibration on training data\")\n\nplot(cal_validate,\n lwd = 2, lty = 1, errbar.col = \"seagreen3\",\n xlab = \"Predicted survival probability\", ylab = \"Actual survival probability\",\n xlim = c(0, 1), ylim = c(0, 1), col = \"seagreen3\", subtitles = FALSE\n)\ntitle(main = \"Calibration on validation data\")\n\n\n\n\nNomogram model calibration curves for TCGA’s BRAC patients at 5-year evaluation time-point.\n\n\n\n\n\n\nModel evaluation (mlr3)\n\n\n\n\n\n\nWe will use the mlr3 ML framework and the mlr3proba R library to:\n\n\n\n\nCreate a survival task from a dataset and split it to training and test (validation) sets\nDefine a Lasso Cox model that can output both linear predictors and survival distribution predictions and train/tune it on the training set\nMake predictions using the trained Lasso Cox model on the separate test set\nMeasure the performance of our model (discrimination and calibration) using several evaluation metrics\nUsing resampling techniques, we will assess our model’s capacity for generalization (prediction on unseen data) and the stability of the model’s selected features\n\n\n\nFor the rest of the analysis, we will borrow the terminology from the mlr3 ecosystem of machine learning packages (e.g. task is a dataset, learner is a model, etc.). See mlr3 book for more details.\nFirst, we load the necessary mlr3 libraries (Lang et al. 2019; Sonabend et al. 2021) and some other useful ones:\n\nlibrary(\"mlr3verse\") # mlr3, mlr3pipeplines, mlr3learners, mlr3tuning, paradox, etc.\nlibrary(\"mlr3proba\") # probabilistic learning and survival analysis\nlibrary(\"mlr3extralearners\") # for lrn('surv.glmnet')\n\n\n\nWorkflow\nWe construct an mlr3 survival task (TCGA BRCA dataset essentially, with normalized PAM50 gene expression features and two clinical/demographic variables) and split it into training and test sets (\\(80\\%/20\\%\\)):\n\n# From 'Penalized Cox models' section:\n# x => gene expression matrix (50 PAM50 genes) + 2 clinical variables\n# y => (time, status) target matrix\n\ndata <- cbind.data.frame(x, y)\n# data = readRDS(file = 'data.rds')\ntask <- mlr3proba::as_task_surv(\n x = data,\n time = \"time\", event = \"status\", id = \"BRCA-TCGA\"\n)\ntask # see useful info about the dataset (#features, #samples, target variables)\n\n# split to train and test sets\nset.seed(42)\nsplit <- mlr3::partition(task, ratio = 0.8)\n# split$train # train indices\n# split$test # test indices\n\n<TaskSurv:BRCA-TCGA> (1047 x 54)\n* Target: time, status\n* Properties: -\n* Features (52):\n - dbl (52): ACTR3B, ANLN, BAG1, BCL2, BIRC5, BLVRA, CCNB1, CCNE1, CDC20, CDC6, CDH3,\n CENPF, CEP55, CXXC5, EGFR, ERBB2, ESR1, EXO1, FGFR4, FOXA1, FOXC1, GPR160, GRB7,\n KIF2C, KRT14, KRT17, KRT5, MAPT, MDM2, MELK, MIA, MKI67, MLPH, MMP11, MYBL2, MYC,\n NAT1, NDC80, NUF2, ORC6, PGR, PHGDH, PTTG1, RRM2, SFRP1, SLC39A6, TMEM45B, TYMS,\n UBE2C, UBE2T, age, ethnicity\nWe create a Lasso Cox mlr3 learner (a wrapper around the glmnet::cv.glmnet() function with the capacity to provide survival predictions), where we specify the two clinical variables to be mandatory (i.e. no penalization) and the \\(s\\) value (\\(\\lambda\\) parameter used for prediction) equal to lambda.min:\n\n# tail(task$feature_names) # age, ethnicity are the 2 last features\npf <- c(rep(1, length(task$feature_names) - 2), rep(0, 2))\n\n# define model\ncoxlasso <- lrn(\"surv.cv_glmnet\",\n alpha = 1, nfolds = 5, s = \"lambda.min\",\n penalty.factor = pf\n)\n# coxlasso # see learner details\n# coxlasso$help() # see learner documentation\n\n\n\n\n\n\n\nNote\n\n\n\nA Cox proportional hazards model (and Lasso Cox as a consequence) is a semi-parametric model, which means that it does not produce survival distribution predictions by default. The mlr3 glmnet survival learner uses internally the function glmnet::survfit.coxnet() to transform the cv.glmnet’s linear predictors (lp) to survival distribution predictions. This transformation uses the Breslow estimator for the cumulative baseline hazard (see stype parameter in survival::survfit.coxph()).\n\n\n\nTrain the Lasso Cox model:\n\nset.seed(3)\ncoxlasso$train(task, row_ids = split$train)\n# view `cv.glmnet` fit\ncoxlasso$model$model\n\nCall: (if (cv) glmnet::cv.glmnet else glmnet::glmnet)(x = data, y = target, nfolds = 5L, alpha = 1, penalty.factor = c(1, 1,..., 1, 0, 0), family = \"cox\")\n\nMeasure: Partial Likelihood Deviance \n\n Lambda Index Measure SE Nonzero\nmin 0.01082 14 12.31 0.2743 15\n1se 0.03626 1 12.35 0.2564 2\nGet the survival distribution predictions (\\(distr\\)) along with the linear predictors (\\(lp\\)):\n\npred <- coxlasso$predict(task, row_ids = split$test)\nhead(as.data.table(pred))\n\n row_ids time status crank lp distr\n1: 5 0.9527721 FALSE -2.346574 -2.346574 <list[1]>\n2: 6 4.0438056 FALSE -2.806708 -2.806708 <list[1]>\n3: 15 1.7385352 FALSE -1.845042 -1.845042 <list[1]>\n4: 45 4.5804244 FALSE -1.715041 -1.715041 <list[1]>\n5: 50 5.1279945 FALSE -2.790122 -2.790122 <list[1]>\n6: 54 6.6858316 FALSE -2.466360 -2.466360 <list[1]>\nSo for every patient in the test set, the Lasso Cox model prediction is a linear predictor of the form \\(lp = \\hat{\\beta} X_{new}\\). \\(crank\\) stands for continuous ranking score and it’s the same as \\(lp\\) for the Lasso Cox model. The \\(distr\\) predictions are the per-patient survival distribution predictions, implemented by the R package distr6 which the mlr3proba imports. See respective documentation on the different prediction types supported.\nAn example of using the distr predictions would be to request for the survival probability at e.g. \\(1,5,10,20\\) years for the first two patients in the test set:\n\ntimes <- c(1, 5, 10, 20)\nids <- c(1,2)\npred$distr[ids]$survival(times)\n\n# same logic for the cumulative hazard\n# pred$distr[ids]$cumHazard(times)\n\n [,1] [,2]\n1 0.9889822 0.9930314\n5 0.8711262 0.9165987\n10 0.6708848 0.7772850\n20 0.4540815 0.6075516\n\n\n\nDiscrimination metrics\nWe want to test our Lasso Cox model and see how well it was able to discriminate the patients in the test set. For this we can use the \\(lp\\) predictions of Lasso Cox model and metrics such as the (time-dependent) C-index and (time-dependent) AUC. \n Harrell’s C-index (Harrell et al. 1982): \n\nharrell_c <- msr(\"surv.cindex\")\nharrell_c$id <- \"surv.cindex.harrell\"\n\n# harrell_c # get some details about the measure\n# harrell_c$minimize # FALSE => higher C-index is better\n# harrell_c$range # [0, 1] => [min, max]\n# harrell_c$predict_type # uses the $crank$ predictions (equal to $lp$ for Lasso Cox\n\npred$score(harrell_c)\n\nsurv.cindex.harrell \n 0.6224306\n\n Uno’s C-index (Uno et al. 2011): (across all time points of the test set): \n\nuno_c <- msr(\"surv.cindex\", weight_meth = \"G2\")\nuno_c$id <- \"surv.cindex.uno\"\n\n# Uno's C needs the train data\npred$score(uno_c, task = task, train_set = split$train)\n\nsurv.cindex.uno \n 0.5932426\n\n Uno’s Integrated AUC (Uno et al. 2007) (across all time points of the test set): \n\nuno_iauc <- msr(\"surv.uno_auc\")\nuno_iauc$id <- \"surv.uno_iauc\"\n# uno_iauc$param_set$values$integrated # integrated = TRUE by default\n# sort(unique(pred$truth[,1])) # time points used\n\n# uno_iauc$properties # needs the train data\npred$score(uno_iauc, task = task, train_set = split$train)\n\nsurv.uno_iauc \n 0.6585791\n\n Uno’s AUC at a specific time point, e.g. \\(10\\) years: \n\nuno_auc <- msr(\"surv.uno_auc\", integrated = FALSE, times = 10)\nuno_auc$id <- \"surv.uno_auc.10\"\n\n# needs the train data\npred$score(uno_auc, task = task, train_set = split$train)\n\nsurv.uno_auc.10 \n 0.667014\n\n\n\nCalibration metrics\nWe want to test how well our Lasso Cox model was calibrated. Andres et al. (2018) and Haider et al. (2020) suggested the distributional (D)-calibration accounting survival probabilities across all times. This can be useful when assessing the entire post-treatment survival prognosis, for example, assessing the post liver transplantation survival utility in Andres et al. (2018).\n D-calibration \n\ndcal <- msr(\"surv.dcalib\")\npred$score(dcal)\n\nsurv.dcalib\n 0.8342449\n\n\n\nOverall metrics\nUsually we derive an estimation of the error between the survival distributions (\\(distr\\) predictions) of the patients in the test set and their actual survival outcomes (corresponding to the survival task’s time and status variables). The most frequently used metric is the Brier Score (Graf et al. 1999): \n Integrated Brier Score (IBS) (across all time points of the test set): \n\nibrier <- msr(\"surv.brier\", proper = TRUE)\n# ibrier$help() # see documentation\n# ibrier$predict_type # uses the `distr` predictions\n\n# better to use the train data for the Kaplan-Meier estimation of the censoring distribution, but can use the test set as well\npred$score(ibrier, task = task, train_set = split$train)\n\nsurv.graf \n0.1157316\nWe can also get the standard error of IBS (the above result is the mean across all the test set’s patients) as follows:\n\nibrier_se <- msr(\"surv.brier\", proper = TRUE, se = TRUE)\npred$score(ibrier_se, task = task, train_set = split$train)\n\n surv.graf \n0.01795861\n\n Brier Score at a specific time point, e.g. \\(10\\) years: \n\nbrier10 <- msr(\"surv.brier\", proper = TRUE, integrated = FALSE, times = 10)\nbrier10$id <- \"surv.graf.10\"\n\n# better to use the train data for the Kaplan-Meier estimation of the censoring distribution, but can use the test set as well\npred$score(brier10, task = task, train_set = split$train)\n\nsurv.graf.10 \n 0.1755639\n\n Right-censored Logarithmic Loss score (RCLL) (Avati et al. 2020; Sonabend 2022): \n\nrcll <- msr(\"surv.rcll\")\npred$score(rcll)\n\nsurv.rcll \n 4.781513\n\n\n\n\n\n\n\nNote\n\n\n\nView all evaluation metrics for survival data implemented in mlr3proba here\n\n\n\n\n\nUncertainty Quantification\nSimilar procedure as followed in a previous section.\nWe will perform a stratified split of the BRCA-TCGA survival task to training and test sets (with a \\(80\\%/20\\%\\) ratio as before). Stratification on the censoring indicator status is important because we want our training and test sets to have the same censoring distribution as the initial dataset. Thus we can avoid measuring performance on test sets with severely different censoring distributions that might influence the performance scores.\nStratify survival task by status:\n\ncoxlasso$reset() # un-train model\n\ntask$col_roles$stratum <- \"status\"\n# task\n\nNext, we define the type of resampling (?mlr_resamplings_subsampling), train the Lasso Cox model on all training sets (\\(100\\)) and store the fitted models for feature selection and evaluation:\n\n# 100 times train/test split (80% for training, 20% for validation)\nsubsampling <- rsmp(\"subsampling\", repeats = 100, ratio = 0.8)\n\nset.seed(42)\nrr <- mlr3::resample(\n task = task, learner = coxlasso,\n resampling = subsampling, store_models = TRUE, store_backends = TRUE\n)\n\nWe can use all the aforementioned evaluation metrics to measure the performance of the Lasso Cox models on the \\(100\\) different test sets. Note that if a metric needs the training dataset it is automatically provided by the ResampleResult object (rr):\n\nmeasures <- list(harrell_c, uno_c, uno_iauc, uno_auc, ibrier, brier10, rcll, dcal)\n\nres <- rr$score(measures = measures)\nhead(res)\n\n task_id learner_id resampling_id iteration surv.cindex.harrell surv.cindex.uno \n1: BRCA-TCGA Lasso Cox subsampling 1 0.5679167 0.6090304\n2: BRCA-TCGA Lasso Cox subsampling 2 0.5422131 0.4884603\n3: BRCA-TCGA Lasso Cox subsampling 3 0.7604049 0.5740556\n4: BRCA-TCGA Lasso Cox subsampling 4 0.6610169 0.5277736\n5: BRCA-TCGA Lasso Cox subsampling 5 0.5800073 0.5655076\n6: BRCA-TCGA Lasso Cox subsampling 6 0.5427837 0.6975740\n surv.uno_iauc surv.uno_auc.10 surv.graf surv.graf.10 surv.rcll surv.dcalib\n1: 0.6628350 0.4719335 0.1736450 0.2582078 4.914347 0.5890987\n2: 0.4023684 0.5652588 0.1855191 0.2185453 4.914264 0.6552037\n3: 0.5941948 0.5235439 0.1014616 0.1329289 5.041827 0.8480314\n4: 0.5360690 0.5110032 0.1122070 0.1157937 5.045144 0.4981035\n5: 0.6160743 0.5388393 0.1099091 0.1940644 5.060352 0.2036418\n6: 0.6494779 0.6400328 0.1377626 0.2445218 5.059521 0.5175357\nHidden columns: task, learner, resampling, prediction\nWe extract and visualize the discrimination and calibration (resampled) performance of our Lasso Cox model using several evaluation metrics:\n\nset.seed(42)\n\n# C-indexes, AUCs (integrated and at t = 10 years)\nres[, .(surv.cindex.harrell, surv.cindex.uno, surv.uno_iauc, surv.uno_auc.10)] %>%\n tidyr::pivot_longer(\n cols = tidyselect::everything(),\n names_to = \"Measure\", values_to = \"Value\"\n ) %>%\n mutate(Measure = case_when(\n Measure == \"surv.cindex.harrell\" ~ \"Harrell's C-index\",\n Measure == \"surv.cindex.uno\" ~ \"Uno's C-index\",\n Measure == \"surv.uno_iauc\" ~ \"Uno's Integrated AUC\",\n Measure == \"surv.uno_auc.10\" ~ \"Uno's AUC (t = 10 years)\",\n )) %>%\n mutate(Measure = factor(Measure, levels = c(\n \"Harrell's C-index\",\n \"Uno's C-index\",\n \"Uno's Integrated AUC\",\n \"Uno's AUC (t = 10 years)\"\n ))) %>%\n ggplot(aes(x = Measure, y = Value, fill = Measure)) +\n geom_boxplot() +\n ylim(c(0.2, 0.8)) +\n geom_hline(yintercept = 0.5, color = \"red\", linetype = \"dashed\") +\n theme_bw(base_size = 14) +\n labs(title = \"Discrimination Measures\") +\n theme(axis.text.x = element_blank())\n\n\n\n\nDiscrimination performance of Lasso Cox on the TCGA-BRCA dataset (expression data of the PAM50 genes and the variables age and ethnicity). Performance metrics used are Harrell’s C-index, Uno’s C-index, Uno’s Integrated AUC and Uno’s AUC at 10 years. The dataset was split to training/validation sets 100 times to allow for the quantification of uncertainty in the different performance estimates.\n\n\n\n# different scales for each measure, so we separate the plots\nset.seed(42)\n\n# Integrated Brier Score and Brier Score at t = 10 years\nres[, .(surv.graf, surv.graf.10)] %>%\n tidyr::pivot_longer(\n cols = tidyselect::everything(),\n names_to = \"Measure\", values_to = \"Value\"\n ) %>%\n mutate(Measure = case_when(\n Measure == \"surv.graf\" ~ \"IBS\",\n Measure == \"surv.graf.10\" ~ \"BS(t=10)\"\n )) %>%\n ggplot(aes(x = Measure, y = Value, fill = Measure)) +\n geom_boxplot(show.legend = FALSE) +\n geom_jitter(color = \"blue\", size = 0.5, alpha = 0.5, show.legend = FALSE) +\n labs(title = \"Integrated Brier Score vs Brier Score (t = 10 years)\") +\n theme_bw(base_size = 14) +\n theme(axis.title.x = element_blank())\n\n# RCLL\nres[, .(surv.rcll)] %>%\n tidyr::pivot_longer(\n cols = tidyselect::everything(),\n names_to = \"Measure\", values_to = \"Value\"\n ) %>%\n mutate(Measure = case_when(\n Measure == \"surv.rcll\" ~ \"RCLL\"\n )) %>%\n ggplot(aes(x = Measure, y = Value)) +\n geom_boxplot(show.legend = FALSE) +\n geom_jitter(color = \"blue\", size = 0.5, alpha = 0.5, show.legend = FALSE) +\n labs(title = \"Right-censored Log Loss\") +\n theme_bw(base_size = 14) +\n theme(axis.title.x = element_blank())\n\n\n Calibration performance of Lasso Cox on the TCGA-BRCA dataset (expression data of the PAM50 genes and the variables age and ethnicity). Performance metrics used are the Integrated Brier Score (IBS), the Brier Score at 10 years and the Right-Censored Logarithmic Loss (RCLL). The dataset was split to training/validation sets 100 times to allow for the quantification of uncertainty in the different performance estimates.\n\n\n\nFeature stability analysis\nWe can extract the selected features from all \\(100\\) trained Lasso Cox models and create a frequency selection table:\n\n# get selected features from all trained learners in a list\nsf_list <- lapply(rr$learners, function(learner) {\n learner$selected_features()\n})\n\n# make frequency selection table\nn <- length(sf_list)\nfs_res <- sort(table(unlist(sf_list)), decreasing = TRUE)\ntimes <- as.vector(unname(fs_res))\ntibble::tibble(feat_name = names(fs_res), times = times, freq = times / n)\n\n# A tibble: 35 × 3\n feat_name times freq\n <chr> <int> <dbl>\n1 age 100 1\n2 ethnicity 100 1\n3 ANLN 43 0.43\n4 BLVRA 41 0.41\n5 BAG1 37 0.37\n6 MIA 34 0.34\n7 TYMS 30 0.3\n8 KRT5 27 0.27\n9 MMP11 27 0.27\n10 BCL2 26 0.26\n# ℹ 25 more rows\nAs age and ethnicity were not penalized, they have non-zero coefficients in all Lasso Cox models and therefore are included in all selected feature sets.\nLastly, we can use the R package stabm (Bommert and Lang 2021) to assess how similar the \\(100\\) selected feature sets were. We will demonstrate the use of three metrics which measure the stability of the Lasso Cox’s feature selection on the TCGA-BRCA dataset:\n\nThe Jaccard index\nNogueira’s metric (corrected for chance, i.e. independent of the number of features; Nogueira, Sechidis, and Brown (2018))\nZucknick’s metric (extension of Jaccard index that considers the correlation between the features; Zucknick, Richardson, and Stronach (2008)):\n\n\nset.seed(42)\njac <- stabm::stabilityJaccard(features = sf_list, correction.for.chance = \"none\")\nnog <- stabm::stabilityNogueira(features = sf_list, p = length(task$feature_names))\n\n# Similarity of each pair of features using Pearson correlation\nsim.mat <- abs(stats::cor(x = task$data(cols = task$feature_names), method = \"p\"))\nzuck <- stabm::stabilityZucknick(\n features = sf_list, sim.mat = sim.mat,\n threshold = 0.9, correction.for.chance = \"estimate\", N = 100\n)\n\ntibble::tibble(jaccard = jac, nogueira = nog, zucknick = zuck)\n\n# A tibble: 1 × 3\n jaccard nogueira zucknick\n <dbl> <dbl> <dbl>\n1 0.474 0.412 0.442\nFrom the above values we conclude that the stability of Lasso Cox’s feature selection is neither poor nor excellent but somewhere in between." + }, + { + "objectID": "survomics.html#multi-omics-integrative-modeling", + "href": "survomics.html#multi-omics-integrative-modeling", + "title": "Supplemental information for ‘Tutorial on survival modeling with applications to omics data’", + "section": "Multi-omics integrative modeling", + "text": "Multi-omics integrative modeling\nIntegration of multi-omics data from various types of omics profiling technologies can improve our understanding of complex disease mechanisms holistically, and hence improve prediction of disease progression and survival (Hasin, Seldin, and Lusis 2017; Subramanian et al. 2020; Vandereyken et al. 2023). One obvious benefit of multi-omics data is the availability of biological information from multiple, partly redundant levels, i.e., different modalities, which allows for the assessment of associations within and between different omics datasets.\nA naive strategy for multi-omics data analysis is to treat all the omics features equally and independently, which may lead to a worse performance in comparison to methods that take into account the group structure. A recent large-scale benchmarking study for survival prediction using TCGA multi-omics data demonstrated a slightly improved survival prediction performance when taking into account the group structure of the multi-omics data sets (Herrmann et al. 2021). However, Herrmann et al. (2021) did not evaluate the survival model performance in terms of feature selection, i.e., selection of the most prognostic omics variables or features. To capture associations within modalities of omics features, penalized regressions with a \\(\\ell_2\\)-norm penalty or Bayesian methods with a group lasso prior can become useful. Focusing on feature selection or parsimonious effects is often beneficial for the purpose of clinical implementation, where one can use the sparse Group-Lasso Cox (N. Simon et al. 2013) and Bayesian Cox with elastic priors (Lee, Chakraborty, and Sun 2015).\nTo capture associations between modalities of omics features, including their overlapping or nested and hierarchical relationships, the use of biological network structures can become useful, as for example in the Bayesian Cox model with spike-and-slab and MRF priors (Madjar et al. 2021). There is an increasing knowledge of biological interconnections across various molecular profiles, and systems biology approaches are being developed that try to capture these deep insights when modeling or predicting disease mechanisms (Yan et al. 2018; Karimi et al. 2022). There is a number of review and benchmarking studies of the methodologies for multi-omics integration, although these reviews do not focus specifically on survival applications; for example, Herrmann et al. (2021) for penalized regression-based methods, boosting-based methods and random forest-based methods; Agamah et al. (2022) for network analyses, Ickstadt, Schäfer, and Zucknick (2018) and Chu et al. (2022) for Bayesian approaches; and Kang, Ko, and Mersha (2022) for deep learning methods.\nHowever, to achieve a comprehensive and biologically meaningful integration of high-dimensional multi-omics data, there is a need for continued development of computational and statistical approaches that consider both technical and biological intricacies of the data and technologies, respectively (Wissel, Rowson, and Boeva 2023). This is currently a very active research field, and we expect to see many improved multi-omics methods for survival prediction in the future." + }, + { + "objectID": "survomics.html#r-session-info", + "href": "survomics.html#r-session-info", + "title": "Supplemental information for ‘Tutorial on survival modeling with applications to omics data’", + "section": "R session info", + "text": "R session info\n\nsessionInfo()\n\nR version 4.3.1 (2023-06-16)\nPlatform: x86_64-apple-darwin20 (64-bit)\nRunning under: macOS Monterey 12.7\n\nMatrix products: default\nBLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib \nLAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0\n\nlocale:\n[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8\n\ntime zone: Europe/Oslo\ntzcode source: internal\n\nattached base packages:\n[1] stats4 stats graphics grDevices utils datasets methods base \n\nother attached packages:\n [1] stabm_1.2.2 mlr3extralearners_0.7.1-9000 mlr3proba_0.5.3 \n [4] mlr3verse_0.2.8 mlr3_0.16.1 regplot_1.1 \n [7] survAUC_1.2-0 rms_6.7-0 Hmisc_5.1-0 \n[10] c060_0.3-0 peperr_1.5 snowfall_1.84-6.2 \n[13] snow_0.4-4 riskRegression_2023.03.22 risksetROC_1.0.4.1 \n[16] MASS_7.3-60 BhGLM_1.1.0 GGally_2.1.2 \n[19] psbcGroup_1.5 mvtnorm_1.2-2 SuppDists_1.1-9.7 \n[22] LearnBayes_2.15.1 SGL_1.3 grpreg_3.4.0 \n[25] plotmo_3.6.2 TeachingDemos_2.12 plotrix_3.8-2 \n[28] Formula_1.2-5 glmnet_4.1-7 Matrix_1.5-4.1 \n[31] M3C_1.22.0 survminer_0.4.9 ggpubr_0.6.0 \n[34] survival_3.5-5 ggplot2_3.4.2 dplyr_1.1.2 \n[37] DESeq2_1.40.2 SummarizedExperiment_1.30.2 Biobase_2.60.0 \n[40] GenomicRanges_1.52.0 GenomeInfoDb_1.36.1 IRanges_2.34.1 \n[43] S4Vectors_0.38.1 BiocGenerics_0.46.0 MatrixGenerics_1.12.2 \n[46] matrixStats_1.0.0 TCGAbiolinks_2.29.6 \n\nloaded via a namespace (and not attached):\n [1] tgp_2.4-21 progress_1.2.2 mlr3hyperband_0.4.5 \n [4] penalized_0.9-52 nnet_7.3-19 Biostrings_2.68.1 \n [7] TH.data_1.1-2 vctrs_0.6.3 digest_0.6.32 \n [10] png_0.1-8 corpcor_1.6.10 shape_1.4.6 \n [13] proxy_0.4-27 parallelly_1.36.0 reshape_0.8.9 \n [16] foreach_1.5.2 withr_2.5.0 param6_0.2.4 \n [19] xfun_0.39 memoise_2.0.1 diptest_0.76-0 \n [22] MatrixModels_0.5-1 zoo_1.8-12 DEoptimR_1.1-1 \n [25] distr6_1.8.0 prettyunits_1.1.1 prabclus_2.3-2 \n [28] KEGGREST_1.40.0 httr_1.4.6 downloader_0.4 \n [31] maptree_1.4-8 rstatix_0.7.2 globals_0.16.2 \n [34] fpc_2.2-10 rstudioapi_0.14 generics_0.1.3 \n [37] base64enc_0.1-3 curl_5.0.1 zlibbioc_1.46.0 \n [40] doSNOW_1.0.20 GenomeInfoDbData_1.2.10 lgr_0.4.4 \n [43] xtable_1.8-4 stringr_1.5.0 doParallel_1.0.17 \n [46] evaluate_0.21 S4Arrays_1.0.4 BiocFileCache_2.8.0 \n [49] hms_1.1.3 colorspace_2.1-0 filelock_1.0.2 \n [52] cmprsk_2.2-11 reticulate_1.30 flexmix_2.3-19 \n [55] magrittr_2.0.3 readr_2.1.4 modeltools_0.2-23 \n [58] lattice_0.21-8 palmerpenguins_0.1.1 future.apply_1.11.0 \n [61] robustbase_0.99-0 SparseM_1.81 XML_3.99-0.14 \n [64] class_7.3-22 pillar_1.9.0 nlme_3.1-162 \n [67] iterators_1.0.14 compiler_4.3.1 RSpectra_0.16-1 \n [70] stringi_1.7.12 paradox_0.11.1 minqa_1.2.5 \n [73] dictionar6_0.1.3 plyr_1.8.8 crayon_1.5.2 \n [76] abind_1.4-5 sm_2.2-5.7.1 locfit_1.5-9.8 \n [79] bit_4.0.5 sandwich_3.0-2 mlr3mbo_0.2.1 \n [82] codetools_0.2-19 multcomp_1.4-25 matrixcalc_1.0-6 \n [85] openssl_2.0.6 e1071_1.7-13 splines_4.3.1 \n [88] Rcpp_1.0.11 quantreg_5.95 dbplyr_2.3.2 \n [91] TCGAbiolinksGUI.data_1.20.0 knitr_1.43 blob_1.2.4 \n [94] utf8_1.2.3 clue_0.3-64 lme4_1.1-34 \n [97] listenv_0.9.0 checkmate_2.2.0 ggsignif_0.6.4 \n[100] tibble_3.2.1 mlr3tuningspaces_0.4.0 statmod_1.5.0 \n[103] tzdb_0.4.0 pkgconfig_2.0.3 tools_4.3.1 \n[106] cachem_1.0.8 RSQLite_2.3.1 rvest_1.0.3 \n[109] DBI_1.1.3 numDeriv_2016.8-1.1 mlr3filters_0.7.1 \n[112] fastmap_1.1.1 rmarkdown_2.22 scales_1.2.1 \n[115] mlegp_3.1.9 grid_4.3.1 mets_1.3.2 \n[118] broom_1.0.5 carData_3.0-5 rpart_4.1.19 \n[121] yaml_2.3.7 foreign_0.8-84 cli_3.6.1 \n[124] purrr_1.0.1 lifecycle_1.0.3 askpass_1.1 \n[127] bbotk_0.7.2 lava_1.7.2.1 kernlab_0.9-32 \n[130] backports_1.4.1 mlr3tuning_0.19.0 BiocParallel_1.34.2 \n[133] gtable_0.3.3 umap_0.2.10.0 parallel_4.3.1 \n[136] mlr3cluster_0.1.8 jsonlite_1.8.7 bitops_1.0-7 \n[139] bit64_4.0.5 Rtsne_0.16 mlr3learners_0.5.6 \n[142] polspline_1.1.23 survMisc_0.5.6 spacefillr_0.3.2 \n[145] htmltools_0.5.5 KMsurv_0.1-5 set6_0.2.6 \n[148] rappdirs_0.3.3 mlr3pipelines_0.5.0-1 glue_1.6.2 \n[151] penalizedSVM_1.1.4 mlr3viz_0.6.1 timereg_2.0.5 \n[154] XVector_0.40.0 RCurl_1.98-1.12 mclust_6.0.0 \n[157] gridExtra_2.3 boot_1.3-28.1 R6_2.5.1 \n[160] tidyr_1.3.0 km.ci_0.5-6 ooplah_0.2.0 \n[163] cluster_2.1.4 beanplot_1.3.1 nloptr_2.0.3 \n[166] mlr3misc_0.13.0 vioplot_0.4.0 DelayedArray_0.26.3 \n[169] tidyselect_1.2.0 htmlTable_2.4.1 xml2_1.3.4 \n[172] mlr3fselect_0.11.0 car_3.1-2 AnnotationDbi_1.62.1 \n[175] future_1.33.0 munsell_0.5.0 data.table_1.14.8 \n[178] htmlwidgets_1.6.2 mlr3data_0.7.0 RColorBrewer_1.1-3 \n[181] biomaRt_2.56.1 rlang_1.1.1 uuid_1.1-1 \n[184] fansi_1.0.4 prodlim_2023.03.31 psbcSpeedUp_2.0.4" + }, + { + "objectID": "survomics.html#references", + "href": "survomics.html#references", + "title": "Supplemental information for ‘Tutorial on survival modeling with applications to omics data’", + "section": "References", + "text": "References\n\n\n\n\nAgamah, Francis E., Jumamurat R. Bayjanov, Anna Niehues, Kelechi F. Njoku, Michelle Skelton, Gaston K. Mazandu, Thomas H. A. Ederveen, Nicola Mulder, Emile R. Chimusa, and Peter A. C. ’t Hoen. 2022. “Computational Approaches for Network-Based Integrative Multi-Omics Analysis.” Frontiers in Molecular Biosciences 9 (November): 967205. https://doi.org/10.3389/fmolb.2022.967205.\n\n\nAndres, Axel, Aldo Montano-Loza, Russell Greiner, Max Uhlich, Ping Jin, Bret Hoehn, David Bigam, James Andrew Mark Shapiro, and Norman Mark Kneteman. 2018. “A Novel Learning Algorithm to Predict Individual Survival After Liver Transplantation for Primary Sclerosing Cholangitis.” Edited by Yinglin Xia. PLOS ONE 13 (3): e0193523. https://doi.org/10.1371/journal.pone.0193523.\n\n\nAntolini, Laura, Patrizia Boracchi, Elia Biganzoli, L Antolini, P Boracchi, and E Biganzoli. 2005. “A time-dependent discrimination index for survival data.” Statistics in Medicine 24 (24): 3927–44. https://doi.org/10.1002/SIM.2427.\n\n\nAvati, Anand, Tony Duan, Sharon Zhou, Kenneth Jung, Nigam H Shah, and Andrew Y Ng. 2020. “Countdown Regression: Sharp and Calibrated Survival Predictions.” In Proceedings of the 35th Uncertainty in Artificial Intelligence Conference, edited by Ryan P Adams and Vibhav Gogate, 115:145–55. Proceedings of Machine Learning Research. PMLR. https://proceedings.mlr.press/v115/avati20a.html.\n\n\nBinder, Harald, and Martin Schumacher. 2008. “Adapting Prediction Error Estimates for Biased Complexity Selection in High-Dimensional Bootstrap Samples.” Statistical Applications in Genetics and Molecular Biology 7 (1). https://doi.org/10.2202/1544-6115.1346.\n\n\nBlanche, Paul, Michael W Kattan, and Thomas A Gerds. 2019. “The c-index is not proper for the evaluation of t-year predicted risks.” Biostatistics 20 (2): 347–57. https://doi.org/10.1093/biostatistics/kxy006.\n\n\nBommert, Andrea, and Michel Lang. 2021. “stabm: Stability Measures for Feature Selection.” Journal of Open Source Software 6 (59): 3010. https://doi.org/10.21105/joss.03010.\n\n\nBreheny, Patrick, and Jian Huang. 2015. “Group Descent Algorithms for Nonconvex Penalized Linear and Logistic Regression Models with Grouped Predictors.” Statistics and Computing 25: 173–87.\n\n\nBreheny, Patrick, Yaohui Zeng, and Ryan Kurth. 2021. “Grpreg: Regularization Paths for Regression Models with Grouped Covariates.” R Package Version 3.4.0. https://CRAN.R-project.org/package=grpreg (23 January 2023, date last accessed).\n\n\nChu, Jiadong et al. 2022. “The Application of Bayesian Methods in Cancer Prognosis and Prediction.” Cancer Genom. Proteom. 19 (1): 1–11. https://doi.org/10.21873/cgp.20298.\n\n\nFriedman, Jerome, Trevor Hastie, and Robert Tibshirani. 2010. “Regularization Paths for Generalized Linear Models via Coordinate Descent.” Journal of Statistical Software 33 (1): 1–22. https://doi.org/10.18637/jss.v033.i01.\n\n\nGraf, Erika, Claudia Schmoor, Willi Sauerbrei, and Martin Schumacher. 1999. “Assessment and comparison of prognostic classification schemes for survival data.” Statistics in Medicine 18: 2529–45. https://doi.org/10.1002/(SICI)1097-0258(19990915/30)18:17/18.\n\n\nHaider, Humza, Bret Hoehn, Sarah Davis, and Russell Greiner. 2020. “Effective Ways to Build and Evaluate Individual Survival Distributions.” Journal of Machine Learning Research 21 (85): 1–63. http://jmlr.org/papers/v21/18-772.html.\n\n\nHarrell, Frank E., Robert M. Califf, David B. Pryor, Kerry L. Lee, and Robert A. Rosati. 1982. “Evaluating the Yield of Medical Tests.” JAMA: The Journal of the American Medical Association 247 (18): 2543–46. https://doi.org/10.1001/jama.1982.03320430047030.\n\n\nHasin, Yehudit, Marcus Seldin, and Aldons Lusis. 2017. “Multi-Omics Approaches to Disease.” Genome Biol. 18 (August): 83. https://doi.org/10.1186/s13059-017-1215-1.\n\n\nHeagerty, Patrick J., and Yingye Zheng. 2005. “Survival Model Predictive Accuracy and ROC Curves.” Biometrics 61 (1): 92–105. https://doi.org/10.1111/J.0006-341X.2005.030814.X.\n\n\nHerrmann, Moritz et al. 2021. “Large-Scale Benchmark Study of Survival Prediction Methods Using Multi-Omics Data.” Brief. Bioinform. 22 (3): bbaa167. https://doi.org/10.1093/bib/bbaa167.\n\n\nIckstadt, Katja, Martin Schäfer, and Manuela Zucknick. 2018. “Toward Integrative Bayesian Analysis in Molecular Biology.” Annual Review of Statistics and Its Application 5 (1): 141–67. https://doi.org/10.1146/annurev-statistics-031017-100438.\n\n\nJohn, Christopher R, David Watson, Dominic Russ, Katriona Goldmann, Michael Ehrenstein, Costantino Pitzalis, Myles Lewis, and Michael Barnes. 2020. “M3C: Monte Carlo Reference-Based Consensus Clustering.” Scientific Reports 10 (1): 1816. https://doi.org/10.1038/s41598-020-58766-1.\n\n\nKang, Mingon, Euiseong Ko, and Tesfaye B Mersha. 2022. “A Roadmap for Multi-Omics Data Integration Using Deep Learning.” Briefings in Bioinformatics 23 (1): bbab454. https://doi.org/10.1093/bib/bbab454.\n\n\nKarimi, Mohammad Reza, Amir Hossein Karimi, Shamsozoha Abolmaali, Mehdi Sadeghi, and Ulf Schmitz. 2022. “Prospects and Challenges of Cancer Systems Medicine: From Genes to Disease Networks.” Briefings in Bioinformatics 23 (1): bbab343. https://doi.org/10.1093/bib/bbab343.\n\n\nLang, Michel, Martin Binder, Jakob Richter, Patrick Schratz, Florian Pfisterer, Stefan Coors, Quay Au, Giuseppe Casalicchio, Lars Kotthoff, and Bernd Bischl. 2019. “mlr3: A modern object-oriented machine learning framework in R.” Journal of Open Source Software 4 (44): 1903. https://doi.org/10.21105/JOSS.01903.\n\n\nLee, Kyu Ha, Sounak Chakraborty, and Jianguo Sun. 2015. “Survival Prediction and Variable Selection with Simultaneous Shrinkage and Grouping Priors.” Stat. Anal. Data Min. 8 (2): 114–27. https://doi.org/10.1002/sam.11266.\n\n\n———. 2021. “psbcGroup: Penalized Parametric and Semiparametric Bayesian Survival Models with Shrinkage and Grouping Priors.” R Package Version 1.5. https://CRAN.R-project.org/package=psbcGroup (23 January 2023, date last accessed).\n\n\nLove, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2.” Genome Biology 15: 550. https://doi.org/10.1186/s13059-014-0550-8.\n\n\nMadjar, Katrin et al. 2021. “Combining Heterogeneous Subgroups with Graph-Structured Variable Selection Priors for Cox Regression.” BMC Bioinformatics 22: 586. https://doi.org/10.1186/s12859-021-04483-z.\n\n\nMounir, Marta AND Silva, Mohamed AND Lucchetta. 2019. “New Functionalities in the TCGAbiolinks Package for the Study and Integration of Cancer Data from GDC and GTEx.” PLOS Computational Biology 15 (3): 1–18. https://doi.org/10.1371/journal.pcbi.1006701.\n\n\nNogueira, Sarah, Konstantinos Sechidis, and Gavin Brown. 2018. “On the Stability of Feature Selection Algorithms.” Journal of Machine Learning Research 18 (174): 1–54. http://jmlr.org/papers/v18/17-514.html.\n\n\nParker, Joel S., Michael Mullins, Maggie C. U. Cheang, Samuel Leung, David Voduc, Tammi Vickery, Sherri Davies, et al. 2009. “Supervised Risk Predictor of Breast Cancer Based on Intrinsic Subtypes.” Journal of Clinical Oncology 27 (8): 1160–67. https://doi.org/10.1200/JCO.2008.18.1370.\n\n\nSill, Martin, Thomas Hielscher, Natalia Becker, and Manuela Zucknick. 2014. “c060: Extended Inference with Lasso and Elastic-Net Regularized Cox and Generalized Linear Models.” Journal of Statistical Software 62 (5): 1–22. https://doi.org/10.18637/jss.v062.i05.\n\n\nSimon, Noah et al. 2013. “A Sparse-Group Lasso.” J. Comput. Graph Stat. 22 (2): 231–45. https://doi.org/10.1080/10618600.2012.681250.\n\n\nSimon, Noah, Jerome Friedman, Trevor Hastie, and Rob Tibshirani. 2019. “SGL: Fit a GLM (or Cox Model) with a Combination of Lasso and Group Lasso Regularization.” R Package Version 1.3. https://CRAN.R-project.org/package=SGL (23 January 2023, date last accessed).\n\n\nSimon, Richard M., Jyothi Subramanian, Ming-Chung Li, and Supriya Menezes. 2011. “Using cross-validation to evaluate predictive accuracy of survival risk classifiers based on high-dimensional data.” Briefings in Bioinformatics 12 (3): 203–14. https://doi.org/10.1093/bib/bbr001.\n\n\nSonabend, Raphael. 2022. “Scoring Rules in Survival Analysis.” https://arxiv.org/abs/2212.05260.\n\n\nSonabend, Raphael, Franz J. Király, Andreas Bender, Bernd Bischl, and Michel Lang. 2021. “mlr3proba: an R package for machine learning in survival analysis.” Bioinformatics 37 (17): 2789–91. https://doi.org/10.1093/BIOINFORMATICS/BTAB039.\n\n\nSubramanian, Indhupriya, Srikant Verma, Shiva Kumar, Abhay Jere, and Krishanpal Anamika. 2020. “Multi-Omics Data Integration, Interpretation, and Its Application.” Bioinformatics and Biology Insights 14 (January): 117793221989905. https://doi.org/10.1177/1177932219899051.\n\n\nUno, Hajime, Tianxi Cai, Michael J. Pencina, Ralph B. D’Agostino, and L. J. Wei. 2011. “On the C-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data.” Statistics in Medicine 30 (10): 1105–17. https://doi.org/10.1002/SIM.4154.\n\n\nUno, Hajime, Tianxi Cai, Lu Tian, and L. J. Wei. 2007. “Evaluating Prediction Rules for t-Year Survivors With Censored Regression Models.” Journal of the American Statistical Association 102 (478): 527–37. https://doi.org/10.1198/016214507000000149.\n\n\nVandereyken, Katy, Alejandro Sifrim, Bernard Thienpont, and Thierry Voet. 2023. “Methods and Applications for Single-Cell and Spatial Multi-Omics.” Nature Reviews Genetics 24 (8): 494–515. https://doi.org/10.1038/s41576-023-00580-2.\n\n\nWissel, David, Daniel Rowson, and Valentina Boeva. 2023. “Systematic Comparison of Multi-Omics Survival Models Reveals a Widespread Lack of Noise Resistance.” Cell Reports Methods 3 (4): 100461. https://doi.org/10.1016/j.crmeth.2023.100461.\n\n\nYan, Jingwen, Shannon L Risacher, Li Shen, and Andrew J. Saykin. 2018. “Network Approaches to Systems Biology Analysis of Complex Disease: Integrative Methods for Multi-Omics Data.” Briefings in Bioinformatics 19 (6): 1370–81. https://doi.org/10.1093/bib/bbx066.\n\n\nYi, Nengjun, Zaixiang Tang, Xinyan Zhang, and Boyi Guo. 2019. “BhGLM: Bayesian hierarchical GLMs and survival models, with applications to genomics and epidemiology.” Bioinformatics 35 (8): 1419–21. https://doi.org/10.1093/bioinformatics/bty803.\n\n\nZhao, Yingdong, Ming-Chung Li, Mariam M. Konate, Li Chen, Biswajit Das, Chris Karlovich, P. Mickey Williams, Yvonne A. Evrard, James H. Doroshow, and Lisa M. McShane. 2021. “TPM, FPKM, or Normalized Counts? A Comparative Study of Quantification Measures for the Analysis of RNA-Seq Data from the NCI Patient-Derived Models Repository.” Journal of Translational Medicine 19 (1): 269. https://doi.org/10.1186/s12967-021-02936-w.\n\n\nZhao, Zhi, John Zobolas, Manuela Zucknick, and Tero Aittokallio. 2024. “Tutorial on survival modeling with applications to omics data.” Bioinformatics, March. https://doi.org/10.1093/bioinformatics/btae132.\n\n\nZhao, Zhi, Manuela Zucknick, Maral Saadati, and Axel Benner. 2023. “Penalized Semiparametric Bayesian Survival Models.” R Package Version 2.0.4. https://CRAN.R-project.org/package=psbcSpeedUp.\n\n\nZucknick, Manuela, Sylvia Richardson, and Euan A Stronach. 2008. “Comparing the Characteristics of Gene Expression Profiles Derived by Univariate and Multivariate Classification Methods.” Statistical Applications in Genetics and Molecular Biology 7 (1). https://doi.org/doi:10.2202/1544-6115.1307.\n\n\nZucknick, Manuela, Maral Saadati, and Axel Benner. 2015. “Nonidentical twins: Comparison of frequentist and Bayesian lasso for Cox models.” Biometrical Journal 57 (6): 959–81. https://doi.org/10.1002/bimj.201400160." } ] \ No newline at end of file diff --git a/docs/survomics.html b/docs/survomics.html index 9bb65ce..f0c9655 100644 --- a/docs/survomics.html +++ b/docs/survomics.html @@ -465,7 +465,7 @@

TCGA omics data
    -

  • Bioconductor might provide an outdated version of TCGAbiolinks. Here, we use the GitHub version TCGAbiolinks_2.29.6. If you encounter some issues when using this tutorial, please check your installed TCGAbiolinks version. If necessary, you can re-install the package from its GitHub repository. Otherwise, download the data from DOI and load the dat object with: load("TCGA_data.rda").

    • +

  • Bioconductor might provide an outdated version of TCGAbiolinks. Here, we use the GitHub version TCGAbiolinks_2.29.6. If you encounter some issues when using this tutorial, please check your installed TCGAbiolinks version. If necessary, you can re-install the package from its GitHub repository. Otherwise, download the data from here and load the dat object with: load("TCGA_data.rda").

  • The package TCGAbiolinks cannot retrieve any proteomics or metabolomics data. It is always useful to look at your data first, in particular the data type and dimensions (i.e. numbers of rows and columns for a data frame or matrix).

@@ -2050,15 +2050,15 @@

Fe -
-

Multi-omics integrative modeling

+
+

Multi-omics integrative modeling

Integration of multi-omics data from various types of omics profiling technologies can improve our understanding of complex disease mechanisms holistically, and hence improve prediction of disease progression and survival (Hasin, Seldin, and Lusis 2017; Subramanian et al. 2020; Vandereyken et al. 2023). One obvious benefit of multi-omics data is the availability of biological information from multiple, partly redundant levels, i.e., different modalities, which allows for the assessment of associations within and between different omics datasets.

A naive strategy for multi-omics data analysis is to treat all the omics features equally and independently, which may lead to a worse performance in comparison to methods that take into account the group structure. A recent large-scale benchmarking study for survival prediction using TCGA multi-omics data demonstrated a slightly improved survival prediction performance when taking into account the group structure of the multi-omics data sets (Herrmann et al. 2021). However, Herrmann et al. (2021) did not evaluate the survival model performance in terms of feature selection, i.e., selection of the most prognostic omics variables or features. To capture associations within modalities of omics features, penalized regressions with a \(\ell_2\)-norm penalty or Bayesian methods with a group lasso prior can become useful. Focusing on feature selection or parsimonious effects is often beneficial for the purpose of clinical implementation, where one can use the sparse Group-Lasso Cox (N. Simon et al. 2013) and Bayesian Cox with elastic priors (Lee, Chakraborty, and Sun 2015).

To capture associations between modalities of omics features, including their overlapping or nested and hierarchical relationships, the use of biological network structures can become useful, as for example in the Bayesian Cox model with spike-and-slab and MRF priors (Madjar et al. 2021). There is an increasing knowledge of biological interconnections across various molecular profiles, and systems biology approaches are being developed that try to capture these deep insights when modeling or predicting disease mechanisms (Yan et al. 2018; Karimi et al. 2022). There is a number of review and benchmarking studies of the methodologies for multi-omics integration, although these reviews do not focus specifically on survival applications; for example, Herrmann et al. (2021) for penalized regression-based methods, boosting-based methods and random forest-based methods; Agamah et al. (2022) for network analyses, Ickstadt, Schäfer, and Zucknick (2018) and Chu et al. (2022) for Bayesian approaches; and Kang, Ko, and Mersha (2022) for deep learning methods.

However, to achieve a comprehensive and biologically meaningful integration of high-dimensional multi-omics data, there is a need for continued development of computational and statistical approaches that consider both technical and biological intricacies of the data and technologies, respectively (Wissel, Rowson, and Boeva 2023). This is currently a very active research field, and we expect to see many improved multi-omics methods for survival prediction in the future.

-
-

R session info

+
+

R session info

sessionInfo()
@@ -2161,8 +2161,8 @@

R session info

[181] biomaRt_2.56.1 rlang_1.1.1 uuid_1.1-1 [184] fansi_1.0.4 prodlim_2023.03.31 psbcSpeedUp_2.0.4
-
-

References

+
+

References

diff --git a/survomics.qmd b/survomics.qmd index 2b9cbf0..4d1ea0c 100644 --- a/survomics.qmd +++ b/survomics.qmd @@ -214,7 +214,7 @@ RNA_count <- RNA_count[, rownames(clin)] Here, we use the GitHub version TCGAbiolinks_2.29.6. If you encounter some issues when using this tutorial, please check your installed **TCGAbiolinks** version. If necessary, you can re-install the package from its [GitHub repository](https://github.com/BioinformaticsFMRP/TCGAbiolinks.git). -Otherwise, download the data from [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.10044741.svg)](https://doi.org/10.5281/zenodo.10044741) and load the `dat` object with: `load("TCGA_data.rda")`. +Otherwise, download the data from [here](https://doi.org/10.5281/zenodo.10044741) and load the `dat` object with: `load("TCGA_data.rda")`. - The package **TCGAbiolinks** cannot retrieve any proteomics or metabolomics data. It is always useful to look at your data first, in particular the data type and dimensions (i.e. numbers of rows and columns for a data frame or matrix). @@ -2110,7 +2110,7 @@ tibble::tibble(jaccard = jac, nogueira = nog, zucknick = zuck) From the above values we conclude that the stability of Lasso Cox's feature selection is neither poor nor excellent but somewhere in between. -# Multi-omics integrative modeling {-} +## Multi-omics integrative modeling {-} Integration of multi-omics data from various types of omics profiling technologies can improve our understanding of complex disease mechanisms holistically, and hence improve prediction of disease progression and survival [@Hasin2017;@Subramanian2020;@Vandereyken2023]. One obvious benefit of multi-omics data is the availability of biological information from multiple, partly redundant levels, i.e., different modalities, which allows for the assessment of associations within and between different omics datasets. @@ -2128,7 +2128,7 @@ There is a number of review and benchmarking studies of the methodologies for m However, to achieve a comprehensive and biologically meaningful integration of high-dimensional multi-omics data, there is a need for continued development of computational and statistical approaches that consider both technical and biological intricacies of the data and technologies, respectively [@Wissel2023]. This is currently a very active research field, and we expect to see many improved multi-omics methods for survival prediction in the future. -# R session info {.unlisted .unnumbered} +## R session info {.unlisted .unnumbered} ```{r} sessionInfo() @@ -2238,5 +2238,5 @@ loaded via a namespace (and not attached): knitr::write_bib("knitr", "references.bib") ``` -# References +## References