{"cells":[{"metadata":{"_uuid":"2b7cba33-2a68-47ef-9b01-888c16af8c90","_cell_guid":"08cbc198-09e4-4b39-9417-11753ce548b1","trusted":true},"cell_type":"code","source":"#!pip install torch==1.6.0+cu101 torchvision==0.7.0+cu101 -f https://download.pytorch.org/whl/torch_stable.html\n#!pip install fastai==2.0.9\n\n!pip install /kaggle/input/fast-v2-offline/torch-1.6.0-cp37-cp37m-manylinux1_x86_64.whl -q\n!pip install /kaggle/input/fast-v2-offline/torchvision-0.7.0-cp37-cp37m-manylinux1_x86_64.whl -q\n!pip install /kaggle/input/fast-v2-offline/dataclasses-0.6-py3-none-any.whl -q\n!pip install /kaggle/input/fast-v2-offline/fastprogress-1.0.0-py3-none-any.whl -q\n!pip install /kaggle/input/fast-v2-offline/fastcore-1.0.1-py3-none-any.whl -q\n!pip install /kaggle/input/fast-v2-offline/fastai-2.0.8-py3-none-any.whl -q","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8cbe6a26-7921-4be7-bbf7-3a6f637deb74","_cell_guid":"0a2c2524-1e17-4505-98cf-1d40e7f61d01","trusted":true},"cell_type":"code","source":"#pip install fastai --upgrade","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"51d8e7d5-1c94-4941-b8de-9814cbbcc0d0","_cell_guid":"02a4f583-404e-447f-b0e2-d7c98d4b73c6","trusted":true},"cell_type":"code","source":"import os, random\nimport numpy as np\nimport pandas as pd\nfrom pandas_profiling import ProfileReport\n\nfrom sklearn.preprocessing import OneHotEncoder\nfrom sklearn import metrics\nfrom sklearn import ensemble\nfrom sklearn.ensemble import RandomForestRegressor\n\nfrom fastai.tabular.all import *\nfrom fastai.medical.imaging import *\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"32325aaf-f75a-4936-b66a-c8ea2a6e032b","_cell_guid":"d16292b0-6f8d-4b02-a6be-20855c4cc453","trusted":true},"cell_type":"code","source":"def seed_everything(seed): \n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\nseed_everything(123)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6d45153f-c28d-43ef-9518-2535cdc3a75a","_cell_guid":"12ea0299-915e-46e0-87c3-f049c72214de","trusted":true},"cell_type":"code","source":"path = Path(\"../input/osic-pulmonary-fibrosis-progression\")\npath.ls()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2de79ac8-acd9-40ef-b7ec-59006545c9de","_cell_guid":"abdcfa25-abef-4564-a4e3-368f7a6dbe18","trusted":true},"cell_type":"markdown","source":"# What Are We Working To Predict?\nWe will be predicting, given an initial FVC measurement and the CT scan images for a patient, what the final 3 FVC measurements are for the patient. Since that initial FVC could be any week, we must create predictions for Weeks of value [-12, 133]. Our final score for this model will be the accuracy of the final 3 week predictions given by the Laplace Log Likelihood.\n\nLet's explore our data first."},{"metadata":{"_uuid":"522d7f85-cc63-42e9-81b3-b9d38a8d11b7","_cell_guid":"b1f416cb-ac33-499e-a079-2b0d5f1fa603","trusted":true},"cell_type":"code","source":"train_df = pd.read_csv(path/'train.csv')\ntest_df = pd.read_csv(path/'test.csv')\nsubmission_df = pd.read_csv(path/'sample_submission.csv')\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.info()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ed5037ac-99f7-4300-affe-4c8fb437c0eb","_cell_guid":"a77cba5f-97e0-4720-b299-da3eb9eedb6b","trusted":true},"cell_type":"code","source":"profile = ProfileReport(train_df, title = \"Profiling Report\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3c2e1745-2c86-4881-a90e-f702ea1a6233","_cell_guid":"41b3ad69-4536-4c8f-83bd-23f495f3cb4f","trusted":true},"cell_type":"code","source":"profile.to_notebook_iframe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"208ca187-c43a-455f-8117-8defacfc1e04","_cell_guid":"bb9bcf75-146f-44a7-b21d-ec6c177a957b","trusted":true},"cell_type":"code","source":"train_df.describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5e5edbf0-ab75-4903-a762-e9bfd0aa1ef6","_cell_guid":"f3832e22-3260-46bd-b4ef-8a268b412a8e","trusted":true},"cell_type":"markdown","source":"# Data Cleaning\n## Duplicates\n\nFirst, let's identify if we have duplicates in our training set and correct for them, before anything else:"},{"metadata":{"_uuid":"82c296d7-d4a0-43e7-906d-e67aed0e409a","_cell_guid":"06b43ded-3b60-4c2c-b6a9-560fd09f77ac","trusted":true},"cell_type":"code","source":"duplicates = train_df[train_df.duplicated(subset = ['Patient', 'Weeks'], keep = False)]\nduplicates","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c8fe60b5-ebd0-4ccb-b30a-ce08c21e1936","_cell_guid":"8c0344a7-24fe-47a0-a442-3457cc3d5554","trusted":true},"cell_type":"markdown","source":"Since we don't have much data to start with, let's average the FVC and Percent values and keep the rest intact."},{"metadata":{"_uuid":"f4371d09-632f-42b9-b61b-60c440b28012","_cell_guid":"a48b295b-1401-4af4-8537-2dec84a70a23","trusted":true},"cell_type":"code","source":"train_df = train_df.groupby(['Patient', 'Weeks']).agg({\n    'FVC': 'mean',\n    'Percent' : 'mean',\n    'Age' : 'mean',\n    'Sex': lambda x: pd.unique(x), \n    'SmokingStatus' : lambda y: pd.unique(y)}).reset_index() #give us back train if test is in train\ntrain_df","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2cf4eeeb-0a66-4d95-ae95-38f064c4b6f8","_cell_guid":"0c738177-2e9e-4a16-a51b-0e5db8b887b2","trusted":true},"cell_type":"code","source":"duplicates = train_df[train_df.duplicated(subset = ['Patient', 'Weeks'], keep = False)]\nduplicates #none please!","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Outliers\nThrough data exploration (and evaluating largest losses in our models), we can identify some data points that don't look particularly correct. By removing these, our hope is our model is more robust."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df[train_df['Patient'] == 'ID00298637202280361773446']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = train_df[train_df['Patient'] != 'ID00298637202280361773446']\ntrain_df[train_df['Patient'] == 'ID00298637202280361773446']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"outliers = {'ID00093637202205278167493' : 16,\n           'ID00207637202252526380974' : 33,\n           'ID00061637202188184085559' : 24,\n           'ID00319637202283897208687' : 16,\n           'ID00355637202295106567614' : 28}\ntrain_df[train_df['Patient'] == 'ID00093637202205278167493']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for key, value in outliers.items():\n    train_df = train_df[(train_df['Patient'] != key) | (train_df['Weeks'] != value)]\ntrain_df[train_df['Patient'] == 'ID00093637202205278167493'] #Week 16 is gone","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"22e1b1d4-3a80-444e-b13c-2513c0513d9a","_cell_guid":"de8bb067-f3d4-45ea-9109-7799de487524","trusted":true},"cell_type":"markdown","source":"# Initial EDA\n## Groups\nWe initial intuition with this problem is to see if we can identify groups of individuals that may trend differently according to the tabular data provided. Let's do some initial EDA to preliminarily identify any trends for FVC over time in regards to sex, age, and smoking status. We will analyze both FVC and FVC Percent over time.\n\nQuestions to answer:\n\nHow do FVC and FVC Percent change over time? And how does this differ for each major group?"},{"metadata":{"_uuid":"c080d760-6c9c-41aa-88bd-966a38d92bb3","_cell_guid":"afd9f3f3-b330-474b-a92f-49356a34af1b","trusted":true},"cell_type":"code","source":"sns.lmplot(x=\"Weeks\", y=\"FVC\", hue=\"Sex\", data = train_df).fig.set_size_inches(10,7)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9693ad71-741d-46d2-ab74-d4f1a239eefa","_cell_guid":"18759a9e-0cf6-4cf9-9641-754cf10eb575","trusted":true},"cell_type":"markdown","source":"Males are much higher FVC - however the general trend actually shows Males decrease over time while Females increase over time. Note there are many more Males than Females. What if instead of FVC, we looked at Percent?"},{"metadata":{"_uuid":"18e51748-2701-4438-9ee7-aceaeea29bb2","_cell_guid":"4a6ac473-08f8-445b-bd69-53ef6d16323d","trusted":true},"cell_type":"code","source":"sns.lmplot(x=\"Weeks\", y=\"Percent\", hue=\"Sex\", data = train_df).fig.set_size_inches(10,7)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5863d210-1c2a-446b-add9-2926c61dc90a","_cell_guid":"d9fa7d23-66b1-4c1c-a9c2-e3ef277d9f5c","trusted":true},"cell_type":"markdown","source":"When we use Percent as opposed to FVC, the male vs female absolute differential goes away. This is expected as FVC Percent is a representation of the % of FVC for a 'typical' patient, given some unknown demographic and health factors (believe it is sex, height, weight?). However, the small trend of males decreasing more than females remains.\n\nIt may be useful to look at how a few individual patients trend over time."},{"metadata":{"_uuid":"2d83ebe0-c633-4be7-95c0-cc01f25518b2","_cell_guid":"971a41ec-f2a9-4da0-adf3-ae1ad8611465","trusted":true},"cell_type":"code","source":"sns.lmplot(x=\"Weeks\", y=\"Percent\", hue=\"Patient\", data = train_df.head(100)).fig.set_size_inches(15,7)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fb1ec981-dd9a-4b10-a72d-dcabc8bf5083","_cell_guid":"b29dc47d-69b2-49d8-828c-0290168171d6","trusted":true},"cell_type":"code","source":"sns.lmplot(x=\"Weeks\", y=\"Percent\", hue=\"Patient\", data = train_df.tail(95)).fig.set_size_inches(15,7)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We're seeing general declines on the whole, but nothing huge. What about Age?"},{"metadata":{"_uuid":"98583c4c-08fb-477d-ab10-4e5de5d6437f","_cell_guid":"e9eb249c-87d4-468f-834b-f16d16ab8b5a","trusted":true},"cell_type":"code","source":"sns.lmplot(x=\"Weeks\", y=\"Percent\", hue=\"Age\", data = train_df).fig.set_size_inches(15,12)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b8048b52-8bde-4b75-a955-647e584ad941","_cell_guid":"8134de29-03f2-428f-8997-b0d41b4aaf11","trusted":true},"cell_type":"markdown","source":"Let's inspect a \"bottom-level\" group, and ascertain if they are similar to one another."},{"metadata":{"_uuid":"3e5adbb1-e6cb-4e7a-87b7-acb40ecd5e5f","_cell_guid":"f4e6dfb6-cf3c-49c6-83f2-557d897e7c81","trusted":true},"cell_type":"code","source":"tr_Group = train_df[(train_df['Age'] == 66) & (train_df['Sex'] == \"Male\") & (train_df['SmokingStatus'] == \"Ex-smoker\")]\nsns.lmplot(x=\"Weeks\", y=\"Percent\", hue=\"Patient\", data = tr_Group).fig.set_size_inches(15,7)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"44c455dc-29c3-470c-a6f8-049cd46cc213","_cell_guid":"81f7b9b9-1890-4b4b-abc9-fda043528d8f","trusted":true},"cell_type":"markdown","source":"Even within our (small) example group where Age, Sex, and Smoking Status are the same, we have wildly different Percent intercepts for Weeks_Adj between 0 and 10, from 90s to high 40s. Therefore the intercept of the initial FVC will be extremely important. What we really hope to learn is how that may change as the disease progresses given the data we have."},{"metadata":{},"cell_type":"markdown","source":"# Data Wrangling and Transformations\n## Combining Datasets\nFirst we'll combine the train and test sets such that our transformations are applied the same to both train and test data."},{"metadata":{"_uuid":"ea143415-bc29-44b5-aea0-88a5a74b72ab","_cell_guid":"3d68da71-425c-4c8f-a927-182a449c2745","trusted":true},"cell_type":"code","source":"submission_df[['Patient','Weeks']] = submission_df['Patient_Week'].str.split(\"_\",expand = True) #split Patient_Week\nsubmission_df = submission_df.drop('FVC', axis=1)\nsubmission_df = submission_df.merge(test_df.drop(['Weeks'], axis = 1), on = \"Patient\")\n\n#introduce a column to indicate the source (train/test) for the data\ntrain_df['source'] = 'train'\ntest_df['source'] = 'test'\nsubmission_df['source'] = 'submission'\n\nsubmission_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_df = train_df.append([test_df, submission_df])\ndata_df.reset_index(inplace = True)\ndata_df = data_df.drop('index', axis=1)\ndata_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And for a single patient from Test:"},{"metadata":{"_uuid":"5b215ff7-b2ef-4b91-91bf-fb456ee0172d","_cell_guid":"2d7ec248-be26-4e16-a96a-b113bc4174de","trusted":true},"cell_type":"code","source":"data_df[data_df['Patient'] == 'ID00419637202311204720264']","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"10126e32-c091-4ce1-bd21-63e94b33697c","_cell_guid":"b499d655-c572-47bb-98d1-3595b83120eb","trusted":true},"cell_type":"markdown","source":"## Feature Engineering\n### Weeks \nLet's add the minimum week for each Patient as a feature as well as a \"Baseline\" Week feature, which is Week relative to this first week. We need to do this carefully so it makes sense for both the training data and the submission data, which includes ALL weeks for the test Patients.\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"data_df['Weeks'] = data_df['Weeks'].astype(int)\ndata_df['min_week'] = data_df['Weeks'] #first create column\ndata_df.loc[data_df['source']=='submission','min_week'] = np.nan #next get rid of submission's \"false\" weeks\ndata_df['min_week'] = data_df.groupby('Patient')['min_week'].transform('min') #now we can groupby Patient and use min of min_week\n\ndata_df['baseline_week'] = data_df['Weeks'] - data_df['min_week'] #add column that represents offset from min Week\n\ndata_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As a dataset with a time component, it is very important to ensure our time measurements are in good shape for prediction. In this case we don't have a proper time feature that represents the progression of the dependent variable (FVC, a proxy for progression of the disease); instead we have Week relative to CT scan. As they are relative to the CT scan, they can be negative, meaning that an FVC test was taken some weeks before the CT scan.\n\nHowever, week relative to CT scan does not matter. What does matter is progression of the disease relative to each patient (the assumption being made is the initial FVC measurement is somewhere near the initial onset and subsequent diagnosis of the disease, while timing of CT scan dependent on other factors like availability/scheduling). To adjust for this, we will create a new feature called Weeks_Adj to adjust the Weeks field so the progression is more comparable. \n\nThe Weeks_Adj will alter negative weeks so as to have the first FVC measure Weeks_Adj = 0 (equivalent to Week 0), and offset all further Week_Adj data points for the patient. Then, for all other patients, their CT Scan is already at Week 0, so no adjustment is needed, assuming that a CT Scan would signal that the disease is present and a scan is needed in order to diagnose. We assume that if they had the CT scan of their chest, it was related to degradation of lung function, and for whatever reason early FVC measurements were not taken or available.\n\nThe hope is that the Week_Adj is now more comparable on a patient-by-patient basis and more predictive of the specific future Weeks we are asked to predict than Weeks relative to CT scan."},{"metadata":{"_uuid":"90bb4fc4-6738-41ab-92c9-f3de979ab03d","_cell_guid":"f496ebea-8316-47c5-8c48-79ca4a8d4391","trusted":true},"cell_type":"code","source":"data_df= data_df.assign(weeks_adj=[y if x < 0 else z for x, y, z in zip(\n                                       data_df['min_week'], \n                                       data_df['baseline_week'], \n                                       data_df['Weeks'])])\ndata_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Note how these fields will appear for our submission data:"},{"metadata":{"_uuid":"dc1d895e-c345-4909-84b0-2e83fa53daa2","_cell_guid":"312cb053-4760-4cad-a159-ff1866bc7898","trusted":true},"cell_type":"code","source":"data_df[data_df['source'] == 'submission']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can quickly confirm only when min_week is less than 0 is when weeks_adj is different than Weeks:"},{"metadata":{"trusted":true},"cell_type":"code","source":"weeksadj_df = data_df[data_df['Weeks'] != data_df['weeks_adj']] #all data points where Weeks is changed\nweeksadj_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weeksadj_df[weeksadj_df['min_week'] > -1] #confirming it is only patients whose min_week is negative","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### FVC Related Features\nNow let's create features for Initial FVC and Initial Percent for each row in our data. We'll use our new Weeks features to make this fairly simple to do!"},{"metadata":{"_uuid":"1baccf98-8a7b-4343-b866-194fb8bbdfe0","_cell_guid":"930765a6-299d-412c-a38e-3b9ffe9fedea","trusted":true},"cell_type":"code","source":"#get initial % and FVC for each patient as features for every row\nfor field in ['FVC', 'Percent']:\n    baseline_df = data_df[data_df['baseline_week'] == 0].groupby('Patient')[field].mean().to_frame()\n    data_df = data_df.merge(baseline_df, on = 'Patient')\n    data_df = data_df.rename(columns={field+\"_x\": field, field+\"_y\": \"initial_\"+field.lower()})\ndata_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#you can run these lines to check that for Submission rows initial fvc and initial percent are correct (the resulting dfs are empty)\n#data_df[(data_df['source'] == 'submission') & (data_df['FVC'] != data_df['initial_fvc'])]\n#data_df[(data_df['source'] == 'submission') & (data_df['Percent'] != data_df['initial_percent'])]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since we are given Percent, we can also assume that there is a \"normal\" FVC for each patient. Not all factors that determine it are provided in the dataset, but I believe age, sex, and height are the major factors. We can confirm by calculating this feature:"},{"metadata":{"_uuid":"49ede80e-6622-4265-858f-5eed25f14945","_cell_guid":"be536e24-6eee-4db0-ad3a-d5ce1b7e70fb","trusted":true},"cell_type":"code","source":"data_df['fvc_norm'] = data_df.FVC / data_df.Percent * 100\ndata_df[(data_df['baseline_week'] == 0) & (data_df['source'] == 'train')].tail()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Indeed we can see Patients with similar Sex and Age can have very different normal FVC levels!"},{"metadata":{"trusted":true},"cell_type":"code","source":"data_df","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fde34199-6bf7-42fc-be72-a2b8f0a2b60c","_cell_guid":"8b3a2633-afc4-48be-ad92-9a36ecd59097","trusted":true},"cell_type":"markdown","source":"### One-Hot Encoding\nNext, we'll one-hot encode our categorial variables - just a few and with only a couple categories."},{"metadata":{"trusted":true},"cell_type":"code","source":"data_df = pd.get_dummies(data_df, columns=['Sex', 'SmokingStatus'], drop_first=True)\ndata_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_df.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Our data is looking good!"},{"metadata":{},"cell_type":"markdown","source":"## Split Data into Train and Test\nNow that our transformations are done, let's split our data for training. Our Source feature makes this quite simple! We also should reset the index of our train_df so the indices match up, as we'll need to use the indices for modeling."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = data_df.loc[data_df['source'] == 'train'].copy()\ntest_df = data_df.loc[data_df['source'] == 'test'].copy()\nsubmission_df = data_df.loc[data_df['source'] == 'submission'].copy()\n\ntrain_df = train_df.reset_index().drop('index', axis=1)\ntrain_df.describe()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ce01c874-197a-4ae1-9fea-c3b4d7194e26","_cell_guid":"23b4f5f3-8ec4-41c6-9668-75bf5c39b288","trusted":true},"cell_type":"markdown","source":"Note min_week - has an average of 14, which means on average the first FVC intake is 14 weeks after the CT scan, with a std of 15. Max is Week 79 after the CT scan. Fairly interesting! Let's look at things a little deeper..."},{"metadata":{"_uuid":"9697e7fd-5206-4441-83d0-8d88a58b2943","_cell_guid":"d3d3db88-eb20-4f70-9def-27b73d35bf85","trusted":true},"cell_type":"markdown","source":"# Data Analysis\n## Final 3 FVC\nRemember that for scoring, what matters is the final 3 measurements. For the training data, we know what these weeks are, but for the test data we don't, so we can't train specifically to predict that - and thus, must submit predictions for all possible Weeks. Regardless, let's create a feature in our training DataFrame for analysis purposes that marks whether the row is a final 3 FVC measurement or not. We may be able to gain some intuition to see if there is anything we can do to train a model that would predict these better."},{"metadata":{"_uuid":"a35336e2-291d-43b0-8e38-700e9e5b21a8","_cell_guid":"65104800-7faa-4368-9955-7bf939723fc4","trusted":true},"cell_type":"code","source":"train_df['final_3'] = train_df['Patient'] != train_df['Patient'].shift(-3)\ntrain_df.head(25)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"What do those final 3 look like? First we could look at the distribution of final 3 FVC measurements over FVC itself:"},{"metadata":{"trusted":true},"cell_type":"code","source":"bins = range(800, 6800, 500)\ntrain_df[train_df['final_3'] == True]['FVC'].plot.hist(bins=bins, alpha=0.5, label='Final 3')\ntrain_df[(train_df['final_3'] == False)]['FVC'].plot.hist(bins=bins, alpha=0.4, label='not Final 3')\nplt.legend(prop={'size': 12})\nplt.title('Distribution of Final 3 FVC Measurements Over FVC')\nplt.xlabel('FVC')\nplt.ylabel('Frequency')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It seems that the Final 3 measurements are more likely to be distributed in FVC < 2000, as we may expect, but not hugely so. The other brackets above 2000 show >75% of the measurements outside of the final 3, but still a good amount in FVCs above 3500+ compared to total measurements.\n\nWhat about Final 3 in relation to Percent?"},{"metadata":{"trusted":true},"cell_type":"code","source":"bins = range(20, 100, 5)\ntrain_df[train_df['final_3'] == True]['Percent'].plot.hist(bins=bins, alpha=0.5, label='Final 3')\ntrain_df[(train_df['final_3'] == False)]['Percent'].plot.hist(bins=bins, alpha=0.4, label='not Final 3')\nplt.legend(prop={'size': 12})\nplt.title('Distribution of Final 3 FVC Measurements Over Percent')\nplt.xlabel('Percent')\nplt.ylabel('Frequency')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Seems like Final 3 % have a fair representation of 25-40% of the measurmements of FVC Percent above 50%. Less than than that though, and the Final 3 are 75%+ of the FVC measurements where FVC Percent is less than 50%.\n\nNow, what about Weeks?"},{"metadata":{"trusted":true},"cell_type":"code","source":"bins = range(-10,140,5)\ntrain_df[train_df['final_3'] == True]['Weeks'].plot.hist(bins=bins, alpha=0.5, label='Final 3')\ntrain_df[train_df['final_3'] == False]['Weeks'].plot.hist(bins=bins, alpha=0.5, label='not Final 3')\nplt.legend(prop={'size': 12})\nplt.title('Distribution of Final 3 FVC Measurements Over Weeks')\nplt.xlabel('Weeks')\nplt.ylabel('Frequency')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Somewhat surprisingly, we have most of our final 3 FVC measurements in the 25-70 weeks range. In fact these are not necessarily backloaded in our Weeks distribution, at least not as much as I would have thought, with the largest bins being 30, 40, and 55 Weeks.\n\nThis graph illustrates a key point. The majority of our training data for our model is helping it learn about early weeks which may not be very helpful to our model. Instead, we would prefer the model learn accordingly with the blue distribution - we want a model that could predict Weeks 25+ very well, at the expense of not being able to predict Weeks 0-25 well whatsoever, let's say.\n\nWith this more stark split between Final 3 and not, we could possibly tune our model to predict certain Weeks better than others as a proxy for Final 3 performance."},{"metadata":{"trusted":true},"cell_type":"code","source":"#for i in range(40, 58):\n  #  cond_train_df = train_df[(train_df['Weeks'] > i)]\n   # final_3_rows = cond_train_df[cond_train_df['final_3'] == True].shape[0]\n   # total_rows = cond_train_df.shape[0] \n   # print(f'Weeks above {i}, have {round(final_3_rows/total_rows*100,2)}% of their data in Final 3 FVC, or {final_3_rows} \\\n#out of {total_rows} total rows.')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Modeling\n## Prep Training DataFrame\nFirst we need to do any last adjustments to our training data before we feed into the Dataloaders. We have some NaN columns that we can simply remove from training."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = train_df.drop(['Patient_Week', 'Confidence'], axis=1)\ntrain_df","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3c5fa0b3-125a-4d5e-8109-dccac9b6f0aa","_cell_guid":"ac547659-f9ae-46e4-8bd1-199e8d6eda0f","trusted":true},"cell_type":"markdown","source":"\n## Validation Set and Features Lists\nNow let's decide on our validation set. This is an extremely important part of the modeling process!\n\nWe need to select a subset of training data that includes rows for a few specific patients so the model will be validated against patients it has never seen before. This mirrors the test set - new patients the model has never seen before. \n\nNext - how much of our small dataset should we set aside for validation?"},{"metadata":{"_uuid":"e3194fc7-08cc-42ae-8e89-c93b0920833f","_cell_guid":"5243af5f-2a12-40c9-b4b4-22b2f0e700ef","trusted":true},"cell_type":"code","source":"valid_set_patient_num = 30\ncond = train_df['Patient'].isin(np.random.choice(train_df.Patient.unique(),valid_set_patient_num,replace=False))\n\nvalid_idx = np.where(cond)[0]\ntrain_idx = np.where(~cond)[0]\nsplits = (list(train_idx),list(valid_idx))\ntrain_idx.shape, valid_idx.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And with our list of validation indices, we prepare the indices in the dataframe for our fastai TabularPandas class, which creates new training and validation objects for us to work with:"},{"metadata":{"_uuid":"04756d33-ee06-4472-adbf-8c35425604d1","_cell_guid":"dbe8a38a-f239-441d-bed3-d13378cef0f9","trusted":true},"cell_type":"code","source":"cat_features = ['Sex_Male', 'SmokingStatus_Ex-smoker', 'SmokingStatus_Never smoked']\ncont_features = ['Weeks','Age', 'min_week', 'baseline_week', 'weeks_adj', 'initial_percent', 'initial_fvc', 'fvc_norm']\nprocs = Categorify\ntab_obj = TabularPandas(train_df, procs, cat_features, cont_features, y_names='FVC', splits=splits)\nlen(tab_obj.train), len(tab_obj.valid)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Note these are the only features we'll have at test time. We won't have Percent, or final_3, so we can't include them in the model.\n\n"},{"metadata":{"_uuid":"7e3abf92-0888-4d8c-a9b5-6724b2b3a2a9","_cell_guid":"531b6af6-e36c-49d7-b23c-b72d9da74616","trusted":true},"cell_type":"code","source":"tab_obj.show(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"eaa6ca01-c70b-45ab-825c-8e34a4187d7a","_cell_guid":"1630b0f5-bd2f-4624-a0d5-87a119dd86a0","trusted":true},"cell_type":"code","source":"tab_obj.items.head(10) #all columns from train_df are still here!","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a9d5cfe5-8ba3-4e83-9618-315ada3ac61d","_cell_guid":"f8c7ace5-7925-4a24-ad43-f681d6c80be0","trusted":true},"cell_type":"markdown","source":"## First Model - Random Forest\n\nNow let's create our Random Forest Function. We'll evaluate using root mean squared error at first, just to get a sense for things."},{"metadata":{"_uuid":"3a06c1ec-f636-4ee2-aaac-19c3f4b0c8ce","_cell_guid":"e6b87d90-786e-4379-a884-e4049f841cdf","trusted":true},"cell_type":"code","source":"def create_rf(xs, y, n_estimators=200, max_features=0.7, min_samples_leaf=7, **kwargs):\n    return RandomForestRegressor(n_jobs=-1, n_estimators=n_estimators, max_features=max_features, \n                                 min_samples_leaf=min_samples_leaf, oob_score=True).fit(xs,y)\n\ndef rmse(pred,y): return round(math.sqrt(((pred-y)**2).mean()), 6)\ndef model_rmse(m, xs, y): return rmse(m.predict(xs),y)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can reassign the train, valid x and y columns easily accessible via Fastai's tabular object into easy to access variables..."},{"metadata":{"_uuid":"a4b21832-63b6-4708-b3a4-fee55d2402df","_cell_guid":"f03c9b4b-6a4e-4288-bf71-b5cd8016a350","trusted":true},"cell_type":"code","source":"xs, y = tab_obj.train.xs, tab_obj.train.y\nvalid_xs, valid_y = tab_obj.valid.xs, tab_obj.valid.y","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Create a forest and test it on RMSE!"},{"metadata":{"_uuid":"80e4dced-128b-42d6-8211-7e2eed81f8b1","_cell_guid":"53e40d25-bf1e-4ca1-99f2-87fdc39f1980","trusted":true},"cell_type":"code","source":"n_trees = 300\nmax_features=.7\nmin_samples_leaf=7\nm = create_rf(xs, y, n_trees, max_features, min_samples_leaf)\nmodel_rmse(m, xs, y), model_rmse(m, valid_xs, valid_y), rmse(m.oob_prediction_, y)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can use the m.estimators_ attribute of our model to access the predictions for each tree and combine them into a matrix called preds, where each column is a set of predictions, and each row is each individual tree in our random forest. The mean of this for each column evaluated against the validation array is the same as our RMSE!"},{"metadata":{"_uuid":"87169681-0538-4174-bc5f-8eb066b1f467","_cell_guid":"8836494f-efc7-4922-bd2f-77cae307cea9","trusted":true},"cell_type":"code","source":"preds = np.stack(t.predict(valid_xs) for t in m.estimators_)\npd.DataFrame(data=preds)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rmse(preds.mean(0), valid_y)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"33ae8b64-585b-4f9e-9fa0-67bc7f8db21b","_cell_guid":"c16a0e0c-ab4d-4190-9449-fc711817c8d2","trusted":true},"cell_type":"markdown","source":"## Evaluation with Laplace Log Likelihood\nNow let's define the competition metric to use as a score: Laplace Log Likelihood. This metric will be negative, and a higher score is better. A perfect model would score about -4.6."},{"metadata":{"_uuid":"0c5680da-1ba9-44b9-833f-7fbf21bf226e","_cell_guid":"90656dd9-cacd-404d-9a21-e299c2f0aba3","trusted":true},"cell_type":"code","source":"def LaplaceLogLikelihood(pred, y, sigma):\n    \n    sigma_clip = np.maximum(sigma, 70)\n    delta = np.minimum(np.absolute(y - pred), 1000.)\n    sq2 = math.sqrt(2.)\n    metric = ((delta / sigma_clip) * -sq2) - np.log(sq2 * sigma_clip)\n    return np.mean(metric)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"But remember, our metric is only applied to the final 3 FVC measurements. We'll need to add functionality to score ourselves on only these data points. With our final_3 feature we created earlier, we easily access the indices of the final 3 FVC measurements, which we can use for scoring."},{"metadata":{"trusted":true},"cell_type":"code","source":"final3_idx = np.where(train_df['final_3'] == True)[0]\nvalid_y_final3 = valid_y.loc[valid_y.index.intersection(final3_idx)]\nfinal3_idx.shape, valid_y_final3.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We create a function that we can pass our model, xs, y, and sigma, and will calculate the LLL score for only the xs and ys that appear in the final 3 FVC measurements."},{"metadata":{"trusted":true},"cell_type":"code","source":"def model_score(m, xs, y, sigma): \n    final3_xs = xs.loc[xs.index.intersection(final3_idx)] #intersection of final3_idx and provided rows\n    final3_y = y.loc[y.index.intersection(final3_idx)]\n    \n    sigma_df = pd.Series(data=sigma, index=xs.index) #transform sigma into Series with indices to match xs\n    final3_sigma = sigma_df.loc[sigma_df.index.intersection(final3_idx)] #filter those indices for final 3s\n    \n    return LaplaceLogLikelihood(m.predict(final3_xs), final3_y, final3_sigma) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For this scoring, we also need a confidence measure, which is related to standard deviation. We can use the preds array we used before to calculate the standard deviation of each prediction across all trees, and use this as a measure of confidence."},{"metadata":{"_uuid":"1b5a941d-1f6e-4eed-8e09-64805e498ddf","_cell_guid":"d03bc41b-34ca-492b-9be1-f65cf160d006","trusted":true},"cell_type":"code","source":"preds_std = preds.std(axis=0)\npreds_std.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"67ed1b94-6a67-4999-a238-12950d03f6fa","_cell_guid":"2ea173b2-b5fe-426a-a4f6-da9b10838119","trusted":true},"cell_type":"code","source":"preds_std[:25]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's try it out!"},{"metadata":{"_uuid":"7d5ada56-e93b-414d-abb4-3b3362cccb4e","_cell_guid":"dd1f0128-b801-4faf-9899-cdd1a697eb49","trusted":true},"cell_type":"code","source":"model_score(m, xs, y, 231), model_score(m, valid_xs, valid_y, preds_std)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And we have a LLL Score for our RF Model! Compared to the variance between trees, what's the best a constant Confidence could give us?"},{"metadata":{"trusted":true},"cell_type":"code","source":"def constant_sigma_choice(m, xs, y):\n    sc = -100.\n    sigma = 70\n    for i in range(71, 500, 5):\n        temp = model_score(m, xs ,y, i)\n        if temp > sc: \n            sc = temp\n            sigma = i\n    return sc, sigma\nscore, sigma = constant_sigma_choice(m, valid_xs, valid_y)\nscore, sigma","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model_score(m, valid_xs, valid_y, sigma)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Final Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"n_trees = 1000\nmax_features= .6\nmin_samples_leaf= 3\nm = create_rf(xs, y, n_trees, max_features, min_samples_leaf)\n#get new preds matrix for std\npreds = np.stack(t.predict(valid_xs) for t in m.estimators_)\npreds_std = preds.std(axis=0)\n[model_score(m, xs, y, 231), \n model_score(m, valid_xs, valid_y, preds_std), \n model_score(m, valid_xs, valid_y, constant_sigma_choice(m, valid_xs, valid_y)[1])]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"eval_df = valid_xs.copy().reset_index().drop('index', axis=1)\npreds = pd.Series(m.predict(valid_xs))\neval_df['FVC_pred'] = preds\neval_df['FVC_true'] = valid_y.reset_index().drop('index', axis=1)\neval_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"eval_df['FVC_error'] = np.absolute(eval_df['FVC_true'] - eval_df['FVC_pred'])\neval_df.groupby('initial_percent').transform('mean').drop_duplicates().sort_values(by=['FVC_error'], ascending=False)[0:25]\n#eval_df[eval_df['initial_fvc'] == 1690]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df[train_df['initial_fvc'] == (6399)]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a120da65-3cac-45f3-a2b5-048dd0c1126d","_cell_guid":"de8599cc-7c82-412a-8559-389c6b608dd1","trusted":true},"cell_type":"markdown","source":"# Feature Importance"},{"metadata":{"_uuid":"b67ab8ca-dc1b-479b-87c3-2db49a2aa87c","_cell_guid":"917c4a33-0893-4554-9927-b7377b4d5f1b","trusted":true},"cell_type":"code","source":"def rf_feat_importance(m, df):\n    return pd.DataFrame({'cols':df.columns, 'imp':m.feature_importances_}).sort_values('imp', ascending=False)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f343f5c-42b2-43a7-88b7-868c819d124f","_cell_guid":"5e0f35f7-e6dd-4e66-b71a-23e7791a61c0","trusted":true},"cell_type":"code","source":"fi = rf_feat_importance(m, xs)\nfi","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7fb455a3-dc28-4016-af23-bf77c48d955e","_cell_guid":"f3e77f7b-0e1a-4235-967d-259ae8b2ef54","trusted":true},"cell_type":"code","source":"from sklearn.inspection import plot_partial_dependence\n\nfig,ax = plt.subplots(figsize=(12, 4))\nplot_partial_dependence(m, valid_xs, ['Sex_Male', 'Age'],\n                        grid_resolution=20, ax=ax);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2bfd06a7-881e-48bc-86ae-25d1107025ff","_cell_guid":"fb19becd-c591-4246-8058-89bdc087ded0","trusted":true},"cell_type":"markdown","source":"# Submission\nLet's use our model to predict the FVC on the submission dataframe and submit!"},{"metadata":{"trusted":true},"cell_type":"code","source":"xs_submit = submission_df[cat_features + cont_features]\nsubmission_df['FVC'] = m.predict(xs_submit)\n#use standard deviations of each Tree as confidence\npreds_submit = np.stack(t.predict(xs_submit) for t in m.estimators_)\nsubmission_df['Confidence'] = preds_submit.std(axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_df = submission_df[['Patient_Week', 'FVC', 'Confidence']]\nsubmission_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_df.to_csv('submission.csv', index=False)\nsubmission_df.describe()","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}