Skip to contents

Abstract

We develop the score and Hessian for Cox’s partial likelihood with weights as well as ties (including both “Breslow” and “Efron” tie-breaking).


There are two terms in Cox’s partial log-likelihood i=1ndiwi[ηilog(j:jiwjeηj)]. \sum_{i=1}^n d_i w_i \left[ \eta_i - \log \left(\sum_{j:j \geq i} w_j e^{\eta_j} \right) \right].

We’ll focus on the second, as the first is simple to deal with, even with ties.

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 tit_i where a failure occurs (di=1d_i=1), we assume there are no other individuals who have the same event time nor are there any start times equal to tit_i. This makes ordering the tit_i’s unambiguous and the distinction between a start time equalling tit_i or not irrelevant.

We’ll assume that we’ve sorted things by event time. That is, j>itj>tij > i \implies t_j > t_i. 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 j>itjtij > i \implies t_j \geq t_i 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 ηk\eta_k) i=1ndiwilog(jiwjeηj) \sum_{i = 1}^n d_i w_i \log \left(\sum_{j \geq i} w_j e^{\eta_j}\right) Below, di=statusid_i=\text{status}_i so it will always be binary even in the presence of ties. Ties will be dealt with shortly.

ηki=1ndiwilog(jiwjeηj)=i=1ndiwijiwjδjkeηkliwleηl=i=1ndiwi1{ki}wkeηkliwleηl=wkeηki=1kdiwi1liwleηl \begin{aligned} \frac{\partial}{\partial \eta_k} \sum_{i = 1}^n d_i w_i \log \left(\sum_{j \geq i} w_j e^{\eta_j}\right) &= \sum_{i=1}^n d_i w_i \frac{ \sum_{j \geq i} w_j \delta_{jk} e^{\eta_k}} {\sum_{l \geq i} w_l e^{\eta_l}} \\ &= \sum_{i=1}^n d_i w_i \frac{ 1_{\{k \geq i\}} w_k e^{\eta_k}} {\sum_{l \geq i} w_l e^{\eta_l}} \\ &= w_k e^{\eta_k} \sum_{i=1}^k d_i w_i \frac{ 1} {\sum_{l \geq i} w_l e^{\eta_l}} \\ \end{aligned}

Second derivative: 2ηk2i=1ndiwilog(jiwjeηj)=wk[wkeηki=1kdiwi1liwleηl]=wkeηk[i=1kdiwi1liwleηli=1kdiwiliwlδlkeηl(liwleηl)2]=wkeηk[i=1kdiwi1liwleηli=1kdiwi1{ki}wkeηk(liwleηl)2]=wkeηk[i=1kdiwi1liwleηli=1kdiwiwkeηk(liwleηl)2] \begin{aligned} \frac{\partial^2}{\partial \eta_k^2} \sum_{i = 1}^n d_i w_i \log \left(\sum_{j \geq i} w_j e^{\eta_j}\right) &= \frac{\partial}{\partial w_k} \left[w_k e^{\eta_k} \sum_{i=1}^k d_i w_i \frac{ 1} {\sum_{l \geq i} w_l e^{\eta_l}} \right] \\ &= w_k e^{\eta_k} \left[\sum_{i=1}^k d_i w_i \frac{ 1} {\sum_{l \geq i} w_l e^{\eta_l}} - \sum_{i=1}^k d_i w_i \frac{ \sum_{l \geq i} w_l \delta_{lk} e^{\eta_l}} {(\sum_{l \geq i} w_l e^{\eta_l})^2} \right] \\ &= w_k e^{\eta_k} \left[\sum_{i=1}^k d_i w_i \frac{ 1} {\sum_{l \geq i} w_l e^{\eta_l}} - \sum_{i=1}^k d_i w_i \frac{ 1_{\{k \geq i\}} w_k e^{\eta_k}} {(\sum_{l \geq i} w_l e^{\eta_l})^2} \right] \\ &= w_k e^{\eta_k} \left[\sum_{i=1}^k d_i w_i \frac{ 1} {\sum_{l \geq i} w_l e^{\eta_l}} - \sum_{i=1}^k d_i w_i \frac{ w_k e^{\eta_k}} {(\sum_{l \geq i} w_l e^{\eta_l})^2} \right] \\ \end{aligned}

Start / stop

Let the start times be (sm)1mn(s_m)_{1 \leq m \leq n}.

We now make the basic observation: ηk(j:jiwjeηjm:smtiwmeηm)=wkeηk[(j:jiδjkm:smtiδmk)]=wkeηk[1{ki}1{skti}] \begin{aligned} \frac{\partial}{\partial \eta_k} \left(\sum_{j:j \geq i} w_j e^{\eta_j} - \sum_{m:s_m \geq t_i} w_m e^{\eta_m} \right) &= w_k e^{\eta_k}\left[ \left(\sum_{j:j \geq i} \delta_{jk} - \sum_{m: s_m \geq t_i} \delta_{mk} \right) \right] \\ &= w_k e^{\eta_k} \left[ 1_{\{k \geq i\}} - 1_{\{s_k \geq t_i\}} \right] \end{aligned}

We need to differentiate (by ηk\eta_k) i=1ndiwilog(jiwjeηjm:smtiwmeηm) \sum_{i = 1}^n d_i w_i \log \left(\sum_{j \geq i} w_j e^{\eta_j} - \sum_{m:s_m \geq t_i} w_m e^{\eta_m} \right)

Proceeding ηki=1ndiwilog(jiwjeηjm:smtiwmeηm)=wkeηki=1ndiwi1{ki}1{skti}jiwjeηjm:smtiwmeηm=wkeηk[i=1kdiwiliwleηlm:smtiwmeηmi:tiskdiwiliwleηlm:smtiwmeηm] \begin{aligned} \frac{\partial}{\partial \eta_k} \sum_{i = 1}^n d_i w_i \log \left(\sum_{j \geq i} w_j e^{\eta_j}- \sum_{m:s_m \geq t_i} w_m e^{\eta_m} \right) &= w_k e^{\eta_k} \sum_{i=1}^n d_i w_i \frac{1_{\{k \geq i\}} - 1_{\{s_k \geq t_i\}} } {\sum_{j \geq i} w_j e^{\eta_j}- \sum_{m:s_m \geq t_i} w_m e^{\eta_m} } \\ &= w_k e^{\eta_k}\biggl[ \sum_{i=1}^k \frac{ d_i w_i} {\sum_{l \geq i} w_l e^{\eta_l}- \sum_{m:s_m \geq t_i} w_m e^{\eta_m} } \\ & \qquad \quad - \sum_{i:t_i \leq s_k} \frac{ d_i w_i} {\sum_{l \geq i} w_l e^{\eta_l}- \sum_{m:s_m \geq t_i} w_m e^{\eta_m} }\biggr] \\ \end{aligned}

This can be computed by computing the cumsum of the sequence (diwiliwleηlm:smtiwmeηm)1in. \left( \frac{ d_i w_i} {\sum_{l \geq i} w_l e^{\eta_l}- \sum_{m:s_m \geq t_i} w_m e^{\eta_m} } \right)_{1 \leq i \leq n}.

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 Ei(ξ)=ji:wj>0ξj=j:tjti,wj>0ξj E_i(\xi) = \sum_{j \geq i: w_j > 0} \xi_j = \sum_{j: t_j \geq t_i, w_j > 0} \xi_j and Si(ξ)=j:sjti,wj>0ξ. S_i(\xi) = \sum_{j: s_j \geq t_i, w_j > 0} \xi. The quantity SiS_i can be expressed in terms of the cumsum of the wjeηjw_j e^{\eta_j} ordered in start time order and then finding the index of tit_i within the start times (si)1in(s_i)_{1 \leq i \leq n}.

The difference j:jiwjeηjm:smtiwmeηm \sum_{j: j \geq i} w_j e^{\eta_j} - \sum_{m:s_m \geq t_i} w_m e^{\eta_m} can be expressed as the risk sum Ri(weη)=Ei(weη)Si(weη). R_i(we^{\eta}) = E_i(we^{\eta}) - S_i(we^{\eta}).

The objective can then be expressed as i=1ndiwilog(Ri(weη)) \sum_{i=1}^n d_i w_i \log(R_i(we^{\eta})) with gradient wkeηki=1ndiwi[1{ki}1{skti}]Ei(weη)Si(weη)=wkeηki=1ndiwi[1{ki}1{skti}]Ri(weη) w_k e^{\eta_k} \sum_{i=1}^n \frac{d_i w_i [1_{\{k \geq i\}} - 1_{\{s_k \geq t_i\}}]}{E_i(we^{\eta})-S_i(we^{\eta})} = w_k e^{\eta_k} \sum_{i=1}^n \frac{d_i w_i [1_{\{k \geq i\}} - 1_{\{s_k \geq t_i\}}]}{R_i(we^{\eta})}

This can be written as a difference of cumsums of the sequence (diwiRi(weη))1in \left(\frac{d_i w_i}{R_i(we^{\eta})}\right)_{1 \leq i \leq n} The cumsums will be ordered with respect to event / stop time, we will need to find the index of the start time sks_k within the event / stop times (ti)1in(t_i)_{1 \leq i \leq n}.

Second derivative

The diagonal of the Hessian is then wkeηki=1ndiwi[1{ki}1{skti}]Ei(weη)Si(weη)wk2e2ηki=1ndiwi[1{ki}1{skti}](Ei(weη)Si(weη))2 w_k e^{\eta_k} \sum_{i=1}^n \frac{d_i w_i [1_{\{k \geq i\}} - 1_{\{s_k \geq t_i\}}]}{E_i(we^{\eta})-S_i(we^{\eta})} - w_k^2 e^{2 \eta_k} \sum_{i=1}^n \frac{d_i w_i [1_{\{k \geq i\}} - 1_{\{s_k \geq t_i\}}]}{(E_i(we^{\eta})-S_i(we^{\eta}))^2} (Used 1{ki}1{skti}{0,1}1_{\{k \geq i\}}-1_{\{s_k \geq t_i\}} \in \{0,1\} and 12=1,02=01^2=1, 0^2=0 so (1{ki}1{skti})2=1{ki}1{skti}(1_{\{k \geq i\}}-1_{\{s_k \geq t_i\}})^2=1_{\{k \geq i\}}-1_{\{s_k \geq t_i\}}.) This can again be expressed in terms of cumsums. The sequence above as well as the sequence (diwi(Ei(weη)Si(weη))2)1in \left(\frac{d_i w_i}{(E_i(we^{\eta})-S_i(we^{\eta}))^2}\right)_{1 \leq i \leq n}

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 i=1ndiwilog(Ri(weη)/Ri(w))=defi=1ndiwiΛi(η). \sum_{i=1}^n d_i w_i \log(R_i(we^{\eta})/R_i(w)) \overset{def}{=} \sum_{i=1}^n d_i w_i \Lambda_i(\eta). The map ηlog(Ri(weη)/Ri(w))\eta \to \log(R_i(we^{\eta})/R_i(w)) 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, ViV_i changes with time, in that it depends on the risk set at tit_i. The random variable has density proportional to wjw_j on {j:tjti,sj<ti}.\{j:t_j \geq t_i, s_j < t_i\}. Let’s call its law PiP_i and the corresponding law under η\etaPi,ηP_{i,\eta}.

This CGF interpretation gives clean interpretation to the kk-th coordinate of the gradient as i=1ndiwiEi,η[1{k}(Vi)]=i=1ndiwiwkeηk1{ksupp(Pi)}Ei[eη]=i=1ndiwiπi,η[k] \begin{aligned} \sum_{i=1}^n d_i w_i E_{i,\eta}[1_{\{k\}}(V_i)] = \sum_{i=1}^n d_i w_i \frac{w_k e^{\eta_k} 1_{\{k \in \text{supp}(P_i)\}}}{E_i[e^\eta]} = \sum_{i=1}^n d_i w_i \pi_{i,\eta}[k] \end{aligned}

The second derivative has (r,c)(r,c) entry i=1ndiwi(πi,η[r]δrcπi,η[r]πi,η[c]). \sum_{i=1}^n d_i w_i\left( \pi_{i,\eta}[r] \delta_{rc} - \pi_{i,\eta}[r] \pi_{i,\eta}[c] \right).

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, jitjtij \geq i \implies t_j \geq t_i. We can recover “disk order” by inverting this particular sort. So series indexed by ii below should be interpreted as indexed by a fixed but arbitrary sort of our event times. Also note that any summation over cases jj will be restricted to the sum over entries wj>0w_j > 0.

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 i=1,wi>0ndi[wiηiwilog(Ri(weη))]. \sum_{i=1, w_i > 0}^n d_i \left[w_i \eta_i - \bar{w}_i \log \left( \bar{R}_i(we^{\eta}) \right) \right].

As above, we recall that we restricted di{0,1}d_i \in \{0,1\}. Again, we note that we only need to correctly compute the risk sums Ri\bar{R}_i for the individuals who failed di=1d_i=1.

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 w\bar{w} and R\bar{R} that are invariant to this choice of ordering.

Specifically, we group individuals into clusters in such a way that for any individual ii who failed (di=1d_i=1) at time tit_i, their cluster is all other individuals who failed at that same time. Given such a clustering CC, set wi={j:Cj=Ci,wj>0wjj:Cj=Ci,wj>01=j:Cj=Ci,wj>0wjKjwi>00otherwise. \bar{w}_i = \begin{cases} \frac{\sum_{j:C_j=C_i, w_j>0} w_j}{\sum_{j:C_j=C_i, w_j>0} 1} = \frac{\sum_{j:C_j=C_i, w_j>0} w_j}{K_j} & w_i > 0 \\ 0 & \text{otherwise.} \end{cases} with Kj=#{j:Cj=Ci,wj>0}. K_j = \# \{j:C_j=C_i, w_j>0\}.

We have again applied our principle of summing only over {j:wj>0}\{j:w_j>0\} in computing the average weight for the cluster.

Again, by the form of the likelihood it only matters that wi\bar{w}_i is corrected “correctly” for individuals who fail (di=1)(d_i=1).

We will choose a clustering that assigns all individuals to an individual cluster unless di=1d_i=1. 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 ii, there is a 𝚏𝚒𝚛𝚜𝚝(i)\mathtt{first}(i) which is the first index in our arbitrary ordering which is in the same cluster as ii and a 𝚕𝚊𝚜𝚝(i)\mathtt{last}(i). The size of each cluster is 𝚕𝚊𝚜𝚝(i)+1𝚏𝚒𝚛𝚜𝚝(i)\mathtt{last}(i) + 1 - \mathtt{first}(i). This will be 1 unless tit_i is a failure (di=1)(d_i=1) and there is more than 1 failure at this tit_i.

From this 𝚏𝚒𝚛𝚜𝚝\mathtt{first} and 𝚕𝚊𝚜𝚝\mathtt{last} we can compute a “centering” weight σ\sigma 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 KK (so a set of KK failures at some particular tt) the corresponding entries of σ\sigma are {0,1/K,2/K,,(K1)/K}.\{0,1/K,2/K,\dots,(K-1)/K\}. Here, KK 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 σi=0\sigma_i=0. That is, suppose we observe a cluster of size 6 with corresponding weight w=[1/2,0,0,3/4,0,2]w=[1/2,0,0,3/4,0,2] then the corresponding σ=[0,0,0,1/3,0,2/3]\sigma=[0,0,0,1/3,0,2/3]. This is the σ\sigma 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 wi>0w_i>0 in the beginning of that segment and into a unique cluster.

We can now define our risk sums Δi(ξ)=j:Cj=Ci,wj>0ξjRi(ξ)=j:tjti,wj>0ξjm:smti,wm>0ξmδσiΔi(ξ) \begin{aligned} \Delta_i(\xi) &= \sum_{j:C_j=C_i, w_j>0} \xi_j \\ \bar{R}_i(\xi) &= \sum_{j: t_j \geq t_{i}, w_j>0} \xi_j - \sum_{m: s_m \geq t_i, w_m>0} \xi_m - \delta \cdot \sigma_i \cdot \Delta_i(\xi) \end{aligned} It is clear from the above that R\bar{R} is linear, i.e. Ri(ξ)=riξ\bar{R}_i(\xi)=\bar{r}_i'\xi for some vector ri\bar{r}_i.

Following through our choice of clustering, we’ll see that Ri(weη)=Ri(weη)R_i(we^{\eta})=\bar{R}_i(we^{\eta}) unless di=1d_i=1. It should also be clear that computing the cumulative sum of Ri(weη)\bar{R}_i(we^{\eta}) over ii doesn’t depend on our initial ordering based on the tit_i’s. Similarly computing the cumsum of ζiRi(weη)\zeta_i \bar{R}_i(we^{\eta}) with ζi\zeta_i constant within clusters (such as the vector w\bar{w}) will also not depend on our initial ordering. Finally, it should be clear that for such ζ\zeta the corresponding Δi\Delta_i can be expressed in terms of the cumsums of (wjζj)1jn(w_j \zeta_j)_{1 \leq j \leq n}.

The Breslow correction uses δ=0\delta=0. The Efron correction uses δ=1\delta=1.

Probabilistic interpretation

In terms of differentiating with respect to η\eta we could again consider differentiating i=1ndiwilog(Ri(weη)/Ri(w)). \sum_{i=1}^n d_i \bar{w}_i \log(\bar{R}_i(we^{\eta})/\bar{R}_i(w)).

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 ii 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 ii cycles through the cluster of event times, downweighting all in the cluster multiplicatively by 1/K1/K as ii cycles through the cluster of event times.

First derivative

Proceeding as previously we want to differentiate i=1ndiwilog(Ri(weη)/Ri(w)). \sum_{i=1}^n d_i \bar{w}_i \log \left( \bar{R}_i(we^{\eta}) / \bar{R}_i(w) \right).

Similarly to when there are no ties, we see this is wkeηki=1ndiwi[1{tkti}1{skti}δσi(1{i𝚕𝚊𝚜𝚝(k)}1{i𝚏𝚒𝚛𝚜𝚝(k)1})]Ri(weη)=i=1ndiwiπi,η[k]. w_k e^{\eta_k}\sum_{i=1}^n \frac{d_i \bar{w}_i [1_{\{t_k \geq t_i\}} - 1_{\{s_k \geq t_i\}} - \delta \cdot \sigma_i \cdot( 1_{\{i \leq \mathtt{last}(k)\}} - 1_{\{i \leq \mathtt{first}(k)-1\}}) ]}{\bar{R}_i(we^{\eta})} = \sum_{i=1}^n d_i \bar{w}_i \pi_{i,\eta}[k].

Second derivative

For the kk-th entry of the diagonal there will be the first term above and we will subtract a second term wk2e2ηki=1ndiwi[1{tkti}1{skti}δσi(1{i𝚕𝚊𝚜𝚝(k)}1{i𝚏𝚒𝚛𝚜𝚝(k)1})]2Ri(weη)2. w^2_k e^{2\eta_k}\sum_{i=1}^n \frac{d_i \bar{w}_i [1_{\{t_k \geq t_i\}} - 1_{\{s_k \geq t_i\}} - \delta \cdot \sigma_i\cdot( 1_{\{i \leq \mathtt{last}(k)\}} - 1_{\{i \leq \mathtt{first}(k)-1\}})]^2}{\bar{R}_i(we^{\eta})^2}.

We can expand this, and we’ll see that this can be expressed in terms of differences of cumsums of the sequences (diwiRi(weη)2,diwiσiRi(weη)2,diwiσi2Ri(weη)2)1in. \left(\frac{d_i \bar{w}_i}{\bar{R}_i(we^{\eta})^2}, \frac{d_i \bar{w}_i \sigma_i}{\bar{R}_i(we^{\eta})^2}, \frac{d_i \bar{w}_i \sigma^2_i}{\bar{R}_i(we^{\eta})^2}\right)_{1 \leq i \leq n}.

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 (Zi)1in(Z_i)_{1 \leq i \leq n} let E(Z)i,S(Z)iE(Z)_i, S(Z)_i denote reversed cumsums of ZZ (padded with a 0 as the n+1n+1-st element) with respect to arbitrary, but otherwise valid, ordering of the tit_i’s and sis_i’s. Let 𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)\mathtt{eventmap}(i) be a sequence such that, whenever di=1d_i=1 we have 𝚎𝚟𝚎𝚗𝚝(i)=#{j:sj<ti}.\mathtt{event}(i) = \#\{j:s_j < t_i\}. Then, di(j:tjtiZjm:smtiZm)=di(E(Z)𝚏𝚒𝚛𝚜𝚝(i)S(Z)𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)). d_i \cdot \left( \sum_{j: t_j \geq t_i} Z_j - \sum_{m: s_m \geq t_i} Z_m \right) = d_i \cdot (E(Z)_{\mathtt{first}(i)} - S(Z)_{\mathtt{eventmap}(i)}). In particular, this implies for ii such that di=1d_i=1 riZ=E(Z)𝚏𝚒𝚛𝚜𝚝(i)S(Z)𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)δσi(E𝚏𝚒𝚛𝚜𝚝(i)(Z)E𝚕𝚊𝚜𝚝(i)+1(Z)). \bar{r}_i'Z = E(Z)_{\mathtt{first}(i)} - S(Z)_{\mathtt{eventmap}(i)} - \delta \cdot \sigma_i \cdot (E_{\mathtt{first}(i)}(Z) - E_{\mathtt{last}(i)+1}(Z)).

Proof: If di=1d_i=1, then E𝚏𝚒𝚛𝚜𝚝(i)E_{\mathtt{first}(i)} sums all entries of ZZ for which tjtit_j \geq t_i. Similarly, under our assumptions S(Z)𝚎𝚟𝚎𝚗𝚝(i)S(Z)_{\mathtt{event}(i)} will be the sum over start times greater than or equal to tit_i.

Note: When we don’t have start times 𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)=n\mathtt{eventmap}(i)=n (0-based indexing) and S(Z)𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)0.S(Z)_{\mathtt{eventmap}(i)} \equiv 0.

Note: This says our risk sum calculation is correct for the failure times. This seems false without the did_i, 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 (Zi)1in(Z_i)_{1 \leq i \leq n} let E(dZ)i\bar{E}(d \cdot Z)_i denote the cumsums of dZd \cdot Z ordered with respect to a valid, but otherwise arbitrary order, (with a 0 padded in as the first element). Then, i=1ndiZi(1{tkti}1{skti})=E(dZ)𝚕𝚊𝚜𝚝(k)+1E(dZ)𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)=E(dZ)𝚕𝚊𝚜𝚝(k)+1E(dZ)𝚏𝚒𝚛𝚜𝚝𝚜𝚝𝚊𝚛𝚝(k) \begin{aligned} \sum_{i=1}^n d_i Z_i (1_{\{t_k \geq t_i\}} - 1_{\{s_k \geq t_i\}}) &= \bar{E}(d \cdot Z)_{\mathtt{last}(k)+1} - \bar{E}(d \cdot Z)_{\mathtt{startmap}(k)} \\ &= \bar{E}(d \cdot Z)_{\mathtt{last}(k)+1} - \bar{E}(d \cdot Z)_{\mathtt{firststart}(k)} \end{aligned} Above, 𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)=#{j:tjsk}\mathtt{startmap}(k) = \# \{j: t_j \leq s_k\}. Further 𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(i)=𝚏𝚒𝚛𝚜𝚝(𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k))𝚏𝚒𝚛𝚜𝚝(k).\mathtt{startmap}(i) = \mathtt{first}(\mathtt{startmap}(k)) \leq \mathtt{first}(k).

Note: With no start times, 𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙0\mathtt{startmap} \equiv 0 (0-based indexing) and E(Z)𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(i)0\bar{E}(Z)_{\mathtt{startmap}(i)} \equiv 0.

Note: The elementwise multiplication by dd is also crucial here.

Derivatives revisited

In the light of the lemmata above, we can rewrite this first derivative now as wkeηki=1ndiwi[1{i𝚕𝚊𝚜𝚝(k)}1{i<𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)}δσi(1{i𝚕𝚊𝚜𝚝(k)}1{i𝚏𝚒𝚛𝚜𝚝(k)1})]Ri(weη) w_k e^{\eta_k}\sum_{i=1}^n \frac{d_i \bar{w}_i [1_{\{i \leq \mathtt{last}(k)\}} - 1_{\{i < \mathtt{startmap}(k)\}} - \delta \cdot \sigma_i \cdot( 1_{\{i \leq \mathtt{last}(k)\}} - 1_{\{i \leq \mathtt{first}(k)-1\}}) ]}{\bar{R}_i(we^{\eta})}

Similarly, that term in the second derivative takes the form wk2e2ηki=1ndiwi[1{i𝚕𝚊𝚜𝚝(k)}1{i<𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)}δσi(1{i𝚕𝚊𝚜𝚝(k)}1{i𝚏𝚒𝚛𝚜𝚝(k)1})]2Ri(weη)2 w_k^2 e^{2 \eta_k}\sum_{i=1}^n \frac{d_i \bar{w}_i [1_{\{i \leq \mathtt{last}(k)\}} - 1_{\{i < \mathtt{startmap}(k)\}} - \delta \cdot \sigma_i \cdot( 1_{\{i \leq \mathtt{last}(k)\}} - 1_{\{i \leq \mathtt{first}(k)-1\}}) ]^2}{\bar{R}_i(we^{\eta})^2}

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 1{i<𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)1}1{i𝚏𝚒𝚛𝚜𝚝(k)1}=1{i<𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)1}. 1_{\{i < \mathtt{startmap}(k)-1\}} 1_{\{i \leq \mathtt{first}(k)-1\}} = 1_{\{i < \mathtt{startmap}(k)-1\}}.

Full Hessian

Suppose we want to multiply the Hessian by a vector (ξi)1in(\xi_i)_{1 \leq i \leq n}. Let’s use (r,c)(r,c) for the rr-th row and cc-th column of the Hessian.

Non-diagonal term

We’ll have a term like wrwceηr+ηci=1ndiwiRi(weη)2rierriec. -w_r w_c e^{ \eta_r + \eta_c} \sum_{i=1}^n \frac{d_i \bar{w}_i }{\bar{R}_i(we^{\eta})^2} \bar{r}_i'e_r \bar{r}_i'e_c.

We are free to rewrite this as: wrwceηr+ηci=1ndiwiRi(weη)2[1{trti}1{srti}δσi(1{i𝚕𝚊𝚜𝚝(r)}1{i𝚏𝚒𝚛𝚜𝚝(r)1})]riec. -w_r w_c e^{ \eta_r + \eta_c} \sum_{i=1}^n \frac{d_i \bar{w}_i }{\bar{R}_i(we^{\eta})^2} [1_{\{t_r \geq t_i\}} - 1_{\{s_r \geq t_i\}} - \delta \cdot \sigma_i \cdot (1_{\{i \leq \mathtt{last}(r)\}} - 1_{\{i \leq \mathtt{first}(r)-1\}})] \bar{r}_i'e_c. 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 ξ\xi yields (with eje_j the elementary basis vectors): wreηri=1ndiwiRi(weη)2[1{trti}1{srti}δσi(1{i𝚕𝚊𝚜𝚝(r)}1{i𝚏𝚒𝚛𝚜𝚝(r)1})](c=1nriecwceηcξc) -w_r e^{\eta_r} \sum_{i=1}^n \frac{d_i \bar{w}_i }{\bar{R}_i(we^{\eta})^2} [1_{\{t_r \geq t_i\}} - 1_{\{s_r \geq t_i\}} - \delta \cdot \sigma_i \cdot (1_{\{i \leq \mathtt{last}(r)\}} - 1_{\{i \leq \mathtt{first}(r)-1\}})] \left(\sum_{c=1}^n \bar{r}_i'e_c w_c e^{\eta_c} \xi_c\right)

The sum over cc is simply computing ri(ξweη)=Ri(ξwη) \bar{r}_i'(\xi w e^{\eta}) = \bar{R}_i(\xi w^{\eta}) Lemma 1 tells us how to do this correctly for the ii with di=1d_i=1. This is done by calling the function _sum_over_risk_set in the code, applied to an argument ξweη\xi w e^{\eta}.

We now have to compute the following sum over iiwreηri=1ndiwiRi(ξweη)Ri(weη)2[1{trti}1{srti}δσi(1{i𝚕𝚊𝚜𝚝(r)}1{i𝚏𝚒𝚛𝚜𝚝(r)1})]. - w_r e^{\eta_r} \sum_{i=1}^n \frac{d_i \bar{w}_i \bar{R}_i(\xi w e^{\eta})}{\bar{R}_i(w e^{\eta})^2} [1_{\{t_r \geq t_i\}} - 1_{\{s_r \geq t_i\}} - \delta \cdot \sigma_i \cdot (1_{\{i \leq \mathtt{last}(r)\}} - 1_{\{i \leq \mathtt{first}(r)-1\}})] .

Up to the scaling by wreηr-w_r e^{\eta_r} computation can be handled by Lemma 2 applied to the sequences wiRi(ξweη)Ri(ξweη)2,wiσiRi(ξweη)Ri(ξweη)2. \frac{\bar{w}_i \bar{R}_i(\xi w e^{\eta})}{\bar{R}_i(\xi w e^{\eta})^2}, \frac{\bar{w}_i \sigma_i \bar{R}_i(\xi w e^{\eta})}{\bar{R}_i(\xi w e^{\eta})^2}. This is handled by the function _sum_over_events in the code.

Diagonal term

The diagonal term will have (r,c)(r,c) entry δrcwreηri=1ndiwi[1{trti}1{srti}δσi(1{i𝚕𝚊𝚜𝚝(r)}1{i𝚏𝚒𝚛𝚜𝚝(r)1})]Ri(weη). \delta_{rc} w_r e^{\eta_r} \sum_{i=1}^n \frac{d_i \bar{w}_i [1_{\{t_r \geq t_i\}} - 1_{\{s_r \geq t_i\}} - \delta \cdot \sigma_i \cdot (1_{\{i \leq \mathtt{last}(r)\}} - 1_{\{i \leq \mathtt{first}(r)-1\}})]}{\bar{R}_i(we^{\eta})}.

Summing over cc in a multiplication by ξ\xi on the right yields c=1nξcδrcwreηri=1ndiwi[1{trti}1{srti}δσi(1{i𝚕𝚊𝚜𝚝(r)}1{i𝚏𝚒𝚛𝚜𝚝(r)1})]Ri(weη)=ξrwreηri=1ndiwi[1{trti}1{srti}δσi(1{i𝚕𝚊𝚜𝚝(r)}1{i𝚏𝚒𝚛𝚜𝚝(r)1})]Ri(weη). \begin{aligned} &\sum_{c=1}^n \xi_c \delta_{rc} w_r e^{\eta_r} \sum_{i=1}^n \frac{d_i \bar{w}_i [1_{\{t_r \geq t_i\}} - 1_{\{s_r \geq t_i\}} - \delta \cdot \sigma_i \cdot (1_{\{i \leq \mathtt{last}(r)\}} - 1_{\{i \leq \mathtt{first}(r)-1\}})]}{\bar{R}_i(we^{\eta})} \\ &= \xi_r w_r e^{\eta_r} \sum_{i=1}^n \frac{d_i \bar{w}_i [1_{\{t_r \geq t_i\}} - 1_{\{s_r \geq t_i\}} - \delta \cdot \sigma_i \cdot (1_{\{i \leq \mathtt{last}(r)\}} - 1_{\{i \leq \mathtt{first}(r)-1\}})]}{\bar{R}_i(we^{\eta})}. \end{aligned} We again apply Lemma 2 to the sequences wiRi(weη),wiσiRi(weη). \frac{\bar{w}_i}{\bar{R}_i(we^{\eta})}, \frac{\bar{w}_i \sigma_i}{\bar{R}_i(we^{\eta})}.

Probabilistic interpretation

The Hessian has (r,c)(r,c) entry i=1ndiwi(δrcπi,η[r]πi,η[r]πi,η[c]). \sum_{i=1}^n d_i \bar{w}_i \left( \delta_{rc} \cdot \pi_{i,\eta}[r] - \pi_{i,\eta}[r]\pi_{i,\eta}[c] \right).

Multiplying on the right by ξ\xi yields i=1ndiwi(ξrπi,η[r]πi,η[r]Ei,η[ξ])=i=1ndiwi(ξrEi,η[ξ])πi,η[r] \sum_{i=1}^n d_i \bar{w}_i \left( \xi_r \pi_{i,\eta}[r] - \pi_{i,\eta}[r] E_{i,\eta}[\xi] \right) = \sum_{i=1}^n d_i \bar{w}_i (\xi_r - E_{i,\eta}[\xi]) \pi_{i,\eta}[r]

Appendix: Proofs of Lemmas

Proof of Lemma 1 (Forward cumsum representation)

Consider the evaluation of the risk sum for an individual ii that experiences an event (di=1d_i = 1). The sum can be decomposed into two parts based on the definition of the risk set at time tit_i, which includes all individuals jj such that sj<titjs_j < t_i \leq t_j: j:sj<titjZj=j:tjtiZjm:smti,tmtiZm. \sum_{j: s_j < t_i \leq t_j} Z_j = \sum_{j: t_j \geq t_i} Z_j - \sum_{m: s_m \geq t_i, t_m \geq t_i} Z_m. Because the start time sms_m is always strictly less than the event time tmt_m (sm<tms_m < t_m), the condition smtis_m \geq t_i naturally implies tmtit_m \geq t_i. Thus, the second sum simplifies: j:sj<titjZj=j:tjtiZjm:smtiZm. \sum_{j: s_j < t_i \leq t_j} Z_j = \sum_{j: t_j \geq t_i} Z_j - \sum_{m: s_m \geq t_i} Z_m.

Let us assume 0-based indexing for the sequences. The sequence ZZ is sorted in ascending order of event times tt. The index 𝚏𝚒𝚛𝚜𝚝(i)\mathtt{first}(i) represents the first occurrence of the event time tit_i in this sorted sequence. Therefore, summing ZjZ_j for all jj where tjtit_j \geq t_i is equivalent to summing the sequence from index 𝚏𝚒𝚛𝚜𝚝(i)\mathtt{first}(i) to the end (n1n-1). By definition, the reversed cumsum E(Z)𝚏𝚒𝚛𝚜𝚝(i)E(Z)_{\mathtt{first}(i)} computes exactly this quantity.

Similarly, let S(Z)S(Z) be the reversed cumsum of ZZ sorted by start times ss. The quantity 𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)=#{j:sj<ti}\mathtt{eventmap}(i) = \#\{j : s_j < t_i\} gives the exact number of start times strictly less than tit_i. In a 0-based indexed array sorted by start times, the indices 0,1,,𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)10, 1, \dots, \mathtt{eventmap}(i)-1 correspond to sm<tis_m < t_i, and the indices from 𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)\mathtt{eventmap}(i) onwards correspond to smtis_m \geq t_i. Thus, the reversed cumsum S(Z)𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)S(Z)_{\mathtt{eventmap}(i)} computes the sum of ZmZ_m for smtis_m \geq t_i.

Multiplying by did_i (which trivially yields 00 on both sides if di=0d_i=0) establishes the first equality: di(j:tjtiZjm:smtiZm)=di(E(Z)𝚏𝚒𝚛𝚜𝚝(i)S(Z)𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)). d_i \cdot \left( \sum_{j: t_j \geq t_i} Z_j - \sum_{m: s_m \geq t_i} Z_m \right) = d_i \cdot \left( E(Z)_{\mathtt{first}(i)} - S(Z)_{\mathtt{eventmap}(i)} \right).

For the second part of the Lemma, recall the definition of Ri(ξ)\bar{R}_i(\xi): Ri(Z)=j:tjtiZjm:smtiZmδσiΔi(Z) \bar{R}_i(Z) = \sum_{j: t_j \geq t_i} Z_j - \sum_{m: s_m \geq t_i} Z_m - \delta \cdot \sigma_i \cdot \Delta_i(Z) where Δi(Z)=j:Cj=CiZj\Delta_i(Z) = \sum_{j: C_j=C_i} Z_j is the sum over the cluster of tied event times at tit_i. The cluster spans indices from 𝚏𝚒𝚛𝚜𝚝(i)\mathtt{first}(i) to 𝚕𝚊𝚜𝚝(i)\mathtt{last}(i). In terms of the reversed cumsum E(Z)E(Z), the sum over this finite slice is simply E(Z)𝚏𝚒𝚛𝚜𝚝(i)E(Z)𝚕𝚊𝚜𝚝(i)+1E(Z)_{\mathtt{first}(i)} - E(Z)_{\mathtt{last}(i)+1}. Substituting this difference into the expression for Ri(Z)\bar{R}_i(Z) yields: riZ=E(Z)𝚏𝚒𝚛𝚜𝚝(i)S(Z)𝚎𝚟𝚎𝚗𝚝𝚖𝚊𝚙(i)δσi(E(Z)𝚏𝚒𝚛𝚜𝚝(i)E(Z)𝚕𝚊𝚜𝚝(i)+1), \bar{r}_i'Z = E(Z)_{\mathtt{first}(i)} - S(Z)_{\mathtt{eventmap}(i)} - \delta \cdot \sigma_i \cdot \left( E(Z)_{\mathtt{first}(i)} - E(Z)_{\mathtt{last}(i)+1} \right), which concludes the proof. \square

Proof of Lemma 2 (Adjoint cumsum representation)

We want to evaluate the sum i=1ndiZi(1{tkti}1{skti})\sum_{i=1}^n d_i Z_i (1_{\{t_k \geq t_i\}} - 1_{\{s_k \geq t_i\}}). The indicator (1{tkti}1{skti})(1_{\{t_k \geq t_i\}} - 1_{\{s_k \geq t_i\}}) equals 11 if and only if sk<titks_k < t_i \leq t_k, and is 00 otherwise. Thus, the sum effectively accumulates diZid_i Z_i over the interval (sk,tk](s_k, t_k].

Assuming 0-based indexing for the arrays sorted by event time tt, let us identify the bounds of this interval in terms of the sequence indices. The upper bound is defined by titkt_i \leq t_k. In the sorted array, the elements satisfying this condition occupy indices 00 through 𝚕𝚊𝚜𝚝(k)\mathtt{last}(k). The total number of such elements is 𝚕𝚊𝚜𝚝(k)+1\mathtt{last}(k) + 1. The lower bound is defined by sk<tis_k < t_i, which means we exclude elements where tiskt_i \leq s_k. The number of elements with tiskt_i \leq s_k is precisely given by 𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)=#{j:tjsk}\mathtt{startmap}(k) = \#\{j : t_j \leq s_k\}. Thus, the excluded indices are 00 through 𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)1\mathtt{startmap}(k) - 1.

Consequently, the interval sk<titks_k < t_i \leq t_k corresponds exactly to the index range i{𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k),,𝚕𝚊𝚜𝚝(k)}i \in \{\mathtt{startmap}(k), \dots, \mathtt{last}(k)\}.

Let E(dZ)\bar{E}(d \cdot Z) denote the forward cumsum sequence of dZd \cdot Z, padded with a 00 at the 0-th index. That is, E(dZ)m=i=0m1diZi\bar{E}(d \cdot Z)_m = \sum_{i=0}^{m-1} d_i Z_i. Using this forward cumsum, the sum over our target index range can be computed as the difference between the sum up to 𝚕𝚊𝚜𝚝(k)\mathtt{last}(k) and the sum up to 𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)1\mathtt{startmap}(k)-1: i=𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)𝚕𝚊𝚜𝚝(k)diZi=E(dZ)𝚕𝚊𝚜𝚝(k)+1E(dZ)𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k). \sum_{i=\mathtt{startmap}(k)}^{\mathtt{last}(k)} d_i Z_i = \bar{E}(d \cdot Z)_{\mathtt{last}(k)+1} - \bar{E}(d \cdot Z)_{\mathtt{startmap}(k)}. This establishes the first equality.

The second equality, 𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k)=𝚏𝚒𝚛𝚜𝚝(𝚜𝚝𝚊𝚛𝚝𝚖𝚊𝚙(k))\mathtt{startmap}(k) = \mathtt{first}(\mathtt{startmap}(k)), holds because all individuals jj that share the exact same event time tj=skt_j = s_k must either entirely fall within or outside the threshold boundary. Hence, the index count maps cleanly to the 𝚏𝚒𝚛𝚜𝚝\mathtt{first} function of that start boundary, ensuring numerical stability in implementations with ties. \square