I am doing several survival analyses on TCGA (The Cancer Genome Atlas) data. As this is my first time doing this kind of analysis, I have a question about it.
To study the influence of the gene expression of several genes on the patients’ survival rate, I have conducted Cox regressions for each gene individually. To do so, I have previously log-transformed, scaled, and centered gene expression values. Then, I wonder if I could assess the prediction power of those genes in a multivariate cox regression model. However, for a particular case, I obtained 30 genes significantly associated with the patients’ survival rate. So, I wonder if I could do some kind of variable selection step with Lasso cox regression using the glmnet R package. After lambda cross-validation, I obtained the genes with non-zero coefficients and now I would like to know the performance of those genes for prediction or a way to assess the quality of the subset of genes obtained. To do this, I selected coefficients reported by glmnet at lambda.min and manually introduced them in coxph function from survival R package through init argument and then prevent any additional fitting with an iter.max=0 argument, just as reported in this previous post Cox regression with lasso regression. However, although each one of the genes was significantly associated with the survival rate on its own, after this step, the p-value for the log-rank test is near 1. Am I doing something wrong? The number of events ranged from 8 to 160 in several analyses. It would be better to do other approximations to assess this problem? If so, could you please provide me with some suggestions?
Edit: Added outputs from glmnet and coxph.
genes = c("HRK.8739", "NKAIN3.286183", "KLK15.55554", "NKX2.4.644524", "ADIG.149685", "DMRTC1.63947", "HTR3B.9177",
"TMPRSS7.344805", "PSKH2.85481", "TH.7054", "PSG3.5671", "ADAM2.2515")
final_clin = final_clin[!final_clin$new_death <= 0,]
b = final_clin[, c(genes)]
x = data.matrix(b)
y = Surv(final_clin$new_death, final_clin$death_event)
coxlasso <- glmnet(x = x, y = y, family = 'cox')
plot(coxlasso, label = T)
cv.lasso_OS <- cv.glmnet(x, y, family="cox", standardize=T, alpha=1, nfolds=10, parallel=T)
cv.lasso_OS
Call: cv.glmnet(x = x, y = y, nfolds = 10, parallel = T, family = "cox", standardize = T, alpha = 1)
Measure: Partial Likelihood Deviance
Lambda Index Measure SE Nonzero
min 0.01255 20 10.93 0.4086 4
1se 0.07351 1 11.33 0.4002 0
(coefs = coef(cv.lasso_OS, s = cv.lasso_OS$lambda.min))
12 x 1 sparse Matrix of class "dgCMatrix"
1
HRK.8739 .
NKAIN3.286183 0.03159785
KLK15.55554 0.27252634
NKX2.4.644524 .
ADIG.149685 .
DMRTC1.63947 -0.19942332
HTR3B.9177 -0.18328332
TMPRSS7.344805 .
PSKH2.85481 0.25614089
TH.7054 .
PSG3.5671 .
ADAM2.2515 -0.09113397
coefs = coefs[coefs != 0]
control = coxph.control(iter.max = 0)
surv_model = coxph(Surv(new_death, death_event) ~ NKAIN3.286183 + KLK15.55554 + DMRTC1.63947
+ HTR3B.9177 + PSKH2.85481 + ADAM2.2515, data = final_clin,
control = control, init = coefs)
summary(surv_model)
Call:
coxph(formula = Surv(new_death, death_event) ~ NKAIN3.286183 +
KLK15.55554 + DMRTC1.63947 + HTR3B.9177 + PSKH2.85481 + ADAM2.2515,
data = final_clin, init = coefs, control = control)
n= 531, number of events= 160
coef exp(coef) se(coef) z Pr(>|z|)
NKAIN3.286183 0.03160 1.03210 0.09904 0.319 0.74970
KLK15.55554 0.27253 1.31328 0.09161 2.975 0.00293 **
DMRTC1.63947 -0.19942 0.81920 0.10511 -1.897 0.05780 .
HTR3B.9177 -0.18328 0.83253 0.10148 -1.806 0.07090 .
PSKH2.85481 0.25614 1.29193 0.09451 2.710 0.00673 **
ADAM2.2515 -0.09113 0.91290 0.10356 -0.880 0.37884
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
NKAIN3.286183 1.0321 0.9689 0.8500 1.253
KLK15.55554 1.3133 0.7615 1.0974 1.572
DMRTC1.63947 0.8192 1.2207 0.6667 1.007
HTR3B.9177 0.8325 1.2012 0.6824 1.016
PSKH2.85481 1.2919 0.7740 1.0735 1.555
ADAM2.2515 0.9129 1.0954 0.7452 1.118
Concordance= 0.69 (se = 0.024 )
Likelihood ratio test= 0 on 6 df, p=1
Wald test = 0 on 6 df, p=1
Score (logrank) test = 2.64 on 6 df, p=0.9