The easy case: no ties
We’ll first compute things when there are no ties in the data.
Individuals may have tied start time. However, for any event time
where a failure occurs
(),
we assume there are no other individuals who have the same event time
nor are there any start times equal to
.
This makes ordering the
’s
unambiguous and the distinction between a start time equalling
or not irrelevant.
We’ll assume that we’ve sorted things by event time. That is,
.
As we have no ties in this section, we can recover “disk order” by
reversing the sort. More generally, with ties we’ll assume that
and that tied event and start times will be sorted in an arbitrary but
fixed order. Therefore, we can unambiguously map between this fixed
order and “disk order”.
Right censored
Ignoring the linear term, and sign, we need to differentiate (by
)
Below,
so it will always be binary even in the presence of ties. Ties will be
dealt with shortly.
Second derivative:
Start / stop
Let the start times be
.
We now make the basic observation:
We need to differentiate (by
)
Proceeding
This can be computed by computing the cumsum of the sequence
Note that this depends only on the risk sums at failure
times! Hence, we only need produce an algorithm that correctly computes
the risk sums at failure times! We will return to this below with
ties.
Difference of cumsums
Let
and
The quantity
can be expressed in terms of the cumsum of the
ordered in start time order and then finding the index of
within the start times
.
The difference
can be expressed as the risk sum
The objective can then be expressed as
with gradient
This can be written as a difference of cumsums of the sequence
The cumsums will be ordered with respect to event / stop time, we will
need to find the index of the start time
within the event / stop times
.
Second derivative
The diagonal of the Hessian is then
(Used
and
so
.)
This can again be expressed in terms of cumsums. The sequence above as
well as the sequence
Again, we will need to cumsum this sequence with respect to event /
stop order.
Probabilistic interpretation
In terms of computing derivatives, it is perhaps helpful to instead
think of differentiating
The map
is the cumulant generating function of an exponential family built on
top of a base measure taking values in the risk set. The law of this
base random variable,
changes with time, in that it depends on the risk set at
.
The random variable has density proportional to
on
Let’s call its law
and the corresponding law under
.
This CGF interpretation gives clean interpretation to the
-th
coordinate of the gradient as
The second derivative has
entry
This also tells us about saturated log-likelihood by
considering multinomials of successive risk sets.
Ties and zero weights
As mentioned previously, we’ll assume that we’ve picked an
arbitrary sort of event times by event time. That is,
.
We can recover “disk order” by inverting this particular sort. So series
indexed by
below should be interpreted as indexed by a fixed but arbitrary sort of
our event times. Also note that any summation over cases
will be restricted to the sum over entries
.
For ties, the issue of determining “which” of a set events occurred
“first” must be disambiguated in some sense. Once a hypothetical
ordering of the events were given, then the risk sums for each failure
at some time be computed.
This is usually not done exactly. The Breslow and Efron methods
“approximate” this computation of risk sums in a certain sense. The two
methods proceed essentially as above, i.e. computing a risk sum for each
individual and then compute the log of the risk sum.
Any reasonable approximation based on computing a risk sum for each
individual (with ties) should have the property that it is exchangable
with respect to a reordering of the failures at a given observed failure
time.
The linear term clearly has no issue, it is the sum of log risk sums
that needs to be addressed. For this, we rewrite the log-likelihood as
As above, we recall that we restricted
.
Again, we note that we only need to correctly compute the risk sums
for the individuals who failed
.
At this point, there are several orderings of the individuals by
event time. The Efron and Breslow, expressed as sums over individuals,
can be viewed as picking an arbitrary such order and proposing
and
that are invariant to this choice of ordering.
Specifically, we group individuals into clusters in such a way that
for any individual
who failed
()
at time
,
their cluster is all other individuals who failed at that same time.
Given such a clustering
,
set
with
We have again applied our principle of summing only over
in computing the average weight for the cluster.
Again, by the form of the likelihood it only matters that
is corrected “correctly” for individuals who fail
.
We will choose a clustering that assigns all individuals to an
individual cluster unless
.
If it is an individual that failed, we find all other individuals who
failed at that time and assign them to the same cluster. In this way,
for each
,
there is a
which is the first index in our arbitrary ordering which is in the same
cluster as
and a
.
The size of each cluster is
.
This will be 1 unless
is a failure
and there is more than 1 failure at this
.
From this
and
we can compute a “centering” weight
to be used for the Efron or Breslow correction. It is easier to describe
in words than in formulae: for a set of indices which is a cluster of
size
(so a set of
failures at some particular
)
the corresponding entries of
are
Here,
is one of the entries of the equation above. When the cluster of tied
event times contains points of weight 0, then these should have
.
That is, suppose we observe a cluster of size 6 with corresponding
weight
then the corresponding
.
This is the
we would get if we dropped those entries with 0 weight from the
(start, stop, event) arrays to begin with. This is
accomplished by a particular lexical sort which takes all points with
equal event times and places those failures with
in the beginning of that segment and into a unique cluster.
We can now define our risk sums
It is clear from the above
that
is linear,
i.e.
for some vector
.
Following through our choice of clustering, we’ll see that
unless
.
It should also be clear that computing the cumulative sum of
over
doesn’t depend on our initial ordering based on the
’s.
Similarly computing the cumsum of
with
constant within clusters (such as the vector
)
will also not depend on our initial ordering. Finally, it should be
clear that for such
the corresponding
can be expressed in terms of the cumsums of
.
The Breslow correction uses
.
The Efron correction uses
.
Probabilistic interpretation
In terms of differentiating with respect to
we could again consider differentiating
This again yields a probabilistic interpretation as in the case of no
ties. In this case, the base measure is handled slightly differently as
to when there are no ties. Breslow’s method fixes the risk set and
weights to be the same for all
in the current cluster of failures. (Recall we have some valid, but
arbitrary ordering by event time.) This is clearly exchangeable. Efron’s
method on the other hand changes the weights as
cycles through the cluster of event times, downweighting all in the
cluster multiplicatively by
as
cycles through the cluster of event times.
First derivative
Proceeding as previously we want to differentiate
Similarly to when there are no ties, we see this is
Second derivative
For the
-th
entry of the diagonal there will be the first term above and we will
subtract a second term
We can expand this, and we’ll see that this can be expressed in terms
of differences of cumsums of the sequences
A lemma or two
To ensure we get the right results, we’ll use a few lemmas.
Lemma 1 (Forward cumsum representation). Given a
sequence
let
denote reversed cumsums of
(padded with a 0 as the
-st
element) with respect to arbitrary, but otherwise valid, ordering of the
’s
and
’s.
Let
be a sequence such that, whenever
we have
Then,
In particular, this implies for
such that
Proof: If
,
then
sums all entries of
for which
.
Similarly, under our assumptions
will be the sum over start times greater than or equal to
.
Note: When we don’t have start times
(0-based indexing) and
Note: This says our risk sum calculation is correct
for the failure times. This seems false without the
,
i.e. we only have the correct sum at failure times. But that’s all we
need thankfully.
Note: Recall that to get this expression into the
“disk order”, we invert by our fixed sort of event times.
Lemma 2 (Adjoint cumsum representation). Given a
sequence
let
denote the cumsums of
ordered with respect to a valid, but otherwise arbitrary order, (with a
0 padded in as the first element). Then,
Above,
.
Further
Note: With no start times,
(0-based indexing) and
.
Note: The elementwise multiplication by
is also crucial here.
Derivatives revisited
In the light of the lemmata above, we can rewrite this first
derivative now as
Similarly, that term in the second derivative takes the form
The last conclusion of Lemma 2 allows us now to resolve products of
all the indicators above. The only that requires any thought is the
product
Full Hessian
Suppose we want to multiply the Hessian by a vector
.
Let’s use
for the
-th
row and
-th
column of the Hessian.
Non-diagonal term
We’ll have a term like
We are free to rewrite this as:
We’ll write it this way as it will allow us to get rid of the two sums
conveniently.
Multiplying this on the right by
yields (with
the elementary basis vectors):
The sum over
is simply computing
Lemma 1 tells us how to do this correctly for the
with
.
This is done by calling the function _sum_over_risk_set in
the code, applied to an argument
.
We now have to compute the following sum over
Up to the scaling by
computation can be handled by Lemma 2 applied to the sequences
This is handled by the function _sum_over_events in the
code.
Diagonal term
The diagonal term will have
entry
Summing over
in a multiplication by
on the right yields
We again apply Lemma 2 to
the sequences
Probabilistic interpretation
The Hessian has
entry
Multiplying on the right by
yields
Appendix: Proofs of Lemmas
Proof of Lemma 1 (Forward cumsum representation)
Consider the evaluation of the risk sum for an individual
that experiences an event
().
The sum can be decomposed into two parts based on the definition of the
risk set at time
,
which includes all individuals
such that
:
Because the start time
is always strictly less than the event time
(),
the condition
naturally implies
.
Thus, the second sum simplifies:
Let us assume 0-based indexing for the sequences. The sequence
is sorted in ascending order of event times
.
The index
represents the first occurrence of the event time
in this sorted sequence. Therefore, summing
for all
where
is equivalent to summing the sequence from index
to the end
().
By definition, the reversed cumsum
computes exactly this quantity.
Similarly, let
be the reversed cumsum of
sorted by start times
.
The quantity
gives the exact number of start times strictly less than
.
In a 0-based indexed array sorted by start times, the indices
correspond to
,
and the indices from
onwards correspond to
.
Thus, the reversed cumsum
computes the sum of
for
.
Multiplying by
(which trivially yields
on both sides if
)
establishes the first equality:
For the second part of the Lemma, recall the definition of
:
where
is the sum over the cluster of tied event times at
.
The cluster spans indices from
to
.
In terms of the reversed cumsum
,
the sum over this finite slice is simply
.
Substituting this difference into the expression for
yields:
which concludes the proof.
Proof of Lemma 2 (Adjoint cumsum representation)
We want to evaluate the sum
.
The indicator
equals
if and only if
,
and is
otherwise. Thus, the sum effectively accumulates
over the interval
.
Assuming 0-based indexing for the arrays sorted by event time
,
let us identify the bounds of this interval in terms of the sequence
indices. The upper bound is defined by
.
In the sorted array, the elements satisfying this condition occupy
indices
through
.
The total number of such elements is
.
The lower bound is defined by
,
which means we exclude elements where
.
The number of elements with
is precisely given by
.
Thus, the excluded indices are
through
.
Consequently, the interval
corresponds exactly to the index range
.
Let
denote the forward cumsum sequence of
,
padded with a
at the 0-th index. That is,
.
Using this forward cumsum, the sum over our target index range can be
computed as the difference between the sum up to
and the sum up to
:
This establishes the first equality.
The second equality,
,
holds because all individuals
that share the exact same event time
must either entirely fall within or outside the threshold boundary.
Hence, the index count maps cleanly to the
function of that start boundary, ensuring numerical stability in
implementations with ties.