Recurrent events are the natural outcome in many medical and epidemiology studies. To assess covariate effects on the gaps between consecutive recurrent events, the Cox proportional hazards model is frequently employed in data analysis. The validity of statistical inference, however, depends on the appropriateness of the Cox model. In this paper, we propose a class of graphical techniques and formal tests for checking the Cox model with recurrent gap time data. The building block of our model checking method is an averaged martingale-like process, based on which a class of multiparameter stochastic processes is proposed. This maneuver is very general and can be used to assess different aspects of model fit. Numerical simulations are conducted to examine finite-sample performance, and the proposed model checking techniques are illustrated with data from the Danish Psychiatric Central Register.
Keywords: Correlated failure times, Induced-dependent censoring, Kaplan–Meier estimator, Renewal processes
. The observation of recurrent events is subject to right censoring. We denote by Ci the time between the initial event and the end of follow-up, and assume that Ci is independent of γi and (Ti1,Ti2,…) conditional on Zi. Let mi ≥ 1 denote the index satisfying 
and 
, with ∑10·![[equivalent]](/corehtml/pmc/pmcents/equiv.gif)
0. Thus, mi is the number of observed gap times, where the first mi − 1 gap times are completely observed and the last gap time is censored by 
. Denote the observed data for subject i by {Yi1,…,Yimi;mi,Zi}, where Yij = Tij for j = 1,…mi − 1 and 
. The observed data across the n subjects are assumed to be iid.
and τ satisfies pr(Yi1 ≥ τ) > 0. Note that we define a![[multiply sign in circle]](/corehtml/pmc/pmcents/otimes.gif)
0 = a,a
is defined as the solution to U*(b) = 0. As pointed out by 
,n} converge uniformly to the same limit, ensuring the consistency of 
is assumed to be positive definite at b = β and 
(b,t) is the limit of S(1)*(b,t)/S(0)*(b,t) as n→∞.
is given by the Breslow-type estimator 
with
behaves similarly to the ordinary martingale residual process in that, for any t,
and 
for large n.
. Under the assumed Cox proportional hazards model, the processes are expected to fluctuate randomly around zero. We show in Appendix A of the 
is asymptotically equivalent to 
, where
are the n-limits of
converges to a zero-mean Gaussian process with covariance function E{Ψi*(t1,z1)Ψi*(t2,z2)T} for any pairs (t1,z1) and (t2,z2). Replacing the unknown quantities in Ψi* by the corresponding sample estimates, we define the function
.
. Conditional on the data, and thus conditioned on 
is a sum of n independent, zero-mean random vectors for every (t,z). Arguing as in section 2.4 of 
converge weakly to the same limiting process. Hence, to approximate the null distribution of W*(t,z), we can simulate a number of realizations from 
by repeatedly generating standard normal random samples (G1,…,Gn). To assess how unusual the observed residual pattern is under the assumed model, we plot the observed multiparameter process along with a few, say 20, processes simulated from 
. An approximate P value for the supremum test 
can be estimated by the empirical probability that the observed value of 
exceeds the observed value of 
, and the P value can be obtained by using the supremum test 
.
and log(Z), and it is easy to see that the pattern is similar to 
where 
is the estimated variance–covariance matrix for 
given by
, where 
is the kth component of 
is the kth diagonal element of 
, where
![[set membership]](/corehtml/pmc/pmcents/x2208.gif)
[0,1]. Let Tij(0) = − log{1 − Φ(Ai + Bij)}, where Φ(·) is the cumulative density function of the standard normal distribution. It can be verified that Ti1(0),Ti2(0),… are realizations of correlated exponential random variables with unit means. Thus by applying a transformation to Tij(0), we can easily impose covariate effects and/or change the shape of the hazard function. As an example, 
has a marginal cumulative baseline hazard function Λ(t) = a×tκ.
is applied to the simulated data and is compared with the LWY test using the first gap times. 
against the set of unique values of 
global test can be constructed by considering 
for all possible values of z = (z1,…,zp) and t. Arguing as in Section 3.1, the distributions of Wl*(z) and Wg*(z,t) also can be approximated through simulating the corresponding zero-mean Gaussian process. Finally, a smoothed plot of 
versus a covariate omitted from the fitted model provides approximately the correct functional form of that covariate.
Model checking methods for the AFT model can be developed based on the corresponding averaged martingale-like process. Plots of the multiparameter stochastic processes with simulated trajectories under the AFT model will help assess model fit. If the AFT model does not fit well, we can repeat this exercise for the Cox model and a linear transformation model. At the end, we can select the model that has the best-fitting residual process.