{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Introduction and Acknowledgement\n\nThis notebook is inspired by [Carlos Souza's Bayesian notebook](https://www.kaggle.com/carlossouza/bayesian-experiments). Go check it out! The core ideas in the code here draw heavily from this notebook.\n\nSome of the differences I implemented are:\n* The values for FVC are in log. Rationale: I'm no expert in pulmonology, but the way I see it (or more precisely, the mental model I have of FVC), FVC is a multiplicative function of various factors. For example, starting from a baseline FVC, whether or not a person smokes would decrease her FVC by X%, a person's age would further decrease by Y%, long-term breathing from low quality air could decrease by another Z%, etc. When we consider all of these factors, the FVC would be FVC_{baseline} X (1 - X%) X (1 - Y%) X (1 - Z%) etc. This means that the variable is likely to have a log-normal distribution and therefore its log would be an appropriate transformation to the data. Another indicator that this could be appropriate is the fact that you can only have positive FVC values, never negative (which would make no sense!). More inform on log-normal distribution in the [usual place](https://en.wikipedia.org/wiki/Log-normal_distribution).\n* I calculate parameters for smoking status and age, not for each individual. Rationale: we won't have the data to estimate posterior distributions for new patients, and since the goal is to forecast out-of-sample, say, when a patient first goes to a doctor. So from my perspective, we should estimate an unpooled model on population information that the doctor has at hand at the time of first consultation, so they can share a prognosis with the patient.\n* TO-DO: incorporate the image data to the tabular predictions. Rationale: potentially very useful information. One idea I will try to incorporate before the timeline for this competition elapses is to use the convolutional layers of one of the EfficientNet models to extract features, and then use dense layers together with max pooling until we get to, say 8 vectors that summarise the whole CT scan for each patient, and then incorporate these vectors in the Bayesian regression, potentially with a hierarchical hyperparameter controlling how tight or loose are my priors for the coefficients on these 8 vectors. (Side note: Ideally the model should be informative without the CT scan, as this would could help the medical community benefit patients without access to CT in the first place.)\n\n# Training a Bayesian log linear model"},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom sklearn.preprocessing import LabelEncoder\nimport pymc3 as pm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"exclude_test_patient_data_from_trainset = True\n\ntrain = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\n\nif exclude_test_patient_data_from_trainset:\n    train = train[~train['Patient'].isin(test['Patient'].unique())]\n\ntrain = pd.concat([train, test], axis=0, ignore_index=True)\\\n    .drop_duplicates()\n\nle_id = LabelEncoder()\ntrain['PatientID'] = le_id.fit_transform(train['Patient'])\n\n# The lines below are my innovations\ntrain['FVC_original'] = train['FVC']\ntrain['FVC'] = np.log(train['FVC'])\n\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['Male'] = train['Sex'].apply(lambda x: 1 if x == 'Male' else 0)\n\ntrain[\"SmokingStatus\"] = train[\"SmokingStatus\"].astype(\n    pd.CategoricalDtype(['Ex-smoker', 'Never smoked', 'Currently smokes'])\n)\naux = pd.get_dummies(train[\"SmokingStatus\"], prefix='ss')\naux.columns = ['ExSmoker', 'NeverSmoked', 'CurrentlySmokes']\ntrain['CurrentlySmokes'] = aux['CurrentlySmokes']\n\naux = train[['Patient', 'Weeks', 'Percent']].sort_values(by=['Patient', 'Weeks'])\naux = train.groupby('Patient').head(1)\naux = aux.rename(columns={'Percent': 'Percent_base'})\ntrain = pd.merge(train, aux[['Patient', 'Percent_base']], how='left',\n                 on='Patient')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.loc[train['Age'] < 50, 'YoungPatient'] = 1\ntrain.loc[train['Age'] >= 50, 'YoungPatient'] = 0\ntrain['YoungPatient'] = pd.get_dummies(train['YoungPatient'])\n\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will build two models, one for smoker patients and the other for non-smokers. Let's see first what the usual time period between the first FVC and last are for each one of these sub-samples."},{"metadata":{"trusted":true},"cell_type":"code","source":"smokers_first_FVC_week = train.loc[train['SmokingStatus'] == 'Currently smokes', ['Patient', 'Weeks']].sort_values(by = ['Patient', 'Weeks']).groupby('Patient').head(1)\nsmokers_last_FVC_week = train.loc[train['SmokingStatus'] == 'Currently smokes', ['Patient', 'Weeks']].sort_values(by = ['Patient', 'Weeks']).groupby('Patient').tail(1)\nsmokers_weeks_range = smokers_last_FVC_week['Weeks'].values - smokers_first_FVC_week['Weeks'].values\nsmokers_weeks_range.mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nonsmokers_first_FVC_week = train.loc[train['SmokingStatus'] != 'Currently smokes', ['Patient', 'Weeks']].sort_values(by = ['Patient', 'Weeks']).groupby('Patient').head(1)\nnonsmokers_last_FVC_week = train.loc[train['SmokingStatus'] != 'Currently smokes', ['Patient', 'Weeks']].sort_values(by = ['Patient', 'Weeks']).groupby('Patient').tail(1)\nnonsmokers_weeks_range = nonsmokers_last_FVC_week['Weeks'].values - nonsmokers_first_FVC_week['Weeks'].values\nnonsmokers_weeks_range.mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Hyperparameter priors for both models (one for current smokers and the other for everyone else)\n\n# For the smokers' model\nFVC_smokers_mean = train.loc[train['SmokingStatus'] == 'Currently smokes', 'FVC'].mean()\nFVC_smokers_stdev = train.loc[train['SmokingStatus'] == 'Currently smokes', 'FVC'].std() / 2\n\nFVC_smokers_beta_mean = train.loc[(train['SmokingStatus'] == 'Currently smokes') & (train['Weeks'] > smokers_last_FVC_week['Weeks'].values.mean()), 'FVC'].mean() - train.loc[(train['SmokingStatus'] == 'Currently smokes') & (train['Weeks'] < smokers_first_FVC_week['Weeks'].values.mean()), 'FVC'].mean()\nFVC_smokers_beta_weeks_mean = FVC_smokers_beta_mean / smokers_weeks_range.mean()\nFVC_smokers_beta_others_sd = 1 # Must be > 0. Values close to zero indicate tighter priors, make it more difficult for posterior distribution of parameters to not include zero.\n\n# And for the non-smokers' model\nFVC_nonsmokers_mean = train.loc[train['SmokingStatus'] != 'Currently smokes', 'FVC'].mean()\nFVC_nonsmokers_stdev = train.loc[train['SmokingStatus'] != 'Currently smokes', 'FVC'].std() / 2\n\nFVC_nonsmokers_beta_mean = train.loc[(train['SmokingStatus'] != 'Currently smokes') & (train['Weeks'] > nonsmokers_last_FVC_week['Weeks'].values.mean()), 'FVC'].mean() - train.loc[(train['SmokingStatus'] != 'Currently smokes') & (train['Weeks'] < nonsmokers_first_FVC_week['Weeks'].values.mean()), 'FVC'].mean()\nFVC_nonsmokers_beta_weeks_mean = FVC_smokers_beta_mean / nonsmokers_weeks_range.mean()\nFVC_nonsmokers_beta_others_sd = 1 # Must be > 0. Values close to zero indicate tighter priors, make it more difficult for posterior distribution of parameters to not include zero.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Data used to update priors\nsmokers_n_patients = train.loc[train['CurrentlySmokes'] == 1, 'Patient'].nunique()\nsmokers_FVC_obs = train.loc[train['CurrentlySmokes'] == 1, 'FVC'].values\nsmokers_PatientID = train.loc[train['CurrentlySmokes'] == 1, 'PatientID'].values\n\nsmokers_Weeks = train.loc[train['CurrentlySmokes'] == 1, 'Weeks'].values\nsmokers_Male = train.loc[train['CurrentlySmokes'] == 1, 'Male'].values\nsmokers_Young = train.loc[train['CurrentlySmokes'] == 1, 'YoungPatient'].values\nsmokers_Percent = train.loc[train['CurrentlySmokes'] == 1, 'Percent_base'].values\n\nnonsmokers_n_patients = train.loc[train['CurrentlySmokes'] == 0, 'Patient'].nunique()\nnonsmokers_FVC_obs = train.loc[train['CurrentlySmokes'] == 0, 'FVC'].values\nnonsmokers_PatientID = train.loc[train['CurrentlySmokes'] == 0, 'PatientID'].values\n\nnonsmokers_Weeks = train.loc[train['CurrentlySmokes'] == 0, 'Weeks'].values\nnonsmokers_Male = train.loc[train['CurrentlySmokes'] == 0, 'Male'].values\nnonsmokers_Young = train.loc[train['CurrentlySmokes'] == 0, 'YoungPatient'].values\nnonsmokers_Percent = train.loc[train['CurrentlySmokes'] == 0, 'Percent_base'].values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with pm.Model() as model_smokers:\n    # create shared variables that can be changed later on\n    FVC_obs_shared = pm.Data(\"FVC_obs_shared\", smokers_FVC_obs)\n    PatientID_shared = pm.Data('PatientID_shared', smokers_PatientID)\n    Weeks_shared = pm.Data('Weeks_shared', smokers_Weeks)\n    Male_shared = pm.Data('Male_shared', smokers_Male)\n    YoungPatient_shared = pm.Data('YoungPatient_shared', smokers_Young)\n    Percent_shared = pm.Data('Percent_shared', smokers_Percent)\n    # Model parameter priors\n    a = pm.Normal('a', mu = FVC_smokers_mean, sigma = FVC_smokers_stdev)\n    b_week = pm.Normal('b_week', mu = FVC_smokers_beta_weeks_mean, sigma = 1.)\n    b_young = pm.Normal('b_young', mu = 0, sigma = FVC_smokers_beta_others_sd)\n    b_weeks_young = pm.Normal('b_weeks_young', mu = 0, sigma = FVC_smokers_beta_others_sd)\n    b_gender = pm.Normal('b_gender', mu = 0, sigma = FVC_smokers_beta_others_sd)\n    b_weeks_gender = pm.Normal('b_weeks_gender', mu = 0, sigma = FVC_smokers_beta_others_sd)\n    b_weeks_young_gender = pm.Normal('b_weeks_young_gender', mu = 0, sigma = FVC_smokers_beta_others_sd)\n    b_percent = pm.Normal('b_percent', mu = 0, sigma = FVC_smokers_beta_others_sd)\n    # Model error\n    sigma = pm.HalfNormal('sigma', np.log(70.))\n    # Model formula\n    FVC_est = a + b_week * Weeks_shared + b_young * YoungPatient_shared + b_weeks_young * Weeks_shared * YoungPatient_shared + b_gender * Male_shared + b_weeks_gender * Weeks_shared * Male_shared + b_weeks_young_gender * Weeks_shared * YoungPatient_shared * Male_shared + b_percent * Percent_shared\n    # Data likelihood\n    FVC_like = pm.Normal('FVC_like', mu = FVC_est, sigma = sigma, observed = FVC_obs_shared)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with pm.Model() as model_nonsmokers:\n    # create shared variables that can be changed later on\n    FVC_obs_shared = pm.Data(\"FVC_obs_shared\", nonsmokers_FVC_obs)\n    PatientID_shared = pm.Data('PatientID_shared', nonsmokers_PatientID)\n    Weeks_shared = pm.Data('Weeks_shared', nonsmokers_Weeks)\n    Male_shared = pm.Data('Male_shared', nonsmokers_Male)\n    YoungPatient_shared = pm.Data('YoungPatient_shared', nonsmokers_Young)\n    Percent_shared = pm.Data('Percent_shared', nonsmokers_Percent) \n    # Model parameter priors\n    a = pm.Normal('a', mu = FVC_nonsmokers_mean, sigma = FVC_nonsmokers_stdev)\n    b_week = pm.Normal('b_week', mu = FVC_nonsmokers_beta_weeks_mean, sigma = 1.)\n    b_young = pm.Normal('b_young', mu = 0, sigma = FVC_nonsmokers_beta_others_sd)\n    b_weeks_young = pm.Normal('b_weeks_young', mu = 0, sigma = FVC_nonsmokers_beta_others_sd)\n    b_gender = pm.Normal('b_gender', mu = 0, sigma = FVC_nonsmokers_beta_others_sd)\n    b_weeks_gender = pm.Normal('b_weeks_gender', mu = 0, sigma = FVC_nonsmokers_beta_others_sd)\n    b_weeks_young_gender = pm.Normal('b_weeks_young_gender', mu = 0, sigma = FVC_nonsmokers_beta_others_sd)\n    b_percent = pm.Normal('b_percent', mu = 0, sigma = FVC_nonsmokers_beta_others_sd)\n    # Model error\n    sigma = pm.HalfNormal('sigma', np.log(70.))\n    # Model formula\n    FVC_est = a + b_week * Weeks_shared + b_young * YoungPatient_shared + b_weeks_young * Weeks_shared * YoungPatient_shared + b_gender * Male_shared + b_weeks_gender * Weeks_shared * Male_shared + b_weeks_young_gender * Weeks_shared * YoungPatient_shared * Male_shared + b_percent * Percent_shared\n    # Data likelihood\n    FVC_like = pm.Normal('FVC_like', mu = FVC_est, sigma = sigma, observed = FVC_obs_shared)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The snippet below samples the posterior distribution of the parameters for the smokers' model.\n\nBeware: it may take a while to run."},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_smokers:\n    trace_smokers = pm.sample(2000, tune = 25000, chains = 1, target_accept = .95, init = \"adapt_diag\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_smokers:\n    pm.traceplot(trace_smokers);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Similarly, the snippet below samples the posterior distribution of the parameters from the model for non-smokers.\n\nBeware: it may take a while to run (more than the smokers' model)."},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_nonsmokers:\n    trace_nonsmokers = pm.sample(2000, tune = 25000, chains = 1, target_accept = .95, init = \"adapt_diag\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_nonsmokers:\n    pm.traceplot(trace_nonsmokers);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The snippet below compiles the training data in a way to facilitate prediction by both models."},{"metadata":{"trusted":true},"cell_type":"code","source":"aux_train = train.groupby('Patient').first().reset_index()\npred_template_train = []\nfor i in range(train['Patient'].nunique()):\n    df_train = pd.DataFrame(columns=['PatientID', 'Weeks'])\n    df_train['Weeks'] = np.arange(-12, 134)\n    df_train['PatientID'] = i\n    df_train['Male'] = aux_train[aux_train['PatientID'] == i]['Male'].values[0]\n    df_train['CurrentlySmokes'] = aux_train[aux_train['PatientID'] == i]['CurrentlySmokes'].values[0]\n    df_train['YoungPatient'] = aux_train[aux_train['PatientID'] == i]['YoungPatient'].values[0]\n    df_train['Percent_base'] = aux_train[aux_train['PatientID'] == i]['Percent_base'].values[0]\n    pred_template_train.append(df_train)\npred_template_train = pd.concat(pred_template_train, ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we will generate predictions for smokers..."},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_smokers:\n    pm.set_data({\n        \"PatientID_shared\": pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 1, 'PatientID'].values.astype(int),\n        \"Weeks_shared\": pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 1, 'Weeks'].values.astype(int),\n        \"Male_shared\": pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 1, 'Male'].values.astype(int),\n        \"YoungPatient_shared\": pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 1, 'YoungPatient'].values.astype(int),\n        \"Percent_shared\": pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 1, 'Percent_base'].values,\n        \"FVC_obs_shared\": np.zeros(len(pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 1, :])).astype(int),\n    })\n    post_pred_smokers_train = pm.sample_posterior_predictive(trace_smokers)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"... and prediction for non-smokers."},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_nonsmokers:\n    pm.set_data({\n        \"PatientID_shared\": pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 0, 'PatientID'].values.astype(int),\n        \"Weeks_shared\": pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 0, 'Weeks'].values.astype(int),\n        \"Male_shared\": pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 0, 'Male'].values.astype(int),\n        \"YoungPatient_shared\": pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 0, 'YoungPatient'].values.astype(int),\n        \"Percent_shared\": pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 0, 'Percent_base'].values,\n        \"FVC_obs_shared\": np.zeros(len(pred_template_train.loc[pred_template_train['CurrentlySmokes'] == 0, :])).astype(int),\n    })\n    post_pred_nonsmokers_train = pm.sample_posterior_predictive(trace_nonsmokers)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we need to (1) identify the patients according to smoking status, since that will determine which model will be used, and (2) use this mapping to fill out the data frame with the predictions accordingly. The first cell implements (1), and the next cell implements (2)."},{"metadata":{"trusted":true},"cell_type":"code","source":"smoker_patients_train = train.loc[train['SmokingStatus'] == 'Currently smokes', 'Patient'].unique().tolist()\nnonsmoker_patients_train = train.loc[train['SmokingStatus'] != 'Currently smokes', 'Patient'].unique().tolist()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = pd.DataFrame(columns=['Patient', 'Weeks', 'FVC_pred', 'sigma'])\ndf_train['Patient'] = le_id.inverse_transform(pred_template_train['PatientID'])\ndf_train['Weeks'] = pred_template_train['Weeks']\n\n# Now we have to populate the forecasts for the smokers and for the nonsmokers separately\ndf_train.loc[df_train['Patient'].isin(smoker_patients_train), 'FVC_pred'] = np.exp(post_pred_smokers_train['FVC_like'].T.mean(axis=1))\ndf_train.loc[df_train['Patient'].isin(nonsmoker_patients_train), 'FVC_pred'] = np.exp(post_pred_nonsmokers_train['FVC_like'].T.mean(axis=1))\n\ndf_train.loc[df_train['Patient'].isin(smoker_patients_train), 'sigma'] = np.exp(post_pred_smokers_train['FVC_like'].T.std(axis=1))\ndf_train.loc[df_train['Patient'].isin(nonsmoker_patients_train), 'sigma'] = np.exp(post_pred_nonsmokers_train['FVC_like'].T.std(axis=1))\n                                                                \ndf_train['FVC_inf'] = df_train['FVC_pred'] / df_train['sigma']\ndf_train['FVC_sup'] = df_train['FVC_pred'] * df_train['sigma']\n\ndf_train = pd.merge(df_train, train[['Patient', 'Weeks', 'FVC_original']], how='left', on=['Patient', 'Weeks'])\ndf_train = df_train.rename(columns={'FVC_original': 'FVC_true'})\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"OK. Let's review the steps we have taken so far:\nwe loaded the training dataset;\nthen we transformed the dataset according to how we are going to model the data (log linear; use of dummies);\nwe set out prior distributions for the parameters in the model;\nthen we used the likelihood from observed data in our model to sample the posterior distributions of our parameters of interest;\nwith this posterior distribution of parameters, we have applied them to the data to generate predictions (think of it as generating the famous \"y hat\" vector.\nWe can now calculate in the following cell the modified Laplace log likelihood according to the definition used in the challenge."},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ndf = pd.merge(df_train, tr[['Patient', 'Weeks']], how='left', on=['Patient', 'Weeks'])\n\nuse_only_last_3_measures = True\n\nif use_only_last_3_measures:\n    y = df.dropna().groupby('Patient').tail(3)\nelse:\n    y = df.dropna()\n\nrmse = ((y['FVC_pred'] - y['FVC_true']) ** 2).mean() ** (1/2)\nprint(f'Root mean squared error (RMSE): {rmse:.1f} ml')\n\nsigma_c = y['sigma'].values.astype('float64')\nsigma_c[sigma_c < 70] = 70\ndelta = (y['FVC_pred'] - y['FVC_true']).abs()\ndelta[delta > 1000] = 1000\nlll = - np.sqrt(2) * delta / sigma_c - np.log(np.sqrt(2) * sigma_c)\nprint(f'Laplace Log Likelihood: {lll.mean():.4f}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Running the model on test data\n\nNow it is time to download the test data, use the parameters as sampled from the posterior distribution to create FVC estimates (and their confidence), and compile all of this to submit to the competition.\n\nThe first snippet downloads the test data and transforms FVC to log. Subsequently the gender and smoking status transformations that were applied to training data are applied to test data. Finally, the same is done for the age < 50 years dummy."},{"metadata":{"trusted":true},"cell_type":"code","source":"test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\ntest['PatientID'] = le_id.fit_transform(test['Patient'])\ntest['FVC_original'] = test['FVC']\ntest['FVC'] = np.log(test['FVC'])\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test['Male'] = test['Sex'].apply(lambda x: 1 if x == 'Male' else 0)\n\ntest[\"SmokingStatus\"] = test[\"SmokingStatus\"].astype(\n    pd.CategoricalDtype(['Ex-smoker', 'Never smoked', 'Currently smokes'])\n)\naux = pd.get_dummies(test[\"SmokingStatus\"], prefix='ss')\naux.columns = ['ExSmoker', 'NeverSmoked', 'CurrentlySmokes']\ntest['CurrentlySmokes'] = aux['CurrentlySmokes']\n\naux = test[['Patient', 'Weeks', 'Percent']].sort_values(by=['Patient', 'Weeks'])\naux = test.groupby('Patient').head(1)\naux = aux.rename(columns={'Percent': 'Percent_base'})\ntest = pd.merge(test, aux[['Patient', 'Percent_base']], how='left',\n                 on='Patient')\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.loc[train['Age'] < 50, 'YoungPatient'] = 1\ntest.loc[train['Age'] >= 50, 'YoungPatient'] = 0\ntest['YoungPatient'] = pd.get_dummies(test['YoungPatient'])\n\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Because we will not train the model in test data, we can now jump directly to the step where we use this pre-processed data (ie, log FVC instead of original FVC, gender dummy, current smoker dummy, young patient dummy) to generate predictions for the test patients, using the parameters sampled from the posterior distribution that we estimated using training data likelihood and our priors."},{"metadata":{"trusted":true},"cell_type":"code","source":"aux_test = test.groupby('Patient').first().reset_index()\npred_template_test = []\nfor i in range(test['Patient'].nunique()):\n    df_test = pd.DataFrame(columns=['PatientID', 'Weeks'])\n    df_test['Weeks'] = np.arange(-12, 134)\n    df_test['PatientID'] = i\n    df_test['Male'] = aux_test[aux_test['PatientID'] == i]['Male'].values[0]\n    df_test['CurrentlySmokes'] = aux_test[aux_test['PatientID'] == i]['CurrentlySmokes'].values[0]\n    df_test['YoungPatient'] = aux_test[aux_test['PatientID'] == i]['YoungPatient'].values[0]\n    df_test['Percent_base'] = aux_test[aux_test['PatientID'] == i]['Percent_base'].values[0]\n    pred_template_test.append(df_test)\npred_template_test = pd.concat(pred_template_test, ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_smokers:\n    if len(pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 1, :]) > 0:\n        pm.set_data({\n            \"PatientID_shared\": pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 1, 'PatientID'].values.astype(int),\n            \"Weeks_shared\": pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 1, 'Weeks'].values.astype(int),\n            \"Male_shared\": pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 1, 'Male'].values.astype(int),\n            \"YoungPatient_shared\": pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 1, 'YoungPatient'].values.astype(int),\n            \"Percent_shared\": pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 1, 'Percent_base'].values,\n            \"FVC_obs_shared\": np.zeros(len(pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 1, :])).astype(int),\n        })\n        post_pred_smokers_test = pm.sample_posterior_predictive(trace_smokers)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_nonsmokers:\n    if len(pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 0, :]) > 0:\n        pm.set_data({\n            \"PatientID_shared\": pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 0, 'PatientID'].values.astype(int),\n            \"Weeks_shared\": pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 0, 'Weeks'].values.astype(int),\n            \"Male_shared\": pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 0, 'Male'].values.astype(int),\n            \"YoungPatient_shared\": pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 0, 'YoungPatient'].values.astype(int),\n            \"Percent_shared\": pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 0, 'Percent_base'].values,\n            \"FVC_obs_shared\": np.zeros(len(pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 0, :])).astype(int),\n        })\n        post_pred_nonsmokers_test = pm.sample_posterior_predictive(trace_nonsmokers)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"smoker_patients_test = test.loc[test['SmokingStatus'] == 'Currently smokes', 'Patient'].unique().tolist()\nnonsmoker_patients_test = test.loc[test['SmokingStatus'] != 'Currently smokes', 'Patient'].unique().tolist()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_test = pd.DataFrame(columns=['Patient', 'Weeks', 'FVC_pred', 'sigma'])\ndf_test['Patient'] = le_id.inverse_transform(pred_template_test['PatientID'])\ndf_test['Weeks'] = pred_template_test['Weeks']\n\n# Now we have to populate the forecasts for the smokers and for the nonsmokers separately\nif len(pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 1, :]) > 0:\n    df_test.loc[df_test['Patient'].isin(smoker_patients_test), 'FVC_pred'] = np.exp(post_pred_smokers_test['FVC_like'].T.mean(axis=1))\nif len(pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 0, :]) > 0:    \n    df_test.loc[df_test['Patient'].isin(nonsmoker_patients_test), 'FVC_pred'] = np.exp(post_pred_nonsmokers_test['FVC_like'].T.mean(axis=1))\n\nif len(pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 1, :]) > 0:\n    df_test.loc[df_test['Patient'].isin(smoker_patients_test), 'sigma'] = np.exp(post_pred_smokers_train['FVC_like'].T.std(axis=1))\nif len(pred_template_test.loc[pred_template_test['CurrentlySmokes'] == 0, :]) > 0:    \n    df_test.loc[df_test['Patient'].isin(nonsmoker_patients_test), 'sigma'] = np.exp(post_pred_nonsmokers_test['FVC_like'].T.std(axis=1))\n                                                                \ndf_test['FVC_inf'] = df_test['FVC_pred'] / df_test['sigma']\ndf_test['FVC_sup'] = df_test['FVC_pred'] * df_test['sigma']\n\ndf_test = pd.merge(df_test, test[['Patient', 'Weeks', 'FVC_original']], how='left', on=['Patient', 'Weeks'])\ndf_test = df_test.rename(columns={'FVC_original': 'FVC_true'})\ndf_test.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The final step is to get the predictions we generated for the test dataset and write a CSV file according to the request format."},{"metadata":{"trusted":true},"cell_type":"code","source":"df = df_test\ndf['Patient_Week'] = df['Patient'] + '_' + df['Weeks'].astype(str)\ndf = df.rename(columns={'FVC_pred': 'FVC'})\ndf['Confidence'] = (df['FVC_sup'] - df['FVC_inf']) / 2\nfinal = df[['Patient_Week', 'FVC', 'Confidence']]\nfinal.to_csv('submission.csv', index = False)\nprint(final.shape)\nfinal.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Done! That's it! Now we have to work on continuing to improve the model."}],"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}