The most known semi-parametric model for analysis of failure time regression data is the proportional hazards (PH) model. There are many tests for the PH model from right censored data given by Cox [5], Moreau et al. [11], Lin [10], Nagelkerke et al. [12], Grambsch and Therneau [6], Quantin et al. [13], Bagdonavičius et al. [4]. Tests for specified covariates are given in Kvaloy and Neef [9], Kraus [8].

We consider tests for proportional hazards assumption concerning specified covariates or groups of covariates. The class of alternatives is wide: log-hazard rates under different values of covariates may cross, approach, go away. The data may be right censored. The limit distribution of the test statistics is derived. Power of the test under approaching alternatives is given. Real examples are considered.

Let S(t|z) , λ(t|z) and Λ(t|z) be the survival, hazard rate and cumulative hazard functions under a m-dimensional explanatory variable (or covariate) z=(z1,…,zm)T .

Let us consider the proportional hazards (PH) model: λ(t|z)=eβTzλ(t), where λ(t) is unknown baseline hazard function, the parameter β=(β1,…,βm)T is m-dimensional.

Under the PH model the ratios of hazard functions under any two different explanatory variables are constant over time, i.e. they are proportional.

Suppose that proportionality of hazard functions with respect to a specified covariate zj or a group of covariates zj1,…,zjk may be violated. Our purpose is to detect such violations. So we seek a wider model which includes not only possibility of constant hazard rates ratios but also time-varying ratios.

We propose a model of the form (1) λ(t|z)=g(z,Λ(t),β,γj)λ(t),Λ(t)=∫0tλ(u)du, where g(z,Λ(t),β,γj)=eβTz+Λ(t)eγjzj1+eγjzj[eΛ(t)eγjzj−1]. β=(β1,…,βm)T and γj are unknown regression parameters, λ(t) and Λ(t) are unknown baseline hazard and cumulative hazard, respectively.

The PH model is a particular case of this model when γj=0 . The model (1) is very wide as compared to the PH model.

Indeed, suppose that two different values z(1) and z(2) of the covariate vector differ only by the value of the covariate zj and denote by c(t) the ratio of hazard rates. For the model (1) c(0)=eβj(zj(2)−zj(1)),c(∞)=e(βj−γj)(zj(2)−zj(1)). So this model gives a large choice of alternatives: the ratio of hazard rates c(t) can vary from any a>0 }0$]]> to any b>0 }0$]]>, whereas c(t) is constant under the PH model. The difference of logarithms of hazard rates under different values of the covariate zj is constant under the PH model, so the logarithms of hazard rates as time functions are parallel, whereas in dependence on the values of its parameters the model (1) includes also the possibilities for these functions to approach, to go away, to intersect. So the model (1) may help to detect non-proportionality of hazard rates (or, equivalently, non-parallelism of log-hazard functions) in above mentioned directions. Other directions are very rare in real data.

We do not discuss application of the model for analysis of survival regression data (which is a subject of another article). Such analysis could be done if the PH model would be rejected. Here the model is used only as a generator of a wide class of alternatives to the PH model.

Let us consider right censored failure time regression data: (X1,δ1,z(1)),…,(Xn,δn,z(n)), where Xi=Ti∧Ci,δi=1{Ti≤Ci}, Ti are failure times and Ci are censoring times. Xi is observation time of the ith unit, the event δi=1 indicates that Xi is the failure time Ti , and the event δi=0 indicates that Xi is the censoring time Ci .

Set Ni(t)=1{Xi≤t,δi=1},Yi(t)=1{Xi≥t},N(t)= ∑i=1nNi(t),Y(t)= ∑i=1nYi(t); here 1A denotes the indicator of the event A. The processes Ni(t) and N(t) are right-continuous counting process showing the numbers of observed failures in the time interval [0,t] for the ith unit and for all n units, respectively. Yi(t) and Y(t) are decreasing (not strongly) left-continuous stochastic processes showing the numbers of units which are still not-failed and not-censored just prior to t for the ith and for all units, respectively.

Suppose that the survival distribution of the i object given z(i) is absolutely continuous with the survival functions Si(t) and the hazard rates λi(t) .

Suppose that the multiplicative intensities model is verified: the compensators of the counting processes Ni with respect to the history of the observed processes are ∫Yiλidu . It is equivalent to the assumption that for any i and t: P(Xi>t)>0 }t)\mathrm{>}0$]]> with almost all s∈[0,t] limh↓01hP{Ti∈[s,s+h)|Xi≥s,x(i)}=limh↓01hP{Ti∈[s,s+h)|Ti≥s,x(i)}=λi(s). It means that for almost all s∈[0,t] the risk to fail just after the time s given that units were non-censored and did not fail to time s is equal to the risk to fail just after the time s given that units did not fail to time s when censoring does not exist. So censoring has no influence to the risk of failure.

Information about time-to-failure distribution contains the points, where the counting processes Ni have jumps. A jump at the point t is possible if Yi(t)=1 . Partial information give the points where the counting processes Ni have not jumps but the predictable processes have jumps.

The multiplicative intensities model is verified in the case of type I, type II, independent random censoring.

In the parametric case with known λ the unknown finite-dimensional parameter consists of the parameters β and γj . The component of the parametric score function U corresponding to γj has the form (2) Uj(β,γj)= ∑i=1n∫0∞{wj(i)(u,Λ,β,γj)−E˜j(u,Λ,β,γj)}dNi(u), where wj(i)(v,Λ,β,γj)=∂∂γjlog{g(z(i),Λ(v),β,γj)},E˜(v,Λ,β,γj)=S˜(1)(v,Λ,β,γj)S˜(0)(v,Λ,β,γj),S˜(0)(v,Λ,β,γj)= ∑i=1nYi(v)g{z(i),Λ(v),β,γj},S˜(1)(v,Λ,β,γj)= ∑i=1nYi(v)∂∂γjg{z(i),Λ(v),β,γj}. We consider semiparametric case when the baseline hazard rate λ is unknown. The test statistic is constructed in the following way. In the expression of Uj the parameter β is replaced by its partial likelihood estimator βˆ under the Cox model, the parameter γj is replaced by 0, and the baseline cumulative intensity Λ is replaced by the Breslow estimator (see Andersen et al. [2]) Λˆ(t)=∫0tdN(u)S(0)(u,βˆ). The estimator βˆ verifies the equation (3) ∑i=1n∫0∞{z(i)−E(t,βˆ)}dNi(t)=0, where E(t,β)=S(1)(t,β)S(0)(t,β),S(0)(t,β)= ∑i=1nYi(t)eβTz(i),S(1)(t,β)= ∑i=1nz(i)Yi(t)eβTz(i). Set Fˆ(t)=1−e−Λˆ(t) . It is an estimator of the baseline cumulative distribution function. After the above mentioned replacement the following statistic is obtained: (4) Uˆγj=− ∑i=1n∫0∞Fˆ(t){zj(i)−Ej(t,βˆ)}dNi(t), where Ej is the component of E corresponding to the covariate zj .

Let us consider asymptotic distribution of the statistic (4) under the PH model. Set S(2)(u,β)= ∑i=1n(z(i))⊗2Yi(u)eβTz(i),V(u,β)=S(2)(u,β)/S(0)(u,β)−E⊗2(u,β), where A⊗2=AAT for any matrix A. Denote by Vj(u,β) the jth column and by Vjj(u,β) the jth diagonal element of the matrix V(u,β) . Set Σˆjj(t)=n−1∫0tFˆ(u)Vjj(u,βˆ)dN(u),Σˆj(t)=n−1∫0tFˆ(u)Vj(u,βˆ)dN(u),Σˆ(t)=n−1∫0tV(u,βˆ)dN(u).

Assumption a) can be weakened considerably but we avoid writing complicated formulas for easier reading. Assumption b) simply means that at some finite time moment (perhaps very remote) observations are stopped. This is a usual assumption for asymptotic results to hold in survival analysis. Assumption c) also can be weakened. Assumption d) means that if censoring would be absent then units might survive after the moment τ with positive probability. Assumption e) is needed to have non-degenerated asymptotic distribution of the test statistic. This assumption is the usual assumption needed for asymptotic normality of the maximum partial likelihood estimator βˆ of the regression parameter β.

Under Assumptions A [3] there exists a neighborhood Θ of β such that (5) supb∈Θ,u∈[0,τ]‖1nS(i)(u,b)−s(i)(u,b)‖→P0,supb∈Θ,u∈[0,τ]‖E(u,b)−e(u,b)‖→P0,supb∈Θ,u∈[0,τ]‖V(u,b)−v(u,b)‖→P0, as n→∞ .

Assumption A a) may be weakened assuming that there exist non-random s(i) such that the convergences (5) hold.

Convergences (5) and Assumption A d) imply that supb∈Θ,t∈[0,τ]‖∫0t(1nS(i)(u,b)−s(i)(u,b))dΛ(u)‖=oP(1),supb∈Θ,t∈[0,τ]‖∫0t(E(u,b)−e(u,b))dΛ(u)‖=oP(1),supb∈Θ,t∈[0,τ]‖∫0t(V(u,b)−v(u,b))dΛ(u)‖=oP(1).

The test: the null hypothesis H0:γj=0 is rejected with approximate significance level α if |T|>zα/2 }{z_{\alpha /2}}$]]>; here zα/2 is the α/2 critical value of the standard normal distribution.

Let us consider the alternative model (1) and suppose that alternatives are approaching: γj=cn , where c≠0 is a fixed constant.

So the model (8) λ(t|z)=exp{βTz+Λ(t)exp{cnzj}}1+exp{cnzj}[exp{Λ(t)exp{cnzj}}−1]λ(t),Λ(t)=∫0tλ(u)du, is considered. Denote by β0 the true value of β under the model (8). The counting process Ni(t) has the form Ni(t)=∫0tg{z(i),Λ(u),β0,c/n}dΛ(u)+Mi(t), where Mi are local martingales with respect to the history generated by the data. Note that g{z(i),Λ(u),β0,c/n}=eβTz(i)(1−zj(i)F(u)cn+oP(1n)).

Let us consider the stochastic process Qn(t,β)=1n[ℓ(t,β)−ℓ(t,β0)], where ℓ(t,β)= ∑i=1n∫0t{βTz(i)−lnS(0)(u,β)}dNi(u).

Note that the derivative with respect to β is Q˙n(t,β)=1nℓ˙(t,β)=1n ∑i=1n∫0t{z(i)−E(u,β)}dNi(u). Constructing the test statistic we defined βˆ as a random vector satisfying the estimating equation ℓ˙(τ,β)=0 . Let us show that βˆ is a consistent estimator of β0 .