Robust Inference For WLW Models

Why WLW Inference Must Respect Subject Clustering

The WLW model represents event orders with separate marginal event-time outcomes. A subject may contribute information to more than one event order, so the resulting records are not independent. In particular, at a fixed follow-up horizon, having experienced a second event implies having experienced a first. More generally, the timing and occurrence of a subject’s later events are connected to that subject’s earlier event history. Treating each event-order record as an independent observation can therefore give an incorrect variance for the estimated coefficients and their contrasts.

The earlier material on the WLW marginal model and WLW risk sets defines the event-order-specific marginal hazards and their risk sets. Here the focus is inference: the estimating equations can be combined across event orders, but their uncertainty must be calculated using the subject—not the event-order record—as the independent sampling unit.

Definition: A subject-clustered sandwich variance uses the sum of a subject’s score contributions across all event orders as that subject’s contribution to the empirical variability of the estimating equations. It allows arbitrary dependence among a subject’s event-order outcomes, while requiring independence between subjects for the usual large-sample justification.

Event-Order Scores and the Combined Equation

Let \(i=1,\ldots,n\) index independent subjects and \(k=1,\ldots,K\) index event order. For a WLW marginal Cox model, write the event-order-specific hazard as

\lambda_{ik}(t\mid Z_i)=\lambda_{0k}(t)\exp(\beta_k^T Z_i),

where \(Z_i\) is the covariate vector, \(\lambda_{0k}(t)\) is the event-order-specific baseline hazard, and \(\beta_k\) is the event-order-specific coefficient vector. The separate marginal models do not assert that the event-order outcomes within a subject are independent.

For event order \(k\), let \(N_{ik}(t)\) count observed events of that order by time \(t\), and let \(Y_{ik}(t)\) indicate that subject \(i\) is in that event order’s WLW risk set just before \(t\). Define the risk-set sums

S_k^{(r)}(\beta_k,t)=\sum_{j=1}^n Y_{jk}(t)\exp(\beta_k^T Z_j)Z_j^{\otimes r}, \qquad r=0,1,

where \(Z_j^{\otimes 0}=1\) and \(Z_j^{\otimes 1}=Z_j\). The event-order score is

U_k(\beta_k)=\sum_{i=1}^n \int_0^\tau \left\{Z_i-\frac{S_k^{(1)}(\beta_k,t)}{S_k^{(0)}(\beta_k,t)}\right\}\,dN_{ik}(t).

The upper limit \(\tau\) denotes the follow-up range used for the analysis. Stack the event-order coefficients and scores as \(\beta=(\beta_1^T,\ldots,\beta_K^T)^T\) and \(U(\beta)=(U_1(\beta_1)^T,\ldots,U_K(\beta_K)^T)^T\). The fitted estimate \(\hat\beta\) solves the combined equations \(U(\hat\beta)=0\). If a model constrains some coefficients to be shared across orders, the equations are instead written in terms of that shared parameter vector; the cluster-based variance principle is unchanged.

For each subject, form the stacked score contribution \(U_i(\beta)\), whose \(k\)th block is that subject’s contribution to \(U_k(\beta_k)\). In the displayed Cox score, that block is the integral involving \(dN_{ik}(t)\). The essential step is to keep all of a subject’s blocks together when estimating variability.

The Subject-Clustered Sandwich

Let \(A\) be the sensitivity matrix, the negative derivative of the combined estimating equations, evaluated at the fitted coefficients. Let \(B\) be the empirical covariance “meat” formed from subject-level scores:

A=-\left.\frac{\partial U(\beta)}{\partial\beta^T}\right|_{\hat\beta}, \qquad B=\sum_{i=1}^n U_i(\hat\beta)U_i(\hat\beta)^T.

The large-sample subject-clustered covariance estimate is

\widehat{\operatorname{Var}}(\hat\beta)=A^{-1}BA^{-T}.

The transpose on the final inverse matters when the sensitivity matrix is not symmetric. For ordinary Cox partial-likelihood equations with separate event-order coefficient blocks, the sensitivity is often block diagonal, because each event-order score depends on its own coefficient block. The meat need not be block diagonal: a subject’s scores for different event orders can covary.

Interpretation: The diagonal entries of the sandwich covariance estimate the variances of the fitted coefficients. Its off-diagonal entries estimate covariances between coefficients. Standard errors are the square roots of the corresponding diagonal entries. A robust variance changes estimated uncertainty; it does not change the fitted coefficients or the marginal model’s estimand.

Why the Sandwich Formula Arises

The justification is a first-order expansion of the estimating equations. Subjects are the independent sampling units, but the displayed Cox score contributions share risk-set quantities calculated from all subjects and are not literally independent in finite samples. The sandwich justification relies on the relevant Cox counting-process and asymptotic theory under suitable regularity conditions, while allowing dependence among event-order components within a subject. Around the true parameter \(\beta_0\), the fitted root approximately satisfies

0=U(\hat\beta)\ \approx\ U(\beta_0)-A(\hat\beta-\beta_0), \qquad \hat\beta-\beta_0\ \approx\ A^{-1}\sum_{i=1}^n U_i(\beta_0).

Consequently, the variability of \(\hat\beta\) depends on the covariance of the whole subject score vector, not on separate event-order records considered as independent. Replacing that covariance by its empirical sum of outer products gives \(B\), and substituting the fitted sensitivity gives the sandwich estimator. The approximation is justified under regularity conditions, including independent subjects, a correctly specified estimating target with a consistent root, a nonsingular sensitivity matrix, and enough independent subjects for the large-sample approximation to be useful.

The calculation also relies on suitable observation and censoring assumptions for the marginal event-order models. In particular, censoring must be noninformative for the modeled event-order outcomes conditional on the covariates included in the analysis, or otherwise handled by a valid method. A sandwich variance does not repair an incorrect risk set, informative censoring, or a misspecified marginal hazard structure.

Fit the combined model.
Estimate the event-order coefficients using the specified WLW risk sets and estimating equations.

Retain subject identities.
For every subject, collect the score contributions from all event orders represented in the fit.

Build the meat.
Sum the outer products of the subject-level stacked scores, \(U_iU_i^T\), not the outer products of individual event-order records.

Combine with sensitivity.
Calculate \(A^{-1}BA^{-T}\), then obtain standard errors and contrast variances from this matrix.

Worked Example: Nested Event-Order Outcomes at a Fixed Horizon

Suppose a hypothetical cohort has 100 subjects, with outcomes recorded by a fixed follow-up horizon. There are 30 subjects with both a first and second event, 30 with a first event but no second event, and 40 with neither event. These counts respect the event-order relationship: no subject has a second event without a first.

Let \(Y_1\) and \(Y_2\) indicate whether a subject has had the first and second event by the horizon. Then \(P(Y_1=1)\) is estimated by \(60/100=0.60\), \(P(Y_2=1)\) by \(30/100=0.30\), and \(P(Y_1=1,Y_2=1)\) by \(30/100=0.30\).

The empirical covariance is

\widehat{\operatorname{Cov}}(Y_1,Y_2) =0.30-(0.60)(0.30)=0.12.

The empirical variances are \(0.60(0.40)=0.24\) and \(0.30(0.70)=0.21\). Thus the correlation is approximately

\frac{0.12}{\sqrt{(0.24)(0.21)}}=0.5345.

This simple horizon-based calculation is not a WLW coefficient analysis. It illustrates why independence across event orders is implausible: the two indicators are positively associated, and the second-event indicator is nested within the first. In a time-to-event analysis, the cluster-score meat carries the corresponding within-subject dependence into coefficient inference.

Calculating a Clustered Covariance Matrix

Worked Example: Sandwich Standard Errors for Two Event-Order Coefficients

Consider a hypothetical fitted WLW analysis with two event orders and a separate coefficient for each order. Suppose its evaluated sensitivity matrix is \(A=\begin{pmatrix}4&0\\0&5\end{pmatrix}\). For illustration, assume the four subject-level stacked scores at the fitted estimate are as follows. Each row is the pair of that subject’s first- and second-order score contributions, and the rows sum to zero in each coordinate as expected at the fitted estimating-equation root.

Subject	First-order score	Second-order score
1	1	0
2	-1	1
3	0	-1
4	0	0

The subject-level meat is the sum of the four outer products:

B= \begin{pmatrix}1&0\\0&0\end{pmatrix} +\begin{pmatrix}1&-1\\-1&1\end{pmatrix} +\begin{pmatrix}0&0\\0&1\end{pmatrix} +\begin{pmatrix}0&0\\0&0\end{pmatrix} = \begin{pmatrix}2&-1\\-1&2\end{pmatrix}.

Since \(A^{-1}=\begin{pmatrix}1/4&0\\0&1/5\end{pmatrix}\), the sandwich covariance is

\widehat{\operatorname{Var}}(\hat\beta) = \begin{pmatrix}1/4&0\\0&1/5\end{pmatrix} \begin{pmatrix}2&-1\\-1&2\end{pmatrix} \begin{pmatrix}1/4&0\\0&1/5\end{pmatrix} = \begin{pmatrix}0.125&-0.050\\-0.050&0.080\end{pmatrix}.

The robust standard errors are \(\sqrt{0.125}=0.3536\) for the first-order coefficient and \(\sqrt{0.080}=0.2828\) for the second-order coefficient. The estimated covariance is negative, so the joint uncertainty is not fully described by the two standard errors alone. This four-subject calculation is deliberately small to show the arithmetic; it is not a setting in which large-sample sandwich inference should be trusted.

Worked Example: A Contrast Uses the Cross-Order Covariance

Using the covariance matrix from the previous calculation, consider the contrast \(\beta_1-\beta_2\). For any coefficient vector, the variance of this contrast is the variance of the first coefficient plus the variance of the second, minus twice their covariance:

\widehat{\operatorname{Var}}(\hat\beta_1-\hat\beta_2) =0.125+0.080-2(-0.050)=0.305, \qquad SE(\hat\beta_1-\hat\beta_2)=\sqrt{0.305}=0.5523.

If the fitted contrast were hypothetically \(0.60\), its Wald statistic for a zero contrast would be \(0.60/0.5523=1.086\). Using a standard normal reference distribution gives a two-sided \(p\)-value of approximately \(0.2773\). This is not strong evidence against a zero contrast at conventional significance levels. The result is an illustration of the calculation only: with four independent subjects, the asymptotic normal approximation is especially doubtful.

If the cross-order covariance had incorrectly been set to zero, the estimated standard error would instead be \(\sqrt{0.125+0.080}=0.4528\). That changes the contrast’s uncertainty, and hence its test statistic and \(p\)-value. It is not valid to discard the covariance simply because the coefficients refer to different event orders.

Assumptions, Diagnostics, and Alternatives

Robust inference is not assumption-free. Its main dependence advantage is that it does not require independent event-order outcomes within a subject. The usual sandwich justification still requires independent sampling clusters, a sufficiently large number of clusters, valid estimating equations, and a well-behaved sensitivity matrix. If subjects are additionally grouped in independent higher-level units—such as households or clinics—the clustering unit must reflect the actual independent sampling structure.

Useful checks include confirming that the software’s cluster identifier is the subject identifier; verifying that each subject’s score contributions are summed across event orders before forming \(B\); checking that the fitted equations converge and that \(A\) is nonsingular; and reviewing whether censoring, risk-set construction, or missing event-order information could be informative. Very large standard errors, unstable estimates, or a nearly singular covariance matrix can indicate sparse information, separation-like behavior, or excessive model complexity.

The model-based covariance usually relies on the assumed event-time model and its within-model variance calculation. It can be more efficient when that model is correct, but can be too small or otherwise misleading if dependence across event orders is ignored. The subject-clustered sandwich allows within-subject score dependence, but it does not fix bias in the coefficient estimates caused by misspecified marginal hazards or invalid censoring assumptions. Bootstrap resampling of whole subjects is another way to reflect within-subject dependence, but it also requires enough independent subjects and appropriate resampling of the study design. Resampling event-order records separately would break the within-subject structure and is not an appropriate alternative.

A further limitation is that the sandwich estimates large-sample uncertainty. With few independent subjects, the empirical meat may be noisy and normal or chi-square reference approximations may perform poorly; small-sample corrections or carefully designed subject-level resampling may be considered, with their assumptions stated. Robust standard errors also do not replace the need to check proportional hazards, functional form, or the scientific meaning of each event-order coefficient.

Common Mistakes

Counting records as independent clusters. The event-order rows for a subject are dependent. Sum the subject’s score contributions first, then form their outer product.
Assuming “robust” means immune to misspecification. The sandwich addresses variance under within-subject dependence; it does not automatically correct an incorrect estimand, hazard structure, censoring assumption, or risk set.
Dropping off-diagonal covariance terms. Cross-order covariance affects inference on contrasts and joint hypotheses, even when separate event-order coefficients are reported.
Interpreting a robust standard error as a changed effect estimate. The sandwich alters uncertainty calculations, not the fitted WLW coefficients.
Using a small number of clusters as if it were a large sample. A robust formula alone does not ensure accurate finite-sample inference.

Key takeaway: In WLW inference, form one stacked score vector per subject, with contributions from every event order that subject informs. The sandwich covariance \(A^{-1}BA^{-T}\) preserves the dependence between those contributions and is the basis for robust standard errors and cross-order contrasts.

Check Your Understanding

Answer each question, then compare with the solutions below.

At a fixed horizon, 25 subjects have both a first and second event, 35 have only a first event, and 40 have neither. How many have a second event but no first event?
For subject-level scores \(U_1=(1,2)^T\) and \(U_2=(-1,-2)^T\), what is their contribution to the meat \(B\)?
Suppose \(A=\begin{pmatrix}2&0\\0&4\end{pmatrix}\) and \(B=\begin{pmatrix}8&0\\0&16\end{pmatrix}\). Find the sandwich covariance and the two standard errors.
Why can a robust sandwich estimate remain invalid even when the software correctly clusters by subject?

Solutions. (1) Zero: a second event by the horizon implies a first event, so the counts are consistent and no subject can have only the second. (2) The contribution is \(U_1U_1^T+U_2U_2^T=\begin{pmatrix}1&2\\2&4\end{pmatrix}+\begin{pmatrix}1&2\\2&4\end{pmatrix}=\begin{pmatrix}2&4\\4&8\end{pmatrix}\). (3) \(A^{-1}=\operatorname{diag}(1/2,1/4)\), so \(A^{-1}BA^{-T}=\operatorname{diag}(2,1)\). The standard errors are \(\sqrt{2}=1.414\) and \(1.000\). (4) Clustering handles within-subject score dependence, but the estimate can still be unreliable with too few independent subjects or invalid estimating equations, informative censoring, an incorrect risk-set construction, or a misspecified marginal model.