{
"cells": [
{
"cell_type": "markdown",
"id": "a2a5c918-abc4-437f-99e4-7bcdd98a0a45",
"metadata": {},
"source": [
"#### Survival inference with pymc"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "dd5eb764-240a-4529-a5c0-f3a2d652fa86",
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"import numpy as np\n",
"import pymc as pm\n",
"import arviz as az\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from sklearn.preprocessing import StandardScaler\n",
"\n",
"from survivalpredict.estimators import CoxProportionalHazard\n",
"from survivalpredict.datasets import load_iranian_telecom_churn\n",
"\n",
"az.style.use(\"arviz-darkgrid\")"
]
},
{
"cell_type": "markdown",
"id": "57ef9356-49a4-41ce-93b8-44cde10999d3",
"metadata": {},
"source": [
"In cases of inference or prediction on a trained models, it might be valuable to gauge uncertainty. PyMC is a Bayesian/‘Markov chain Monte Carlo’ framework. Simply put, PyMC allows us to build models that can capture the uncertainty and possible variability of our model, given the observed data and model assumptions. There is much literature on PyMC/MCMC; we will assume some familiarity with the topic.\n",
"\n",
"Curretly surivivalpredict supports generating equivalent PyMC models for two of its estimator classes. The `CoxProportionalHazard`s `get_pymc_model` method will return an equivalent pymc model. Models with either ‘efron’ or ‘breslow’ ties, stratification, L2 regularization, and left-censorship are all supported. The `ParametricDiscreteTimePH` actually uses PyMC under the hood, calling the `get_pymc_model` will return the underlying pymc model."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "36b01f7b-e117-4546-bdf3-77e38ce7aa67",
"metadata": {},
"outputs": [],
"source": [
"#loading some sample data\n",
"iranian_telecom_churn_dict = load_iranian_telecom_churn()\n",
"X = iranian_telecom_churn_dict['X']\n",
"#in this case is important to scale our data\n",
"X = StandardScaler().fit_transform(X)\n",
"times = iranian_telecom_churn_dict['times'].astype(np.int64)\n",
"events = iranian_telecom_churn_dict['events'].astype(np.bool_)\n",
"column_names = iranian_telecom_churn_dict['column_names']"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "e0a77edb-b6f0-4407-b39f-98bc3d7ca112",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 0.50777218 0.48238747 -0.44392135 0.42167482 -2.31404596 -2.99677386\n",
" -0.37073225 -0.417258 0.26621988 -0.03414636 0.30532322 1.86236737]\n"
]
}
],
"source": [
"cox = CoxProportionalHazard(ties='efron').fit(X,times,events)\n",
"print(cox.coef_)"
]
},
{
"cell_type": "markdown",
"id": "8141e36c-734e-46b6-a2f6-7caddb9a776c",
"metadata": {},
"source": [
"The `CoxProportionalHazard` class comes with a pymc model generator that builds pymc models that are equivalent to the survivalpredict model. It should be noted that adding time_start/'left-censorship' will drastically reduce the samples per second of the pymc class, as it triggers more computation. Stratification will also somewhat slow down sampling.\n",
"\n",
"\n",
"When calling `get_pymc_model`, if the `empirical_bayes` argument is set to True, and our estimator class has already been trained, the initial values of our PyMC model coefficients will be the coefficients trained by our estimator class."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "55ab5e09-f348-4392-904a-e4f10cd35c75",
"metadata": {},
"outputs": [],
"source": [
"#Getting our pymc model, setting empirical_bayes to False as part of our demonstration.\n",
"pymc_model = cox.get_pymc_model(X, times, events,labes_names=column_names,empirical_bayes=False)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "ab54336f-055e-44ab-ade4-4f9f3e8dad05",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pm.model_to_graphviz(pymc_model)"
]
},
{
"cell_type": "markdown",
"id": "3a1a4f44-e6ef-4daf-bc58-57bf67fc4869",
"metadata": {},
"source": [
"A quick demonstration to show that this PyMC model is equivalent to our survivalpredict model would be to generate the ‘maximum likelihood estimation’/’maximum a posteriori’ estimates and compare them to our coefficients. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "f6c06259-5267-49c6-ad0e-904407fad073",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "dc42e01b26d14eae9090c264ea0b26d0",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Output()"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n"
],
"text/plain": []
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"with pymc_model:\n",
" mle = pm.find_MAP()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "da814133-ab0b-4331-b490-01ed8c2eb657",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'coefs': array([ 0.50733508, 0.48238074, -0.44546389, 0.43406345, -2.31074651,\n",
" -2.9503162 , -0.37137833, -0.4172489 , 0.26588864, -0.03518892,\n",
" 0.30293834, 1.81416004]),\n",
" 'n_log_likelihood': array(3052.59817273)}"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mle"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "93d9cc09-6392-45a2-a4b6-6cb2aa628ea9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"coefs: [ 0.50777218 0.48238747 -0.44392135 0.42167482 -2.31404596 -2.99677386\n",
" -0.37073225 -0.417258 0.26621988 -0.03414636 0.30532322 1.86236737]\n",
"n_log_likelihood: 3052.597200555113\n"
]
}
],
"source": [
"print('coefs: ',cox.coef_)\n",
"print('n_log_likelihood: ',cox.n_log_likelihood)\n",
"#looks pretty close to me"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "7eef1673-1481-46fc-9498-016a6c1bf713",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"NUTS[nutpie]: [coefs]\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "6a2367362e08477e9da8286d61a75989",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Output()"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n"
],
"text/plain": []
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#generating the trace/sampling\n",
"#this might take a bit, it can be faster with a gpu.\n",
"with pymc_model:\n",
" idata = pm.sample()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "414efe75-a9a2-476d-bfa0-ca1857e976e6",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" mean | \n",
" sd | \n",
" eti89_lb | \n",
" eti89_ub | \n",
" ess_bulk | \n",
" ess_tail | \n",
" r_hat | \n",
" mcse_mean | \n",
" mcse_sd | \n",
"
\n",
" \n",
" \n",
" \n",
" | coefs[call_failure] | \n",
" 0.51 | \n",
" 0.078 | \n",
" 0.38 | \n",
" 0.64 | \n",
" 2108 | \n",
" 2427 | \n",
" 1.00 | \n",
" 0.0017 | \n",
" 0.0012 | \n",
"
\n",
" \n",
" | coefs[complains] | \n",
" 0.482 | \n",
" 0.0289 | \n",
" 0.43 | \n",
" 0.53 | \n",
" 5043 | \n",
" 3338 | \n",
" 1.00 | \n",
" 0.00041 | \n",
" 0.00029 | \n",
"
\n",
" \n",
" | coefs[charge_amount] | \n",
" -0.456 | \n",
" 0.142 | \n",
" -0.69 | \n",
" -0.23 | \n",
" 2603 | \n",
" 2461 | \n",
" 1.00 | \n",
" 0.0028 | \n",
" 0.002 | \n",
"
\n",
" \n",
" | coefs[seconds_of_use] | \n",
" 0.41 | \n",
" 0.48 | \n",
" -0.35 | \n",
" 1.2 | \n",
" 888 | \n",
" 1283 | \n",
" 1.01 | \n",
" 0.016 | \n",
" 0.012 | \n",
"
\n",
" \n",
" | coefs[frequency_of_use] | \n",
" -2.32 | \n",
" 0.34 | \n",
" -2.9 | \n",
" -1.8 | \n",
" 1650 | \n",
" 2063 | \n",
" 1.00 | \n",
" 0.0084 | \n",
" 0.0062 | \n",
"
\n",
" \n",
" | coefs[frequency_of_sms] | \n",
" -3.04 | \n",
" 1.07 | \n",
" -4.8 | \n",
" -1.4 | \n",
" 828 | \n",
" 1112 | \n",
" 1.01 | \n",
" 0.037 | \n",
" 0.025 | \n",
"
\n",
" \n",
" | coefs[distinct_called_numbers] | \n",
" -0.366 | \n",
" 0.14 | \n",
" -0.59 | \n",
" -0.14 | \n",
" 3186 | \n",
" 2430 | \n",
" 1.00 | \n",
" 0.0025 | \n",
" 0.0017 | \n",
"
\n",
" \n",
" | coefs[age_group] | \n",
" -0.41 | \n",
" 0.18 | \n",
" -0.69 | \n",
" -0.12 | \n",
" 1411 | \n",
" 1802 | \n",
" 1.00 | \n",
" 0.0048 | \n",
" 0.0034 | \n",
"
\n",
" \n",
" | coefs[tariff_plan] | \n",
" 0.24 | \n",
" 0.126 | \n",
" 0.031 | \n",
" 0.43 | \n",
" 5427 | \n",
" 3313 | \n",
" 1.00 | \n",
" 0.0017 | \n",
" 0.0013 | \n",
"
\n",
" \n",
" | coefs[status] | \n",
" -0.034 | \n",
" 0.068 | \n",
" -0.14 | \n",
" 0.074 | \n",
" 2733 | \n",
" 2441 | \n",
" 1.00 | \n",
" 0.0013 | \n",
" 0.00091 | \n",
"
\n",
" \n",
" | coefs[age] | \n",
" 0.3 | \n",
" 0.197 | \n",
" -0.011 | \n",
" 0.61 | \n",
" 1221 | \n",
" 1672 | \n",
" 1.00 | \n",
" 0.0056 | \n",
" 0.004 | \n",
"
\n",
" \n",
" | coefs[customer_value] | \n",
" 1.86 | \n",
" 1.14 | \n",
" 0.1 | \n",
" 3.7 | \n",
" 814 | \n",
" 1001 | \n",
" 1.01 | \n",
" 0.04 | \n",
" 0.027 | \n",
"
\n",
" \n",
" | n_log_likelihood | \n",
" 3058.6 | \n",
" 2.5 | \n",
" 3100 | \n",
" 3100 | \n",
" 1529 | \n",
" 1939 | \n",
" 1.00 | \n",
" 0.065 | \n",
" 0.062 | \n",
"
\n",
" \n",
"
\n",
"