{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 1. Bayesian Experiments\nIn our [previous notebook](https://www.kaggle.com/carlossouza/probabilistic-machine-learning-a-diff-approach), we publishe a different approach to the competition: training a **Bayesian Linear Regression model** that could not only predict the decline in FVC for each patient but also the confidence. The key benefit of **Probabilistic Machine Learning models** is that they are very good at predicting uncertainty, which is required in several applications such as personalized medicine.\n\nHowever, the first model was not that simple. We used all tabular data and engineered features, creating a hierarquical model with several parameters. The model was **very unstable**, and its submissions only occasionally generated results without **errors**. The purpose of this notebook is to fix that. Our goals are:\n1. Build a stable **error-free** simplest possible **Bayesian Hierarchical Linear Regression model**\n2. Extend this model to use **more tabular data**\n3. **Compare a visualize** both models\n4. Submit the best model solution\n\n# 2. A simpler way to look at the data\nSeveral notebooks, including our first notebook, add a column feature of baseline FVC and Percent. Then, they go on and calculate the time difference in weeks between the baseline FVC and the target FVC to be predicted. We gave a lot of thought on these engineered features, and we decided to **drop them and experiment using a simpler dataset**.\n\nLet us explain why. We believe adding the baseline FVC and Percent and calculating the time difference in weeks between the baseline FVC and the target FVC to be predicted makes sense **only when we are not informing the model about who the patient is**. However, **that's completely unnecessary when we inform the model who the patients are**, identifying the owner of each data point.\n\nBut informing the model about who the patients are wouldn't be a form of data leakage? Again, after giving a lot of thought, we believe that's not the case. **We must inform the model who the patients are**. That's because **the data points are NOT independent**. FVC measures from a given patient are **strongly correlated**. We should not discard that information, it is too valuable.\n\nThe insight that sparked this line of thought was to look at this problem as a **matrix completion** task. It resembles a recommender system challenge or a missing value imputation challenge, such as 2009's [Netflix Prize](https://www.netflixprize.com/). The image below reframes how we are proposing looking at the problem:\n\n<img src=\"https://i.ibb.co/0Z9kW8H/matrix-completion.jpg\" alt=\"drawing\" width=\"800\"/>\n\nUsually in a recommendation engine, the rows are the customers, the columns are the products (e.g. Amazon) or movies (e.g. Netflix), and the cells are ratings from 1 to 5. Collaborative filtering techniques, such as Probabilistic Matrix Factorization, are used to **fill in this (very sparse) table**, recommending new products that users might like on the basis of reactions by **similar users**. The data points (ratings) from a given user are **strongly correlated**; and **by informing who the users are, the model learns to group users with similar interests**.\n\nThe same way of thinking can be applied in this competition. A model might **learn to appropriately group patients** by inferring the stage of the disease, only by being fed with the FVC decline curves from each patient (with several missing values!). Of course there are differences: in a recommender AI system, the order of the columns doesn't matter, and the missing values are constrained in the interval [1, 5]. Here, **the order of the columns matter** (time growing from left to right), and the FVC values are not constrained in a particular interval. Nevertheless, we believe we can benefit from seeing the similarities between these 2 worlds and adjusting the algorithms accordingly where they diverge.\n\nTo prove that this approach works (and actually produces good results!), in the first model we will use only:\n- FVC\n- Week\n- Patient ID\n\nYou might argue it will be \"easy\" for a model to predict future FVC measures once it is trained with several past FVC measures from a given patient. So, **we will make it harder: we will remove all data points from the 5 test patients from the training set**, except the first baseline FVC measurement. This means that **our model will have only a single point** from these 5 patients to use in order to generate future FVC predictions.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","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","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","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\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And for now, that's all there is to data pre-processing! We can again just inspect the FVC declines:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def chart(patient_id, ax):\n    data = train[train['Patient'] == patient_id]\n    x = data['Weeks']\n    y = data['FVC']\n    ax.set_title(patient_id)\n    ax = sns.regplot(x, y, ax=ax, ci=None, line_kws={'color':'red'})\n    \n\nf, axes = plt.subplots(1, 3, figsize=(15, 5))\nchart('ID00007637202177411956430', axes[0])\nchart('ID00009637202177434476278', axes[1])\nchart('ID00010637202177584971671', axes[2])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3. Model A: Bayesian Hierarchical Linear Regression with Partial Pooling\nOur first model (A) is the simplest possible Bayesian Hierarchical Linear Regression model. It assumes **every patient has a personalized FVC decline curve**, with a single intercept ($\\alpha$) and slope ($\\beta$).\n\nActually the simplest possible linear regression, not hierarchical, would assume all FVC decline curves have the same $\\alpha$ and $\\beta$. That's the **pooled model**. In the other extreme, we could assume a model where each patient has a personalized FVC decline curve, **and these curves are completely unrelated**. That's the **unpooled model**, where each patient has completely separate regressions.\n\nHere, we will use the middle ground: **Partial pooling**. Specifically, we will assume that while $\\alpha$'s and $\\beta$'s are different for each patient as in the unpooled case, the **coefficients all share similarity**. We can model this by assuming that each individual coefficient comes from a common group distribution. For more details about these 3 different ways of modelling, check [GLM: Hierarchical Linear Regression](https://docs.pymc.io/notebooks/GLM-hierarchical.html) great tutorial. The graphic model below represents this model:\n\n<img src=\"https://i.ibb.co/H7NgBfR/Artboard-2-2x-100.jpg\" alt=\"drawing\" width=\"600\"/>\n\nMathematically, the model is described by the following equations:\n$$\n\\mu_{\\alpha} \\sim \\mathcal{N}(1700, 400) \\\\\n\\sigma_{\\alpha} \\sim |\\mathcal{N}(0, 1000)| \\\\\n\\mu_{\\beta} \\sim \\mathcal{N}(-4, 1) \\\\\n\\sigma_{\\beta} \\sim |\\mathcal{N}(0, 5)| \\\\\n\\alpha_i \\sim \\mathcal{N}(\\mu_{\\alpha}, \\sigma_{\\alpha}) \\\\\n\\beta_i \\sim \\mathcal{N}(\\mu_{\\beta}, \\sigma_{\\beta}) \\\\\n\\sigma \\sim |\\mathcal{N}(0, 150)| \\\\\nFVC_{ij} \\sim \\mathcal{N}(\\alpha_i + t \\beta_i, \\sigma)\n$$\n\nwhere *t* is the time in weeks. You might ask why we have chosen these specific priors. We did some iterations, at first with pretty vague/uninformative priors. After observing the first results, we converged to these distributions described above. Although PyMC3 samplers can work pretty well with uninformative priors, other PPLs such as Pyro requires more specific priors (otherwise their samplers take a really long time to converge).\n\nImplementing it in PyMC3 is very straightforward:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import pymc3 as pm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_patients = train['Patient'].nunique()\nFVC_obs = train['FVC'].values\nWeeks = train['Weeks'].values\nPatientID = train['PatientID'].values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with pm.Model() as model_a:\n    # create shared variables that can be changed later on\n    FVC_obs_shared = pm.Data(\"FVC_obs_shared\", FVC_obs)\n    Weeks_shared = pm.Data('Weeks_shared', Weeks)\n    PatientID_shared = pm.Data('PatientID_shared', PatientID)\n    \n    mu_a = pm.Normal('mu_a', mu=1700., sigma=400)\n    sigma_a = pm.HalfNormal('sigma_a', 1000.)\n    mu_b = pm.Normal('mu_b', mu=-4., sigma=1)\n    sigma_b = pm.HalfNormal('sigma_b', 5.)\n\n    a = pm.Normal('a', mu=mu_a, sigma=sigma_a, shape=n_patients)\n    b = pm.Normal('b', mu=mu_b, sigma=sigma_b, shape=n_patients)\n\n    # Model error\n    sigma = pm.HalfNormal('sigma', 150.)\n\n    FVC_est = a[PatientID_shared] + b[PatientID_shared] * Weeks_shared\n\n    # Data likelihood\n    FVC_like = pm.Normal('FVC_like', mu=FVC_est,\n                         sigma=sigma, observed=FVC_obs_shared)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 4. Fiting model A\nJust press the inference button (TM)! :)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Inference button (TM)!\nwith model_a:\n    trace_a = pm.sample(2000, tune=2000, target_accept=.9, init=\"adapt_diag\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 5. Checking model A\n## 5.1. Inspecting the learned parameters\nFirst, let's inspect the parameters learned:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_a:\n    pm.traceplot(trace_a);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Great! It looks like our model learned different $\\alpha$'s and $\\beta$'s for each patient, pooled from the same source distribution.\n\n## 5.2. Visualizing FVC decline curves for some patients\nNow, let's visually inspect FVC decline curves predicted by our model. We will completely fill in the FVC table, predicting all missing values. The first step is to create a table to fill:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_template = []\nfor i in range(train['Patient'].nunique()):\n    df = pd.DataFrame(columns=['PatientID', 'Weeks'])\n    df['Weeks'] = np.arange(-12, 134)\n    df['PatientID'] = i\n    pred_template.append(df)\npred_template = pd.concat(pred_template, ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, let's generate all predictions and fill in the table:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_a:\n    pm.set_data({\n        \"PatientID_shared\": pred_template['PatientID'].values.astype(int),\n        \"Weeks_shared\": pred_template['Weeks'].values.astype(int),\n        \"FVC_obs_shared\": np.zeros(len(pred_template)).astype(int),\n    })\n    post_pred = pm.sample_posterior_predictive(trace_a)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.DataFrame(columns=['Patient', 'Weeks', 'FVC_pred', 'sigma'])\ndf['Patient'] = le_id.inverse_transform(pred_template['PatientID'])\ndf['Weeks'] = pred_template['Weeks']\ndf['FVC_pred'] = post_pred['FVC_like'].T.mean(axis=1)\ndf['sigma'] = post_pred['FVC_like'].T.std(axis=1)\ndf['FVC_inf'] = df['FVC_pred'] - df['sigma']\ndf['FVC_sup'] = df['FVC_pred'] + df['sigma']\ndf = pd.merge(df, train[['Patient', 'Weeks', 'FVC']], how='left', on=['Patient', 'Weeks'])\ndf = df.rename(columns={'FVC': 'FVC_true'})\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, let's visualize the predictions for 6 patients:\n- In the first row, we can see predictions for 3 patients for which our model had several points to learn personalized $\\alpha$'s and $\\beta$'s\n- In the second row, we can see predictions for 3 patients for which our model had only a single point to learn personalized $\\alpha$'s and $\\beta$'s","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def chart(patient_id, ax):\n    data = df[df['Patient'] == patient_id]\n    x = data['Weeks']\n    ax.set_title(patient_id)\n    ax.plot(x, data['FVC_true'], 'o')\n    ax.plot(x, data['FVC_pred'])\n    ax = sns.regplot(x, data['FVC_true'], ax=ax, ci=None, \n                     line_kws={'color':'red'})\n    ax.fill_between(x, data[\"FVC_inf\"], data[\"FVC_sup\"],\n                    alpha=0.5, color='#ffcd3c')\n    ax.set_ylabel('FVC')\n\nf, axes = plt.subplots(2, 3, figsize=(15, 10))\nchart('ID00007637202177411956430', axes[0, 0])\nchart('ID00009637202177434476278', axes[0, 1])\nchart('ID00011637202177653955184', axes[0, 2])\nchart('ID00419637202311204720264', axes[1, 0])\nchart('ID00421637202311550012437', axes[1, 1])\nchart('ID00422637202311677017371', axes[1, 2])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The results are exactly what we expected to see! Highlight observations:\n- The model adequately learned Bayesian Linear Regressions! The orange line (learned predicted FVC mean) is very inline with the red line (deterministic linear regression). But most important: it learned to predict uncertainty, showed in the light orange region (one sigma above and below the mean FVC line)\n- In the first row, where the model had several points for each patient, we can see the model predicts a higher uncertainty where the data points are more disperse (1st and 3rd patients). Conversely, where the points are more grouped together (2nd patient, first row), the model predicts a higher confidence (narrower light orange region)\n- In the 2nd row, we can see our model does a good job even when we supply it with only a single data point per patient to use as base to make inferences!\n- Comparing the first 3 patients in the first row with the last 3 patients on the 2nd row, we can see that our model correctly estimates a higher confidence for the first 3 patients, for which it had more data points to make inferences!\n- Finally, in all 6 patients, we can see that the uncertainty grows as the look more into the future: the light orange region widens as the # of weeks grow! That makes a lot of sense!\n\n## 5.3. Computing the modified Laplace Log Likelihood and RMSE\nThese are straightforward to compute:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Add the test data points back to calculate the metrics\ntr = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ndf = pd.merge(df, tr[['Patient', 'Weeks', 'FVC']], 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']) ** 2).mean() ** (1/2)\nprint(f'RMSE: {rmse:.1f} ml')\n\nsigma_c = y['sigma'].values\nsigma_c[sigma_c < 70] = 70\ndelta = (y['FVC_pred'] - y['FVC']).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":"These look like very encouraging numbers, especially taking into consideration that **we only used 3 columns: FVC, weeks and Patient ID**, discarding all the rest. Now, let's go to the 2nd objective.\n\n# 6. Extending the model to use more tabular data\nWe will implement very few changes to make our model use all possible tabular data. The first is simply one-hot encode the categorical features **Sex** and **SmokingStatus**. In order to avoid multicollinearity issues, we will drop the last column after transforming both. (For more about multicollinearity issues, check [here](https://www.analyticsvidhya.com/blog/2020/03/what-is-multicollinearity/)). The second is adding **percent at baseline** feature.","execution_count":null},{"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['ExSmoker'] = aux['ExSmoker']\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')\n\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And that's all there is to it.\n\n# 7. Model B: Hierarchical Model with Partial Pooling using more tabular data\nThe model is essentially the same as model A. There is, though, a slight change. The equation that predicts the FVC is changed to:\n\n$$\nFVC_{ij} \\sim \\mathcal{N}(\\alpha_i + X_{ij} \\beta_i, \\sigma)\n$$\n\n\nNow, multiplying $\\beta_i$ we have $X_{ij}$ instead of *t*. $X_{ij}$ is a vector of the patient *i* at timestep *j* that contains:\n- Male binary code from patient *i* at timestep *j* (which is constant at all timesteps)\n- Ex-smoker binary code from patient *i* at timestep *j* (which is constant at all timesteps)\n- Currently-smokes binary code from patient *i* at timestep *j* (which is constant at all timesteps)\n- Percent from patient *i* at baseline moment\n- The week (timestep *j*)\n\nSome good questions and our best answers:\n- Why aren't you including **age**? Empirically, when we add age to the features, it significantly degrades the model performance. Don't know exactly why this happens. We hypothesize that this is due to multicollinearity. Anyways, in a model trained to predict personalized curves of FVC decline, in which we inform the Patient ID, in theory everything that qualifies the patient should be already encoded in its ID.\n- Why are you considering **Percent at baseline moment** only, discarding the other values? That's a very good question :) The percent feature varies in time, and it is strongly correlated with FVC (the variable we are interested in predicting). However, at inference time, we don't have percent feature in every week from [-12, 133]; we only have the percent at baseline moment. At inference time, if we consider percent flat from week -12 to week 133 (as we saw in some notebooks), that would be conceptually wrong IMHO. The only number we have is the percent at baseline moment. So, that's the number we feed our model.\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"n_patients = train['Patient'].nunique()\nFVC_obs = train['FVC'].values\nX = train[['Weeks', 'Male', 'ExSmoker', 'CurrentlySmokes', \n           'Percent_base']].values\nPatientID = train['PatientID'].values","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with pm.Model() as model_b:\n    # create shared variables that can be changed later on\n    FVC_obs_shared = pm.Data(\"FVC_obs_shared\", FVC_obs)\n    X_shared = pm.Data('X_shared', X)\n    PatientID_shared = pm.Data('PatientID_shared', PatientID)\n    \n    mu_a = pm.Normal('mu_a', mu=1700, sigma=400)\n    sigma_a = pm.HalfNormal('sigma_a', 1000.)\n    mu_b = pm.Normal('mu_b', mu=-4., sigma=1., shape=X.shape[1])\n    sigma_b = pm.HalfNormal('sigma_b', 5.)\n\n    a = pm.Normal('a', mu=mu_a, sigma=sigma_a, shape=n_patients)\n    b = pm.Normal('b', mu=mu_b, sigma=sigma_b, shape=(n_patients, X.shape[1]))\n\n    # Model error\n    sigma = pm.HalfNormal('sigma', 150.)\n\n    FVC_est = a[PatientID_shared] + (b[PatientID_shared] * X_shared).sum(axis=1)\n\n    # Data likelihood\n    FVC_like = pm.Normal('FVC_like', mu=FVC_est,\n                         sigma=sigma, observed=FVC_obs_shared)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 8. Fiting model B\nJust press the inference button (TM)! :)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Inference button (TM)!\nwith model_b:\n    trace_b = pm.sample(2000, tune=2000, target_accept=.9, init=\"adapt_diag\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 9. Checking model B\n## 9.1. Inspecting the learned parameters\nFirst, let's inspect the parameters learned:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_b:\n    pm.traceplot(trace_b);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Great! Again, it looks like our model learned different $\\alpha$'s and $\\beta$'s for each patient, pooled from the same source distribution. Note: this model is unstable, and sometimes do not converge.\n\n## 9.2. Visualizing FVC decline curves for some patients\nNow, let's visually inspect FVC decline curves predicted by our model. We will completely fill in the FVC table, predicting all missing values. The first step is to create a table to fill:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"aux = train.groupby('Patient').first().reset_index()\npred_template = []\nfor i in range(train['Patient'].nunique()):\n    df = pd.DataFrame(columns=['PatientID', 'Weeks'])\n    df['Weeks'] = np.arange(-12, 134)\n    df['PatientID'] = i\n    df['Male'] = aux[aux['PatientID'] == i]['Male'].values[0]\n    df['ExSmoker'] = aux[aux['PatientID'] == i]['ExSmoker'].values[0]\n    df['CurrentlySmokes'] = aux[aux['PatientID'] == i]['CurrentlySmokes'].values[0]\n    df['Percent_base'] = aux[aux['PatientID'] == i]['Percent_base'].values[0]\n    pred_template.append(df)\npred_template = pd.concat(pred_template, ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, let's fill the table:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"with model_b:\n    pm.set_data({\n        \"PatientID_shared\": pred_template['PatientID'].values.astype(int),\n        \"X_shared\": pred_template[['Weeks', 'Male', 'ExSmoker', \n                                   'CurrentlySmokes', \n                                   'Percent_base']].values.astype(int),\n        \"FVC_obs_shared\": np.zeros(len(pred_template)).astype(int),\n    })\n    post_pred = pm.sample_posterior_predictive(trace_b)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.DataFrame(columns=['Patient', 'Weeks', 'FVC_pred', 'sigma'])\ndf['Patient'] = le_id.inverse_transform(pred_template['PatientID'])\ndf['Weeks'] = pred_template['Weeks']\ndf['FVC_pred'] = post_pred['FVC_like'].T.mean(axis=1)\ndf['sigma'] = post_pred['FVC_like'].T.std(axis=1)\ndf['FVC_inf'] = df['FVC_pred'] - df['sigma']\ndf['FVC_sup'] = df['FVC_pred'] + df['sigma']\ndf = pd.merge(df, train[['Patient', 'Weeks', 'FVC']], how='left', on=['Patient', 'Weeks'])\ndf = df.rename(columns={'FVC': 'FVC_true'})\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And finally, let's visualize again 6 patients, the same 6 we saw with Model A:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"def chart(patient_id, ax):\n    data = df[df['Patient'] == patient_id]\n    x = data['Weeks']\n    ax.set_title(patient_id)\n    ax.plot(x, data['FVC_true'], 'o')\n    ax.plot(x, data['FVC_pred'])\n    ax = sns.regplot(x, data['FVC_true'], ax=ax, ci=None, \n                     line_kws={'color':'red'})\n    ax.fill_between(x, data[\"FVC_inf\"], data[\"FVC_sup\"],\n                    alpha=0.5, color='#ffcd3c')\n    ax.set_ylabel('FVC')\n\nf, axes = plt.subplots(2, 3, figsize=(15, 10))\nchart('ID00007637202177411956430', axes[0, 0])\nchart('ID00009637202177434476278', axes[0, 1])\nchart('ID00011637202177653955184', axes[0, 2])\nchart('ID00419637202311204720264', axes[1, 0])\nchart('ID00421637202311550012437', axes[1, 1])\nchart('ID00422637202311677017371', axes[1, 2])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It looks like Model B generates the same results as Model A. Note: this model is unstable, and sometimes do not converge.\n\n## 9.3. Computing the modified Laplace Log Likelihood and RMSE\nFinally, let's check the metrics for model B:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Add the test data points back to calculate the metrics\ntr = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ndf = pd.merge(df, tr[['Patient', 'Weeks', 'FVC']], 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'RMSE: {rmse:.1f} ml')\n\nsigma_c = y['sigma'].values\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":"Again, we see good encouraging numbers (Note: this model is unstable, and sometimes do not converge.). However, these are remarkably close to the numbers we obtained when we used only FVC, weeks and Patient ID - actually a little bit worse. Why is that the case? **Why aren't these new features helping the model?** We discuss this question and other learnings in the next section.\n\n# 10. Discussion\nOur first objective was to **build a stable error-free simplest possible Bayesian Hierarchical Linear Regression model**. We achieved it, implementing a model with partial pooling, where **each patient has a personalized FVC decline curve with personalized parameters**, but **all these parameters share similarities**, being drawn from the same shared distributions.\n\nThe results observed from the experiment with model A met our expectations. The model adequately learned not only to predict good FVC declines (low RMSE) but also **learned to properly gauge confidence** (low Laplace Log Likelihood). We saw that our model effectively made predictions **even in the challenging situation where it had only a single FVC measure** to base its future FVC forecasts. We visually saw **higher predicted confidence when the FVC data points were grouped together** in a more clear line, and a **lower predicted confidence when the FVC measurements were more disperse**. Our model was also able to properly predict **decreasing confidence as we looked more and more into the future**.\n\nWhile experimenting with model B, our 2nd objective, we had to **deal with multicollinearity**, a common problem of linear regression models in high dimensions. After dropping some features to avoid it, we were able to observe the results, which were very similar than the ones obtained with model A. So, with this modelling, **new features did not improve our predictive error**.\n\nWe hypothesise that this happens because **all the information that qualifies a patient is already encoded in the Patient ID**. Once we feed the model with this information, nothing else matters: it is the maximum of personalisation. Any other feature that groups patients are, by definition, more general than the Patient ID (the maximum of personalisation). This refers back to the original idea that sparked this line of thought: **collaborative filtering**. There are basically 2 approaches to mass-customize recommendations (adopted by giants such as Amazon, Netflix, Facebook, etc): collaborative filtering and content-base filtering. The former learns to group similar users based on their ratings (nothing else matters), while the latter recommends based on user & item’s features. So, the fact that model B produce similar scores of model A, even with more features, is **consistent with the adopted approach**.\n\nMaybe with a larger base of patients, more features would help the model to group them together. However, only 176 patients is too few. \n\n## Occam’s Razor\n\nNow, let’s re-train the model using all data, and submit a solution. **But which model, A or B?**\n\nThe principle of simplicity in scientific models and theories is commonly called Ockham's razor, or Occham's razor. It is popularly attributed to 1400s English friar and philosopher William of Ockham, also known as William of Occham. The razor alludes to the shaving away of unneeded detail. A common paraphrase of Ockham's principle, originally written in Latin, is **\"All things being equal, the simplest solution tends to be the best one.\"**\n\nThis is such a popular concept that we will leave [Lisa Simpson’s explanation of Occam’s Razor principle](https://www.youtube.com/watch?v=Ly0YzGpi63M) as our final remark :)\n\n# 11. Final submission","execution_count":null},{"metadata":{"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')\ntrain = pd.concat([train, test], axis=0, ignore_index=True)\\\n    .drop_duplicates()\nle_id = LabelEncoder()\ntrain['PatientID'] = le_id.fit_transform(train['Patient'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_patients = train['Patient'].nunique()\nFVC_obs = train['FVC'].values\nWeeks = train['Weeks'].values\nPatientID = train['PatientID'].values\n\nwith pm.Model() as model_a:\n    # create shared variables that can be changed later on\n    FVC_obs_shared = pm.Data(\"FVC_obs_shared\", FVC_obs)\n    Weeks_shared = pm.Data('Weeks_shared', Weeks)\n    PatientID_shared = pm.Data('PatientID_shared', PatientID)\n    \n    mu_a = pm.Normal('mu_a', mu=1700., sigma=400)\n    sigma_a = pm.HalfNormal('sigma_a', 1000.)\n    mu_b = pm.Normal('mu_b', mu=-4., sigma=1)\n    sigma_b = pm.HalfNormal('sigma_b', 5.)\n\n    a = pm.Normal('a', mu=mu_a, sigma=sigma_a, shape=n_patients)\n    b = pm.Normal('b', mu=mu_b, sigma=sigma_b, shape=n_patients)\n\n    # Model error\n    sigma = pm.HalfNormal('sigma', 150.)\n\n    FVC_est = a[PatientID_shared] + b[PatientID_shared] * Weeks_shared\n\n    # Data likelihood\n    FVC_like = pm.Normal('FVC_like', mu=FVC_est,\n                         sigma=sigma, observed=FVC_obs_shared)\n    \n    # Fitting the model\n    trace_a = pm.sample(2000, tune=2000, target_accept=.9, init=\"adapt_diag\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_template = []\nfor p in test['Patient'].unique():\n    df = pd.DataFrame(columns=['PatientID', 'Weeks'])\n    df['Weeks'] = np.arange(-12, 134)\n    df['Patient'] = p\n    pred_template.append(df)\npred_template = pd.concat(pred_template, ignore_index=True)\npred_template['PatientID'] = le_id.transform(pred_template['Patient'])\n\nwith model_a:\n    pm.set_data({\n        \"PatientID_shared\": pred_template['PatientID'].values.astype(int),\n        \"Weeks_shared\": pred_template['Weeks'].values.astype(int),\n        \"FVC_obs_shared\": np.zeros(len(pred_template)).astype(int),\n    })\n    post_pred = pm.sample_posterior_predictive(trace_a)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.DataFrame(columns=['Patient', 'Weeks', 'Patient_Week', 'FVC', 'Confidence'])\ndf['Patient'] = pred_template['Patient']\ndf['Weeks'] = pred_template['Weeks']\ndf['Patient_Week'] = df['Patient'] + '_' + df['Weeks'].astype(str)\ndf['FVC'] = post_pred['FVC_like'].T.mean(axis=1)\ndf['Confidence'] = post_pred['FVC_like'].T.std(axis=1)\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":"**Next step:**\n- Replace the linear model by a non-linear model... isn't that obvious?! :)\n\n<img src=\"https://upload.wikimedia.org/wikipedia/commons/1/18/Bayes%27_Theorem_MMB_01.jpg\" alt=\"drawing\" width=\"800\"/>","execution_count":null},{"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}