To fit Cox Proportional Hazards we have to maximize the partial likihood of beta. Thankfully, we don’t have to face up to this somewhat menacing optimization ourselves because we can use the great package lifelines written by Cam Davidson-Pilon.

I recommend the tutorial on the lifelines website here on how to do it, but I have included a code snippet below inspired by that tutorial because I want to add on loglik, which will be important for Bayesian model averaging

```
%matplotlib inline
import matplotlib.pyplot as plt
from lifelines.datasets import load_rossi
from lifelines import CoxPHFitter
cf = lifelines.CoxPHFitter()
cf.fit(load_rossi(), 'week', event_col = 'arrest')
#Plot the results
fig, axes = plt.subplots(nrows=1, ncols=2, sharex=True)
cf.baseline_cumulative_hazard_.plot(ax = axes[0], title = "Baseline Cumulative Hazard")
cf.baseline_survival_.plot(ax = axes[1], title = "Baseline Survival")
#Pedict a value
import numpy as np
coefficients = np.array([[0,0,0,0,0,0,0]])
cf.predict_survival_function(coefficients)
fig, axis = plt.subplots(nrows=1, ncols=1, sharex=True)
survivalZeros.plot(ax = axis, title = "Baseline Cumulative Hazard for coefficients")
#Get the log likihood
loglik = cf._log_likelihood
```