In general, the hazard ratio may be a function of time, and estimation of the hazard ratio function may provide useful insights into temporal aspects of treatment effects. For example, Gilbert and others (2002) develop a nonparametric estimation procedure for the log-hazard ratio function with simultaneous confidence bands, for use as an exploratory data analytic tool. Naturally, confidence bands may be wide with such a nonparametric estimator, particularly at longer follow-up times where data may be sparse. See also Gray (1992), Kooperberg and others (1995), Cai and Sun (2003), Tian and others (2005), Abrahamowicz and Mackenzie (2007), and Peng and Huang (2007), and references therein, for additional related work.

Parametric or semiparametric hazard ratio models have potential to contribute valuably to treatment effect assessment. Hazard ratio models having parameters of useful interpretation, and that embrace a range of hazard ratio shapes, may be particularly valuable. The Cox model allows time-varying covariates to be defined that can, for example, allow separate hazard ratios for the elements of a partition of the time axis or allow the hazard ratio to be a parametric function of follow-up time more generally. Various other semiparametric regression models have been proposed for failure time data analyses, including accelerated failure time models, proportional odds models, and linear transformation models, many of which are embraced by the broad class of models for which Zeng and Lin (2007) develop maximum likelihood estimation procedures. Some more semiparametric models can be found in Vaupel and others (1979), Hsieh (1996), Chen and Wang (2000), Tsodikov (2002), Yang and Prentice (2005), and Chen and Cheng (2006). Many of these models induce a semiparametric class of models for the hazard ratio function that includes proportional hazards as a special case. Hazard ratio estimators under such semiparametric models can avoid the instability that may attend nonparametric hazard ratio function estimators.

One of these, proposed by Yang and Prentice (2005), involves short-term and long-term hazard ratio parameters, and a hazard ratio function that depends also on the control group survivor function. Assume absolutely continuous failure times and label the 2 groups control and treatment, with hazard functions λ_C(t) and λ_T(t), respectively. Let h(t) = λ_T(t)/λ_C(t) be the hazard ratio function and S_C(t) the survivor function of the control group. The model postulates thatwhere β₁ and β₂ are scalar parameters and

This model includes the proportional hazards model and the proportional odds model as special cases. It has a monotone h(t) with a variety of patterns, including proportional hazards, no initial effect, disappearing effect, and crossing hazards, among others. Thus, the model presumably entails sufficient flexibility for many applications. It has also been studied for current status data in Tong and others (2007).

In comparison, for many commonly used special cases of the accelerated failure time model either lim_t↓0h(t) = 1 or lim_t↑τ0h(t){0,1,∞} and the hazard ratio stays above or below one when λ_C is increasing. This is less flexible than desired. For the class of linear transformation models, with the logarithmic transformation, the hazard ratio also inherits some of these restrictions at many common baseline distributions. Similar properties hold as well for many other semiparametric models.

Under model (1.1), estimation procedures to date have focused on the finite-dimensional parameters, as has mostly been the case also for estimation under other semiparametric models. Here, we extend the estimation to pointwise and simultaneous inference on the hazard ratio function itself. First, consistency and asymptotic normality of the estimate at a fixed time point are established. Then procedures for constructing pointwise confidence intervals and simultaneous confidence bands for the hazard ratio function are developed, and some modifications are implemented to improve moderate sample size performance. For additional display of the treatment effect, simultaneous confidence bands are also obtained for the average hazard ratio function over a time interval. The average hazard ratio gives a summary measure of treatment comparison and provides a picture of the cumulative treatment effect to augment display of the temporal pattern of the hazard ratio. These hazard ratio estimation procedures are applied to data from the Women's Health Initiative (WHI) estrogen plus progestin clinical trial (Writing Group For the Women's Health Initiative Investigators, 2002; Manson and others, 2003), which yielded a hazard ratio function for the primary coronary heart disease outcome that was decidedly nonproportional. Understanding the hazard ratio function shape in this setting was important to integrating the clinical trial data with a large body of preceding observational literature that had failed to identify an early hazard ratio increase (e.g. Manson and others, 2003; Prentice and others, 2005).

Definewhere b = (b₁,b₂)^T. Let τ < τ₀ be such thatwith probability 1. For t ≤ τ, letwhere ΔH₂(s;b) denotes the jump of H₂(s;b) in s and denotes the left continuous (in s) version of Define the martingale residuals

Yang and Prentice (2005) proposed a pseudo maximum likelihood estimator of β, which is the zero of Q(b), where

with f_i = (f_1i,f_2i)^T, where

Once is obtained, R(t) can be estimated by , and the hazard ratio h(t) can be estimated by

In Appendix A of the Supplementary Material available at Biostatistics online, we show that is strongly consistent for h(t) under model (2.1).

To study the distributional properties of , let

For the asymptotic distribution of , define

From Theorem A2 of Yang and Prentice (2005) and some algebra,

where

Now for , from Lemma A3 in Yang and Prentice (2005) and some algebra,whereLet

For t ≤ τ, define the process

With the representations for Q(β) and , in Appendix B of the Supplementary Material available at Biostatistics online, we show that W_n is asymptotically equivalent to which converges weakly to a zero-mean Gaussian process W^*. The weak convergence of W_n thus follows. The limiting covariance function σ(s,t) of W^* involves the derivative D(t;β) and the derivative matrix in U. Although analytic forms of these derivatives are available, they are quite complicated and cumbersome. Here, we approximate them by numerical derivatives. For the functions B(t),C(t), μ₁(t), μ₂(t), ν₁(t), and ν₂(t), define corresponding by replacing β with , R(t) with and D(t;β) with the numerical derivatives. Similarly, let be the numerical approximation of U. Simulation studies show that the results are fairly stable with respect to the choice of the jump size in the numerical derivatives, and that the choice of n ^{− 1/2} works well. With these approximations, we can estimate σ(s,t), s ≤ t ≤ τ, by

where s(t) converges in probability, uniformly in t over I, to a bounded function s^*(t) > 0. From the weak convergence of W_n to W^* and the functional delta method, we have the weak convergence of V_n to W^*/s^*. Thus, if c_α is the upper αth percentile of sup_tI|W^*/s^*|, an asymptotic 100(1 − α)% simultaneous confidence band for h(t),tI, can be obtained as

It is difficult to obtain c_α analytically. One obvious alternative would be the bootstrapping method. However, it is very time-consuming and results in lower than nominal coverage probabilities in some simulation studies. Lin and others (1993) used a normal resampling approximation to simulate the asymptotic distribution of sums of martingale residuals for checking the Cox regression model. The normal resampling approach reduces computing time significantly and has become a standard method. It has been used in many works, including Lin and others (1994), Cheng and others (1997), Gilbert and others (2002), Tian and others (2005), and Peng and Huang (2007). We will modify this approach for our problem here.

For t ≤ τ, define the process

where ϵ_i,i = 1,…,n, are independent variables that are also independent from the data. Furthermore, these variables have mean zero and variance converging to one as n→∞. In the normal resampling approach mentioned above, the ϵ_i's are the standard normal variables. However, the standard normal variables often result in lower coverage probabilities in various simulation studies. Thus, with moderate sized samples, we need to make some adjustment.

Conditional on (X_i,δ_i,Z_i),i = 1,…,n, is a sum of n independent variables at each time point. In Appendix B of the Supplementary Material available at Biostatistics online, we show that given the data converges weakly to W^*. It follows that given the data converges in distribution to . Therefore, c_α can be estimated empirically from a large number of realizations of the conditional distribution of given the data.

Several choices of the weight s arise from recommendations in the literature for confidence bands of the survivor function and the cumulative hazard function in the one sample case. The choice results in equal precision bands (Nair, 1984), which differ from pointwise confidence intervals in that c_α replaces z_α/2. The choice results in the Hall–Wellner type bands recommended by Bie and others (1987), which often have narrower widths in the middle of data range and wider widths near the extremes of data range (Lin and others, 1994). One could also choose This choice does not involve and hence is easier to implement. It may be adequate when only varies mildly over time.

Let a(0,τ) and define the average hazard ratio, over [a,t],

Note that the average hazard ratio involves an integral of the hazard ratio rather than a ratio of integrated hazards. It provides a measure for the cumulative treatment effect over a time interval to augment the temporal effect display from the hazard ratio estimates. It can be estimated by

To obtain simultaneous confidence bands for the average hazard ratio, let

In Appendix B of the Supplementary Material available at Biostatistics online, we show that converges weakly to the zero-mean Gaussian process . Also, behaves more stably than and its covariance function is comparatively insensitive to the choice of weight function. Hence, for simplicity, we consider only the process

From the functional delta method, it follows that converges weakly to . Thus, if is the upper αth percentile of an asymptotic 100(1 − α)% simultaneous confidence band for can be obtained as

Next, we report the results from some representative simulation studies. Here and for the real data application in Section 5 later, τ was set to exclude the last-order statistic. All numerical computations were done in “Matlab.“ First, under the model of Yang and Prentice (2005), lifetime variables were generated with R(t) chosen to yield the standard exponential distribution for the control group. The values of β were (log(0.9),log(1.2)) and (log(1.2),log(0.8)), representing 1/3 increase or decrease over time from the initial hazard ratio, respectively. The censoring variables were independent and identically distributed with the log-normal distribution, where the normal distribution had mean c and standard deviation 0.5, with c chosen to achieve various censoring rates. The empirical coverage probabilities were obtained from 1000 repetitions, and for each repetition, the critical values c_α and were calculated empirically from 1000 realizations of relevant conditional distributions. The results of these simulations are summarized in Table 1, where the equal precision bands, Hall–Wellner type bands and unweighted bands for the hazard ratio are denoted by EP, HW, and UW, respectively. Results for simultaneous confidence bands of the average hazard ratio are also included with the column header . From Table 1, the empirical coverage probabilities for the hazard ratio were mostly close to the nominal level. The empirical coverage probabilities for the average hazard ratio were mostly conservative. The conservative results were partially due to the finite-sample modifications intended for the hazard ratio. Those modifications improved the performance of the hazard ratio estimation procedure under some scenarios, while yielding conservatism in others, particularly for the more stable average hazard ratio estimator. The coverage probabilities for the equal precision bands overall were closer to the nominal level than other types of bands.

To check the robustness of the proposed procedures, we carried out various simulation studies with monotone hazard ratio not satisfying the model of Yang and Prentice (2005). For Table 2, the control group lifetime variables were standard exponential. The hazard ratio was linear from 0 to the 99th percentile of the standard exponential and continuous and constant afterward. The initial and end hazard ratios again were (0.9,1.2) and (1.2,0.8), respectively, and the censoring variables were the same as before. From Table 2, the results are similar to those in Table 1, with slight undercoverage under some scenarios.

To examine model adequacy, we can use a residual plot that is similar to the method for the Cox regression model (Cox and Snell, 1968). Let Λ_C and Λ_T be the cumulative hazard functions of the 2 groups, respectively. Then Λ_C(T_i),i ≤ n₁,Λ_T(T_i),i > n₁ are i.i.d. from the standard exponential distribution. Let _C and _T be the model-based estimator of Λ_C and Λ_T, respectively, and define the residuals If model (2.1) is correct, the residuals should behave like a censored sample from the standard exponential distribution. Thus, the Aalen–Nelson cumulative hazard estimator based on them should be close to the identity function. If there is noticeable deviation, then model (2.1) is questionable. Similarly, the residual plot can be obtained for the piecewise constant hazards ratio model used in Prentice and others (2005). Both residual plots, not shown here, suggest that the 2 models fit the data adequately, with similar residual behaviors.

The 95% pointwise confidence intervals and simultaneous confidence bands for the hazard ratio function are given in Figure 1. For comparison, the 95% confidence intervals for 0–2, 2–5, and > 5 years from Prentice and others (2005) are included, over the median of uncensored data in each time interval. Compared with the piecewise constant hazards ratio model, the confidence bands do not depend on partitioning of the data range and provide more continuously changing display of the treatment effect. The confidence bands are generally in agreement with the results from Prentice and others (2005). The UW band is wider than the other 2 bands most of the time. The HW band is the narrowest in the middle section but is quite wide at the beginning. Both the EP band and the HW band give narrower intervals for the middle portion of the data range than the piecewise Cox model. Near the end of the data range, all 3 bands have about the same width as the confidence interval from Prentice and others (2005). Overall the EP band matches most closely with the results for the piecewise constant hazards ratio model. The width of the EP band is less than or equal to the piecewise model–based confidence intervals for most of the data range, except at the beginning. Note that the constant function 1 is not excluded in the HW and UW bands. In comparison, the EP band stays above 1 for about the first 600 days. From Prentice and others (2005), the confidence interval for 0 − 2yr excludes 1, indicating an elevation in coronary heart disease risk for the treatment early on. For this data set, the standard error of the estimated hazard ratio begins at 0.43, quickly comes down to below 0.20 at 600 days and stays below 0.20 for the rest of data range. Since the UW band does not take the variance into account and the HW band emphasizes the middle range, the elevated standard error at early follow-up times likely explains the discrepancy among the results. Compared with the original analysis that showed an overall difference between the 2 groups, the results here and those from Prentice and others (2005) give more detailed analysis on the dependence of the hazard ratio on time and help explaining the discrepancy between the results of the WHI clinical trial and preceding observational research, much of which involved cohorts where women could be enrolled some years after initiating hormone therapy.

For the average hazard ratio function, the estimator and the 95% simultaneous confidence band are given in Figure 2. The standard error of the estimated average hazard ratio varies more mildly over time, and both the estimated average hazard ratio and the confidence band are changing much more smoothly compared with the results for the hazard ratio in Figure 1. Note that the confidence band stays above 1 for t < 700 days. This is in agreement with the results of Prentice and others (2005).

To compare with the nonparametric approach, Figure 3 gives the estimated hazard ratio, the 95% pointwise confidence intervals and simultaneous confidence band of Gilbert and others (2002), based on the R programs from the author's site. The same scale as that in Figure 1 is used for comparison and results in truncation of some portion of the plot. The estimated hazard ratio suggests that the hazard ratio is reasonably monotonic. The nonparametric hazard ratio estimate is somewhat lower than the hazard ratio estimates in Figure 1 under either model (2.1) or the piecewise constant hazards ratio model. The confidence band is wider than those in Figure 1 for the beginning and later parts of the data range, reflecting the difficulty in making nonparametric inference on the hazard functions, especially with heavy censoring and in the tail region.

To assess the cumulative treatment effect, we have worked with the average hazard ratio function here, partly due to its close connection with the hazard ratio and its corresponding ready interpretation. Alternatively, the ratios S_T(t)/S_C(t) and (1 − S_T(t))/(1 − S_C(t)) or the difference S_T(t) − S_C(t), could be considered.