<jupyter_start><jupyter_text>We use two datasets for validation of our Cox model: one synthetic and one real  Validation against closed-form MLE fit A synthetic dataset for which we can fit the Cox partial likelihood Maximum Likelihood Estimate by hand, to validate numerical PyTorch implementation against.<jupyter_code>valdf = pd.DataFrame({'id':['Bob','Sally','James','Ann'], 'time':[1,3,6,10], 'status':[1,1,0,1], 'smoke':[1,0,0,1]})
valdf<jupyter_output><empty_output><jupyter_text>Specify model Cox(t,d ~ smoke) for times $t$ and status indicator $d$.  $L(T, D, X) = \prod_{t_i^e} \frac{h(t_i^e,X_i)}{\sum_{j:\, t_j \geq t_i^e} h(t_j, X_j)}$  where $t_i^e$ denote event times (i.e. times where $d=1$), and $t_i$ denotes all times (whether event or censored).  In this way the ratio $\frac{h(t_i^e,X_i)}{\sum_{j:\, t_j \geq t_i} h(t_j, X_j)}$ is the ratio of the hazard for the event observed at time $t_i^e$ divided by the sum of hazards for all subjects at risk at that time (so with $t_j \geq t_i^e$, i.e. have not had an event yet, have not died, so still at risk, including the subject who had an event at that time since that person must have been at risk too, so inequality includes time $t_i^e$).   We therefore see that for each of the subjects in the data just above the hazard, $h(t, X) = h_0(t)\,\exp(X \beta)$, is:  - Bob: $h(t, X) = h_0(t)\,\exp(\beta)$- Sally: $h(t, X) = h_0(t)\,\exp(0)$- James: $h(t, X) = h_0(t)\,\exp(0)$- Ann: $h(t, X) = h_0(t)\,\exp(\beta)$ The likelihood for this dataset then has three factors, one for each event time, and in each factor is the hazard of the person who had the event divided by the sum of the hazards of all the people who were still at risk then (including the person who had the event:    $\begin{align}L(T, D, X) = &\left[\frac{h_0(t)\exp(\beta)}{h_0(t)\exp(\beta) \,+\, h_0(t)\exp(0) \,+\, h_0(t)\exp(0) \,+\, h_0(t)\exp(\beta)} \right] \\&\times \left[ \frac{h_0(t)\exp(0)}{h_0(t)\exp(0) \,+\, h_0(t)\exp(0) \,+\, h_0(t)\exp(\beta)} \right] \\&\times \left[ \frac{h_0(t)\exp(\beta)}{h_0(t)\exp(\beta)} \right]\end{align}$ Note that the baseline hazard, $h_0(t)$ will always cancel everywhere in the likelihood. We then have:    $\begin{align}L(T, D, X) &= \left[\frac{\exp(\beta)}{2(1 + \exp(\beta))} \right] \times \left[ \frac{1}{2+ \exp(\beta)} \right] \times 1\end{align}$ Take the natural logarithm of this as the logarithm is a monotone transformation and will not change the position of the minimum, but will simplify the computation of the derivative and make the computation more numerically stable. Also we multipy it by $-1$ so that to maximise the likelihood we minimise the quantity below. $-\ln[L(T,D,X)] = -\beta +\ln[2(1+\exp(\beta)] - \ln[1] + \ln[2+\exp(\beta)]$ Taking the derivative of this wrt $\beta$ gives  $\begin{align}-\frac{d}{d\beta} \ln[L(T,D,X)] &= -1 + \frac{d}{d\beta}\ln[(1+\exp(\beta)] + \frac{d}{d\beta}\ln[2+\exp(\beta)] \\&= -1 + \frac{1}{1+\exp(\beta)}\exp(\beta) + \frac{1}{2+\exp(\beta)}\exp(\beta)\end{align}$ Setting this equal to zero, $-\frac{d}{d\beta} \ln[L(T,D,X)] = 0$, gives  $1 = \frac{\exp(\beta)}{1+\exp(\beta)} + \frac{\exp(\beta)}{2+\exp(\beta)}$ So that this is the equation we need to solve for $\beta$ to find the stationary point (hopefully a minimum, we haven't shown it but this score function is convex so this is indeed a minimum), and this will be the Maximum Likelihood Estimate (MLE) for $\beta$. This gives $\beta = \ln(2)/2 \approx 0.34657$ (Note that we chose the solution in the Reals here, the solution is in fact $\beta = \ln(\sqrt{2})$ where the square root has two solutions and only the positive one gives a real value for $\beta$.) Now we fit our Cox model code to verify that we indeed obtain the closed-form answer as expected<jupyter_code>tname = 'time'
Xnames = ['smoke']
dname = 'status'

coxmod1 = TorchCox.TorchCox(lr=1)

coxmod1.fit(valdf, Xnames=Xnames, tname=tname, dname=dname)

coxmod1.beta.detach().numpy()[0]<jupyter_output>[0.34657338]<jupyter_text>And we do! :)<jupyter_code>np.log(2)/2<jupyter_output><empty_output><jupyter_text>Validation against R's `survival` package Now we will compare our result against the R package 'survival' on another dataset.<jupyter_code>df_tied = pd.read_csv('../data/ovarian_deduplicated.csv')
df_tied
coxmod = TorchCox.TorchCox()

tname = 'tyears'
Xnames = ['Karn', 'Ascites']
dname = 'd'
#sudo pip install rpy2
%load_ext rpy2.ipython
%%R
library(readr)
library(survival)
library(dplyr)
library(tidyr)

df2 = read_csv("/home/ilan/Desktop/TorchCox/data/ovarian_deduplicated.csv")

starttime = Sys.time()

rmod = coxph(Surv(tyears, d) ~ Karn + Ascites, df2, ties="breslow")
print(coef(rmod))

endtime = Sys.time()

print(endtime-starttime)
%%time

coxmod.fit(df_tied, Xnames=Xnames, tname=tname, dname=dname, basehaz=False)<jupyter_output>[-0.23638311  0.4170416 ]
CPU times: user 944 ms, sys: 9.92 ms, total: 954 ms
Wall time: 276 ms<jupyter_text>We indeed match that result as well!  We are about 10x slower than the R package (which runs on C code in the background). But simplicity and extensibility of our code compensates for that in my view.  (Also timings here most likely dominated by overhead of loading libraries, comparison on larger dataset required.) To be more specific, R's survival package is a blazingly fast and mature package, however the cost of this is that the code itself is non-trivial to edit or modify, see e.g.https://github.com/therneau/survival/blob/master/src/coxfit6.c Contrast this with the code for this implementation which is comparably readable and can be found in `torchcox/TorchCox.py`. You can easily see how you could modify the likelihood there to insert regularisation, change the optimiser, modify the log-linear dependence on the covariates, etc. Fit the model again, but this time computing the baseline hazard, to ensure that works as well.<jupyter_code>%%time

coxmod.fit(df_tied, Xnames=Xnames, tname=tname, dname=dname, basehaz=True)
coxmod.basehaz<jupyter_output><empty_output><jupyter_text>Predict on the training set to ensure predict_proba() method works as well.<jupyter_code>df_tied['pred'] = coxmod.predict_proba(df_tied, Xnames=Xnames, tname=tname)
df_tied<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14854026307346807, 'fraction_numerical': 0.01940968880333654, 'mean_word_length': 4.882075471698113, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 23, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '294013', 'license': 'MIT', 'n_tokens_mistral': 2257, 'n_tokens_neox': 2160, 'n_words': 753, 'path': 'notebooks/Validation.ipynb', 'repo_name': 'ilanfri/TorchCox', 'size': '6234'}