Survival Analysis for C-MAPSS turbofan jet engines ~ predictive maintenance

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import os
sns.set(rc = {'figure.figsize':(16,9)})

import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)

columns = ["machine_name", "cycle", "operational_setting_1", "operational_setting_2", "operational_setting_3"] + [f'sensor_measurement_{i:02}' for i in range(1,22)]

# Read data
turbofan_df = pd.read_csv("train_FD002.txt", header = None, sep = "\s+", names = columns)
turbofan_df_test = pd.read_csv("test_FD002.txt", header = None, sep = "\s+", names = columns)
turbofan_df.head()

turbofan_df[turbofan_df.machine_name == 9
           ]

# backup timeseries data
timeseries_df = turbofan_df.copy()
timeseries_df_test = turbofan_df_test.copy()
timeseries_df.head()

# Select maximum cycle
max_cycle = turbofan_df.groupby(by = "machine_name")['cycle'].transform(max)
#max_cycle = turbofan_df.groupby(by = "machine_name")['cycle'].max()
max_cycle

turbofan_df = turbofan_df[turbofan_df["cycle"] == max_cycle].set_index('machine_name')
turbofan_df.head()

# Lollipop plot for each machine name
plt.hlines(y=turbofan_df.index, xmin=1, xmax=turbofan_df['cycle'], color='darkgrey')
plt.plot(turbofan_df['cycle'], turbofan_df.index, "o", color='royalblue')
plt.plot([1 for i in range(len(turbofan_df))], turbofan_df.index, "o")

ylim = plt.gca().get_ylim()
plt.ylim(ylim[1], ylim[0])

# Add titles and axis names
plt.title("Cycle runtime of engine until failure event")
plt.xlabel("Cycle")
plt.ylabel('Machine ID')

# Show the plot
plt.show()

# Create status column
turbofan_df['status'] = turbofan_df['cycle'].apply(lambda x: False if x > 215 else True)

turbofan_df.head(5)

# Distribution of each status
turbofan_df['status'].value_counts()

turbofan_df.describe()

turbofan_df.iloc[:,8:].head()

turbofan_df.nunique()

# Change to category
cat = ['operational_setting_3', 'sensor_measurement_16']

turbofan_df[cat] = turbofan_df[cat].astype('category')
turbofan_df.info()

sns.heatmap(turbofan_df.corr()**2, annot = True,)

selected_columns = ['sensor_measurement_04','sensor_measurement_11','sensor_measurement_14']
cleaned_data = turbofan_df.loc[:, selected_columns + cat + ['status', 'cycle']]
sns.heatmap(cleaned_data.corr(), annot = True)

# One Hot Encoding for Categorical Variable
from sksurv.preprocessing import OneHotEncoder

data_x = OneHotEncoder().fit_transform(cleaned_data.iloc[:, :-2])
data_x.head()

# Preprocessing for target variable

data_y = list(cleaned_data.loc[:, ["status", "cycle"]].itertuples(index = None, name = None))
data_y = np.array(data_y, dtype=[('status', bool), ('cycle', float)])
data_y[:5]

# Preprocessing for target variable

target_data = list(cleaned_data.loc[:, ["status", "cycle"]].itertuples(index = None, name = None))
target_data = np.array(target_data, dtype=[('status', bool), ('cycle', float)])
target_data[:5]

from lifelines import KaplanMeierFitter

# Instantiate Kaplan Meier object
kmf = KaplanMeierFitter()

# Fit Kaplan Meier estimator
kmf.fit(durations=target_data['cycle'], event_observed=target_data['status'])

# Plot survival function with confidence intervals
kmf.plot_survival_function()

plt.title("Kaplan-Meier Estimator")
plt.ylabel("est. probability of survival $\hat{S}(t)$")
plt.xlabel("time $t$")

plt.show()

os3_unique = list(cleaned_data['operational_setting_3'].unique())
sm16_unique = list(cleaned_data['sensor_measurement_16'].unique())
print(f'Unique values operational_setting_3: {os3_unique}')
print(f'Unique values sensor_measurement_16: {sm16_unique}')

# plt.xlabel callable bugfix
import matplotlib.pyplot as plt
from importlib import reload
plt=reload(plt)

fig, ax = plt.subplots(1, 1, sharex='col', sharey='row')

# Mask for operational_setting_3
fullos3 = (cleaned_data['operational_setting_3'] == 100)

kmf = KaplanMeierFitter()
kmf.fit(durations=cleaned_data[fullos3]['cycle'], event_observed=cleaned_data[fullos3]['status'])
kmf.plot_survival_function(ax=ax)

kmf = KaplanMeierFitter()
kmf.fit(durations=cleaned_data[~fullos3]['cycle'], event_observed=cleaned_data[~fullos3]['status'])
kmf.plot_survival_function(ax=ax)


plt.title(f"Kaplan-Meier Estimator for both operational settings 3")
plt.legend(loc="best")
plt.xlabel('cycles')
plt.ylabel("est. probability of survival $\hat{S}(t)$")
plt.show()

# plt.xlabel callable bugfix
import matplotlib.pyplot as plt
from importlib import reload
plt=reload(plt)

threshold = 8140

fig, ax = plt.subplots(1, 1, sharex='col', sharey='row')

# Mask
sm14 = (cleaned_data['sensor_measurement_14'] > threshold)

kmf_exp = KaplanMeierFitter()
kmf_exp.fit(durations=cleaned_data[sm14]['cycle'], event_observed=cleaned_data[sm14]['status'], label='exp')
kmf_exp.plot_survival_function(ax=ax, show_censors=True, ci_legend=True, 
                               label='sm14 greater than ' + str(threshold)
                              )

kmf_control = KaplanMeierFitter()
kmf_control.fit(durations=cleaned_data[~sm14]['cycle'], event_observed=cleaned_data[~sm14]['status'], label='control')
kmf_control.plot_survival_function(ax=ax, show_censors=True, ci_legend=True, 
                                   label='sm14 not greater than ' + str(threshold)
                                  )


from lifelines.plotting import add_at_risk_counts
add_at_risk_counts(kmf_exp, kmf_control, ax=ax)

plt.title(f"Kaplan-Meier Estimator for each Category of Sensor Measurement 14")
plt.ylabel("est. probability of survival $\hat{S}(t)$")
plt.xlabel('cycles')
plt.show()

from sksurv.compare import compare_survival

cat.append('sensor_measurement_14')

p_value_list = []
for column in cat:
    p_value = compare_survival(data_y, turbofan_df[column])[1]
    p_value_list.append(p_value)
    
result = pd.DataFrame({'columns': cat, 'p-value': p_value_list}).set_index('columns')
result['conclusion'] = result['p-value'].apply(lambda x: "significant" if x < 0.05 else "not significant")
result

cleaned_data

cleaned_data[sm14]['cycle'].count()

cleaned_data[~sm14]['cycle'].count()

# Import logrank_test
from lifelines.statistics import logrank_test

sm16 = (cleaned_data['sensor_measurement_16'] == 0.02)

# Run log-rank test to compare groups by survival function
mask_dict = {'sensor_measurement_14': sm14, 'sensor_measurement_16': sm16, 'operational_setting_3': fullos3}
p_value_list = []

for item in mask_dict.items():
    lrt = logrank_test(durations_A = cleaned_data[item[1]]['cycle'], 
                       durations_B = cleaned_data[~item[1]]['cycle'], 
                       event_observed_A = cleaned_data[item[1]]['status'], 
                       event_observed_B = cleaned_data[~item[1]]['status'])
    p_value_list.append(lrt.p_value)

test_results = pd.DataFrame({'feature': mask_dict.keys(), 'p-value': p_value_list}).set_index('feature')
test_results['conclusion'] = test_results['p-value'].apply(lambda x: "significant" if x < 0.05 else "not significant")
test_results

# plt.xlabel callable bugfix
import matplotlib.pyplot as plt
from importlib import reload
plt=reload(plt)

# Import WeibullFitter class
from lifelines import WeibullFitter

# Instantiate WeibullFitter
wb = WeibullFitter()

# Fit data to wb
wb.fit(data_y["cycle"], data_y["status"])

# Plot survival function
wb.survival_function_.plot()
plt.title(f"time-to-failure with Weibull model")
plt.ylabel("est. probability of survival $\hat{S}(t)$")
plt.xlabel('cycles')

plt.show()

threshold = 8100

# Mask for operational_setting_3
sm14 = (cleaned_data['sensor_measurement_14'] > threshold)

# Fit to sm14 greater than 8140
wb.fit(#durations=prison[parole]['week'], event_observed=prison[parole]['arrest']
      durations=cleaned_data[sm14]['cycle'], event_observed=cleaned_data[sm14]['status']
       #, label='exp'
      )

# Print rho
print("The rho_ parameter of group sm14 greater than 8140: ", wb.rho_)

# Fit to engines sm14 lower equal than 8140
wb.fit(durations=cleaned_data[~sm14]['cycle'], event_observed=cleaned_data[~sm14]['status'])

# Print rho
print("The rho_ parameter of group sm14 lower equal than 8140 survival function is: ", wb.rho_)

# Import WeibullAFTFitter and instantiate
from lifelines import WeibullAFTFitter
aft = WeibullAFTFitter()

# Fit heart_patients data into aft
aft.fit(cleaned_data, duration_col='cycle', event_col='status')

# Print the summary
print(aft.summary)

# Calculate the exponential of sensor measurement 11 coefficient
exp_sm11 = np.exp(aft.params_.loc['lambda_'].loc['sensor_measurement_11'])
print('When sensor_measurement_11 increases by 1, the average survival duration changes by a factor of', exp_sm11)

# Fit data to aft
aft.fit(cleaned_data, duration_col='cycle', event_col='status')

# Plot partial effects of prio
#aft.plot_partial_effects_on_outcome(  covariates='sensor_measurement_14',   values=np.arange(8090, 8170, 15))
aft.plot_partial_effects_on_outcome(  covariates='sensor_measurement_11',   values=np.arange(42, 45, 0.5))

plt.axvline(x=250,ls="--",c=".5")

plt.title(f"Weibull baseline model vs. varied values of sensor measurements 11")
plt.ylabel("survival probability")
plt.xlabel('cycles')
plt.show()

# Predict median of new data
aft_pred = aft.predict_median(cleaned_data)

# Print average predicted time to arrest
print("On average, the median number of cycles for an engine to fail is: ", np.mean(aft_pred))

from lifelines.fitters.exponential_fitter import ExponentialFitter 
from lifelines.fitters.log_normal_fitter import LogNormalFitter
from lifelines.fitters.log_logistic_fitter import LogLogisticFitter

# Instantiate each fitter
wb = WeibullFitter()
exp = ExponentialFitter()
lognorm = LogNormalFitter()
loglog = LogLogisticFitter()

# Fit to data
for model in [wb, exp, lognorm, loglog]:
      model.fit(durations=cleaned_data[sm14]['cycle'], event_observed=cleaned_data[sm14]['status'] )
      # Print AIC
      print(model.__class__.__name__, model.AIC_)

fig, axes = plt.subplots(2, 2, figsize=(10, 8))
axes = axes.reshape(4,)

for i, model in enumerate([WeibullFitter(), LogNormalFitter(), LogLogisticFitter(), ExponentialFitter()]):
    model.fit(durations=cleaned_data[sm14]['cycle'], event_observed=cleaned_data[sm14]['status'] )
    qq_plot(model, ax=axes[i])

plt.suptitle(f"QQ Plot of Fitted models of different distributions of those engines with Sensor Measurement 14 above threshold {threshold}" )
plt.show()

data_x = data_x[data_x.columns[~data_x.columns.isin(["sensor_measurement_16=0.03"])]]

from sksurv.linear_model import CoxPHSurvivalAnalysis

# Define model
estimator = CoxPHSurvivalAnalysis(verbose = 10)

# Fitting model
estimator.fit(data_x, data_y)

cleaned_data = cleaned_data.drop('sensor_measurement_16', axis=1)

# Import CoxPHFitter class
from lifelines import CoxPHFitter

# Instantiate CoxPHFitter class cph
cph = CoxPHFitter(
    #penalizer=0.001
)

# Fit cph to data
cph.fit(df=cleaned_data, duration_col="cycle", event_col="status")

# Print model summary
print(cph.summary)

# Coef model
pd.Series(estimator.coef_, index = data_x.columns)

# Assign summary to summary_df
summary_df = cph.summary

# Create new column of survival time ratios
summary_df["surv_ratio"] = 1 / summary_df['exp(coef)']

summary_df

# Print surv_ratio for prio
survivalratio = summary_df.loc['sensor_measurement_04', "surv_ratio"]
survivalratio

TimeToFailureChange = survivalratio - 1
TimeToFailureChange

print(f'Time-to-failure changes by {TimeToFailureChange * 100:.2f} %')

covars = list(summary_df["surv_ratio"].index)
covars

sorted(covars)

# Plot partial effects
cph.plot_partial_effects_on_outcome(covariates=sorted(covars),
         values=[[100.0,1140,42,8080], 
                 [100.0,1130,43,8100], 
                 [100.0,1150,45,8120],
                 [60,   1150,40,8080]])

# Show plot
plt.title("Indiviual survival curves of four prototype Turbofan Jet engines with specific covariate properties - Time to Failure Estimation" )
plt.ylabel("est. probability of survival")
#plt.xlabel("cycles")
plt.legend(loc="upper right")
plt.show()

# Check PH assumption
print(cph.check_assumptions(training_df=cleaned_data, p_value_threshold=0.05, show_plots=True))

cleaned_strata_data = cleaned_data.copy()
cleaned_data.sensor_measurement_14.describe()

cleaned_strata_data

cph2 = CoxPHFitter()
cph2.fit(df=cleaned_strata_data, duration_col="cycle", event_col="status", 
         #strata=['sm14_strata'], 
         formula="bs(sensor_measurement_14, df=3, lower_bound=7850, upper_bound=8300) + sensor_measurement_04 + sensor_measurement_11 + operational_setting_3"
        )

cph2.print_summary(3)
sns.set(rc = {'figure.figsize':(7,7)})
cph2.plot()
plt.show()

import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)

cph2.check_assumptions(cleaned_strata_data)

# Select engines that have not failed yet
engines_still_working = cleaned_strata_data.loc[cleaned_strata_data['status'] == False].copy()

engines_still_working.loc[:, 'age_cycles'] = 200
engines_still_working

# Existing durations of employees that have not churned
engines_still_working_last_obs = engines_still_working['age_cycles']

# Predict survival function conditional on existing durations
cph2.predict_survival_function(engines_still_working, 
                              conditional_after=engines_still_working_last_obs)

# Predict median remaining times for current employees
pred = cph2.predict_median(engines_still_working, conditional_after=engines_still_working_last_obs)

# Print the smallest median remaining time
print(pred)

numbers =pred != np.inf 
numbers.sum()

#train_cols = index_names + remaining_sensors + ['start', 'breakdown']
#predict_cols = ['time_cycles'] + remaining_sensors + ['start', 'breakdown']  # breakdown value will be 0

timevary_df = timeseries_df.copy()
timevary_df_test = timeseries_df_test.copy()
timevary_df.head()

timevary_df.groupby(by="machine_name")['cycle'].max()

# RUL
RUL_timevary_df = timevary_df.merge(timevary_df.groupby(by="machine_name")['cycle'].max().to_frame(name='fail_cycle'), left_on='machine_name', right_index=True)
RUL_timevary_df_test = timevary_df_test.merge(timevary_df_test.groupby(by="machine_name")['cycle'].max().to_frame(name='fail_cycle'), left_on='machine_name', right_index=True)
RUL_timevary_df 

RUL_timevary_df['start'] = RUL_timevary_df['cycle'] - 1
RUL_timevary_df_test['start'] = RUL_timevary_df_test['cycle'] - 1

#timevary_df['status'] = timevary_df['cycle'].apply(lambda x: False if x > timevary_df['fail_cycle'] else True)
RUL_timevary_df['status'] = np.where(RUL_timevary_df['cycle']>=RUL_timevary_df ['fail_cycle'], True, False)
RUL_timevary_df_test['status'] = np.where(RUL_timevary_df_test['cycle']>=RUL_timevary_df_test['fail_cycle'], True, False)

# Calculate remaining useful life for each row
RUL_timevary_df["RUL"] = RUL_timevary_df["fail_cycle"] - RUL_timevary_df["cycle"]
RUL_timevary_df_test["RUL"] = RUL_timevary_df_test["fail_cycle"] - RUL_timevary_df_test["cycle"]
RUL_timevary_df_test

from lifelines import CoxTimeVaryingFitter

#RUL_timevary_df_training = RUL_timevary_df.drop(['fail_cycle', 'RUL'], axis=1)
# drop additional 
RUL_timevary_df_training = RUL_timevary_df[['machine_name', 'cycle', 'operational_setting_1',
       'operational_setting_2', 'operational_setting_3',
        'sensor_measurement_02',
       'sensor_measurement_03', 'sensor_measurement_04',
      
       'sensor_measurement_07', 'sensor_measurement_08',
       'sensor_measurement_09', 'sensor_measurement_10',
       'sensor_measurement_11', 'sensor_measurement_12',
       'sensor_measurement_13', 'sensor_measurement_14',
       'sensor_measurement_15', 
       'sensor_measurement_17',  'sensor_measurement_20',
       'sensor_measurement_21', 'start', 'status']]

ctv = CoxTimeVaryingFitter()
ctv.fit(RUL_timevary_df_training, id_col="machine_name", event_col='status', 
        start_col='start', stop_col='cycle', show_progress=True, step_size=1)

ctv.print_summary()
plt.figure(figsize=(10,5))
ctv.plot()
plt.show()
plt.close()

RUL_timevary_df_rcensored = RUL_timevary_df[['machine_name', 'cycle', 'operational_setting_1',
       'operational_setting_2', 'operational_setting_3',
        'sensor_measurement_02',
       'sensor_measurement_03', 'sensor_measurement_04',
      
       'sensor_measurement_07', 'sensor_measurement_08',
       'sensor_measurement_09', 'sensor_measurement_10',
       'sensor_measurement_11', 'sensor_measurement_12',
       'sensor_measurement_13', 'sensor_measurement_14',
       'sensor_measurement_15', 
       'sensor_measurement_17',  'sensor_measurement_20',
       'sensor_measurement_21', 'start', 'status', "RUL"]].copy()[RUL_timevary_df['cycle'] <= 215]

RUL_timevary_df_rcensored

preds_df = RUL_timevary_df_rcensored.groupby("machine_name").last()
preds_df = preds_df[preds_df.status == False].copy()

preds_df_wo_RUL = preds_df.drop(['RUL'], axis=1)

predictions = pd.DataFrame(ctv.predict_log_partial_hazard(preds_df_wo_RUL), columns=['predictions'])
predictions.index = preds_df.index

predictions['RUL'] = preds_df['RUL']
predictions.head(10)

# plt.xlabel callable bugfix
import matplotlib.pyplot as plt
from importlib import reload
plt=reload(plt)

plt.figure(figsize=(15,5))
plt.plot(predictions['RUL'], predictions['predictions'], '.b')
xlim = plt.gca().get_xlim()
plt.xlim(xlim[1], xlim[0])
plt.xlabel('RUL')
plt.ylabel('log_partial_hazard')
plt.title('log_partial_hazard vs. Remaining useful Lifetime of all engines still working after cycle 215')
plt.show()

partial_hazard_df = RUL_timevary_df[['machine_name', 'cycle', 
       'operational_setting_2', 'operational_setting_3',
        'sensor_measurement_02',
       'sensor_measurement_03', 'sensor_measurement_04',
       'sensor_measurement_07', 'sensor_measurement_08',
       'sensor_measurement_09', 'sensor_measurement_10',
       'sensor_measurement_11', 'sensor_measurement_12',
       'sensor_measurement_13', 'sensor_measurement_14',
       'sensor_measurement_15', 
       'sensor_measurement_17',  'sensor_measurement_20',
       'sensor_measurement_21',  'start', 'status', "RUL"]].copy()

partial_hazard_df_test = RUL_timevary_df_test[['machine_name', 'cycle', 
       'operational_setting_2', 'operational_setting_3',
        'sensor_measurement_02',
       'sensor_measurement_03', 'sensor_measurement_04',
       'sensor_measurement_07', 'sensor_measurement_08',
       'sensor_measurement_09', 'sensor_measurement_10',
       'sensor_measurement_11', 'sensor_measurement_12',
       'sensor_measurement_13', 'sensor_measurement_14',
       'sensor_measurement_15', 
       'sensor_measurement_17',  'sensor_measurement_20',
       'sensor_measurement_21',  'start', 'status', "RUL"]].copy()

partial_hazard_df_wo_RUL =  partial_hazard_df.drop(['RUL'], axis=1)
partial_hazard_df_test_wo_RUL =  partial_hazard_df_test.drop(['RUL'], axis=1)
partial_hazard_df_wo_RUL.info()

ctv_all = CoxTimeVaryingFitter()
ctv_all.fit(partial_hazard_df_wo_RUL, id_col="machine_name", event_col='status', 
        start_col='start', stop_col='cycle', show_progress=True, step_size=1)

partial_hazard_df['log_partial_hazard'] = ctv_all.predict_log_partial_hazard(partial_hazard_df_wo_RUL)
partial_hazard_df_test['log_partial_hazard'] = ctv_all.predict_log_partial_hazard(partial_hazard_df_test_wo_RUL)
partial_hazard_df[-220:-1]

rul10ph = partial_hazard_df[partial_hazard_df.RUL < 11][['RUL', 'log_partial_hazard']]
rul10ph

# plt.xlabel callable bugfix
import matplotlib.pyplot as plt
from importlib import reload
plt=reload(plt)

sns.set(rc = {'figure.figsize':(16,6)})
plt.scatter(partial_hazard_df['RUL'], partial_hazard_df['log_partial_hazard'])

plt.ylabel('log_partial_hazard')
plt.xlabel('RUL')
plt.axhline(y=5,ls="--",c=".5")
xlim = plt.gca().get_xlim()
plt.xlim(xlim[1], xlim[0])
plt.title('Train data - RUL vs. log partial hazard')

plt.show()

sns.set(rc = {'figure.figsize':(16,9)})
g = sns.FacetGrid(rul10ph, col="RUL", 
                  #height=12, 
                  #aspect=1.2, 
                  col_wrap=4)
g.map(sns.histplot, "log_partial_hazard")
g.map(plt.axvline,x=10,ls="--",c=".5", color='r',)
plt.show()

%%time
sns.set(rc = {'figure.figsize':(8,5)})

def loghaz_stopheuristic(df, threshold, recurrence_limit, method='consecutive', window=11):
    machinelist = df.machine_name.unique()
    decision_df = pd.DataFrame(index=machinelist)
    df['stop_decision'] = np.nan
    df['stop_decision_when'] = np.nan
    for m in machinelist:
        mask = (df.machine_name == m)
        temp = df[mask & (df.RUL < (window + 1))]
        overcount = 0
        for index, row in enumerate(temp.itertuples(), 1):
            if row.log_partial_hazard > threshold:
                 overcount += 1
            elif method=='consecutive':
                 overcount = 0
            if overcount > recurrence_limit:
                if row.RUL != 0:
                    decision_df.loc[m , 'stop_decision'] = True
                decision_df.loc[m , 'stop_decision_when'] = row.RUL
                #print(f'machine {m} was stopped' )
                break
            
    return df, decision_df
            

result, decisions = loghaz_stopheuristic(partial_hazard_df, threshold=5.5, recurrence_limit=5, method='consecutive', window=16)

print(f'stopped: {decisions.sum()} total: {decisions.shape[0]}')
      
print(f'\nmedian rul: {decisions.stop_decision_when.median()}')
decisions.stop_decision_when.plot(kind='box')
plt.ylabel('RUL')
plt.title('Distribution of engines retired with remaining RUL before failure (RUL = 0)')
plt.show()
decisions.stop_decision_when.plot(kind='hist', bins=35)
#plt.xlabel('RUL')
plt.title('Histogramm of engines retired with remaining RUL')
plt.show()

r, decisions_test = loghaz_stopheuristic(partial_hazard_df_test, threshold=5, recurrence_limit=3, method='consecutive', window=30)

print(f'stopped: {decisions_test.sum()} total: {decisions_test.shape[0]} median rul: {decisions_test.stop_decision_when.median()}')
decisions_test.stop_decision_when.plot(kind='box')
plt.ylabel('RUL')
plt.title('Distribution of engines retired with remaining RUL before failure (RUL = 0)')
plt.show()
decisions_test.stop_decision_when.plot(kind='hist', bins=35)
#plt.xlabel('RUL')
plt.title('Histogramm of engines retired with remaining RUL')
plt.show()