<jupyter_start><jupyter_text>Lambda School Data Science - Survival Analysis![My normal approach is useless here, too.](https://imgs.xkcd.com/comics/probability.png)https://xkcd.com/881/The aim of survival analysis is to analyze the effect of different risk factors and use them to predict the duration of time between one event ("birth") and another ("death").  LectureSurvival analysis was first developed by actuaries and medical professionals to predict (as its name implies) how long individuals would survive. However, it has expanded into include many different applications.* it is referred to as **reliability analysis** in engineering* it can be referred to more generally as **time-to-event analysis**In the general sense, it can be thought of as a way to model anything with a finite duration - retention, churn, completion, etc. The culmination of this duration may have a "good" or "bad" (or "neutral") connotation, depending on the situation. However old habits die hard, so most often it is called survival analysis and the following definitions are still commonly used:* birth: the event that marks the beginning of the time period for observation* death: the event of interest, which then marks the end of the observation period for an individual Examples* Customer churn * birth event: customer subscribes to a service * death event: customer leaves the service* Employee retention * birth event: employee is hired * death event: employee quits* Engineering, part reliability * birth event: part is put in use * death event: part fails* Program completion * birth event: student begins PhD program * death event: student earns PhD* Response time * birth event: 911 call is made * death event: police arrive* Lambda School * birth event: student graduates LambdaSchool * death event: student gets a job! Take a moment and try to come up with your own specific example or two. So... if all we're predicting here is a length of time between two events, why can't we just use regular old Linear Regression?Well... if you have all the data, go for it. In some situations it may be reasonably effective. But, data for survival times are often highly skewed and, more importantly, we don't always get a chance to observe the "death" event. The current time or other factors interfere with our ability to observe the time of the event of interest. These observations are said to be _censored_.Additionally, the occurrence or non-occurrence of an event is binary - so, while the time is continuous, the event itself is in some ways similar to a binary event in logistic regression. Censorship in DataSuppose a new cancer treatment is developed. Researchers select 50 individuals for the study to undergo treatment and participate in post-treatment obsesrvation. Birth Event = Participant begins trial Death Event = Participant dies due to cancer or complications of cancerDuring the study:1. Some participants die during the course of the study--triggering their death event 2. Some participants drop out or the researchers otherwise lose contact with them. The researchers have their data up intil the time they dropped out, but they don't have a death event to record3. Some participants are still be alive at the end of the observation period. So again, researchers have their data up until some point, but there is no death event to recordWe only know the interval between the "birth" event and the "death" event for participants in category 1. All others we only know that they survived _up to_ a certain point. Dealing with Censored DataWithout survival analysis, we could deal with censored data in two ways:* We could just treat the end of the observation period as the time of the death event* (Even worse) We could drop the censored data using the rationale that we have "incomplete data" for those observationsBut... both of these will underestimate survival rates for the purpose of the study. We **know** that all those individuals "survived" the "death event" past a certain point.Luckily, in the 1980s a pair of smarty pants named David (main author Cox and coauthor Oakes) did the hard math work to make it possible to incorporate additional features as predictive measures to survival time probabilities. (Fun fact, the one named Cox also came up with logistic regression with non-David coauthor, <NAME>.) lifelinesIt wasn't until 2014 that some other smart people made an implementation of survival analysis in Python called lifelines. It is built over Pandas and follows the same conventions for usage as scikit-learn._Additional note: scikit pushed out a survival analysis implementation last year (2018) named scikit-survival that is imported by the name `sksurv`. It's super new so it may/may not have a bunch of bugs... but if you're interested you can check it out in the future. (For comparison, scikit originally came out in 2007 and Pandas came out in 2008)._<jupyter_code>!pip install lifelines
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lifelines
# lifelines comes with some datasets to get you started playing around with it.
# Most of the datasets are cleaned-up versions of real datasets. Here we will
# use their Leukemia dataset comparing 2 different treatments taken from
# http://web1.sph.emory.edu/dkleinb/allDatasets/surv2datasets/anderson.dat

from lifelines.datasets import load_leukemia
leukemia = load_leukemia()

leukemia.head()<jupyter_output><empty_output><jupyter_text>You can you any Pandas DataFrame with lifelines.  The only requirement is that the DataFrame includes features that describe:* a duration of time for the observation* a binary column regarding censorship (`1` if the death event was observed, `0` if the death event was not observed)Sometimes, you will have to engineer these features. How might you go about that? What information would you need?<jupyter_code>leukemia.info()
leukemia.describe()
time = leukemia.t.values
event = leukemia.status.values

ax = lifelines.plotting.plot_lifetimes(time, event_observed=event)
ax.set_xlim(0, 40)
ax.grid(axis='x')
ax.set_xlabel("Time in Months")
ax.set_title("Lifelines for Survival of Leukemia Patients");
plt.plot();<jupyter_output><empty_output><jupyter_text>Kaplan-Meier survival estimate The Kaplan-Meier method estimates survival probability from observed survival times. It results in a step function that changes value only at the time of each event, and confidence intervals can be computer for the survival probabilities. The KM survival curve,a plot of KM survival probability against time, provides a useful summary of the data.It can be used to estimate measures such as median survival time.It CANNOT account for risk factors and is NOT regression. It is *non-parametric* (does not involve parameters).However it is a good way to visualize a survival dataset, and can be useful to compare the effects of a single categorical variable.<jupyter_code>kmf = lifelines.KaplanMeierFitter()

kmf.fit(time, event_observed=event)
!pip install -U matplotlib  # Colab has matplotlib 2.2.3, we need >3
kmf.survival_function_.plot()
plt.title('Survival Function Leukemia Patients');
print(f'Median Survival: {kmf.median_} months after treatment')
kmf.survival_function_.plot.line()
ax = plt.subplot(111)

treatment = (leukemia["Rx"] == 1)
kmf.fit(time[treatment], event_observed=event[treatment], label="Treatment 1")
kmf.plot(ax=ax)
print(f'Median survival time with Treatment 1: {kmf.median_} months')

kmf.fit(time[~treatment], event_observed=event[~treatment], label="Treatment 0")
kmf.plot(ax=ax)
print(f'Median survival time with Treatment 0: {kmf.median_} months')

plt.ylim(0, 1);
plt.title("Survival Times for Leukemia Treatments");<jupyter_output>Median survival time with Treatment 1: 8.0 months
Median survival time with Treatment 0: 23.0 months<jupyter_text>Cox Proportional Hazards Model -- Survival RegressionIt assumes the ratio of death event risks (hazard) of two groups remains about the same over time.This ratio is called the hazards ratio or the relative risk.All Cox regression requires is an assumption that ratio of hazards is constant over time across groups.The good news—we don’t need to know anything about overall shape of risk/hazard over timeThe bad news—the proportionality assumption can be restrictive<jupyter_code># Using Cox Proportional Hazards model
cph = lifelines.CoxPHFitter()
cph.fit(leukemia, 't', event_col='status')
cph.print_summary()<jupyter_output><lifelines.CoxPHFitter: fitted with 42 observations, 12 censored>
      duration col = 't'
         event col = 'status'
number of subjects = 42
  number of events = 30
    log-likelihood = -69.59
  time fit was run = 2019-01-20 19:09:48 UTC

---
        coef  exp(coef)  se(coef)    z      p  log(p)  lower 0.95  upper 0.95     
sex     0.31       1.37      0.45 0.69   0.49   -0.72       -0.58        1.21     
logWBC  1.68       5.38      0.34 5.00 <0.005  -14.36        1.02        2.34  ***
Rx      1.50       4.50      0.46 3.26 <0.005   -6.79        0.60        2.41    *
---
Signif. codes: 0 '***' 0.0001 '**' 0.001 '*' 0.01 '.' 0.05 ' ' 1

Concordance = 0.85
Likelihood ratio test = 47.19 on 3 df, log(p)=-21.87<jupyter_text>Interpreting the Results`coef`: usually denoted with $b$, the coefficient`exp(coef)`: $e^{b}$, equals the estimate of the hazard ratio. Here, we can say that participants who received treatment 1 had ~4.5 times the hazard risk (risk of death) compared to those who received treatment 2. And for every unit the `logWBC` increased, the hazard risk increased >5 times.`se(coef)`: standard error of the coefficient (used for calculating z-score and therefore p-value)`z`: z-score $\frac{b}{se(b)}$`p`: p-value. derived from z-score. describes statistical significance. more specifically, it is the likelihood that the variable has no effect on the outcome`log(p)`: natural logarithm of p-value... used to more easily see differences in significance`lower/upper 0.95`: confidence levels for the coefficients. in this case, we can confidently say that the coefficient for `logWBC` is somewhere _between_ 1.02 and 2.34.`Signif. codes`: easily, visually identify significant variables! The more stars, the more solid (simply based on p-value). Here `logWBC` is highly significant, `Rx` is significant, and `sex` has no statistical significance`Concordance`: a measure of predictive power for classification problems (here looking at the `status` column. a value from 0 to 1 where values above 0.6 are considered good fits (the higher the better)`Likelihood ratio (LR) test`: this is a measure of how likely it is that the coefficients are not zero, and can compare the goodness of fit of a model versus an alternative null model. Is often actually calculated as a logarithm, resulting in the log-likelihood ratio statistic and allowing the distribution of the test statistic to be approximated with [Wilks' theorem](https://en.wikipedia.org/wiki/Wilks%27_theorem).<jupyter_code>cph.plot_covariate_groups(covariate='logWBC', groups=np.arange(1.5,5,.5))
cph.plot_covariate_groups(covariate='sex', groups=[0,1])<jupyter_output><empty_output><jupyter_text>Remember how the Cox model assumes the ratio of death events between groups remains constant over time?Well we can check for that.<jupyter_code>cph.check_assumptions(leukemia)
# We can see that the sex variable is not very useful by plotting the coefficients
cph.plot()
# Let's do what the check_assumptions function suggested
cph = lifelines.CoxPHFitter()
cph.fit(leukemia, 't', event_col='status', strata=['sex'])
cph.print_summary()
cph.baseline_cumulative_hazard_.shape<jupyter_output><lifelines.CoxPHFitter: fitted with 42 observations, 12 censored>
      duration col = 't'
         event col = 'status'
            strata = ['sex']
number of subjects = 42
  number of events = 30
    log-likelihood = -55.73
  time fit was run = 2019-01-20 19:09:51 UTC

---
        coef  exp(coef)  se(coef)    z      p  log(p)  lower 0.95  upper 0.95     
logWBC  1.45       4.28      0.34 4.23 <0.005  -10.64        0.78        2.13  ***
Rx      1.00       2.71      0.47 2.11   0.04   -3.35        0.07        1.93    .
---
Signif. codes: 0 '***' 0.0001 '**' 0.001 '*' 0.01 '.' 0.05 ' ' 1

Concordance = 0.85
Likelihood ratio test = 74.90 on 2 df, log(p)=-37.45<jupyter_text>Notice that this regression has `Likelihood ratio test = 74.90 on 2 df, log(p)=-37.45`, while the one that included `sex` had `Likelihood ratio test = 47.19 on 3 df, log(p)=-21.87`. The LRT is higher and log(p) is lower, meaning this is likely a better fitting model.<jupyter_code>cph.plot()
cph.compute_residuals(leukemia, kind='score')
cph.predict_cumulative_hazard(leukemia[:5])
surv_func = cph.predict_survival_function(leukemia[:5])

exp_lifetime = cph.predict_expectation(leukemia[:5])

plt.plot(surv_func)

exp_lifetime
# lifelines comes with some datasets to get you started playing around with it
# The Rossi dataset originally comes from Rossi et al. (1980), 
# and is used as an example in Allison (1995). 

# The data pertain to 432 convicts who were released from Maryland state prisons 
# in the 1970s and who were followed up for one year after release. Half the 
# released convicts were assigned at random to an experimental treatment in 
# which they were given financial aid; half did not receive aid.

from lifelines.datasets import load_rossi
recidivism = load_rossi()

recidivism.head()

# Looking at the Rossi dataset, how long do you think the study lasted?

# All features are coded with numerical values, but which features do you think 
# are actually categorical?
recidivism.info()
recidivism.describe()<jupyter_output><empty_output><jupyter_text>These are the "lifelines" of the study participants as they attempt to avoid recidivism<jupyter_code>recidivism_sample = recidivism.sample(n=25)

duration = recidivism_sample.week.values
arrested = recidivism_sample.arrest.values

ax = lifelines.plotting.plot_lifetimes(duration, event_observed=arrested)
ax.set_xlim(0, 78)
ax.grid(axis='x')
ax.vlines(52, 0, 25, lw=2, linestyles='--')
ax.set_xlabel("Time in Weeks")
ax.set_title("Recidivism Rates");
plt.plot();
kmf = lifelines.KaplanMeierFitter()

duration = recidivism.week
arrested = recidivism.arrest

kmf.fit(duration, arrested)
kmf.survival_function_.plot()
plt.title('Survival Curve:\nRecidivism of Recently Released Prisoners');
kmf.plot()
plt.title('Survival Function of Recidivism Data');
print(f'Median time before recidivism: {kmf.median_} weeks')
kmf_w_aid = lifelines.KaplanMeierFitter()
kmf_no_aid = lifelines.KaplanMeierFitter()

ax = plt.subplot(111)

w_aid = (recidivism['fin']==1)

t = np.linspace(0, 70, 71)
kmf_w_aid.fit(duration[w_aid], event_observed=arrested[w_aid], timeline=t, label="Received Financial Aid")
ax = kmf_w_aid.plot(ax=ax)
#print("Median survival time of democratic:", kmf.median_)

kmf_no_aid.fit(duration[~w_aid], event_observed=arrested[~w_aid], timeline=t, label="No Financial Aid")
ax = kmf_no_aid.plot(ax=ax)
#print("Median survival time of non-democratic:", kmf.median_)

plt.ylim(.5,1)
plt.title("Recidivism for Participants Who Received Financial Aid \vs. Those Who Did Not");
naf = lifelines.NelsonAalenFitter()
naf.fit(duration, arrested)

print(naf.cumulative_hazard_.head())
naf.plot()
naf_w_aid = lifelines.NelsonAalenFitter()
naf_no_aid = lifelines.NelsonAalenFitter()

naf_w_aid.fit(duration[w_aid], event_observed=arrested[w_aid], timeline=t, label="Received Financial Aid")
ax = naf_w_aid.plot(loc=slice(0, 50))
naf_no_aid.fit(duration[~w_aid], event_observed=arrested[~w_aid], timeline=t, label="No Financial Aid")
ax = naf_no_aid.plot(ax=ax, loc=slice(0, 50))
plt.title("Recidivism Cumulative Hazard\nfor Participants Who Received Financial Aid \nvs. Those Who Did Not");
plt.show()
cph = lifelines.CoxPHFitter()
cph.fit(recidivism, duration_col='week', event_col='arrest', show_progress=True)

cph.print_summary()
cph.plot()
cph.plot_covariate_groups('fin', [0, 1])
cph.plot_covariate_groups('prio', [0, 5, 10, 15])
r = cph.compute_residuals(recidivism, 'martingale')
r.head()
cph = lifelines.CoxPHFitter()
cph.fit(recidivism, duration_col='week', event_col='arrest', show_progress=True)

cph.print_summary()
cph.plot();
cph.plot_covariate_groups('prio', [0, 5, 10, 15]);
cph.check_assumptions(recidivism)<jupyter_output><lifelines.StatisticalResult>
         test_name = proportional_hazard_test
 null_distribution = chi squared
degrees_of_freedom = 1

---
               test_statistic      p  log(p)    
age  identity           12.06 <0.005   -7.57  **
     km                 11.03 <0.005   -7.02  **
     log                13.07 <0.005   -8.11  **
     rank               11.09 <0.005   -7.05  **
fin  identity            0.06   0.81   -0.21    
     km                  0.02   0.89   -0.12    
     log                 0.49   0.48   -0.73    
     rank                0.02   0.90   -0.11    
mar  identity            0.75   0.39   -0.95    
     km                  0.60   0.44   -0.82    
     log                 1.03   0.31   -1.17    
     rank                0.67   0.41   -0.88    
paro identity            0.12   0.73   -0.32    
     km                  0.12   0.73   -0.31    
     log                 0.01   0.92   -0.09    
     rank                0.14   0.71   -0.34    
prio identity            0.01[...]<jupyter_text>The Intuition - Hazard and Survival Functions Hazard Function - the dangerous bathtubThe hazard function represents the *instantaneous* likelihood of failure. It can be treated as a PDF (probability density function), and with real-world data comes in three typical shapes.![Different hazard functions](https://upload.wikimedia.org/wikipedia/commons/2/25/Compsyseng17_04.jpg)Increasing and decreasing failure rate are fairly intuitive - the "bathtub" shaped is perhaps the most surprising, but actually models many real-world situations. In fact, life expectancy in general, and most threats to it, assume this shape.What the "bathtub" means is that - threats are highest at youth (e.g. infant mortality), but then decrease and stabilize at maturity, only to eventually re-emerge in old age. Many diseases primarily threaten children and elderly, and middle aged people are also more robust to physical trauma.The "bathtub" is also suitable for many non-human situations - often with reliability analysis, mechanical parts either fail early (due to manufacturing defects), or they survive and have a relatively long lifetime to eventually fail out of age and use. Survival Function (aka reliability function) - it's just a (backwards) CDFSince the hazard function can be treated as a probability density function, it makes sense to think about the corresponding cumulative distribution function (CDF). But because we're modeling time to failure, it's actually more interesting to look at the CDF backwards - this is called the complementary cumulative distribution function.In survival analysis there's a special name for it - the survival function - and it gives the probability that the object being studied will survive beyond a given time.![4 survival functions](https://upload.wikimedia.org/wikipedia/commons/e/e0/Four_survival_functions.svg)As you can see they all start at 1 for time 0 - at the beginning, all things are alive. Then they all move down over time to eventually approach and converge to 0. The different shapes reflect the average/expected retention of a population subject to this function over time, and as such this is a particularly useful visualization when modeling overall retention/churn situations. Ways to estimate/model survival analysis - terms to be aware ofKey Components Necessary for these models - duration, and whether observation is censored.- Kaplan Meier Estimator- Nelson-Aalen Estimator- Proportional Hazards (Cox Model, integrates covariates)- Additive Hazards Model (Aalen's Additive Model, when covariates are time-dependent)As with most statistics, these are all refinements of the general principles, with the math to back them up. Software packages will tend to select reasonable defaults, and allow you to use parameters to tune or select things. The math for these gets varied and deep - but feel free to [dive in](https://en.wikipedia.org/wiki/Survival_analysis) if you're curious!  Live! Let's try modeling heart attack survivalhttps://archive.ics.uci.edu/ml/datasets/echocardiogram<jupyter_code># Loading the data to get you started
churn_data = pd.read_csv(
    'https://raw.githubusercontent.com/treselle-systems/'
    'customer_churn_analysis/master/WA_Fn-UseC_-Telco-Customer-Churn.csv')
churn_data.head()
!pip install -U matplotlib
churn_data.info()
def encode_cat(df):
    for col_name in df.columns:
        if(df[col_name].dtype == 'object'):
            df[col_name]= df[col_name].astype('category')
            df[col_name] = df[col_name].cat.codes
    return df

churn_data = encode_cat(churn_data)
churn_data['TotalCharges'] = pd.to_numeric(churn_data['TotalCharges'], errors = 'coerce')
churn_data.isna().sum()
total_charges_mean = churn_data['TotalCharges'].mean()
churn_data['TotalCharges'].fillna(total_charges_mean, inplace=True)
churn_data['Churn'] = churn_data['Churn'].replace('No', 0)
churn_data['Churn'] = churn_data['Churn'].replace('Yes', 1)
churn_data.tenure.value_counts()
churn_data.tenure.describe()
sample_churn_data = churn_data.sample(n=40, random_state=42)
churn_data.Churn.value_counts()
time = sample_churn_data.tenure.values
event = sample_churn_data.Churn.values

ax = lifelines.plotting.plot_lifetimes(time, event_observed = event)
ax.set_xlim(0, 75)
ax.grid(axis = 'x')
ax.set_xlabel("Time in Months")
ax.set_title("Customer Retention");
plt.plot();
kmf = lifelines.KaplanMeierFitter()

kmf.fit(time, event_observed = event)
kmf.survival_function_.plot()
plt.title('Kaplan-Meier Estimate Function Customer Retention');
print(f'Median: {kmf.median_} months after commitment')
churn_data.dtypes
ax = plt.subplot(111)

treatment = (sample_churn_data["Churn"] == 1)
kmf.fit(time[treatment], event_observed = event[treatment])

kmf.plot(ax = ax)
print(f'Median retention time: {kmf.median_} months')

plt.ylim(0, 1);
plt.title("Retention Time for Customers");
# Using Cox Proportional Hazards model
cph = lifelines.CoxPHFitter()
cph.fit(churn_data, 'tenure', event_col = 'Churn')
cph.print_summary()<jupyter_output><lifelines.CoxPHFitter: fitted with 7043 observations, 5174 censored>
      duration col = 'tenure'
         event col = 'Churn'
number of subjects = 7043
  number of events = 1869
    log-likelihood = -12677.95
  time fit was run = 2019-01-23 00:52:11 UTC

---
                  coef  exp(coef)  se(coef)      z      p  log(p)  lower 0.95  upper 0.95     
customerID       -0.00       1.00      0.00  -1.65   0.10   -2.32       -0.00        0.00     
gender           -0.04       0.96      0.05  -0.94   0.35   -1.06       -0.13        0.05     
SeniorCitizen     0.04       1.04      0.06   0.69   0.49   -0.71       -0.07        0.15     
Partner          -0.18       0.84      0.05  -3.24 <0.005   -6.75       -0.29       -0.07    *
Dependents       -0.09       0.91      0.07  -1.32   0.19   -1.68       -0.23        0.04     
PhoneService     -0.22       0.81      0.12  -1.84   0.07   -2.71       -0.45        0.01     
MultipleLines    -0.08       0.92      0.03  -2.95 <0.005   -5.74       -[...]
{'fraction_non_alphanumeric': 0.0876553793399057, 'fraction_numerical': 0.034376339477068155, 'mean_word_length': 3.9503500954805855, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 39, 'https://': 8, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '26142', 'license': 'MIT', 'n_tokens_mistral': 7632, 'n_tokens_neox': 6836, 'n_words': 2715, 'path': 'module2-survival-analysis/LS_DS_232_Survival_Analysis.ipynb', 'repo_name': 'DanielleRomanoff/DS-Unit-2-Sprint-3-Advanced-Regression', 'size': '23335'}