Statistics review 12: Survival analysis
=======================================

R code accompanying
`paper <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1065034/pdf/cc2955.pdf>`__

Key learning points
-------------------

-  Analysis of time-to-event data
-  Kaplan–Meier methods
-  log rank test
-  Cox’s proportional hazards

.. code:: r

    suppressPackageStartupMessages(library(tidyverse))

.. code:: r

    options(repr.plot.width=4, repr.plot.height=3)

Survival Data
-------------

.. code:: r

    pt <- 1:15
    days <- c(1,1,4,5,6,8,9,9,12,15,22,25,37,55,72)
    outcome <- c('d', 'd', 'd', 'd', NA, 'd', 's', 'd', 'd',
                NA, 'd', 's', 'd', 'd', 's')
    treatment <- c(2,2,2,2,2,1,2,2,1,1,2,1,1,1,1)
    age <- c(75,79,85,76,66,75,72,70,71,73,66,73,68,59,61)
    
    df <- data.frame(pt=pt, days=days, outcome=outcome, treatment=treatment, age=age)
    df



.. raw:: html

    <table>
    <thead><tr><th scope=col>pt</th><th scope=col>days</th><th scope=col>outcome</th><th scope=col>treatment</th><th scope=col>age</th></tr></thead>
    <tbody>
    	<tr><td> 1</td><td> 1</td><td>d </td><td>2 </td><td>75</td></tr>
    	<tr><td> 2</td><td> 1</td><td>d </td><td>2 </td><td>79</td></tr>
    	<tr><td> 3</td><td> 4</td><td>d </td><td>2 </td><td>85</td></tr>
    	<tr><td> 4</td><td> 5</td><td>d </td><td>2 </td><td>76</td></tr>
    	<tr><td> 5</td><td> 6</td><td>NA</td><td>2 </td><td>66</td></tr>
    	<tr><td> 6</td><td> 8</td><td>d </td><td>1 </td><td>75</td></tr>
    	<tr><td> 7</td><td> 9</td><td>s </td><td>2 </td><td>72</td></tr>
    	<tr><td> 8</td><td> 9</td><td>d </td><td>2 </td><td>70</td></tr>
    	<tr><td> 9</td><td>12</td><td>d </td><td>1 </td><td>71</td></tr>
    	<tr><td>10</td><td>15</td><td>NA</td><td>1 </td><td>73</td></tr>
    	<tr><td>11</td><td>22</td><td>d </td><td>2 </td><td>66</td></tr>
    	<tr><td>12</td><td>25</td><td>s </td><td>1 </td><td>73</td></tr>
    	<tr><td>13</td><td>37</td><td>d </td><td>1 </td><td>68</td></tr>
    	<tr><td>14</td><td>55</td><td>d </td><td>1 </td><td>59</td></tr>
    	<tr><td>15</td><td>72</td><td>s </td><td>1 </td><td>61</td></tr>
    </tbody>
    </table>



Estimating the survival curve using the Kaplan–Meier method
-----------------------------------------------------------

.. code:: r

    library(survival)
    library(ggfortify)

.. code:: r

    df$surv <- with(df, Surv(days, outcome=='d'))

.. code:: r

    head(df)



.. raw:: html

    <table>
    <thead><tr><th scope=col>pt</th><th scope=col>days</th><th scope=col>outcome</th><th scope=col>treatment</th><th scope=col>age</th><th scope=col>surv</th></tr></thead>
    <tbody>
    	<tr><td>1 </td><td>1 </td><td>d </td><td>2 </td><td>75</td><td>1 </td></tr>
    	<tr><td>2 </td><td>1 </td><td>d </td><td>2 </td><td>79</td><td>1 </td></tr>
    	<tr><td>3 </td><td>4 </td><td>d </td><td>2 </td><td>85</td><td>4 </td></tr>
    	<tr><td>4 </td><td>5 </td><td>d </td><td>2 </td><td>76</td><td>5 </td></tr>
    	<tr><td>5 </td><td>6 </td><td>NA</td><td>2 </td><td>66</td><td>6?</td></tr>
    	<tr><td>6 </td><td>8 </td><td>d </td><td>1 </td><td>75</td><td>8 </td></tr>
    </tbody>
    </table>



Overall
~~~~~~~

.. code:: r

    km1 <- survfit(surv ~ 1, data=df)

.. code:: r

    autoplot(km1)





.. image:: SR12_Survival_analysis_files/SR12_Survival_analysis_11_1.png


By treatment
~~~~~~~~~~~~

.. code:: r

    km2 <- survfit(surv ~ treatment, data=df)

.. code:: r

    autoplot(km2)





.. image:: SR12_Survival_analysis_files/SR12_Survival_analysis_14_1.png


Plotting cumulative hazard
~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    autoplot(km2, fun='event')





.. image:: SR12_Survival_analysis_files/SR12_Survival_analysis_16_1.png


Comparing survival curves of two groups using the log rank test
---------------------------------------------------------------

.. code:: r

    survdiff(surv ~ treatment, data=df)



.. parsed-literal::

    Call:
    survdiff(formula = surv ~ treatment, data = df)
    
    n=13, 2 observations deleted due to missingness.
    
                N Observed Expected (O-E)^2/E (O-E)^2/V
    treatment=1 6        4     6.99      1.28      5.27
    treatment=2 7        6     3.01      2.98      5.27
    
     Chisq= 5.3  on 1 degrees of freedom, p= 0.0216 


Cox’s proportional hazards model
--------------------------------

Treatment only
~~~~~~~~~~~~~~

.. code:: r

    m1 <- coxph(surv ~ treatment, data=df)

.. code:: r

    summary(m1)



.. parsed-literal::

    Call:
    coxph(formula = surv ~ treatment, data = df)
    
      n= 13, number of events= 10 
       (2 observations deleted due to missingness)
    
                coef exp(coef) se(coef)    z Pr(>|z|)  
    treatment 1.7603    5.8144   0.8465 2.08   0.0376 *
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
              exp(coef) exp(-coef) lower .95 upper .95
    treatment     5.814      0.172     1.107     30.55
    
    Concordance= 0.713  (se = 0.099 )
    Rsquare= 0.327   (max possible= 0.948 )
    Likelihood ratio test= 5.14  on 1 df,   p=0.02333
    Wald test            = 4.32  on 1 df,   p=0.03756
    Score (logrank) test = 5.28  on 1 df,   p=0.0216



Treatment and age
~~~~~~~~~~~~~~~~~

.. code:: r

    m2 <- coxph(surv ~ treatment + age, data=df)

.. code:: r

    summary(m2)



.. parsed-literal::

    Call:
    coxph(formula = surv ~ treatment + age, data = df)
    
      n= 13, number of events= 10 
       (2 observations deleted due to missingness)
    
                 coef exp(coef) se(coef)     z Pr(>|z|)  
    treatment 1.65039   5.20901  0.96397 1.712   0.0869 .
    age       0.21523   1.24015  0.08562 2.514   0.0119 *
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
              exp(coef) exp(-coef) lower .95 upper .95
    treatment     5.209     0.1920    0.7875    34.458
    age           1.240     0.8064    1.0486     1.467
    
    Concordance= 0.868  (se = 0.118 )
    Rsquare= 0.646   (max possible= 0.948 )
    Likelihood ratio test= 13.5  on 2 df,   p=0.001173
    Wald test            = 8.23  on 2 df,   p=0.01632
    Score (logrank) test = 11.43  on 2 df,   p=0.003298



Check for violation of proportional hazard (constant HR over time)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: r

    res <- cox.zph(m2)
    res



.. parsed-literal::

                rho  chisq     p
    treatment 0.110 0.0994 0.753
    age       0.274 0.4489 0.503
    GLOBAL       NA 0.5000 0.779


.. code:: r

    plot(res)



.. image:: SR12_Survival_analysis_files/SR12_Survival_analysis_28_0.png



.. image:: SR12_Survival_analysis_files/SR12_Survival_analysis_28_1.png


Exercise
--------

.. code:: r

    library(KMsurv)

.. code:: r

    data(tongue)

.. code:: r

    help(tongue)

.. code:: r

    head(tongue)



.. raw:: html

    <table>
    <thead><tr><th scope=col>type</th><th scope=col>time</th><th scope=col>delta</th></tr></thead>
    <tbody>
    	<tr><td>1 </td><td> 1</td><td>1 </td></tr>
    	<tr><td>1 </td><td> 3</td><td>1 </td></tr>
    	<tr><td>1 </td><td> 3</td><td>1 </td></tr>
    	<tr><td>1 </td><td> 4</td><td>1 </td></tr>
    	<tr><td>1 </td><td>10</td><td>1 </td></tr>
    	<tr><td>1 </td><td>13</td><td>1 </td></tr>
    </tbody>
    </table>



**1**. Create a KM plot for the tongue data set.


**2**. Evaluate if the type of tumor affects survival using a log-rank
test


**3**. Evaluate if the type of tumor affects survival using Cox PH test.


**4**. Check if the Cox PH assumptions are met.