{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":20604,"databundleVersionId":1357052,"sourceType":"competition"}],"dockerImageVersionId":29994,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-02-15T07:47:03.667579Z","iopub.execute_input":"2024-02-15T07:47:03.667914Z","iopub.status.idle":"2024-02-15T07:47:03.672858Z","shell.execute_reply.started":"2024-02-15T07:47:03.667884Z","shell.execute_reply":"2024-02-15T07:47:03.671423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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')","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","execution":{"iopub.status.busy":"2024-02-15T07:47:06.388030Z","iopub.execute_input":"2024-02-15T07:47:06.388355Z","iopub.status.idle":"2024-02-15T07:47:06.405173Z","shell.execute_reply.started":"2024-02-15T07:47:06.388326Z","shell.execute_reply":"2024-02-15T07:47:06.403833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head(3\n          )","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:08.678247Z","iopub.execute_input":"2024-02-15T07:47:08.678768Z","iopub.status.idle":"2024-02-15T07:47:08.692059Z","shell.execute_reply.started":"2024-02-15T07:47:08.678733Z","shell.execute_reply":"2024-02-15T07:47:08.690940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test.head(3)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:11.293069Z","iopub.execute_input":"2024-02-15T07:47:11.293605Z","iopub.status.idle":"2024-02-15T07:47:11.305473Z","shell.execute_reply.started":"2024-02-15T07:47:11.293560Z","shell.execute_reply":"2024-02-15T07:47:11.304293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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('ID00419637202311204720264', axes[0])\nchart('ID00009637202177434476278', axes[1])\nchart('ID00010637202177584971671', axes[2])","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:13.688140Z","iopub.execute_input":"2024-02-15T07:47:13.688682Z","iopub.status.idle":"2024-02-15T07:47:14.039421Z","shell.execute_reply.started":"2024-02-15T07:47:13.688632Z","shell.execute_reply":"2024-02-15T07:47:14.038741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The decline in lung capacity is very clear. We see, though, they are very different from patient to patient.\n\n\"Many large data sets are in fact large collections of small data sets. For example, in areas such as **personalized medicine** and recommendation systems, there might be a large amount of data, but there is still a relatively **small amount of data for each patient** or client, respectively. To **customize predictions for each person** it becomes necessary to **build a model for each person** — with its inherent **uncertainties** — and to couple these models together in a hierarchy so that **information can be borrowed from other similar people**. We call this the **personalization of models**, and it is naturally implemented using **hierarchical Bayesian approaches** (...).\" (Ghahramani, 2015)\n\nThis is exactly what we will try here. Moving on.","metadata":{}},{"cell_type":"markdown","source":"# 3. Postulate the model\nTime to be creative. There are a multitude of ways we could model this tabular dataset. Some of the tools we could use:\n- Hidden Markov Models\n- Gaussian Processes\n- Variational Auto Encoders\n\nAs I am learning, I will try first the simplest possible model: a **linear regression**.\nHowever, we will sophisticate a little bit. Here are our assumptions:\n- Every patient has unique linear regression parameters ($\\alpha$ and $\\beta$). So, by inferring the right parameters, we will be able to predict the line(s) for each patient, thus being able to predict his FVC in any week\n- However, these parameters are not completely independent. There is an underlying model that governs them for all patients.\n- Both $\\alpha$ and $\\beta$ are normally distributed with different means and variances\n- These means and variances are functions of the baseline measure (baseline week, FVC and Percent), and patient's age, sex and smoking status\n- In the next notebook, we will sophisticate even further, by assuming the parameters are also function of latent variables learned from the CT scans. But that will come later. Baby steps :)\n\nOur model is represented by the following Bayesian Network:\n\n<img src=\"https://i.ibb.co/DCxKbdT/Asset-1-2x-100.jpg\" alt=\"drawing\" width=\"600\"/>\n\nLet me explain the logic behind this model:\n- $FVC_{ij}$ is the observed variable we are interested in. At any week $j$, $-12 \\leq j \\leq 133$, the FVC of the patient $i$ is presumed to be normally distributed with mean $\\alpha_i + \\beta_i i$ and $\\sigma_i^2$ (the confidence asked)\n- $\\alpha_i$, the intercept of the decline function for each patient $i$, logically is a function of $FVC_i^b$ (the baseline measurement for patient $i$) and $w_i^b$ (the week when the baseline FVC was measured). We assume it is normally distributed with mean $FVC_i^b + w_i^b \\beta^{int}$ and variance $\\sigma_{int}^2$ (int is superscript, but I couldn't get latex to behave rs)\n- $\\beta_i$, the slope of the decline function for each patient $i$, logically is a function of $A_i$ (patient's age), sex and smoking status. We assume it is normally distributed with mean $\\alpha^s + A_i \\beta_c^s$, with variance $\\sigma_s^2$ (again, s should be superscript). We considered 6 different $\\beta_c^s$: for women who currently smoke, men who currently smoke, women ex-smokers, men ex-mokers, women who never smoked and men who never smoked.\n- For now, to simplify, we left Percent random variable out. We will include in a second version.\n- Finally, we know nothing about the priors $\\beta^{int}$, $\\alpha^s$, $\\sigma_i$, $\\sigma^{int}$ and $\\sigma^s$. We will model the first 2 as normals, and the last 3 as half-normals.\n\nMathematically, the model specification is\n$$\nFVC_{ij} \\sim \\mathcal{N}(\\alpha_i + j \\beta_i, \\sigma_i) \\\\\n\\sigma_i \\sim |\\mathcal{N}(0, 200)| \\\\\n\\alpha_i \\sim \\mathcal{N}(FVC_i^b + w_i^b \\beta^{int}, \\sigma^{int}) \\\\\n\\beta_i \\sim \\mathcal{N}(\\alpha^s + A_i \\beta_c^s, \\sigma^s)\\\\\n\\beta^{int} \\sim \\mathcal{N}(0, 100) \\\\\n\\sigma^{int} \\sim |\\mathcal{N}(0, 100)| \\\\\n\\beta_c^s \\sim \\mathcal{N}(0, 100) \\\\\n\\alpha^s \\sim \\mathcal{N}(0, 100) \\\\\n\\sigma^s \\sim |\\mathcal{N}(0, 100)|\n$$","metadata":{}},{"cell_type":"code","source":"# Kaggle, please add Pyro/PyTorch support!\nimport pymc3 as pm\nimport theano\nimport arviz as az\nfrom sklearn import preprocessing","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:20.312611Z","iopub.execute_input":"2024-02-15T07:47:20.312988Z","iopub.status.idle":"2024-02-15T07:47:20.317236Z","shell.execute_reply.started":"2024-02-15T07:47:20.312949Z","shell.execute_reply":"2024-02-15T07:47:20.316350Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.1. Simple data prep","metadata":{}},{"cell_type":"code","source":"# Very simple pre-processing: adding patient class\ndef patient_class(row):\n    if row['Sex'] == 'Male':\n        if row['SmokingStatus'] == 'Currently smokes':\n            return 0\n        elif row['SmokingStatus'] == 'Ex-smoker':\n            return 1\n        elif row['SmokingStatus'] == 'Never smoked':\n            return 2\n    else:\n        if row['SmokingStatus'] == 'Currently smokes':\n            return 3\n        elif row['SmokingStatus'] == 'Ex-smoker':\n            return 4\n        elif row['SmokingStatus'] == 'Never smoked':\n            return 5\n\ntrain['Class'] = train.apply(patient_class, axis=1)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:24.838041Z","iopub.execute_input":"2024-02-15T07:47:24.838391Z","iopub.status.idle":"2024-02-15T07:47:24.872928Z","shell.execute_reply.started":"2024-02-15T07:47:24.838357Z","shell.execute_reply":"2024-02-15T07:47:24.871776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train['Patient'].nunique()","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:27.787538Z","iopub.execute_input":"2024-02-15T07:47:27.787901Z","iopub.status.idle":"2024-02-15T07:47:27.793966Z","shell.execute_reply.started":"2024-02-15T07:47:27.787863Z","shell.execute_reply":"2024-02-15T07:47:27.793119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"aux = train[['Patient', 'Weeks']].groupby('Patient')\\\n    .min().reset_index()\naux = pd.merge(aux, train[['Patient', 'Weeks', 'FVC']], how='left', \n               on=['Patient', 'Weeks'])\n\n# aux = aux.groupby('Patient').mean().reset_index()\na =aux.groupby(['Patient'],as_index=False)['Weeks'].count()\naux[aux['Patient']=='ID00048637202185016727717']","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:30.263033Z","iopub.execute_input":"2024-02-15T07:47:30.263386Z","iopub.status.idle":"2024-02-15T07:47:30.289002Z","shell.execute_reply.started":"2024-02-15T07:47:30.263350Z","shell.execute_reply":"2024-02-15T07:47:30.287773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Very simple pre-processing: adding FVC and week baselines\naux = train[['Patient', 'Weeks']].groupby('Patient')\\\n    .min().reset_index()\naux = pd.merge(aux, train[['Patient', 'Weeks', 'FVC']], how='left', \n               on=['Patient', 'Weeks'])\naux = aux.groupby('Patient').mean().reset_index()\naux['Weeks'] = aux['Weeks'].astype(int)\naux['FVC'] = aux['FVC'].astype(int)\naux\ntrain = pd.merge(train, aux, how='left', on='Patient', suffixes=('', '_base'))","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:32.768052Z","iopub.execute_input":"2024-02-15T07:47:32.768380Z","iopub.status.idle":"2024-02-15T07:47:32.789981Z","shell.execute_reply.started":"2024-02-15T07:47:32.768349Z","shell.execute_reply":"2024-02-15T07:47:32.788804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:35.803146Z","iopub.execute_input":"2024-02-15T07:47:35.803519Z","iopub.status.idle":"2024-02-15T07:47:35.827554Z","shell.execute_reply.started":"2024-02-15T07:47:35.803477Z","shell.execute_reply":"2024-02-15T07:47:35.826604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Very simple pre-processing: creating patient indexes\nle = preprocessing.LabelEncoder()\ntrain['PatientID'] = le.fit_transform(train['Patient'])\n\npatients = train[['Patient', 'PatientID', 'Age', 'Class', 'Weeks_base', 'FVC_base']].drop_duplicates()\nfvc_data = train[['Patient', 'PatientID', 'Weeks', 'FVC']]\n\n# patients","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:38.922680Z","iopub.execute_input":"2024-02-15T07:47:38.922988Z","iopub.status.idle":"2024-02-15T07:47:38.935363Z","shell.execute_reply.started":"2024-02-15T07:47:38.922962Z","shell.execute_reply":"2024-02-15T07:47:38.933766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fvc_data.head()","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:41.557614Z","iopub.execute_input":"2024-02-15T07:47:41.557959Z","iopub.status.idle":"2024-02-15T07:47:41.568601Z","shell.execute_reply.started":"2024-02-15T07:47:41.557930Z","shell.execute_reply":"2024-02-15T07:47:41.567807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.2. Modeling in PyMC3\nProbabilistic Programming Languages are very very very cool :)","metadata":{}},{"cell_type":"code","source":"len(fvc_data['Weeks'])\n","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:44.428121Z","iopub.execute_input":"2024-02-15T07:47:44.428487Z","iopub.status.idle":"2024-02-15T07:47:44.435396Z","shell.execute_reply.started":"2024-02-15T07:47:44.428451Z","shell.execute_reply":"2024-02-15T07:47:44.434191Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"FVC_b = patients['FVC_base'].values\nw_b = patients['Weeks_base'].values\nage = patients['Age'].values\npatient_class = patients['Class'].values\n\nt = fvc_data['Weeks'].values\nFVC_obs = fvc_data['FVC'].values\npatient_id = fvc_data['PatientID'].values\n\nwith pm.Model() as hierarchical_model:\n    # Hyperpriors for Alpha\n    beta_int = pm.Normal('beta_int', 0, sigma=100)\n    sigma_int = pm.HalfNormal('sigma_int', 100)\n    \n    # Alpha\n    mu_alpha = FVC_b + beta_int * w_b\n    alpha = pm.Normal('alpha', mu=mu_alpha, sigma=sigma_int, \n                      shape=train['Patient'].nunique())\n    \n    # Hyperpriors for Beta\n    sigma_s = pm.HalfNormal('sigma_s', 100)\n    alpha_s = pm.Normal('alpha_s', 0, sigma=100)\n    beta_cs = pm.Normal('beta_cs', 0, sigma=100, shape=6)\n    \n    # Beta\n    mu_beta = alpha_s + age * beta_cs[patient_class]\n    beta = pm.Normal('beta', mu=mu_beta, sigma=sigma_s,\n                     shape=train['Patient'].nunique())\n    \n    # Model variance\n    sigma = pm.HalfNormal('sigma', 200)\n    \n    # Model estimate\n    FVC_est = alpha[patient_id] + beta[patient_id] * t\n    \n    # Data likelihood\n    FVC_like = pm.Normal('FVC_like', mu=FVC_est,\n                          sigma=sigma, observed=FVC_obs)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:47:46.778196Z","iopub.execute_input":"2024-02-15T07:47:46.778630Z","iopub.status.idle":"2024-02-15T07:47:56.140887Z","shell.execute_reply.started":"2024-02-15T07:47:46.778599Z","shell.execute_reply":"2024-02-15T07:47:56.139497Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. Fit the model\nJust press the inference button (TM)! :)","metadata":{}},{"cell_type":"code","source":"# Inference button (TM)!\nwith hierarchical_model:\n    trace = pm.sample(2000, tune=2000, target_accept=.9)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:48:04.863640Z","iopub.execute_input":"2024-02-15T07:48:04.864000Z","iopub.status.idle":"2024-02-15T07:50:43.990180Z","shell.execute_reply.started":"2024-02-15T07:48:04.863967Z","shell.execute_reply":"2024-02-15T07:50:43.988908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We just sampled 4000 different models that explain the data! Very cool! :)","metadata":{}},{"cell_type":"markdown","source":"# 5. Check the model\nLet's see the generative model we've created.","metadata":{}},{"cell_type":"code","source":"with hierarchical_model:\n    pm.traceplot(trace);","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:51:02.924077Z","iopub.execute_input":"2024-02-15T07:51:02.924459Z","iopub.status.idle":"2024-02-15T07:51:15.435542Z","shell.execute_reply.started":"2024-02-15T07:51:02.924422Z","shell.execute_reply":"2024-02-15T07:51:15.434818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(trace['sigma'])","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:51:40.122725Z","iopub.execute_input":"2024-02-15T07:51:40.123254Z","iopub.status.idle":"2024-02-15T07:51:40.128168Z","shell.execute_reply.started":"2024-02-15T07:51:40.123220Z","shell.execute_reply":"2024-02-15T07:51:40.127359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Very cool!!! Looks like our model learned personalized alphas and betas for each patient!","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"## 5.1. Checking some patients\nPyMC3 comes with a very powerful visualization tool called [ArviZ](https://arviz-devs.github.io/arviz/index.html). However, I didn't figure out how to use yet... Let's use Matplotlib and Seaborn.","metadata":{}},{"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    x2 = np.arange(-12, 133, step=0.1)\n    \n    pid = patients[patients['Patient'] == patient_id]['PatientID'].values[0]\n#     print(pid)\n    for sample in range(100):\n        alpha = trace['alpha'][sample, pid]\n        beta = trace['beta'][sample, pid]\n        sigma = trace['sigma'][sample]\n        y2 = alpha + beta * x2\n        ax.plot(x2, y2, linewidth=0.1, color='green')\n        y2 = alpha + beta * x2 + sigma\n        ax.plot(x2, y2, linewidth=0.1, color='yellow')\n#         y2 = alpha + beta * x2 - sigma\n#         ax.plot(x2, y2, linewidth=0.1, color='yellow')\n\nf, axes = plt.subplots(1, 3, figsize=(15, 5))\nchart('ID00007637202177411956430', axes[0])\nchart('ID00009637202177434476278', axes[1])\nchart('ID00010637202177584971671', axes[2])","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:51:50.072881Z","iopub.execute_input":"2024-02-15T07:51:50.073247Z","iopub.status.idle":"2024-02-15T07:51:51.821681Z","shell.execute_reply.started":"2024-02-15T07:51:50.073203Z","shell.execute_reply":"2024-02-15T07:51:51.820972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here I plotted 100 out of the 4000 personalized models each patient has! In green we can see the fitted regression line, in yellow the standard deviation. Let's ensemble all that!","metadata":{}},{"cell_type":"markdown","source":"# 6. (Iterate and) Use the model\nLet's use our generative model now! Iteration will be left as an exercise to the reader :)","metadata":{}},{"cell_type":"markdown","source":"## 6.1. Simple data prep","metadata":{}},{"cell_type":"code","source":"# Very simple pre-processing: adding patient class\ndef patient_class(row):\n    if row['Sex'] == 'Male':\n        if row['SmokingStatus'] == 'Currently smokes':\n            return 0\n        elif row['SmokingStatus'] == 'Ex-smoker':\n            return 1\n        elif row['SmokingStatus'] == 'Never smoked':\n            return 2\n    else:\n        if row['SmokingStatus'] == 'Currently smokes':\n            return 3\n        elif row['SmokingStatus'] == 'Ex-smoker':\n            return 4\n        elif row['SmokingStatus'] == 'Never smoked':\n            return 5\n\ntest['Class'] = test.apply(patient_class, axis=1)\ntest = test.rename(columns={'FVC': 'FVC_base', 'Weeks': 'Weeks_base'})\ntest.head()","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:51:55.818071Z","iopub.execute_input":"2024-02-15T07:51:55.818622Z","iopub.status.idle":"2024-02-15T07:51:55.836290Z","shell.execute_reply.started":"2024-02-15T07:51:55.818588Z","shell.execute_reply":"2024-02-15T07:51:55.835703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# prepare submission dataset\nsubmission = []\nfor i, patient in enumerate(test['Patient'].unique()):\n    df = pd.DataFrame(columns=['Patient', 'Weeks', 'FVC'])\n    df['Weeks'] = np.arange(-12, 134)\n    df['Patient'] = patient\n    df['PatientID'] = i\n    df['FVC'] = 0\n    submission.append(df)\n    \nsubmission = pd.concat(submission).reset_index(drop=True)\nsubmission.head()","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:51:58.950002Z","iopub.execute_input":"2024-02-15T07:51:58.950328Z","iopub.status.idle":"2024-02-15T07:51:58.985976Z","shell.execute_reply.started":"2024-02-15T07:51:58.950294Z","shell.execute_reply":"2024-02-15T07:51:58.985008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 6.2. Posterior prediction\nThere are 2 ways of generating predictions on unseen held-out data using PyMC3. The first involves using `theano.shared` variables. It's pretty straightforward, 4-5 lines of code and we are done. I tried that, and although it worked perfectly while I was runnning the notebook, when I submitted Kaggle server complained, outputting **Submission CSV Not Found** error msg.\n\nMotivated by that, I will try the 2nd approach. It's a little bit longer than the 4-5 lines of code, but way more educational to me. The idea is outlined by PyMC3 developers in [this answer from Luciano Paz](https://discourse.pymc.io/t/how-do-we-predict-on-new-unseen-groups-in-a-hierarchical-model-in-pymc3/2571). To predict FVCs on hold-out data, we will create a 2nd model, using as priors the distributions for the parameters learned on the 1st model. It's Bayes spirit/philosophy: we keep constantly updating our models as we see more data. :)","metadata":{}},{"cell_type":"code","source":"trace['beta_cs'].mean()","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:52:03.042760Z","iopub.execute_input":"2024-02-15T07:52:03.043105Z","iopub.status.idle":"2024-02-15T07:52:03.049364Z","shell.execute_reply.started":"2024-02-15T07:52:03.043068Z","shell.execute_reply":"2024-02-15T07:52:03.048202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"FVC_b = test['FVC_base'].values\nw_b = test['Weeks_base'].values\nage = test['Age'].values\npatient_class = test['Class'].values\nt = submission['Weeks'].values\npatient_id = submission['PatientID'].values\n            \nwith pm.Model() as new_model:\n    # Hyperpriors for Alpha\n    beta_int = pm.Normal('beta_int', \n                         trace['beta_int'].mean(), \n                         sigma=trace['beta_int'].std())\n    sigma_int = pm.TruncatedNormal('sigma_int', \n                                   trace['sigma_int'].mean(),\n                                   sigma=trace['sigma_int'].std(),\n                                   lower=0)\n    \n    # Alpha\n    mu_alpha = FVC_b + beta_int * w_b\n    alpha = pm.Normal('alpha', mu=mu_alpha, sigma=sigma_int, \n                      shape=test['Patient'].nunique())\n    \n    # Hyperpriors for Beta\n    sigma_s = pm.TruncatedNormal('sigma_s', \n                                 trace['sigma_s'].mean(),\n                                 sigma=trace['sigma_s'].std(),\n                                 lower=0)\n    alpha_s = pm.Normal('alpha_s', \n                        trace['alpha_s'].mean(), \n                        sigma=trace['alpha_s'].std())\n    cov = np.zeros((6, 6))\n    np.fill_diagonal(cov, trace['beta_cs'].var(axis=0))\n    beta_cs = pm.MvNormal('beta_cs',\n                          mu=trace['beta_cs'].mean(axis=0),\n                          cov=cov,\n                          shape=6)\n    \n    # Beta\n    mu_beta = alpha_s + age * beta_cs[patient_class]\n    beta = pm.Normal('beta', mu=mu_beta, sigma=sigma_s,\n                     shape=test['Patient'].nunique())\n    \n    # Model variance\n    sigma = pm.TruncatedNormal('sigma', \n                               trace['sigma'].mean(),\n                               sigma=trace['sigma'].std(),\n                               lower=0)\n    \n    # Model estimate\n    # Here, there are 2 ways of estimating FVC. One is deterministic, the other\n    # stochastic. Assuming FVC is deterministic, we calculate sigma later, by\n    # evaluating std dev over the 4000 different models. This yields a higher \n    # confidence (lower sigmas). Assuming FVC is stochastic (commented code below)\n    # yields irregular lines. The mean FVC values are about the same, but the \n    # confidence is much lower (higher sigmas, about 2x the first case). Let's\n    # try submitting both cases, starting by the 1st assumption.\n    FVC_est = pm.Deterministic('FVC_est', alpha[patient_id] + beta[patient_id] * t)\n    \n    # sigma = pm.HalfNormal('sigma', 200)\n    # FVC_like = pm.Normal('FVC_like', mu=alpha[patient_id] + beta[patient_id] * t, \n    #                      sigma=sigma,\n    #                      shape=submission.shape[0])","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:52:05.182780Z","iopub.execute_input":"2024-02-15T07:52:05.183139Z","iopub.status.idle":"2024-02-15T07:52:16.380826Z","shell.execute_reply.started":"2024-02-15T07:52:05.183103Z","shell.execute_reply":"2024-02-15T07:52:16.379779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with new_model:\n    trace2 = pm.sample(2000, tune=2000, target_accept=.9)","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:52:20.169778Z","iopub.execute_input":"2024-02-15T07:52:20.170162Z","iopub.status.idle":"2024-02-15T07:53:07.574543Z","shell.execute_reply.started":"2024-02-15T07:52:20.170122Z","shell.execute_reply":"2024-02-15T07:53:07.573693Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trace2['FVC_est']","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:53:23.916415Z","iopub.execute_input":"2024-02-15T07:53:23.916855Z","iopub.status.idle":"2024-02-15T07:53:23.924581Z","shell.execute_reply.started":"2024-02-15T07:53:23.916826Z","shell.execute_reply":"2024-02-15T07:53:23.923703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trace2","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:53:26.189460Z","iopub.execute_input":"2024-02-15T07:53:26.190160Z","iopub.status.idle":"2024-02-15T07:53:26.195314Z","shell.execute_reply.started":"2024-02-15T07:53:26.190124Z","shell.execute_reply":"2024-02-15T07:53:26.194471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There we go! 4000 predictions for each point! Now, let's merge with the submission and calculate $FVC_{mean}$ and Confidence ($\\sigma$):","metadata":{}},{"cell_type":"markdown","source":"## 6.3. Generating final predictions","metadata":{}},{"cell_type":"code","source":"preds = pd.DataFrame(data=trace2['FVC_est'].T)\nsubmission = pd.merge(submission, preds, left_index=True, right_index=True)\nsubmission['Patient_Week'] = submission['Patient'] + '_' + submission['Weeks'].astype(str)\nsubmission = submission.drop(columns=['Patient', 'Weeks', 'FVC', 'PatientID'])\n\nFVC = submission.iloc[:, :-1].mean(axis=1)\nconfidence = submission.iloc[:, :-1].std(axis=1)\nsubmission['FVC'] = FVC\nsubmission['Confidence'] = confidence\nsubmission = submission[['Patient_Week', 'FVC', 'Confidence']]\nsubmission.to_csv('submission.csv', index=False)\nsubmission.head()","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:53:30.701593Z","iopub.execute_input":"2024-02-15T07:53:30.701934Z","iopub.status.idle":"2024-02-15T07:53:31.155969Z","shell.execute_reply.started":"2024-02-15T07:53:30.701904Z","shell.execute_reply":"2024-02-15T07:53:31.154962Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Define function to calculate patient class\ndef calculate_patient_class(sex, smoking_status):\n    if sex == 'Male':\n        if smoking_status == 'Currently smokes':\n            return 0\n        elif smoking_status == 'Ex-smoker':\n            return 1\n        elif smoking_status == 'Never smoked':\n            return 2\n    else:\n        if smoking_status == 'Currently smokes':\n            return 3\n        elif smoking_status == 'Ex-smoker':\n            return 4\n        elif smoking_status == 'Never smoked':\n            return 5\n\n# Define function to predict FVC for upcoming weeks\ndef predict_fvc_for_patient(patient_id, age, sex, smoking_status, current_week, upcoming_weeks):\n    # Check if patient exists in the test dataset\n    if patient_id not in test['Patient'].values:\n        print(f\"Patient with ID '{patient_id}' not found in the test dataset.\")\n        return None\n    \n    # Calculate patient class\n    patient_class = calculate_patient_class(sex, smoking_status)\n    \n    # Get patient index\n    patient_index = test[test['Patient'] == patient_id].index[0]\n    \n    # Predict FVC for upcoming weeks\n    predicted_fvc = []\n    for week in upcoming_weeks:\n        FVC_est = trace['alpha'][:, patient_index] + trace['beta'][:, patient_index] * week\n        FVC_est += trace['beta_int'] * current_week  # Adjust FVC for baseline week\n        FVC_est += trace['alpha_s'] + trace['beta_cs'][:, patient_class] * age  # Adjust FVC for age and class\n        predicted_fvc.append(FVC_est.mean())  # Taking mean FVC from posterior samples\n    \n    # Create a DataFrame for results\n    results_df = pd.DataFrame({\n        'Week': upcoming_weeks,\n        'Predicted_FVC': predicted_fvc\n    })\n    \n    return results_df\n\n# Example usage:\npatient_id = 'ID00421637202311550012437'\nage = 50\nsex = 'Male'\nsmoking_status = 'Ex-smoker'\ncurrent_week = 50\nupcoming_weeks = [51, 52, 53]  # Weeks for which we want to predict FVC\n\n# Call the function to predict FVC for upcoming weeks\nprediction_results = predict_fvc_for_patient(patient_id, age, sex, smoking_status, current_week, upcoming_weeks)\n\n# Display the prediction results\nif prediction_results is not None:\n    print(prediction_results)\n\n","metadata":{"execution":{"iopub.status.busy":"2024-02-15T07:53:33.651771Z","iopub.execute_input":"2024-02-15T07:53:33.652125Z","iopub.status.idle":"2024-02-15T07:53:33.674757Z","shell.execute_reply.started":"2024-02-15T07:53:33.652084Z","shell.execute_reply":"2024-02-15T07:53:33.673082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}