{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Why might we want to go Bayesian here?\n\nParticularly with very little data, when you need uncertainty quantification and when you need to optimize a metric that is not your models loss function/likelihood, a Bayesian + decision theoretic approach is very attractive.\n\nThis notbook follows a rather classical Bayesian approach:\n* I define a sensible model in terms of a likelihood function (in this case a linear decline model with both intercept and slope influenced by patient-specific random effects and some covariates/features).\n* I chose to model log-transformed Percent (=FVC percent predicted of normal).\n    * Using Percent is logical, because it already accounts for some covariates like sex, height, and age, rather than modeling FVC directly (you can always get a predicted FVC from the predicted Percent from the ratio of FVC to Percent at baseline). \n    * A log-transformation seemed like a good idea, because residuals should be more normal that way and heteroskedasticity is less of a concern.\n* One nice idea as also pointed out by others before is to have  patient-specific random effects on the intercept and the slope. \n    * This has the nice property to sensibly reflect the correlation structure in the data - you just do not learn as much from multiple observations from the same patient as from multiple independent observations from separate patients. The posterior uncertainty from the Bayesian model nicely reflects this.\n    * Another rationale for why I thought this was a very good idea is that the baseline measurement contains measurement error and we account for this more clearly by having both the baseline observation as a coveriate and a patient specific random effect. I tried experiments where I made the random effects of the intercept and the slope correlated, but I did not seem to gain much from that.\n    * For those more used to neural network, you can think of these random effects as 2D-embeddings for a patient that have very specific meanings (i.e. first one is about where the patient lies on average relative to the population of patients and the second one describes whether the patient has an above or below average slope over time).\n    * You might think, can we not do what you did here with a neural network? That certainly occurred to me, but at least my attempts produced worse results (more my problem than that of the neural networks).\n* I specified some reasonable prior distributions based on what we know a-priori (e.g. we have an approximate idea of the expected decline over time - I work in respiratory drug development, so I had at least some idea of what the priors should be). I experimented with tweaking the prior parameters using cross-validation, but I failed to beat the priors I specified based on my knowledge.\n* Then, we let the MCMC sampler do its magice and give us MCMC samples from the posterior implied by what we defined above.\n* Then we take your posterior belief about the truth as represented by the MCMC samples from the posterior distribution and perform a decision analysis to optimize the metric we care about - i.e. the competition metric. \n* Also note that there are limits beyond which you never want to go with the competition metric due to the way it gets clipped.\n\nI had worked on this notebook even before seeing the [nice notebook](https://www.kaggle.com/carlossouza/bayesian-experiments) by [Carlos Souza](https://www.kaggle.com/carlossouza), who used [PyMC3](https://docs.pymc.io/) and has some very nice diagrams of the model structure (definitely worth visiting and looking at). The main differences to that notebook is that I model log-transformed Percent (rather than FVC directly), and try to let additional covariates/features affect the intercept and slope for each patient. The other - but less important difference - is that I am using [pystan](https://pystan.readthedocs.io/en/latest/) (and you could do the same thing in R with rstan - see the [Stan homepage](https://mc-stan.org/) for more options and information on the Stan project), but the two packages both use the NUTS-MCMC sampler. I used [pystan](https://pystan.readthedocs.io/en/latest/) simply because I know Stan pretty well and felt more like I knew exactly what I was doing.\n\nThis went into our team's final submission as part of a bunch of models that we stacked, but I thought I'd share it, since there's a couple of nice ideas here. I left out the cross-validation (we used a 5-fold that never split a patient across folds and stacked on the last 3 out-of-fold observations per patient) here and just provide code for doing a submission, in part because this code takes pretty long to run. In my day-to-day work, I would try to speed this up with within chain-parallization, but within a Kaggle notebook we do not have the necessary CPUs."},{"metadata":{},"cell_type":"markdown","source":"# Read the training and test data"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nimport os.path\nimport pystan\nfrom scipy.optimize import minimize\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 len(test)==5:\n    test['Patient'] = test['Patient'] + \"_t\" # This ensures in development situation test patients are treated differently than in same record for train patients\n\ntest2 = test[['Patient', 'Age', 'Sex', 'SmokingStatus']].merge(pd.DataFrame([(Patient, np.int(Weeks)) for Weeks in range(-12,134) for Patient in set(test['Patient'])],\n                        columns=['Patient', 'Weeks']), \n           on=\"Patient\", how=\"left\")\ntest2['fold'] = 0\n\ntrain = train.append(test).reset_index(drop=True)\ntrain['fold'] = 1\n\ntrain = train.append(test2).reset_index(drop=True)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Obtain baseline values & number patient\nI.e. first FVC and Percent values that we have even on the test set."},{"metadata":{"trusted":true},"cell_type":"code","source":"base = (train.loc[train[train['fold']==1].groupby('Patient')['Weeks'].idxmin()]\n    .groupby('Patient')\n    .agg(\n        base = (\"FVC\", lambda x: np.mean(x)),\n        logbase = (\"FVC\", lambda x: np.mean(np.log(x))),\n        basepercent = (\"Percent\", lambda x: np.mean(x)),        \n        logbasepercent = (\"Percent\", lambda x: np.mean(np.log(x))),\n        basewk = (\"Weeks\", lambda x: np.mean(x))\n    )\n    .filter(items=['Patient', 'base', 'logbase', 'basewk', 'basepercent', 'logbasepercent'])\n    .reset_index())\n\ntrain2 = (train\n  .merge(base, on='Patient', how='left')\n  .assign(sexm = lambda x: np.select([x.Sex.eq('Male'), x.Sex.eq('Female')], [1, 0]),\n          relweek = lambda x: x.Weeks-x.basewk,\n          logpercent = lambda x: np.log(x.Percent),\n          smoker = lambda x: np.select( [x.SmokingStatus.eq(\"Never smoked\"),\n                                        x.SmokingStatus.eq(\"Ex-smoker\"),\n                                        x.SmokingStatus.eq(\"Currently smokes\")],\n                                      [0, 1, 2])))\n\ntrain2.Patient = pd.Categorical(train2.Patient)\ntrain2['patno'] = train2.Patient.cat.codes\ntrain2['denom'] = 100*train2['base']/train2['basepercent']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Generate training (and test/validation data in the format pystan wants)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Function to generate training and validation data in the format pystan wants it\ndef get_train_data(fold, hyperparams):\n    train2tmp = (train2.copy())[ (train2['fold']!=fold) ]\n    train_dat = {'records': len(train2tmp),\n            'patients': len(np.unique(train2['patno'])), # We want the actual number of patients so that we immediately create random effects even for out-of-fold patients\n            'patno': train2tmp['patno']+1,            \n            'denom': train2tmp['denom'],\n            'fvc': train2tmp['FVC'],\n            'base': train2tmp['base'],\n            'basenormal': 1.0*(train2tmp['basepercent']>=80),\n            'basepercent': train2tmp['basepercent']*(train2tmp['basepercent']<80.0) + (train2tmp['basepercent']>=80)*(80.0 + (train2tmp['basepercent']-80.0)*2.0/3.0),\n            'logpercent': train2tmp['logpercent'],\n            'logbasepercent': train2tmp['logbasepercent'],\n            'relweek': train2tmp['relweek'],\n            'Weeks': train2tmp['Weeks'],\n            'smoker': train2tmp['smoker'],             \n            'basewk': train2tmp['basewk'],\n            'hyperparams': hyperparams}\n    train2tmp2 = (train2.copy())[ (train2['fold']==fold) ]\n    test_dat = {'records': len(train2tmp2),\n                'Patient': train2tmp2['Patient'],\n                'Weeks': train2tmp2['Weeks'],\n            'patients': len(np.unique(train2['patno'])),\n            'patno': train2tmp2['patno']+1,            \n            'denom': train2tmp2['denom'],\n            'fvc': train2tmp2['FVC'],\n            'base': train2tmp2['base'],\n            'basenormal': 1.0*(train2tmp2['basepercent']>=80),\n            'basepercent': train2tmp2['basepercent']*(train2tmp2['basepercent']<80.0) + (train2tmp2['basepercent']>=80)*(80.0 + (train2tmp2['basepercent']-80.0)*2.0/3.0),\n            'logpercent': train2tmp2['logpercent'],\n            'logbasepercent': train2tmp2['logbasepercent'],\n            'relweek': train2tmp2['relweek'],            \n            'Weeks': train2tmp2['Weeks'],            \n            'smoker': train2tmp2['smoker'],             \n            'basewk': train2tmp2['basewk'],\n            'hyperparams': hyperparams}\n    return train_dat, test_dat","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# How do we get point predictions and \"confidence\" from posterior MCMC samples?\nThis is a classic Bayesian approach, you take your posterior belief about the truth as represented by the MCMC samples from the posterior distribution and then you perform decision analysis to optimize the metric you care about - i.e. the competition metric. Also note that there are limits beyond which you never want to go with the competition metric due to the way it gets clipped."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Function to obtain confidence based on posterior predictive distribution samples\n\ndef fct(x, samples):\n    confidence = np.maximum(70.0, x[1])\n    return - np.mean( - np.sqrt(2.0)*np.minimum( np.absolute( x[0] - samples ), 1000.0)/x[1] - np.log(np.sqrt(2.0)*x[1]) )\n#def minimize_metric(samples):\n#samples = np.random.normal(loc=0.0, scale=1.0, size=10000)\n\n# Even when we are wrong by 1000 or more, a confidence of >1414 is just losing on the metric\n# when we are at least wrong by 70, then <99 is pointless, but if it's 35, then that's 50, so let's set the lower limit to 70\ndef minimize_metric(samples):\n    res = minimize(fct, x0=np.array( [np.maximum(101,np.minimum(8999,np.mean(samples))), \n                                      np.minimum(1413, np.maximum(71,np.std(samples)))]), \n                   args=samples, bounds = ((100, 9000), (70, 1414))).x\n    return res[0], res[1]\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Specifying the Stan  model\n\nThis specifies the Stan model. For more information on the Stan modeling language see the [Stan homepage ](https://mc-stan.org/)and the documentation you can finde there. It's a basic linear random effects model."},{"metadata":{"trusted":true},"cell_type":"code","source":"stan_code = \"\"\"\ndata {\n    int<lower=0> records; // Total number of records\n    int<lower=0> patients; // Total number of patients\n    int patno[records]; // Patient number for each record\n    \n    real logpercent[records]; // Log-Percent for each record\n    \n    int relweek[records]; // Week relative to first FVC for each record\n    real logbasepercent[records]; // baseline (first record) for the log-Percent that we are modeling\n    real basepercent[records]; // baseline Percent\n    real basenormal[records]; // Is the baseline observation normal (Percent>= 80%)\n    real smoker[records]; // Is the patient a smoker\n    real basewk[records]; // What week was the baseline observation taken in?\n    \n    real hyperparams[24]; // Hyperparameters that specify the priors - we hand these over this way rather than hard-coding them so that we do not need to re-compile the model each time we want to tweak them.\n}\n\nparameters {\n    real mu[2]; // Prior means for population intercept and slope\n    real beta[7]; // coefficients for covariates\n    \n    real<lower=0> sigma; // residual error SD\n    \n    real patre[patients]; // Patient random effects on intercept\n    real<lower=0> patre_sigma; // SD for those\n    real patslopere[patients]; // Patient random effects on slope\n    real<lower=0> patslopere_sigma; // SD for those\n}\n\nmodel {\n    // Specifying priors\n    mu[1] ~ normal(hyperparams[1], hyperparams[2]);\n    mu[2] ~ normal(hyperparams[3], hyperparams[4]);\n    \n    sigma ~ normal(hyperparams[5], hyperparams[6]);\n        \n    beta[1] ~ normal( hyperparams[7], hyperparams[8]);    \n    beta[2] ~ normal( hyperparams[9], hyperparams[10]);        \n    beta[3] ~ normal( hyperparams[11], hyperparams[12]);    \n    beta[4] ~ normal( hyperparams[13], hyperparams[14]);\n    beta[5] ~ normal( hyperparams[15], hyperparams[16]);\n    beta[6] ~ normal( hyperparams[17], hyperparams[18]);\n    beta[7] ~ normal( hyperparams[19], hyperparams[20]);\n            \n    // It turned out that this particular way of parameterizing sampled better\n    // than specifying patre_sigma ~ normal(0,1); and then multiplying by SD and adding mean later,\n    // which can often be a crucial reparameterizing for getting better MCMC sampling (but not here).\n    patre_sigma ~ normal(hyperparams[21], hyperparams[22]);\n    patslopere_sigma ~ normal(hyperparams[23], hyperparams[24]);\n    \n    patre ~ std_normal();\n    patslopere ~ std_normal();\n\n  // Specifying the likelihood\n  { \n  real tmpmeans[records]; \n  for (r in 1:records){\n      tmpmeans[r] = mu[1] + patre[patno[r]]*patre_sigma + \n                    beta[1]*logbasepercent[r]  + beta[2] * basewk[r]/52.0 + beta[3]*basepercent[r] +\n                    relweek[r]/52.0 * \n                        ( mu[2] + patslopere[patno[r]]*patslopere_sigma + \n                          beta[4]*(smoker[r]>=1) + beta[5]*(smoker[r]==2) + beta[6]*basewk[r]/52.0 + beta[7]*basenormal[r]);\n    }\n  logpercent ~ normal(tmpmeans, sigma);\n  }\n  \n}\n\"\"\"\n\nsm = pystan.StanModel(model_code=stan_code) # Compile Stan model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Fitting the Stan model\n\nYou can still see the left-overs from our cross-validation here, but we just fit to one fold (=training data) and predict the other (=test data)."},{"metadata":{"trusted":true},"cell_type":"code","source":"iter_no = 2500\nhalf_iter_no = 1250\n#iter_no = 200 # Low number for testing, if you just want to iterate\n#half_iter_no = 100\nnum_of_chains = 4\n\ntotal_number_of_patients = len(set(train2['Patient'])) # Should be number in training and validation/test sets (e.g. for developement this is 176, \n    # for submitting higher, but unknown, the baseline observation of all gets always included in training to learn random effect of patient - you can think of it as a super-simple embedding)\n\nparam = {'hyperparam1': 2.65, 'hyperparam2': 0.125, 'hyperparam3': -0.12, 'hyperparam4': 0.05, 'hyperparam5': 0.06, 'hyperparam6': 0.05, 'hyperparam7': 0.25, \n         'hyperparam8': 0.5, 'hyperparam9': 0, 'hyperparam10': 0.2, 'hyperparam11': 0.03, 'hyperparam12': 0.1, 'hyperparam13': 0, 'hyperparam14': 0.1, 'hyperparam15': 0, \n         'hyperparam16': 0.1, 'hyperparam17': 0, 'hyperparam18': 0.05, 'hyperparam19': 0, 'hyperparam20': 0.05, 'hyperparam21': 0.1, 'hyperparam22': 0.25, 'hyperparam23': 0.1, 'hyperparam24': 0.25}\n\nfor fold in range(1):\n    \n    # Get the data for the fold\n    \n    print(\"Fold: \" + str(fold))    \n    train_dat, test_dat = get_train_data(fold=fold, \n                                             hyperparams=[param['hyperparam1'],\n                                                         param['hyperparam2'],\n                                                         param['hyperparam3'],\n                                                         param['hyperparam4'],\n                                                         param['hyperparam5'],\n                                                         param['hyperparam6'],\n                                                         param['hyperparam7'],\n                                                         param['hyperparam8'],\n                                                         param['hyperparam9'],\n                                                         param['hyperparam10'],\n                                                         param['hyperparam11'],\n                                                         param['hyperparam12'],\n                                                         param['hyperparam13'],\n                                                         param['hyperparam14'],\n                                                         param['hyperparam15'],\n                                                         param['hyperparam16'],\n                                                         param['hyperparam17'],\n                                                         param['hyperparam18'],\n                                                         param['hyperparam19'],\n                                                         param['hyperparam20'],\n                                                         param['hyperparam21'],\n                                                         param['hyperparam22'],\n                                                         param['hyperparam23'],\n                                                         param['hyperparam24']])    \n        \n    # Fit the Stan model\n    \n    fit = sm.sampling(data=train_dat, iter=iter_no, chains=num_of_chains, n_jobs=-1, control = {'adapt_delta': 0.99})\n    \n    print( fit.stansummary(pars=['mu', 'beta', 'sigma', 'patre_sigma', 'patslopere_sigma']) )\n    \n    # Get the posterior samples\n    \n    sigma = fit.extract(permuted=False, pars=\"sigma\")['sigma'].reshape(half_iter_no*num_of_chains)\n    mu = fit.extract(permuted=False, pars=\"mu\")['mu'].reshape(half_iter_no*num_of_chains, 2)\n    beta = fit.extract(permuted=False, pars=\"beta\")['beta'].reshape(half_iter_no*num_of_chains, 7)\n    \n    pat_sigma = fit.extract(permuted=False, pars=\"patre_sigma\")['patre_sigma'].reshape(half_iter_no*num_of_chains)\n    slope_sigma = fit.extract(permuted=False, pars=\"patslopere_sigma\")['patslopere_sigma'].reshape(half_iter_no*num_of_chains)\n    \n    patre = fit.extract(permuted=False, pars=\"patre\")['patre'].reshape(half_iter_no*num_of_chains, total_number_of_patients)\n    patslopere = fit.extract(permuted=False, pars=\"patslopere\")['patslopere'].reshape(half_iter_no*num_of_chains, total_number_of_patients)\n    \n    # Create out-of-fold (or here test set) predictions\n    \n    foldresults = {'fold': fold,                   \n                   'true': test_dat['fvc'].copy().values,\n                   'relweek': test_dat['relweek'].copy().values,\n                   'patno': test_dat['patno'].copy().values,\n                   'Patient': test_dat['Patient'].copy().values,\n                   'Weeks': test_dat['Weeks'].copy().values,\n                   'pred':  np.zeros(shape=len(test_dat['fvc'].values)),\n                   'pred_sd': np.zeros(shape=len(test_dat['fvc'].values))\n                  }\n\n    # Lopping through all the records in the validation or test set\n    # I wish I had implemented a better way of doing this than just spelling out the regression equation (apologies).\n    # (clearly, you can re-write this with a design matrix for fixed and random effects and then write as a matrix multiplcation)\n    for r in range(test_dat['records']):      \n\n        avgpred = mu[:, 0] + beta[:,0]*test_dat['logbasepercent'].values[r] + beta[:,1]*test_dat['basewk'].values[r]/52.0 + \\\n            beta[:,2]*test_dat['basepercent'].values[r] + patre[:, (foldresults['patno'][r]-1) ] * pat_sigma + \\\n            test_dat['relweek'].values[r]/52.0 * ( mu[:, 1] + patslopere[:, (foldresults['patno'][r]-1) ]*slope_sigma + \\\n                                                  beta[:,3]*(test_dat['smoker'].values[r]>=1) + beta[:,4]*(test_dat['smoker'].values[r]==2) + \\\n                                                  beta[:,5]*test_dat['basewk'].values[r]/52.0 +  beta[:,6]*test_dat['basenormal'].values[r] )\n        \n        samples = np.exp( np.repeat(avgpred, 100) + np.random.normal(loc=0, scale=1, size=100*len(sigma)) * np.repeat(sigma, 100) ) * test_dat['denom'].values[r]/100.0        \n        foldresults['pred'][r], foldresults['pred_sd'][r] = minimize_metric(samples)        \n        \n    if fold==0:\n        results = pd.DataFrame(foldresults)\n    else:\n        results = results.append(pd.DataFrame(foldresults))\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Create submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"results['Patient_Week'] = results['Patient'].astype(str) + '_' + results['Weeks'].astype(str)\nresults['FVC'] = results['pred']\nresults['Confidence'] = results['pred_sd']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"results[['Patient_Week', 'FVC', 'Confidence']].to_csv('submission.csv', index=False)","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}