## Nonproportional Hazards Model

Non-proportional hazards modeling has been widely studied in the past decades. In this section, the Bagdonavicius and Nikulin's [BN99] generalized proportional hazards model and its variant are introduced, which also can deal with nonconstancy as well as heterogeneity. Before the discussion, several nonproportional hazards models are reviewed.

First consider a parametric model like the Weibull regression taking the cumulative hazard function as (see, for example, Gore, Pocock, and Kerr [GPK84, page 185] for a survivor-function expression):

A{t) = tbeßo+ßizi+---+ßp Zp = tbeßT z, which can be viewed as a Weibull-class with universal shape parameter b, but with different scales indexed by ßTz. By this, the Hsieh model of Section 3 have a parametric (Weibull) regression model as a special case:

X RT rp

A(t) = te eß z, for Y x = yo + Yixi + ... + Yqx q-

There are also several nonproportional hazards model studied in Gore at al. [GPK84] Among them, a Cox-type model with 'varying-proportionality hazards' is of interest:

where 0(z,t) is a smooth function of covariate z and time variable t. Model (19) is capable of modeling heterogeneity plus nonconstancy. Because the time dependence is described simultaneously by Ao(t) and 0( • ,t), the form of 0(z,t) should be 'pre-specified' to make the setting identifiable. For modeling only the heterogeneity over the covariate space, a very flexible partly linear setting can be imposed on the relative risk function: A(t; x, z) = A0(t)exp{3Tx + g(z)}; see Sasieni [SAS92a], Nielsen, Linton, and Bickel [NLB98], and Heller [HEL01]. When z is categorical, the partly linear model reduces to the stratified proportional hazards model, by which the heterogeneity effect over z is further stratified out by putting a number of unknown baseline hazards. Moreover, a continuously stratified Cox model introduced in Sasieni [SAS92b]) and the 'time-varying' coefficients Cox model studied by Murphy and Sen [MS91], Murphy [MUR93] and Martinussen, Scheike, and Skovgaard [MSS01] (among others) are important works on the Cox-type relative risk modeling for heterogeneity and/or nonconstancy.

An alternative Cox-type model (other than the Hsieh model) which can deal with heterogeneity and nonconstancy together is the generalized PH model proposed by Bagdonavicius and Nikulin [BN99], see also Bagdonavicius and Nikulin [BN02] for more related works and models. The generalized proportional hazards model has the following form, in terms of hazard function,

where g(-) is a positive function, z(t) is a set of time-dependent covariate, and Az is the corresponding cumulative hazard. A special form of (20) derived in Bagdonavicius, Hafdi, and Nikulin [BHN04] is

where Ao(t) = f0 Ao(u)du is the baseline cumulative hazard function. Model (21) is called hereafter the Bagdonavicius-Nikulin model. Similar to the Hsieh model, the Bagdonavicius-Nikulin model also gives cross-effect for the cumulative hazards, but not necessarily for the hazard functions. When the logrelative risk is the effect of concern, the main differences between the Hsieh and the Bagdonavicius-Nikulin models are: (i) the former assumes the relative risk between groups to be possibly very large when t approaches 0, the model design of the latter relaxes this assumption; and (ii) relative risk of the Hsieh model is increasing or decreasing according to the relative direction of the heteroscedasticity part yTx, the Bagdonavicius-Nikulin model has more complex situation which depends on the configurations of ¡3 and 7 in (21). Figure 2 gives a similar illustration to Figure 1 and corresponds to part of the configurations listed in Bagdonavicius et al. [BNLZ05]. When 7 is large (0.5 or -0.5), the dependence of log-relative risks (LRR) on time is more significant. Also for larger 7, difference of LRRs in different z's is apparent, showing heterogeneity effect over the covariate-space. In Figure 2, the LRR also has an order in z (for z = -5, ..., 5), so only z = 5 is plotted by a solid line; for all cases, ¡3 =1 and 0 < t < 2. The estimation of the Bagdonavicius-Nikulin model depends on iteratively solving the score equations derived from the 'partial likelihood' and a Breslow-type estimating equation for the baseline Fig. 2(a): Bagdonavicius-Nikulin model, gamma=0.1 Fig. 2(b): Bagdonavicius-Nikulin model, gamma=0.5

Fig. 2(a): Bagdonavicius-Nikulin model, gamma=0.1

Fig. 2(b): Bagdonavicius-Nikulin model, gamma=0.5

SiaiaiïiiiïiHiiiïiiiiiiiaiBiBIBIliaiBIBIBiî

Fig. 2(c): Bagdonavicius-Nikulin model, gamma=-0.1

Fig. 2(d): Bagdonavicius-Nikulin model, gamma=-0.5

Fig. 2. Heterogeneous log-relative risks in the Bagdonavicius-Nikulin model.

cumulative hazard. With contrast to the score-type tests derived from the Hsieh model, the Bagdonavicius-Nikulin model can also be used as the alternative hypothesis and thus testing for homogeneity of survival distributions can be implemented (Bagdonavicius et al. [BNLZ05]).

### 5 Extensions and Brief Discussion

There are numerous papers targeting at heterogeneity and associated statistical inference. This paper is not intended to review the entire development, but only to focus on heterogeneity effect that can be modeled through collected samples.

The early works of Hsieh [HSI95, HSI96a, HSI96b, HSI96c] make statistical inferences on location-shift and scale-change parameters through EPA for two-sample models with complete or right-censored data. In Hsieh [HSI01], however, hazards regression is discussed and estimation and testing problems are solved with slightly different manner, though with similar spirit to the EPA method. The Hsieh model and the later introduced Bagdonavicius-Nikulin model allow for cross-effect modeling. Cross-effect or nonproportional hazards problem still is an important issue for future researches. In actual data analysis of follow-up data, multiple cross-effect may exist and extensions to these two models are of interests. First, consider the Hsieh model equipped with varying coefficients (in terms of hazard function):

X(t; z, x) = Xo(t){Ao(t)}eY(t)TX-1e/(t)Tz+Y(t)Tx, (22)

where ¡(t) and y(t) are two sets of varying coefficients. Model (22) can be viewed as dealing with the 'time-heteroscedasticity interaction'. In practice, the same set of time-partition for a piecewise-constant approximation can be applied to Xo(t), ¡3(t) and 7(t), and estimating equations similar to (15) and (16) while taking account the approximation can be used (Wu and Hsieh [WH04]). By suitably choosing the cutoff points which constitute the piecewise-constant intervals, multiple crossings among cumulative hazards according to different groups can be modeled. Second, a variant of model (21) is recently proposed by Bagdonavicius and Nikulin [BN04] by which at least two crossings can be properly captured:

X(t; z) = Xo(t)efT z(t){1 + yT z(t)Ao(t) + ST z(t)A20(t)}, (23)

where 7 and S are two sets of parameters to be estimated.

Diagnostics for a hazards regression model is important for model validity and goodness-of-fit problem. By simply plotting the estimated (log-) relative risks versus time or the important covariates, time-constancy and effect homogeneity/heterogeneity can be checked. For example, Valsecchi, Silvestri, and Sasieni [VSS96] used the plot to check nonconstancy as well as proportionality for various explanatory variables of ovarian cancer patient's survival. Further applications of plotting the (estimated) relative risk as a model discrimination and diagnostics tool for the Hsieh and the Bagdonavicius-Nikulin models is interesting for future researches.

### ACKNOWLEDGMENT

This research is partly supported by Grant NSC93-2118-M-039-001 of Taiwan's National Science Council. The relevant works of the author were mostly inspired by Professors Fushing Hsieh and Mikhail Nikulin; and also by Professor Marvin Zelen in a conversation during his visit to Taiwan.