Predicting Lung Cancer Survival

Importing the Packages

import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sb
import re # regular expressions, using this to trim the column names 
import csv
import numpy as np
from lifelines import CoxPHFitter

Understanding and Visualizing the training data.

df = pd.read_csv('x train(Q0)/clinical_data.csv')
df.head(20)

Variables in the radiomics Data

Below is the description of the variables in radiomics data

Data Wrangling train data

# Read in the radiomics CSV file
df2 = pd.read_csv('x train(Q0)/radiomics.csv', header=[0, 1, 2])

# Concatenate the three rows of the header
header = ['_'.join(col).strip() for col in df2.columns.values]

# Update the column names with the concatenated header
df2.columns = header

# Output the modified CSV file
df2.to_csv('x train(Q0)/radiomics_concatenated.csv', index=False)

df2 = df2.rename(columns={df2.columns[0]: 'PatientID'})

regex = r"_Unnamed: \d+_level_2"

# Replace the matched string with an empty string in all column headers
fixed_columns = [re.sub(regex, "", col) for col in df2.columns]

# Update the column headers of the dataframe
df2.columns = fixed_columns

# Output the modified CSV file
df2.to_csv('x train(Q0)/radiomics_concatenated.csv', index=False)

df2.head(10)

full_df = pd.merge(df, df2, on=["PatientID"])
print(full_df.shape[0])
full_df.head(20)

Cleaning the train data

full_df = full_df.dropna()
full_df.drop(full_df[full_df['Tstage'] > 4].index, inplace = True)
full_df.drop(full_df[full_df['Mstage'] > 1].index, inplace = True)
full_df.drop(full_df[full_df['Nstage'] > 3].index, inplace = True)
print("Count after dropping missing/incorrect values: " + str(full_df.shape[0]))

no_PatientID = full_df.drop('PatientID', axis=1)
no_PatientID.describe()

# check which values contain "NOS" (case-insensitive)
mask = full_df['Histology'].str.contains('NOS', case=False, na=False)

# replace values with "nos"
full_df.loc[mask, 'Histology'] = 'nos'

mask = full_df['Histology'].str.contains('Squamous cell carcinoma', case=False, na=False)
full_df.loc[mask, 'Histology'] = 'squamous cell carcinoma'

mask = full_df['Histology'].str.contains('Adenocarcinoma', case=False, na=False)
full_df.loc[mask, 'Histology'] = 'adenocarcinoma'

histology_counts = full_df.groupby('Histology').count()

# plot the counts for each histology as a bar chart
plt.bar(histology_counts.index, histology_counts.iloc[:,0])

# set the plot title and axis labels
plt.title('Distribution of Histology')
plt.xlabel('Histology')
plt.ylabel('Count')

plt.xticks(rotation=90)

# show the plot
plt.show()

pairplot_cols = full_df[['Histology', 'Mstage', 'Nstage', 'Tstage', 'age', 'shape_original_shape_Compactness1', 'shape_original_shape_Compactness2', 'shape_original_shape_SurfaceVolumeRatio', 'textural_original_glcm_ClusterTendency']]
sb.pairplot(pairplot_cols, x_vars=['Mstage', 'Nstage', 'Tstage', 'age'], y_vars=['shape_original_shape_Compactness1', 'shape_original_shape_Compactness2', 'shape_original_shape_SurfaceVolumeRatio', 'textural_original_glcm_ClusterTendency'], hue='Histology')
plt.show()

sb.regplot(x=full_df[['age']], y=full_df[['shape_original_shape_Compactness1']], data=full_df);

sb.regplot(x=full_df[['age']], y=full_df[['shape_original_shape_Compactness2']], data=full_df);

sb.regplot(x=full_df[['age']], y=full_df[['textural_original_glcm_ClusterTendency']], data=full_df);

full_df = full_df.drop(['Histology','SourceDataset'], axis=1)
full_df.head(10)

full_df.to_csv('x train(Q0)/train_full_df.csv', index=False)

Random Forest Regression

from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from sklearn.feature_selection import RFE
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV
import random

# our two dataframes for X and y splits
y_csv = pd.read_csv('ytrain.csv')
merged_df = pd.merge(full_df, y_csv, on='PatientID')

# separating features
X = merged_df.drop('SurvivalTime', axis=1)
y = merged_df['SurvivalTime']

# splits the data into train and test sets

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

# feature scaling
scaler = MinMaxScaler()
scaler.fit(X_train)
scaler.transform(X_train)
scaler.transform(X_test)

param_grid = {
    'n_estimators': [100, 500, 1000],
    'max_depth': [60],
    'min_samples_split': [15],
    'min_samples_leaf': [5]
}

rfe = RandomForestRegressor()
rfe.fit(X_train, y_train)

importances = rfe.feature_importances_
indices = np.argsort(importances)[::-1]
print("Feature ranking:")
for f in range(X.shape[1]):
    print("%d. feature %d (%f)" % (f + 1, indices[f], importances[indices[f]]))

top_features = X.columns[indices][:20]

X_train = X_train[top_features]
X_test = X_test[top_features]

grid_search = GridSearchCV(estimator=rfe, param_grid=param_grid, cv=5, n_jobs=-1, scoring='neg_mean_squared_error')
grid_search.fit(X_train, y_train)

best_rfe = grid_search.best_estimator_
print(best_rfe)

y_pred = best_rfe.predict(X_test)
print('Predicted Values: ')
y_pred = np.maximum(y_pred, 0) # changing negative values 
np.set_printoptions(suppress=True, precision=0)
print(y_pred)

mse = mean_squared_error(y_test, y_pred)
print("Mean squared error:", mse)

r2_score(y_test, y_pred)
best_rfe.score(X_test, y_test)

Random Forest Scatterplot

y_test_small = y_test[:81]
y_pred_noneg_small = y_pred[:81]

fig, ax = plt.subplots()
ax.scatter(y_test_small, y_pred_noneg_small, c='b', label='Predicted Values')
ax.scatter(y_test_small, y_test_small, c='r', label='Actual Values')
ax.set_xlabel('Actual Survival Time')
ax.set_ylabel('Predicted Survival Time')
ax.set_title('Random Forest Regression Results')
ax.legend()
plt.show()

Linear Regression Model

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GridSearchCV
from sklearn.utils import resample

y_csv = pd.read_csv('ytrain.csv')
merged_df = pd.merge(full_df, y_csv, on='PatientID')

X = merged_df.drop('SurvivalTime', axis=1)
y = merged_df['SurvivalTime']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

# feature scaling
scaler = StandardScaler()
scaler.fit(X_train)
scaler.transform(X_train)
scaler.transform(X_test)

# model
model = LinearRegression()
model.fit(X_train, y_train)

print('Predicted Values: ')
y_pred = model.predict(X_train)
y_pred = np.maximum(y_pred, 0)
np.set_printoptions(suppress=True, precision=0)
print(y_pred)

r2_score(y_train, y_pred)
model.score(X_test, y_test)

# mse = mean_squared_error(y_test, y_pred)
# print("Mean squared error:", mse)

import matplotlib.pyplot as plt

plt.scatter(y_train, y_pred, color='blue')
plt.plot(y_train, y_train, color='red', linewidth=2)
plt.xlabel('Actual Values')
plt.ylabel('Predicted Values')
plt.title('Linear Regression Model')
plt.show()

Kaplan Meier

import pandas as pd
from lifelines import KaplanMeierFitter
from lifelines.statistics import logrank_test
from sklearn.model_selection import train_test_split


# Load the clinical data, features data, and survival data into separate dataframes

survival_data_train = pd.read_csv('ytrain.csv')
train_full_data = pd.read_csv('x train(Q0)/train_full_df.csv')


# Merge the dataframes based on a common identifier column
merged_data = pd.merge(train_full_data, survival_data_train, on='PatientID')
train_data1, train_data2 = train_test_split(merged_data, test_size=0.5, random_state=40)
#print(merged_data.dtypes)
#print(merged_data)
#print(train_data1)
#print(train_data2)

# Extract the relevant columns for analysis
time_column = 'SurvivalTime'  # column containing the time-to-event data
event_column = 'Event'  # column containing the event indicator data


# columns containing the covariate data (excluding string columns)
#feature_columns = [col for col in merged_data.columns if col not in [time_column, event_column] and merged_data[col].dtype != 'O']  
#test_feature_columns = [col for col in test_merged_data.columns if col not in [time_column, event_column] and test_merged_data[col].dtype != 'O']  

# Create a Kaplan-Meier estimator object and fit the survival data
kmf1 = KaplanMeierFitter()
kmf2 = KaplanMeierFitter()

kmf1.fit(train_data1[time_column], event_observed=train_data1[event_column])
kmf2.fit(train_data2[time_column], event_observed=train_data2[event_column])


# Plot the survival curve
kmf1.plot(color='blue', label='Dataset 1')
kmf2.plot(color='red', label='Dataset 2')

confidence_interval1 = kmf1.confidence_interval_
survival_function1 = kmf1.survival_function_
event_table1 = kmf1.event_table

confidence_interval2 = kmf2.confidence_interval_
survival_function2 = kmf2.survival_function_
event_table2 = kmf2.event_table

# # print confident Interval
print('Confidence Interval for dataset 1: ')
print(confidence_interval1)

print('Confidence Interval for dataset 2: ')
print(confidence_interval2)

# # print the survival function
print('Survival function for dataset 1:')
print(survival_function1)

print('Survival function for dataset 2:')
print(survival_function2)

# # print the event table
print('Event table for dataset 1 :')
print(event_table1)

print('Event table for dataset 2 :')
print(event_table2)

results = logrank_test(train_data1[time_column], train_data2[time_column],train_data1[event_column], train_data2[event_column])

# print the results
print('Log-rank test results:')
print(results.summary)
print('p-value: {:.3f}'.format(results.p_value))
print('test statistic: {:.3f}'.format(results.test_statistic))

Results