I have normalized expression values of some genes in some patients and the how many days they have survived from their diagnosis with cancer like below


I have taken mean of the expression values for each patient and I have how many days each patient survived since diagnosis as below


I done so but nothing I got

> fit <- survfit(Surv(time) ~ gene, data = km)
> print(fit)
Call: survfit(formula = Surv(time) ~ gene, data = km)

                  n events median 0.95LCL 0.95UCL
gene=-0.333140816 1      1    859      NA      NA
gene=-0.307846735 1      1    347      NA      NA
gene=-0.303559694 1      1   1339      NA      NA
gene=-0.290518776 1      1     61      NA      NA

The question is if these genes increase survival time or not but likely I am doing wrong

By your help

> fit <- coxph(Surv(time) ~ gene, data = km)
> print(fit)
coxph(formula = Surv(time) ~ gene, data = km)

       coef exp(coef) se(coef)     z     p
gene 0.8664    2.3783   1.4520 0.597 0.551

Likelihood ratio test=0.35  on 1 df, p=0.5543
n= 69, number of events= 69 

For plotting I guess I do need another information because


> ggsurvplot(fit, data = km)
Error in ggsurvplot(fit, data = a) : object 'ggsurv' not found


I divided patients based on the median of survival days to Up and Down category like below but still I am failing to visualize that

 fit <- survfit(Surv(Time, Status) ~ Gene, data = km)

> ggsurvplot(fit, data = km)
Error in data.frame(..., check.names = FALSE) : 
  arguments imply differing number of rows: 70, 0, 140

Any help please?



1 Answer 1


In general, survival analysis can be said to be composed of two steps; Cox regression, with which you calculate the "hazard ratio" based on your variables, and a "Kaplan-Meier (KM) estimate", which is used to visuazlize the data.

Here is a nice tutorial for doing survival analysis with the survival and survminer packages. The latter includes ggplot2 kind visualization of the KM curves.

For the Cox regression, you can use continuous(*) or discrete variables, in your case these would be the normalized expression value of your gene and categories of your gene respectively. You can categorize the expression level of a gene as low vs high, or low vs medium vs high, based on different quantiles. Another approach to categorize your expression data would be to find an "optimal" cut point with maximally selected rank statistics; from the link:

This is an outcome-oriented methods providing a value of a cutpoint that correspond to the most significant relation with outcome (here, survival)

Basically you search for the cut-point where your data is split into two such that the two resulting parts of the data show the most separation/difference in terms of survival.

* Can a Cox Proportional-Hazards Model be built only with continuous predictors?

  • 1
    $\begingroup$ In a scenario where you have a "gene module", yes. Basically you would be creating a new variable (mean expression of your module genes) and use this in your Cox regression. If the aim would be looking at the synergistic effects of genes to survival, that is a whole different story for which you would need to use interaction terms in your Cox regression. $\endgroup$
    – haci
    Oct 11, 2019 at 9:27
  • 1
    $\begingroup$ Can you try fit <- coxph(Surv(time) ~ gene, data = km) instead of fit <- survfit(Surv(time) ~ gene, data = km)? coxph seems to be way when using continuous variables. I am adding a link to the answer so it does not get lost here int the comments. $\endgroup$
    – haci
    Oct 11, 2019 at 11:21
  • 1
    $\begingroup$ And regarding lack of patient status in your dataset, without that info, the model assumes an event has happened for all, the following is from ?Surv: Although unusual, the event indicator can be omitted, in which case all subjects are assumed to have an event. That is why you are not getting an error but please pay attention to the consequences of this assumption on your results. $\endgroup$
    – haci
    Oct 11, 2019 at 11:28
  • 1
    $\begingroup$ For the KM curves, you should use categorical data, so your options are using i) survfit (with categories of your variable) for Cox and KM or ii) coxph with your continuous variable for Cox and survfit with categories of your data for KM. $\endgroup$
    – haci
    Oct 11, 2019 at 11:43
  • 1
    $\begingroup$ fit <- survfit(Surv(Time, Status) ~ Gene, data = km), here the Status is not about your variable of interest, it is whether or not the event has happened. You can try fit <- survfit(Surv(Time) ~ Status, data = km). With this, you attempt to see the effect of your gene's expression categories, high and low specifically. BTW, it might be misleading to choose Status as a variable name for the categories of your gene. $\endgroup$
    – haci
    Oct 16, 2019 at 14:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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