{"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":"2025-11-05T10:37:16.891369Z","iopub.execute_input":"2025-11-05T10:37:16.891752Z","iopub.status.idle":"2025-11-05T10:37:16.896751Z","shell.execute_reply.started":"2025-11-05T10:37:16.891719Z","shell.execute_reply":"2025-11-05T10:37:16.895426Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T10:37:16.904136Z","iopub.execute_input":"2025-11-05T10:37:16.904539Z","iopub.status.idle":"2025-11-05T10:37:16.923305Z","shell.execute_reply.started":"2025-11-05T10:37:16.904502Z","shell.execute_reply":"2025-11-05T10:37:16.920787Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train.head(3)","metadata":{"execution":{"iopub.status.busy":"2025-11-05T10:37:16.926112Z","iopub.execute_input":"2025-11-05T10:37:16.926533Z","iopub.status.idle":"2025-11-05T10:37:16.951206Z","shell.execute_reply.started":"2025-11-05T10:37:16.926494Z","shell.execute_reply":"2025-11-05T10:37:16.948627Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"test.head(3)","metadata":{"execution":{"iopub.status.busy":"2025-11-05T10:37:16.953473Z","iopub.execute_input":"2025-11-05T10:37:16.954040Z","iopub.status.idle":"2025-11-05T10:37:16.976399Z","shell.execute_reply.started":"2025-11-05T10:37:16.953985Z","shell.execute_reply":"2025-11-05T10:37:16.975081Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T10:42:46.877334Z","iopub.execute_input":"2025-11-05T10:42:46.877731Z","iopub.status.idle":"2025-11-05T10:42:47.761698Z","shell.execute_reply.started":"2025-11-05T10:42:46.877685Z","shell.execute_reply":"2025-11-05T10:42:47.760755Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T10:37:17.725959Z","iopub.execute_input":"2025-11-05T10:37:17.726326Z","iopub.status.idle":"2025-11-05T10:37:23.355581Z","shell.execute_reply.started":"2025-11-05T10:37:17.726291Z","shell.execute_reply":"2025-11-05T10:37:23.353762Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T10:42:47.763128Z","iopub.execute_input":"2025-11-05T10:42:47.763439Z","iopub.status.idle":"2025-11-05T10:42:47.803439Z","shell.execute_reply.started":"2025-11-05T10:42:47.763409Z","shell.execute_reply":"2025-11-05T10:42:47.802387Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train['Patient'].nunique()","metadata":{"execution":{"iopub.status.busy":"2025-11-05T10:37:23.406838Z","iopub.execute_input":"2025-11-05T10:37:23.407928Z","iopub.status.idle":"2025-11-05T10:37:23.422503Z","shell.execute_reply.started":"2025-11-05T10:37:23.407861Z","shell.execute_reply":"2025-11-05T10:37:23.421084Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T10:37:23.424122Z","iopub.execute_input":"2025-11-05T10:37:23.424641Z","iopub.status.idle":"2025-11-05T10:37:23.490527Z","shell.execute_reply.started":"2025-11-05T10:37:23.424585Z","shell.execute_reply":"2025-11-05T10:37:23.488151Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T10:37:23.492635Z","iopub.execute_input":"2025-11-05T10:37:23.493132Z","iopub.status.idle":"2025-11-05T10:37:23.528880Z","shell.execute_reply.started":"2025-11-05T10:37:23.492995Z","shell.execute_reply":"2025-11-05T10:37:23.527611Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train","metadata":{"execution":{"iopub.status.busy":"2025-11-05T10:37:23.530473Z","iopub.execute_input":"2025-11-05T10:37:23.530961Z","iopub.status.idle":"2025-11-05T10:37:23.562290Z","shell.execute_reply.started":"2025-11-05T10:37:23.530912Z","shell.execute_reply":"2025-11-05T10:37:23.560599Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T10:37:23.564228Z","iopub.execute_input":"2025-11-05T10:37:23.564707Z","iopub.status.idle":"2025-11-05T10:37:23.583263Z","shell.execute_reply.started":"2025-11-05T10:37:23.564666Z","shell.execute_reply":"2025-11-05T10:37:23.581274Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fvc_data.head()","metadata":{"execution":{"iopub.status.busy":"2025-11-05T10:37:23.585377Z","iopub.execute_input":"2025-11-05T10:37:23.586105Z","iopub.status.idle":"2025-11-05T10:37:23.614202Z","shell.execute_reply.started":"2025-11-05T10:37:23.586045Z","shell.execute_reply":"2025-11-05T10:37:23.612685Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"len(fvc_data['Weeks'])\n","metadata":{"execution":{"iopub.status.busy":"2025-11-05T10:37:23.615608Z","iopub.execute_input":"2025-11-05T10:37:23.616015Z","iopub.status.idle":"2025-11-05T10:37:23.642879Z","shell.execute_reply.started":"2025-11-05T10:37:23.615929Z","shell.execute_reply":"2025-11-05T10:37:23.641198Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T10:37:23.645355Z","iopub.execute_input":"2025-11-05T10:37:23.647482Z","iopub.status.idle":"2025-11-05T10:37:59.393053Z","shell.execute_reply.started":"2025-11-05T10:37:23.647017Z","shell.execute_reply":"2025-11-05T10:37:59.391360Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T10:45:02.105079Z","iopub.execute_input":"2025-11-05T10:45:02.105448Z","iopub.status.idle":"2025-11-05T10:47:03.740507Z","shell.execute_reply.started":"2025-11-05T10:45:02.105416Z","shell.execute_reply":"2025-11-05T10:47:03.738840Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"We just sampled 4000 different models that explain the data! Very cool! :)","metadata":{}},{"cell_type":"code","source":"len(trace['sigma'])","metadata":{"execution":{"iopub.status.busy":"2025-11-05T10:41:18.539342Z","iopub.execute_input":"2025-11-05T10:41:18.539683Z","iopub.status.idle":"2025-11-05T10:41:18.545813Z","shell.execute_reply.started":"2025-11-05T10:41:18.539649Z","shell.execute_reply":"2025-11-05T10:41:18.544867Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T11:26:33.454848Z","iopub.execute_input":"2025-11-05T11:26:33.455288Z","iopub.status.idle":"2025-11-05T11:26:33.481892Z","shell.execute_reply.started":"2025-11-05T11:26:33.455255Z","shell.execute_reply":"2025-11-05T11:26:33.480620Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T11:26:37.609863Z","iopub.execute_input":"2025-11-05T11:26:37.610248Z","iopub.status.idle":"2025-11-05T11:26:37.654281Z","shell.execute_reply.started":"2025-11-05T11:26:37.610215Z","shell.execute_reply":"2025-11-05T11:26:37.653069Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"trace['beta_cs'].mean()","metadata":{"execution":{"iopub.status.busy":"2025-11-05T11:26:41.856170Z","iopub.execute_input":"2025-11-05T11:26:41.856586Z","iopub.status.idle":"2025-11-05T11:26:41.865182Z","shell.execute_reply.started":"2025-11-05T11:26:41.856531Z","shell.execute_reply":"2025-11-05T11:26:41.864024Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T11:26:45.448623Z","iopub.execute_input":"2025-11-05T11:26:45.449005Z","iopub.status.idle":"2025-11-05T11:26:47.028815Z","shell.execute_reply.started":"2025-11-05T11:26:45.448971Z","shell.execute_reply":"2025-11-05T11:26:47.027177Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"with new_model:\n    trace2 = pm.sample(2000, tune=2000, target_accept=.9)","metadata":{"execution":{"iopub.status.busy":"2025-11-05T10:41:37.831796Z","iopub.execute_input":"2025-11-05T10:41:37.832269Z","iopub.status.idle":"2025-11-05T10:42:44.760922Z","shell.execute_reply.started":"2025-11-05T10:41:37.832220Z","shell.execute_reply":"2025-11-05T10:42:44.759759Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"trace2['FVC_est']","metadata":{"execution":{"iopub.status.busy":"2025-11-05T10:42:44.762957Z","iopub.execute_input":"2025-11-05T10:42:44.763438Z","iopub.status.idle":"2025-11-05T10:42:44.776715Z","shell.execute_reply.started":"2025-11-05T10:42:44.763387Z","shell.execute_reply":"2025-11-05T10:42:44.775693Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"trace2","metadata":{"execution":{"iopub.status.busy":"2025-11-05T10:42:44.778385Z","iopub.execute_input":"2025-11-05T10:42:44.778703Z","iopub.status.idle":"2025-11-05T10:42:44.785479Z","shell.execute_reply.started":"2025-11-05T10:42:44.778671Z","shell.execute_reply":"2025-11-05T10:42:44.784401Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"**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":"2025-11-05T10:42:44.787284Z","iopub.execute_input":"2025-11-05T10:42:44.787663Z","iopub.status.idle":"2025-11-05T10:42:45.376429Z","shell.execute_reply.started":"2025-11-05T10:42:44.787629Z","shell.execute_reply":"2025-11-05T10:42:45.375431Z"},"trusted":true},"outputs":[],"execution_count":null},{"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":"2025-11-05T11:32:50.514241Z","iopub.execute_input":"2025-11-05T11:32:50.514855Z","iopub.status.idle":"2025-11-05T11:32:50.564627Z","shell.execute_reply.started":"2025-11-05T11:32:50.514800Z","shell.execute_reply":"2025-11-05T11:32:50.563444Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(test.columns.tolist())\nprint (train.columns.tolist())\nprint(train.head())\nprint (test.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T11:19:24.959400Z","iopub.execute_input":"2025-11-05T11:19:24.959839Z","iopub.status.idle":"2025-11-05T11:19:24.978720Z","shell.execute_reply.started":"2025-11-05T11:19:24.959801Z","shell.execute_reply":"2025-11-05T11:19:24.977356Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from sklearn import preprocessing\nfrom sklearn.metrics import r2_score\nimport numpy as np\n\n# --- 0️⃣ Encode PatientID in test (same as train) ---\nle = preprocessing.LabelEncoder()\nle.fit(train['Patient'])\ntest['PatientID'] = le.transform(test['Patient'])\n\n# --- 1️⃣ Observed FVC ---\ny_true = test['FVC_base'].values\n\n# --- 2️⃣ Prepare indices ---\npatient_indices = test['PatientID'].values\npatient_classes = test['Class'].values\nweeks = test['Weeks_base'].values  # use Weeks_base if you only have that\nages = test['Age'].values\n\n# --- 3️⃣ Posterior means from PyMC trace ---\nalpha_mean = trace['alpha'].mean(axis=0)\nbeta_mean = trace['beta'].mean(axis=0)\nbeta_int_mean = trace['beta_int'].mean()\nalpha_s_mean = trace['alpha_s'].mean()\nbeta_cs_mean = trace['beta_cs'].mean(axis=0)\n\n# --- 4️⃣ Predicted FVC ---\ny_pred = (\n    alpha_mean[patient_indices] +\n    beta_mean[patient_indices] * weeks +\n    beta_int_mean +\n    alpha_s_mean +\n    beta_cs_mean[patient_classes] * ages\n)\n\n# --- 5️⃣ Compute R² ---\nr2 = r2_score(y_true, y_pred)\nprint(\"R² score:\", r2)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T11:39:32.582516Z","iopub.execute_input":"2025-11-05T11:39:32.582943Z","iopub.status.idle":"2025-11-05T11:39:32.604706Z","shell.execute_reply.started":"2025-11-05T11:39:32.582906Z","shell.execute_reply":"2025-11-05T11:39:32.603416Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# -------------------------------\n# 0️⃣ Imports\n# -------------------------------\nimport pandas as pd\nimport numpy as np\nimport pymc3 as pm\nimport arviz as az\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import r2_score\n\n# -------------------------------\n# 1️⃣ Load data\n# -------------------------------\ntrain = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\n\n# -------------------------------\n# 2️⃣ Preprocessing\n# -------------------------------\n# Add baseline week and FVC per patient\n\n\ntrain['Weeks_base'] = train.groupby('Patient')['Weeks'].transform('min')\ntrain['FVC_base'] = train.groupby('Patient')['FVC'].transform('first')\n\n# For test, create Weeks_base and FVC_base if missing\n\n\nif 'Weeks_base' not in test.columns:\n    test['Weeks_base'] = test.groupby('Patient')['Weeks'].transform('min')\nif 'FVC_base' not in test.columns:\n    test['FVC_base'] = test.groupby('Patient')['FVC'].transform('first')\n\n# Encode PatientID\nle = LabelEncoder()\ntrain['PatientID'] = le.fit_transform(train['Patient'])\ntest['PatientID'] = le.transform(test['Patient'])\n\n# Define patient class\ndef calculate_patient_class(sex, smoking_status):\n    if sex == 'Male':\n        return {'Currently smokes':0,'Ex-smoker':1,'Never smoked':2}[smoking_status]\n    else:\n        return {'Currently smokes':3,'Ex-smoker':4,'Never smoked':5}[smoking_status]\n\ntrain['Class'] = train.apply(lambda row: calculate_patient_class(row['Sex'], row['SmokingStatus']), axis=1)\ntest['Class'] = test.apply(lambda row: calculate_patient_class(row['Sex'], row['SmokingStatus']), axis=1)\n\n# -------------------------------\n# 3️⃣ Map unique patients to indices\n# -------------------------------\nunique_patients = train[['PatientID', 'FVC_base', 'Weeks_base', 'Age', 'Class']].drop_duplicates().sort_values('PatientID')\npatient_idx_map = dict(zip(unique_patients['PatientID'], range(len(unique_patients))))\n\ntrain['pid_idx'] = train['PatientID'].map(patient_idx_map)\ntest['pid_idx'] = test['PatientID'].map(patient_idx_map)\n\n# -------------------------------\n# 4️⃣ Prepare patient-level data\n# -------------------------------\n\nFVC_base_pat = unique_patients['FVC_base'].values\nWeeks_base_pat = unique_patients['Weeks_base'].values\nAge_pat = unique_patients['Age'].values\nClass_pat = unique_patients['Class'].values\n\n# -------------------------------\n# 5️⃣ Build PyMC3 Hierarchical Model\n# -------------------------------\nwith pm.Model() as hierarchical_model:\n    # Hyperpriors for alpha\n    beta_int = pm.Normal('beta_int', 0, sigma=100)\n    sigma_alpha = pm.HalfNormal('sigma_alpha', 100)\n    alpha = pm.Normal('alpha', mu=FVC_base_pat, sigma=sigma_alpha, shape=len(FVC_base_pat))\n    \n    # Hyperpriors for beta\n    sigma_beta = pm.HalfNormal('sigma_beta', 100)\n    alpha_s = pm.Normal('alpha_s', 0, sigma=100)\n    beta_cs = pm.Normal('beta_cs', 0, sigma=100, shape=6)\n    gamma_cs = pm.Normal('gamma_cs', 0, sigma=10, shape=6)  # optional quadratic effect\n    \n    # Beta for each patient\n    beta_mu = alpha_s + Age_pat * beta_cs[Class_pat]  # linear age effect\n    beta = pm.Normal('beta', mu=beta_mu, sigma=sigma_beta, shape=len(FVC_base_pat))\n    \n    # Model variance\n    sigma = pm.HalfNormal('sigma', 200)\n    \n    # FVC estimate per row in train\n    FVC_est = alpha[train['pid_idx'].values] + beta[train['pid_idx'].values] * train['Weeks'].values\n    FVC_like = pm.Normal('FVC_like', mu=FVC_est, sigma=sigma, observed=train['FVC'].values)\n    \n    # Sample from posterior\n    trace = pm.sample(500, tune=500, target_accept=0.95, cores=2)\n\n# -------------------------------\n# 6️⃣ Predict on test set\n# -------------------------------\nalpha_mean = trace['alpha'].mean(axis=0)\nbeta_mean = trace['beta'].mean(axis=0)\n\ny_pred = alpha_mean[test['pid_idx'].values] + beta_mean[test['pid_idx'].values] * test['Weeks_base'].values\ny_true = test['FVC_base'].values\n\n# -------------------------------\n# 7️⃣ Compute R²\n# -------------------------------\nr2 = r2_score(y_true, y_pred)\nprint(\"✅ R² score:\", r2)\n\n# -------------------------------\n# 8️⃣ Save model trace\n# -------------------------------\naz.to_netcdf(trace, 'fvc_model_trace.nc')\nprint(\"✅ Model saved as fvc_model_trace.nc\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T11:51:09.952257Z","iopub.execute_input":"2025-11-05T11:51:09.952711Z","iopub.status.idle":"2025-11-05T11:52:18.597180Z","shell.execute_reply.started":"2025-11-05T11:51:09.952675Z","shell.execute_reply":"2025-11-05T11:52:18.595790Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# -------------------------------\n# 9️⃣ Compute Correlation Metrics\n# -------------------------------\nfrom scipy.stats import pearsonr, spearmanr\n\n# Pearson correlation\npearson_corr, pearson_p = pearsonr(y_true, y_pred)\nprint(f\"✅ Pearson correlation: {pearson_corr:.4f} (p-value: {pearson_p:.4e})\")\n\n# Spearman correlation\nspearman_corr, spearman_p = spearmanr(y_true, y_pred)\nprint(f\"✅ Spearman correlation: {spearman_corr:.4f} (p-value: {spearman_p:.4e})\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T13:00:09.206398Z","iopub.execute_input":"2025-11-05T13:00:09.207078Z","iopub.status.idle":"2025-11-05T13:00:09.236247Z","shell.execute_reply.started":"2025-11-05T13:00:09.207024Z","shell.execute_reply":"2025-11-05T13:00:09.234454Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# -------------------------------\n# 9️⃣ Correlation Metrics Table\n# -------------------------------\nfrom scipy.stats import pearsonr, spearmanr\nimport pandas as pd\n\n# Calculate correlations\npearson_corr, pearson_p = pearsonr(y_true, y_pred)\nspearman_corr, spearman_p = spearmanr(y_true, y_pred)\n\n# Create a results table\nmetrics_table = pd.DataFrame({\n    'Metric': ['R²', 'Pearson', 'Spearman'],\n    'Value': [r2, pearson_corr, spearman_corr],\n    'p-value': [None, pearson_p, spearman_p]\n})\n\nprint(\"📊 Correlation Metrics:\")\nprint(metrics_table)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T13:00:48.369409Z","iopub.execute_input":"2025-11-05T13:00:48.369946Z","iopub.status.idle":"2025-11-05T13:00:48.398177Z","shell.execute_reply.started":"2025-11-05T13:00:48.369907Z","shell.execute_reply":"2025-11-05T13:00:48.396493Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# -------------------------------\n# Correlation Matrix Heatmap\n# -------------------------------\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\n# Select numeric features\nfeatures = ['Weeks', 'FVC', 'Percent', 'Age']\n\n# Compute correlation matrix\ncorr_matrix = train[features].corr()\n\n# Plot heatmap\nplt.figure(figsize=(8,6))\nsns.heatmap(corr_matrix, annot=True, cmap='coolwarm', center=0, linewidths=1, linecolor='white')\nplt.title(\"Correlation Matrix of Features\", fontsize=14)\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T13:05:13.648152Z","iopub.execute_input":"2025-11-05T13:05:13.648626Z","iopub.status.idle":"2025-11-05T13:05:13.995968Z","shell.execute_reply.started":"2025-11-05T13:05:13.648588Z","shell.execute_reply":"2025-11-05T13:05:13.994960Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# -------------------------------\n# 0️⃣ Imports\n# -------------------------------\nimport pandas as pd\nimport numpy as np\nimport pymc3 as pm\nimport arviz as az\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import r2_score\n\n# -------------------------------\n# 1️⃣ Load data\n# -------------------------------\ntrain = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv')\ntest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\n\n# -------------------------------\n# 2️⃣ Preprocessing\n# -------------------------------\n# Add baseline week and FVC per patient\ntrain['Weeks_base'] = train.groupby('Patient')['Weeks'].transform('min')\ntrain['FVC_base'] = train.groupby('Patient')['FVC'].transform('first')\n\n# For test, create Weeks_base and FVC_base if missing\nif 'Weeks_base' not in test.columns:\n    test['Weeks_base'] = test.groupby('Patient')['Weeks'].transform('min')\nif 'FVC_base' not in test.columns:\n    test['FVC_base'] = test.groupby('Patient')['FVC'].transform('first')\n\n# Encode PatientID\nle = LabelEncoder()\ntrain['PatientID'] = le.fit_transform(train['Patient'])\ntest['PatientID'] = le.transform(test['Patient'])\n\n# Define patient class\ndef calculate_patient_class(sex, smoking_status):\n    if sex == 'Male':\n        return {'Currently smokes':0,'Ex-smoker':1,'Never smoked':2}[smoking_status]\n    else:\n        return {'Currently smokes':3,'Ex-smoker':4,'Never smoked':5}[smoking_status]\n\ntrain['Class'] = train.apply(lambda row: calculate_patient_class(row['Sex'], row['SmokingStatus']), axis=1)\ntest['Class'] = test.apply(lambda row: calculate_patient_class(row['Sex'], row['SmokingStatus']), axis=1)\n\n# -------------------------------\n# 3️⃣ Map unique patients to indices\n# -------------------------------\nunique_patients = train[['PatientID', 'FVC_base', 'Weeks_base', 'Age', 'Class']].drop_duplicates().sort_values('PatientID')\npatient_idx_map = dict(zip(unique_patients['PatientID'], range(len(unique_patients))))\n\ntrain['pid_idx'] = train['PatientID'].map(patient_idx_map)\ntest['pid_idx'] = test['PatientID'].map(patient_idx_map)\n\n# -------------------------------\n# 4️⃣ Prepare patient-level data\n# -------------------------------\nFVC_base_pat = unique_patients['FVC_base'].values\nWeeks_base_pat = unique_patients['Weeks_base'].values\nAge_pat = unique_patients['Age'].values\nClass_pat = unique_patients['Class'].values\n\n# -------------------------------\n# 5️⃣ PyMC3 Hierarchical Model with nonlinear effects\n# -------------------------------\nwith pm.Model() as hierarchical_model:\n    # Hyperpriors for alpha (patient-specific intercept)\n    sigma_alpha = pm.HalfNormal('sigma_alpha', 100)\n    alpha = pm.Normal('alpha', mu=FVC_base_pat, sigma=sigma_alpha, shape=len(FVC_base_pat))\n    \n    # Hyperpriors for beta (patient-specific slope)\n    sigma_beta = pm.HalfNormal('sigma_beta', 100)\n    \n    # Age × Class effect\n    beta_cs = pm.Normal('beta_cs', 0, sigma=100, shape=6)\n    \n    # Weeks × Class linear effect\n    beta_week_cs = pm.Normal('beta_week_cs', 0, sigma=10, shape=6)\n    \n    # Weeks² × Class nonlinear effect\n    gamma_week_cs = pm.Normal('gamma_week_cs', 0, sigma=5, shape=6)\n    \n    # Compute beta for each patient\n    beta = (Age_pat * beta_cs[Class_pat] +\n            Weeks_base_pat * beta_week_cs[Class_pat] +\n            (Weeks_base_pat**2) * gamma_week_cs[Class_pat])\n    \n    # Model variance\n    sigma = pm.HalfNormal('sigma', 200)\n    \n    # FVC estimate per training row\n    FVC_est = (alpha[train['pid_idx'].values] +\n               beta[train['pid_idx'].values] * train['Weeks'].values +\n               gamma_week_cs[Class_pat[train['pid_idx'].values]] * (train['Weeks'].values**2))\n    \n    # Likelihood\n    FVC_like = pm.Normal('FVC_like', mu=FVC_est, sigma=sigma, observed=train['FVC'].values)\n    \n    # Sample from posterior\n    trace = pm.sample(500, tune=500, target_accept=0.95, cores=2)\n\n# -------------------------------\n# 6️⃣ Predict on test set\n# -------------------------------\nalpha_mean = trace['alpha'].mean(axis=0)\nbeta_cs_mean = trace['beta_cs'].mean(axis=0)\nbeta_week_cs_mean = trace['beta_week_cs'].mean(axis=0)\ngamma_week_cs_mean = trace['gamma_week_cs'].mean(axis=0)\n\ny_pred = (\n    alpha_mean[test['pid_idx'].values] +\n    (test['Age'].values * beta_cs_mean[test['Class'].values]) +\n    (test['Weeks_base'].values * beta_week_cs_mean[test['Class'].values]) +\n    ((test['Weeks_base'].values**2) * gamma_week_cs_mean[test['Class'].values])\n)\ny_true = test['FVC_base'].values\n\n# -------------------------------\n# 7️⃣ Compute R²\n# -------------------------------\nr2 = r2_score(y_true, y_pred)\nprint(\"✅ R² score:\", r2)\n\n# -------------------------------\n# 8️⃣ Save model trace\n# -------------------------------\naz.to_netcdf(trace, 'fvc_model_trace1.nc')\nprint(\"✅ Model saved as fvc_model_trace1.nc\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-05T12:04:13.341796Z","iopub.execute_input":"2025-11-05T12:04:13.342329Z","iopub.status.idle":"2025-11-05T12:09:06.252601Z","shell.execute_reply.started":"2025-11-05T12:04:13.342190Z","shell.execute_reply":"2025-11-05T12:09:06.251251Z"}},"outputs":[],"execution_count":null}]}