{"cells":[{"metadata":{},"cell_type":"markdown","source":"When I first saw the osic pulmonary fibrosis competition I thought it looked like quite an interesting competition for its mix of tabular and image data and trying to produce a forecast based on that. This provides scope to try all sorts of things such as feature engineering, image embedding models and time series models. To get started I have focused on the tabular data only which as it turns out can get us quite far before even touching the images. I've tried to stay away from spending too much time on model selection and instead opted for feature experimentation.\n\nBut first a quick overview (disclaimer I'm not a medical expert so apologies if I get any of the medical stuff wrong). Pulmonary Fibrosis is a lung disease that causes the lungs to decline over time. The rate of decline can range between very rapid and very slow. We can measure the diseases progress through [FVC](https://lunginstitute.com/blog/fev1-and-fvc/#:~:text=The%20forced%20vital%20capacity%20(FVC,the%20severity%20of%20the%20condition)) (Forced Vital Capacity) which involves the patient taking a deep breath and blowing as hard as they can into a tube. The amount of air they blow out is the FVC measured in ml. \n\nSo the challenge is when given a starting FVC for a patient to forecast the decline in the FVC value for all the upcoming weeks.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom sklearn import linear_model, ensemble\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error\n\nimport tensorflow as tf\n\nfrom tqdm.notebook import tqdm\n\nimport os\nfrom PIL import Image","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load data\n\nLet's begin by loading the data.","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\nsubmission = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.info()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So there's not a huge volume of tabular data. 1.5k of training examples with only 5 columns (weeks, percent, age, sex, smoking status) to construct features from.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Merge datasets\n\nI've learned recently that it is good practice to merge the train, validation and test sets at the start of a notebook. This ensures that exactly the same transformation is applied to every example. First though, I'll make sure there are no duplicates in the training dataset.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Then form the submission dataset. The test dataset needs expanding out across the 146 weeks per patient that the competition requires. This can be achieved by joining it to the sample submission.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"submission['Patient'] = (\n    submission['Patient_Week']\n    .apply(\n        lambda x:x.split('_')[0]\n    )\n)\n\nsubmission['Weeks'] = (\n    submission['Patient_Week']\n    .apply(\n        lambda x: int(x.split('_')[-1])\n    )\n)\n\nsubmission =  submission[['Patient','Weeks', 'Confidence','Patient_Week']]\n\nsubmission = submission.merge(test.drop('Weeks', axis=1), on=\"Patient\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Mark each example in each dataset with the name of the dataset they come from. This enables me to quickly split the dataset back up into the three component pieces at the end of the notebook.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train['Dataset'] = 'train'\ntest['Dataset'] = 'test'\nsubmission['Dataset'] = 'submission'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Merge the datasets into one and reset the index.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data = train.append([test, submission])\n\nall_data = all_data.reset_index()\nall_data = all_data.drop(columns=['index'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Quick data analysis\n\nI think it's worth having a look at how the FVC (label) value declines for a sample of patients in the training dataset. Let's pick the first five in the data and plot the decline of the FVC.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_patients = train.Patient.unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(5, 1, figsize=(10, 20))\n\nfor i in range(5):\n    patient_log = train[train['Patient'] == train_patients[i]]\n\n    ax[i].set_title(train_patients[i])\n    ax[i].plot(patient_log['Weeks'], patient_log['FVC'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So the decline is kinda linear as it does generally trend down over time. There are some spikes back up along the way which could cause a few issues. However this means that a simple linear model could have a good go at producing forecasts on this challenge as the main thing we need to do is predict the rate of decline for a patient. Like a trend line for these charts.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Feature Engineering\n\nThere is some good scope for engineering new features for this model. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### First FVC and First Week\nSome useful features might be the first FVC recorded per patient and the week it was recorded in","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data['FirstWeek'] = all_data['Weeks']\nall_data.loc[all_data.Dataset=='submission','FirstWeek'] = np.nan\nall_data['FirstWeek'] = all_data.groupby('Patient')['FirstWeek'].transform('min')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"first_fvc = (\n    all_data\n    .loc[all_data.Weeks == all_data.FirstWeek][['Patient','FVC']]\n    .rename({'FVC': 'FirstFVC'}, axis=1)\n    .groupby('Patient')\n    .first()\n    .reset_index()\n)\n\nall_data = all_data.merge(first_fvc, on='Patient', how='left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Weeks Passed\nThis feature measures how many weeks have passed since the patients first FVC reading.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data['WeeksPassed'] = all_data['Weeks'] - all_data['FirstWeek']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Patient height\n\nApparently the height of a patient is an important variable when predicting FVC. Maybe tall people have larger lungs and thus more air to exhale. Thanks to Srikanth Potukuchi whose [notebook](https://www.kaggle.com/srikanthpotukuchi/osic-random-forest-new-height-extracted) showed me how to estimate a patients height using their FVC.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def calculate_height(row):\n    if row['Sex'] == 'Male':\n        return row['FirstFVC'] / (27.63 - 0.112 * row['Age'])\n    else:\n        return row['FirstFVC'] / (21.78 - 0.101 * row['Age'])\n\nall_data['Height'] = all_data.apply(calculate_height, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Categorical columns\n\nThe sex and smoking status columns are categorical columns that need some transformation to turn them into numbers. Pandas get dummies makes this easy to achieve.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data = pd.concat([\n    all_data,\n    pd.get_dummies(all_data.Sex),\n    pd.get_dummies(all_data.SmokingStatus)\n], axis=1)\n\nall_data = all_data.drop(columns=['Sex', 'SmokingStatus'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Scale features\n\nNow scale all the features to get them onto the same range of numbers (0-1).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def scale_feature(series):\n    return (series - series.min()) / (series.max() - series.min())\n\nall_data['Weeks'] = scale_feature(all_data['Weeks'])\nall_data['Percent'] = scale_feature(all_data['Percent'])\nall_data['Age'] = scale_feature(all_data['Age'])\nall_data['FirstWeek'] = scale_feature(all_data['FirstWeek'])\nall_data['FirstFVC'] = scale_feature(all_data['FirstFVC'])\nall_data['WeeksPassed'] = scale_feature(all_data['WeeksPassed'])\nall_data['Height'] = scale_feature(all_data['Height'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Specify what columns will be used as features. This is for easy filtering of the datasets later.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"feature_columns = [\n    'Percent',\n    'Age',\n    'FirstWeek',\n    'FirstFVC',\n    'WeeksPassed',\n    'Height',\n    'Female',\n    'Male', \n    'Currently smokes',\n    'Ex-smoker',\n    'Never smoked',\n]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Split dataframe\n\nSplit the data back into the three dataframes they started as.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train = all_data.loc[all_data.Dataset == 'train']\ntest = all_data.loc[all_data.Dataset == 'test']\nsubmission = all_data.loc[all_data.Dataset == 'submission']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And take a look at the features that will be used to train the model.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train[feature_columns].head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model\n\nAs I mentioned at the start of the notebook I didn't intend to experiment too much with model selection in this notebook so I have stuck with a simple Linear Regressor from SKLearn. While these non-deep models don't usually produce as accurate results as deep models they are super fast to train and easy to evaluate. This is great for notebooks focussing primarily on features.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"model = linear_model.HuberRegressor(max_iter=200)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With the model loaded, insert the features and labels for training.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"model.fit(train[feature_columns], train['FVC'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Make predictions","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"predictions = model.predict(train[feature_columns])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Evaluate\n\nLet's begin by having a look at the models weights. This gives us a good indication of what features are driving the models predictions.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.bar(train[feature_columns].columns.values, model.coef_)\nplt.xticks(rotation=90)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"While mean squared error isn't the competition metric it is a simple loss metric to help understand how close the models predictions are to the actual labels. The limitation of this error number though is that it can't be too close to zero as that would indicate over-fitting a model that should only be producing a trend line.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"mse = mean_squared_error(\n    train['FVC'],\n    predictions,\n    squared=False\n)\n\nmae = mean_absolute_error(\n    train['FVC'],\n    predictions\n)\n\nprint('MSE Loss: {0:.2f}'.format(mse))\nprint('MAE Loss: {0:.2f}'.format(mae))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Found code for competition metric [here](https://www.kaggle.com/titericz/tabular-simple-eda-linear-model#Calculate-competition-metric)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def competition_metric(trueFVC, predFVC, predSTD):\n    clipSTD = np.clip(predSTD, 70 , 9e9)  \n    deltaFVC = np.clip(np.abs(trueFVC - predFVC), 0 , 1000)  \n    return np.mean(-1 * (np.sqrt(2) * deltaFVC / clipSTD) - np.log(np.sqrt(2) * clipSTD))\n    \n\nprint(\n    'Competition metric: ', \n    competition_metric(train['FVC'].values, predictions, 285) \n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's also include a scatterplot and histogram to see an overview of how close the predictions are to the labels.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train['prediction'] = predictions","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(predictions, train['FVC'])\n\nplt.xlabel('predictions')\nplt.ylabel('FVC (labels)')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"delta = predictions - train['FVC']\nplt.hist(delta, bins=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally take the first five patients as a sample and compare the true FVC readings against the models predicted FVC readings.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(5, 1, figsize=(10, 20))\n\nfor i in range(5):\n    patient_log = train[train['Patient'] == train_patients[i]]\n\n    ax[i].set_title(train_patients[i])\n    ax[i].plot(patient_log['WeeksPassed'], patient_log['FVC'], label='truth')\n    ax[i].plot(patient_log['WeeksPassed'], patient_log['prediction'], label='prediction')\n    ax[i].legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Submission\n\nWith the model trained we can use it to make predictions. As the submission dataframe went through the same transformation as the train dataset it is ready for the model to infer on.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"submission[feature_columns].head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Add the predictions to the submission dataframe.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sub_predictions = model.predict(submission[feature_columns])\nsubmission['FVC'] = sub_predictions","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plot the forecasted FVC for five patients in the test set.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"test_patients = list(submission.Patient.unique())\nfig, ax = plt.subplots(5, 1, figsize=(10, 20))\n\nfor i in range(5):\n    patient_log = submission[submission['Patient'] == test_patients[i]]\n\n    ax[i].set_title(test_patients[i])\n    ax[i].plot(patient_log['WeeksPassed'], patient_log['FVC'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Drop the surplus columns not needed for submission and add a dummy confidence column.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = submission[['Patient_Week', 'FVC']]\n\nsubmission['Confidence'] = 285","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And write the submission to file.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}