I am performing a survival analysis using a glm and predict function and want to validate the results using k-fold cross validation. I am wondering how I interpret the output of the cross-validation. Could you please look at he example-data and tell me what the plot_kfold_cv and cv.glm really shows?

genes.test <- c("ENSG00000092470","ENSG00000164118")

df.t <- structure(list(ENSG00000092470 = c(1.46425216384238, 4.99996145653544, 
0.955931911035645, 1.30190334510162, 2.64295709016934, 1.769816490456
), ENSG00000164118 = c(2.45086739393274, 4.54384728835685, 1.62555232861294, 
2.26725511511905, 3.33218577119466, 2.46092825618336), OverallStatus.test = c(1L, 
1L, 0L, 0L, 1L, 0L), Comb = c(23.2015732395141, 27.0638350315371, 
-253.137052536687, -22.9959568690356, 121.031215304535, -66.296876247146
)), row.names = c("86", "175", "217", "394", "444", "618"), class = "data.frame")

OverallStatus.test <- c(1L, 1L, 0L, 0L, 1L, 0L)
OverallTime.test <- structure(c(71.0416666666667, 601.958333333333, 1485, 1748, 117.041666666667, 
1863), class = "difftime", units = "days")

FML.mrna <- as.formula(paste0('OverallStatus.test~',paste(genes.test,collapse = '+')))
multifit.mrna <- glm(FML.mrna,family=binomial(), data = df.t)
df.t$Comb <- predict.glm(multifit.mrna,type = "link")
                   marker = df.t$Comb,
                   predict.time =  cutoff, method="KM")
plot(Marker$FP, Marker$TP,type="l",col="#BC3C29FF")

legend(0.2,0.2,paste("AUC =",round(Marker$AUC,3)),
       x.intersp=1, y.intersp=0.8,
       lty= 1 ,lwd= 2,col=c("#BC3C29FF","#0072B5FF"),
       bty = "n",

plot_kfold_cv(df.t, multifit.mrna,k = 5)
cv.mse <- cv.glm(df.t, multifit.mrna) 

Predict using coxph:

coxph <- coxph(Surv(OverallTime.test, OverallStatus.test) ~ ENSG00000092470+ENSG00000164118, data = df.t)

[1]  0.8746219  0.8813691 -3.7266492  0.1066562  2.4872752 -0.6232731
  • 1
    $\begingroup$ Is there some reason why you are analyzing your survival data as a binomial/logistic regression instead of as a survival model that takes time-to-event values into account? $\endgroup$
    – EdM
    Dec 4 '20 at 16:44
  • $\begingroup$ I want to combine multiple genes and found this to be the only way to make the combined signature. If you have suggestion how to make and plot a combined signature of multiple genes (the Comb object) I am happy to hear! Thanks. $\endgroup$ Dec 4 '20 at 18:16

One of the difficulties in using software based on multiple packages, like R, is keeping track of which packages are suitable for use in your research program. Here, you seem to be trying to use gene-expression values as a potential biomarker for survival outcomes. Each of the steps you have taken might make sense in some contexts, but there are better ways to proceed with your data.

First, you use a logistic regression (glm() with binomial family) to develop the biomarker. Your "Comb" values from predict() with type = "link" are what's called the linear predictors: the sum, for each case, of terms from the products of each predictor value times its fitted coefficient, plus the intercept. In the binomial model, that provides the estimated log-odds of the outcome.

That model, however, throws away all the information about survival times. That's not a good idea in general. It would be better to get a linear predictor that takes survival times into account. You can do that by fitting a coxph() model from the survival package, then using predict() on the coxph output object with the default type = "lp".

Second, the survivalROC() function from the eponymous package evaluates how well a particular biomarker is associated with survival outcomes at a particular point in time, estimating a receiver operating characteristic (ROC) curve in a way that incorporates censoring in the survival data. You certainly could use the linear predictor from a Cox model instead of from your binomial model with this function. I've used that myself. It would be better for your "biomarker" also to include information about critical clinical variables, which could be linear predictors from a model that incorporates clinical variables along with the gene-expression values. If a Cox survival model fits well, however, the standard C-index reported for a model that includes the biomarker provides information related to the area under the ROC curve, a reasonable measure of model performance, although without providing the ROC curve itself.

Third, the plot_kfold_cv() function evidently is from the sjPlot package. According to the manual page, however, "Currently, only linear, poisson and negative binomial regression models are supported." Neither your binomial model nor the Cox model that I suggest would thus seem to provide interpretable output from that function.

Fourth, the cv.glm() function, by default as you invoke it, reports the mean-square error over the 5 folds, both raw and corrected for bias. I'm not sure, but I suspect that with a logistic/binomial model it would use the squares of the differences between the predicted probabilities and the observed {0,1} outcomes, related to the Brier score. Check the code to be sure. For a survival model you would want to use some other measure of model performance, like observed versus predicted survival probabilities at some point in time. I prefer bootstrapping to evaluate model performance in terms of validation and calibration, and use cross-validation instead to choose things like hyperparatmeter settings for penalized models with LASSO or ridge regression.

It might make sense for you to step back a bit and engage in an overview of regression techniques and their application to things like survival data. In addition to many pages on this site, take a look at the resources provided by Frank Harrell under "RMS" (regression modeling strategies).

  • $\begingroup$ Thanks for a detailed answer, cleared up a few things. Would you mind sharing the script how to do this on my example data: "You can do that by fitting a coxph() model from the survival package, then using predict() on the coxph output object with the default type = "lp"" $\endgroup$ Dec 4 '20 at 20:28
  • $\begingroup$ I updated the answer with a coxph predict instead. Is this correct? $\endgroup$ Dec 4 '20 at 20:33
  • $\begingroup$ @user2300940 yes, that's how to proceed with a coxph survival model. A particular predict() function is chosen based on the type of object it's working on, and the default with predict.coxph is to return the linear predictor. It's always good to check the defaults on programs that you use, to make sure they are what you want. It would be good to read up on survival models, from the survival package documentation and the Harrell resources. $\endgroup$
    – EdM
    Dec 4 '20 at 20:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.