This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Predicting survival of cancer patients using glmnet cox

I want to use glmnet cox regression approach to predict survival from methylation data for cancer patients. But I couldn't find any descriptive reference from the authors except this one https://cran.r-project.org/web/packages/glmnet/vignettes/Coxnet.pdf, specially discussing its prediction values.

I am looking for answers to some of my basic questions.

  1. can we predict the survival time (number of days for which patient will survive after diagnosis) and vital status of a cancer patient from its gene expression or methylation data using glmnet cox regression? If not: a) What does it actually predict then? If it is hazard ratio, then what is its statistical significance? b) How can I use these predicted values to infer the survival of these patient. So that I can compare these predictions with the actual survival information of the patients to test the accuracy of my model.

  2. I have also used a real patient data and divided it into training (80%) and testing (20%). For training my model, I used the "overall survival" as well as "vital status of the patient" from the training set . Later on, used test data for prediction, where it is giving some negative values. But I have no clue what are these values and what does this negative sign mean, if it has some?

Please accept my apology for my short knowledge in statistics. Any help and suggestions are welcome.

Thanks in advance

r survival glmnet cox cancer
• 11,590 views • link
updated by Kevin Blighe Kevin Blighe • written by Researcher Researcher

can we predict the survival time (number of days for which patient will survive after diagnosis) and vital status of a cancer patient from its gene expression or methylation data using glmnet cox regression?

You can have days / time to death as the outcome variable, in which case the model becomes a linear regression. If you have vital status as the outcome, then it becomes a binary logistic regression model, as people can only either be alive or dead. Researchers do not typically use survival data in this way.

a) What does it actually predict then? If it is hazard ratio, then what is its statistical significance? b) How can I use these predicted values to infer the survival of these patient. So that I can compare these predictions with the actual survival information of the patients to test the accuracy of my model.

To be pedantic about it: models used in this way are not 'predicting' anything. If you set-up your model correctly, then time to death and vital status are taken into account when calculating the hazard ratio. Hazard ratio for gene expression obviously has a different interpretation than a hazard ratio for a binary categorical variable.

Usually the log rank p-value is taken and quoted in manuscripts.

For a better answer on the statistics behind all of this, please post your question at CrossValidated.

I have also used a real patient data and divided it into training (80%) and testing (20%). For training my model, I used the "overall survival" as well as "vital status of the patient" from the training set . Later on, used test data for prediction, where it is giving some negative values. But I have no clue what are these values and what does this negative sign mean, if it has some?

So, you are using TCGA data. In order for me to answer this part, you will have to show the code that you have used.

• 0 views • link
written by Kevin Blighe Kevin Blighe

Hi @Kevin, Thank you for your reply and a nice explanation. Yes I am using TCGA data and need some more clarification for your comments.

Researchers do not typically use survival data in this way.

What do you suggest then, how should I use the survival data of patients to predict their survival?

Usually the log rank p-value is taken and quoted in manuscripts.

In this statement you are referring to which manuscript.

Thanks

• 0 views • link
written by Researcher Researcher

What do you suggest then, how should I use the survival data of patients to predict their survival?

It could just be the way that you are describing it, but it seems like you simply want a model like this:

VitalStatus ~ gene

I was just implying that we typically have something like this:

Surv(time, VitalStatus) ~ gene

I have an entire tutorial, here: Survival analysis with gene expression

----------------------------------------

In this statement you are referring to which manuscript.

It was just a general observation from having worked with many physicians and having read manuscripts. I also just found This, which states: "The standard method, the log-rank test, was used for statistical comparison of survival times."

• 0 views • link
written by Kevin Blighe Kevin Blighe

1 answer

After our comments, I am just posting a fully reproducible example of how anyone can perform a Cox regression analysis with Lasso via glmnet, here using gene expression data. I post this due to the fact that your question title is likely to attract many hits from search engines, and therefore requires a reasonable answer, not a commentary.

You should also check the vignette for Cox models with lasso: https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html#cox

Following my example from here to get some data: Survival analysis via Cox Proportional Hazards regression

Download and set-up data

library(Biobase)
library(GEOquery)

# load series and platform data from GEO
gset <- getGEO('GSE2990', GSEMatrix =TRUE, getGPL=FALSE)
x <- exprs(gset[[1]])

# remove Affymetrix control probes
x <- x[-grep('^AFFX', rownames(x)),]

# transform the expression data to Z scores
x <- t(scale(t(x)))

# extract information of interest from the phenotype data (pdata)
idx <- which(colnames(pData(gset[[1]])) %in%
  c('age:ch1', 'distant rfs:ch1', 'er:ch1',
    'ggi:ch1', 'grade:ch1', 'node:ch1',
    'size:ch1', 'time rfs:ch1'))
metadata <- data.frame(pData(gset[[1]])[,idx],
  row.names = rownames(pData(gset[[1]])))

# remove samples from the pdata that have any NA value
discard <- apply(metadata, 1, function(x) anyis.na(x)))
metadata <- metadata[!discard,]

# filter the Z-scores expression data to match the samples in our pdata
x <- x[,which(colnames(x) %in% rownames(metadata))]

# check that sample names match exactly between pdata and Z-scores 
all((colnames(x) == rownames(metadata)) == TRUE)

# create a merged pdata and Z-scores object
coxdata <- data.frame(metadata, t(x))

# tidy column names
colnames(coxdata)[1:8] <- c('Age', 'Distant.RFS', 'ER',
  'GGI', 'Grade', 'Node',
  'Size', 'Time.RFS')

# prepare certain phenotypes
coxdata$Age <- as.numeric(gsub('^KJ', '', coxdata$Age))
coxdata$Distant.RFS <- as.numeric(coxdata$Distant.RFS)
coxdata$ER <- factor(coxdata$ER, levels = c(0, 1))
coxdata$Grade <- factor(coxdata$Grade, levels = c(1, 2, 3))
coxdata$Time.RFS <- as.numeric(gsub('^KJX|^KJ', '', coxdata$Time.RFS))

Now perform the Cox survival analysis with lasso:

require(glmnet)
require(survival)

coxdata[1:10,1:15]
         Age Distant.RFS ER       GGI Grade Node Size Time.RFS  X1007_s_at
GSM65752  40           0  0   2.48005     3    0  1.2     2280  0.61091071
GSM65753  46           0  1 -0.633592     1    0  1.3     2675  0.66320599
GSM65754  67           1  1  -1.02972     1    0    6      426  0.39016262
GSM65755  41           1  1   1.04395     3    0  3.3      182 -0.01019652
GSM65756  38           1  1  0.683559     3    0  3.2       46  0.71052390
GSM65757  34           0  1   1.05919     2    0  1.6     3952  0.92185361
GSM65758  46           1  1  -1.23306     2    0  2.1     1824  0.01764556
GSM65760  57           1  1  0.679034     3    0  2.2      699  0.23028475
GSM65761  63           1  1   0.31585     2    0  2.8      730  1.19710288
GSM65762  54           1  1  -1.17846     2    0  1.7      882  0.70807685
            X1053_at   X117_at   X121_at  X1255_g_at    X1294_at  X1316_at
GSM65752  0.80321777 0.9576007 0.3236762  0.20730041 -0.57647252 0.5170007
GSM65753  0.12525380 0.4033209 0.1810776  0.08599057  0.43754667 0.5622666
GSM65754 -0.24124193 0.4739672 0.4422982  0.24498162  0.01942667 0.8180497
GSM65755  0.78790284 0.3364908 0.2314423  0.04583503 -0.62393026 0.3874972
GSM65756  0.28652687 0.4308167 0.3995143  0.36486763 -0.71256876 0.9834658
GSM65757  0.35280009 0.5448720 0.2622013 -0.03252980  0.08691080 0.5999085
GSM65758  0.11911061 0.6634477 0.2808262 -0.02828385  0.49225720 0.7275434
GSM65760 -0.07668636 0.3542536 0.4505841  0.33018455  0.63742514 0.6240488
GSM65761  0.31422159 0.3428862 0.3198296  0.17794074  1.01971818 0.7655406
GSM65762 -0.19605026 0.4011981 0.2035890  0.08022851  0.14603593 0.4257209


coxlasso <- glmnet(
  x = data.matrix(coxdata[,10:ncol(coxdata)]),
  y = Surv(coxdata$Time.RFS, coxdata$Distant.RFS),
  family = 'cox')

Kevin

• 0 views • link
written by Kevin Blighe Kevin Blighe

Dear Dr Blighe, After obtaining coxlasso, How should one use predict? (or if I did it correctly what is output?)

b= (coxdata[1,10:ncol(coxdata)])
predict(coxlasso,newx=t(matrix(unlist(b)))) # gives the following

structure(c(0, -0.0573107373872287, -0.111778770449609, -0.162756412874245, 
-0.200548555731299, -0.273323342727612, -0.37012635158656, -0.469871946543072, 
-0.560372678843767, -0.603628666594683, -0.655720887665505, -0.744350261134244, 
-0.745442838491973, -0.714564041893826, -0.717049862434844, -0.70521708276932, 
-0.679319806284637, -0.665634536861567, -0.700412594114427, -0.631192863204038, 
-0.518927033204136, -0.59903960519009, -0.710614096489971, -0.734474495657648, 
-0.689713249193243, -0.655382990137035, -0.651263237099303, -0.73407802999023, 
-0.986691917670182, -1.20827896728222, -1.51349499203079, -1.81573235348754, 
-2.10732385404113, -2.22227393784832, -2.28715605040497, -2.37237736865864, 
-2.39892781588181, -2.36043683952605, -2.30550417286347, -2.20642321062464, 
-2.07954256878492, -1.93909488063093, -1.80323178627091, -1.68509558186175, 
-1.5868039850474, -1.47271332547507, -1.4496698727763, -1.44238444153376, 
-1.42238792669507, -1.40830105312875, -1.42415625831459, -1.49491224976645, 
-1.56582208573533, -1.62935815886695, -1.6888692912116, -1.74952813431818, 
-1.79478410463216, -1.84859609609203, -1.91299001106592, -1.9701781067282, 
-2.00258384919503, -2.02809472579634, -2.01731529287855, -2.01089094538337, 
-2.01409549624958, -2.0841219068267, -2.17383465452615, -2.27144773940283, 
-2.35637473554323, -2.47712835478697, -2.64082729928918, -2.82987980091043, 
-3.04828718845536, -3.19363361314437, -3.3121089786443, -3.44275077136097, 
-3.59026154364923, -3.73271512101116, -3.88856516780995, -4.0598688214655, 
-4.24089644933194, -4.41183710459872, -4.54472144701383, -4.67485738991895, 
-4.81328885214305, -4.95878120189189, -5.12172618371055, -5.28584411091286, 
-5.46974157546711, -5.66872554927844, -5.79469241653852, -5.94584536181941, 
-6.0989671561301, -6.25322278976369, -6.4462496497583, -6.68000499424007, 
-6.91742293950293, -7.14866284337644, -7.40415287117432, -7.66418089026084
), .Dim = c(1L, 100L), .Dimnames = list(NULL, c("s0", "s1", "s2", 
"s3", "s4", "s5", "s6", "s7", "s8", "s9", "s10", "s11", "s12", 
"s13", "s14", "s15", "s16", "s17", "s18", "s19", "s20", "s21", 
"s22", "s23", "s24", "s25", "s26", "s27", "s28", "s29", "s30", 
"s31", "s32", "s33", "s34", "s35", "s36", "s37", "s38", "s39", 
"s40", "s41", "s42", "s43", "s44", "s45", "s46", "s47", "s48", 
"s49", "s50", "s51", "s52", "s53", "s54", "s55", "s56", "s57", 
"s58", "s59", "s60", "s61", "s62", "s63", "s64", "s65", "s66", 
"s67", "s68", "s69", "s70", "s71", "s72", "s73", "s74", "s75", 
"s76", "s77", "s78", "s79", "s80", "s81", "s82", "s83", "s84", 
"s85", "s86", "s87", "s88", "s89", "s90", "s91", "s92", "s93", 
"s94", "s95", "s96", "s97", "s98", "s99")))

Thanks

• 0 views • link
written by english.server english.server

Dear Kevin,

Can you please explain me why you use [,10;ncol(coxdata)] -> why from column 10?

coxlasso <- glmnet( x = data.matrix(coxdata[,10:ncol(coxdata)]), y = Surv(coxdata$Time.RFS, coxdata$Distant.RFS), family = 'cox')

• 0 views • link
written by Nona • Nona

Hi, I think that it relates to the fact that the genes begin only from column 10.

• 0 views • link
written by Kevin Blighe Kevin Blighe

Log in to answer this question.

More posts like this