All p-values rendered non-significant in Cox model after initializing coefficients from glmnet

Question

I have a set of 200 genes that are split into numerical high and low, encoded as (1/2).

I have set this variable this way for linearity of the model. Also, stratified by cancer and normal cases.

I tried to do coxph() for all of them and almost all of them are significant with a very low p-values. But the exp(coef) (HR) were crazy huge like 115.0 or something.

so, I decided to do elastic net selection, so as to not bias my selection of covariates to include in the model.

I am using this code snippet below to fit the model, to do a simple cross-validation, and then to get p-values from coxph().

I got 80 non-zero coefficients from glmnet() and I forced these coefficients into coxph()

# x is a model matrix
# y is a Surv matrix

fit <- glmnet(x, stratifySurv(y, genesX$case) , family = "cox" ,
              type.measure ="deviance", maxit = 3000, alpha = 0.5)

# find lambda for which dev.ratio is max
max.dev.index     <- which.max(fit$dev.ratio)
optimal.lambda <- fit$lambda[max.dev.index]

# take beta for optimal lambda
optimal.beta  <- fit$beta[,max.dev.index]

# find non zero beta coef
nonzero.coef <- abs(optimal.beta)>0
table(nonzero.coef)
## 80 genes

sig_genes  <- names(nonzero.coef[nonzero.coef == T])
selectedBeta <- optimal.beta[nonzero.coef] 


# CROSS Validation
cvfit <- cv.glmnet(x, stratifySurv(y, genesX$case) , family = "cox" ,
                   type.measure ="deviance", maxit = 3000, alpha = 0.5, nfolds = 10)

plot(cvfit)
cvfit$lambda.1se
cvfit$lambda.min


# Cox model
coxfit <-  coxph( as.formula(paste("Surv(Stime, OSS) ~  ", paste(sig_genes,collapse="+"), "+strata(case)" )) ,
                 init = selectedBeta,
                 iter = 0,
                 ties = "breslow",
                 data = genesX)


plot(selectedBeta,coef(coxfit))

I find if I force the model like that, I get zero significant genes in the final model. Otherwise, if I run coxph() like this:

coxfit <-  coxph( as.formula(paste("Surv(Stime, OSS) ~  ", paste(sig_genes,collapse="+"), "+strata(case)" )) ,
                 data = t_v_genesX)

I get totally different results where the coefficients plots are not correlated, some of the selected genes are non-significant and again the sky-high coefficients results.

Answer 1

With 54 events, a model might be able to handle 3 to 5 unpenalized predictors without overfitting. The LASSO part of your elastic net with alpha = 0.5 returned 80 predictors, about half of your 200 candidate genes, at the minimum deviance. So presumably the ridge-regression part of your elastic net (the L2 penalty) had done a good deal of penalization on the coefficients of the retained predictors. That is, the magnitudes of those coefficient values were a good deal lower than they would have been in an unpenalized regression.

So you shouldn't be surprised that the magnitudes of those coefficient values might be so small that you couldn't distinguish them from them from 0, in what you hoped to be a standard frequentist test of your null hypotheses. Furthermore, the test you apparently attempted (forcing a standard Cox model to accept those coefficient values and not optimizing further by setting iter = 0) won't do what you hoped. I believe that the Wald tests used to evaluate coefficient significance assume that you are working at the maximum partial likelihood values, while you most certainly are not.

Then, taking your list of 80 genes and just plugging them into an unpenalized Cox regression is undoing all the work that you did to avoid overfitting with the original elastic-net penalization. A Cox model with 80 predictors and only 54 events is going to be wildly overfit if it converges at all.

You need to think carefully about the hypotheses you want to test, and your approach.

First, simply dichotomizing the continuous predictors at their medians is poor practice, throwing away any information you have about their associations with outcome as continuous predictors. Think about it: you have about 300 cases total, so a 50:50 split puts 150 into each of your 2 categories. But you only have 54 events overall: just how good a fit do you think that can provide?

Second, saying "measurement was at different time points by RNA-seq" raises potential red flags. If some of those measurements were only done after the time = 0 point for your survival function, there's a risk that the RNA-seq values are the result of longer survival rather than predictors of longer survival.

Third, your apparent stratification by cancer/not status means that your associations of gene expression with outcome will be irrespective of cancer status, assumed to be the same in both strata (although the baseline survival will differ between strata). If that's what you want to model, fine. But one might expect a study like this to look instead at interactions of gene expression with cancer status in terms of outcome.