<p>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.</p>

<p>Here is a <a href="http://www.sthda.com/english/wiki/survival-analysis-basics" rel="nofollow noreferrer">nice tutorial</a> for doing survival analysis with the <code>survival</code> and <code>survminer</code> packages. The latter includes <code>ggplot2</code> kind visualization of the KM curves.</p>

<p>For the Cox regression, you can use continuous(*) or discrete variables, in your case these would be <em>the normalized expression value of your gene</em> and <em>categories of your gene</em> respectively. You can categorize the expression level of a gene as <em>low vs high</em>, or <em>low vs medium vs high</em>, based on different quantiles.  Another approach to categorize your expression data would be to find an "optimal" cut point with <a href="https://rpkgs.datanovia.com/survminer/reference/surv_cutpoint.html" rel="nofollow noreferrer">maximally selected rank statistics</a>; from the link: </p>

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

<p>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.</p>

<p>* <a href="https://stats.stackexchange.com/questions/228493/can-a-cox-proportional-hazards-model-be-built-only-with-continuous-predictors">Can a Cox Proportional-Hazards Model be built only with continuous predictors?</a></p>

CC BY-SA 4.0
2019-10-11T09:11:44.107Z
Thank you but I have about to 130 genes here for which I want to know if they together are related to survival or not; So should I take mean of their expression?
Zizogolu
CC BY-SA 4.0
2019-10-11T09:27:45.747Z
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.
haci
CC BY-SA 4.0
2019-10-11T10:30:36.310Z
Sorry I edited my question because I faced a problem
Zizogolu
CC BY-SA 4.0
2019-10-11T10:54:55.653Z
The survival object should be created with time and status/event (i.e. if the patient is dead/alive) like `Surv(time, status)`, you lack the status variable  in your data.
haci
CC BY-SA 4.0
2019-10-11T11:03:11.177Z
Sorry because I don't know if these patients alive/dead; All I know is survival days from GSE19417 accession :(
Zizogolu
CC BY-SA 4.0
2019-10-11T11:21:01.390Z
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.
haci
CC BY-SA 4.0
2019-10-11T11:28:44.033Z
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.
haci
CC BY-SA 4.0
2019-10-11T11:43:18.200Z
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.
haci
CC BY-SA 4.0
2019-10-11T13:41:49.813Z
Sorry I have a categorical column to data for visualization in my main post but I am getting error
Zizogolu
CC BY-SA 4.0
2019-10-16T14:19:43.940Z
`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.
haci
CC BY-SA 4.0
2019-10-10T17:09:23.877Z
2019-10-11T11:22:54.373Z
2019-10-11T11:22:54.373Z
haci
<p>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</p>

<blockquote>
  <p><strong><em>EDITED</em></strong></p>
</blockquote>

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

<pre><code>&gt;
</code></pre>

<p>I done so but nothing I got</p>

<pre><code>&gt; fit &lt;- survfit(Surv(time) ~ gene, data = km)
&gt; 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
</code></pre>

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

<blockquote>
  <p><strong><em>By your help</em></strong></p>
</blockquote>

<pre><code>&gt; fit &lt;- coxph(Surv(time) ~ gene, data = km)
&gt; print(fit)
Call:
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 
</code></pre>

<p>For plotting I guess I do need another information because</p>

<pre><code> library("survminer")

&gt; ggsurvplot(fit, data = km)
Error in ggsurvplot(fit, data = a) : object 'ggsurv' not found
</code></pre>

<blockquote>
  <p><strong><em>EDITED</em></strong></p>
</blockquote>

<p>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</p>

<pre><code> fit &lt;- survfit(Surv(Time, Status) ~ Gene, data = km)


&gt; ggsurvplot(fit, data = km)
Error in data.frame(..., check.names = FALSE) : 
  arguments imply differing number of rows: 70, 0, 140
</code></pre>

<p>Any help please?</p>

<p>Thanks</p>

<r><gene-expression><linear-regression>
Doing plot with this data
# Doing plot with this data

(Asked by: username_0 on 2019-10-10)

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

EDITED

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)
Call:
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
 library("survminer")

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


EDITED

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?
Thanks

**Tags:** r, gene-expression, linear-regression

---

## 1 answers

### Answer by username_1 on 2019-10-10. Score: 3 | Accepted answer

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](http://www.sthda.com/english/wiki/survival-analysis-basics) 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](https://rpkgs.datanovia.com/survminer/reference/surv_cutpoint.html); 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?](https://stats.stackexchange.com/questions/228493/can-a-cox-proportional-hazards-model-be-built-only-with-continuous-predictors)

> **Comments:**
> **username_0**: Thank you but I have about to 130 genes here for which I want to know if they together are related to survival or not; So should I take mean of their expression?  
> **username_1**: 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.  
> **username_0**: Sorry I edited my question because I faced a problem  
> **username_1**: The survival object should be created with time and status/event (i.e. if the patient is dead/alive) like `Surv(time, status)`, you lack the status variable  in your data.  
> **username_0**: Sorry because I don't know if these patients alive/dead; All I know is survival days from GSE19417 accession :(  
> **username_1**: 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.  
> **username_1**: 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.  
> **username_1**: 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.  
> **username_0**: Sorry I have a categorical column to data for visualization in my main post but I am getting error  
> **username_1**: `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.  

---

