- Email: [email protected]
A test for abrupt change in hazard regression models with Weibull baselines
Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
Contents lists available at SciVerse ScienceDirect
Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi
A test for abrupt change in hazard regression models with Weibull baselines Matthew R. Williams a,b,n, Dong-Yun Kim a a b
Department of Statistics, Virginia Tech, Blacksburg, VA, USA Research and Development Division, National Agricultural Statistics Service, United States Department of Agriculture, Fairfax, VA, USA
a r t i c l e i n f o
abstract
Article history: Received 5 June 2012 Received in revised form 18 January 2013 Accepted 18 April 2013
We develop a likelihood ratio test for an abrupt change point in Weibull hazard functions with covariates, including the two-piece constant hazard as a special case. We first define the log-likelihood ratio test statistic as the supremum of the profile log-likelihood ratio process over the interval which may contain an unknown change point. Using local asymptotic normality (LAN) and empirical measure, we show that the profile loglikelihood ratio process converges weakly to a quadratic form of Gaussian processes. We determine the critical values of the test and discuss how the test can be used for model selection. We also illustrate the method using the Chronic Granulomatous Disease (CGD) data. & 2013 Elsevier B.V. All rights reserved.
Keywords: Change point Hazard regression Likelihood ratio test Local asymptotic normality
1. Introduction Due to their flexibility, Weibull distributions are suitable for parametric modeling of time-to-event data and they have been widely used in scientific applications, e.g., in survival analysis, reliability and industrial engineering, and hydrology, among others. To incorporate covariate information into modeling, regression models including covariates with the Weibull baseline are often used for the estimation of the hazard function. There are cases where one or more covariate coefficients may change for some unknown reason, and one may wish to test whether such change exists. In this paper we propose a likelihood-ratio test (LRT) statistic that can be used to test whether there is an abrupt change at an unknown point in the covariate coefficient vector in the Weibull hazard function. Since a change in the covariate parameter vector can be tested using a single test statistic, it is much more convenient than other procedures that test one parameter at a time, and it is useful for model selection purposes. The test allows for staggered entry and type I censoring, as often found in large-scale clinical trials or experiments, and the censoring mechanism is explicitly incorporated into the test statistic. We derive the asymptotic distribution of the test statistic using local asymptotic normality (LAN) and the weak convergence of the empirical measure. An advantage of this approach is that it allows us to approximate the profile log-likelihood ratio test statistic by a quadratic process even when closed-form estimators for model parameters are not available, which is the case even when the change point is assumed to be known.
n Corresponding author at: Research and Development Division, National Agricultural Statistics Service, United States Department of Agriculture, Fairfax, VA, USA. Tel.: +1 703 877 8000x160. E-mail addresses: [email protected] (M.R. Williams), [email protected] (D.-Y. Kim).
0378-3758/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jspi.2013.04.006
iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
2
To be more specific, assume that the hazard function hðt i Þ is of the form ( k−1 t i expðzTi β1 þ wTi αÞ; t i oν hðt i ; zi ; wi Þ ¼ t k−1 expðzTi β2 þ wTi αÞ; t i ≥ν i
ð1Þ
where t 1 ; …; t n are a sequence of independent observations, and β1 , β2 , and α are unknown covariate coefficient vectors. The shape parameter k 4 0 is assumed to be known, but α may be unknown nuisance parameters. The null hypothesis H0: β1 ¼ β2 ¼ β0 leads to an irregular test, since the unknown change point ν vanishes under the null hypothesis, and in this case (1) becomes a Weibull model and the scale parameter is modeled by a regression model via a log link function. For the exponential case (k ¼1), Dupuy (2009) derived non-asymptotic bounds for significance and type II error rates for the LRT, assuming that all coefficients change and data are subject to random right censoring. If semi-parametric approaches using the partial likelihood are used, some computational complexity arising from the baseline hazard λ0 ðt i Þ ¼ t k−1 can be i avoided. In this framework, there are several results for the case when β is a scalar: Using Davies' approximation (Davies, 1987). O'Quigley and Natarajan (2004) obtained a score test for a case equivalent to when α ¼ 0. Liang et al. (1990) proved the weak convergence of the score test process to an Ornstein–Uhlenbeck process. Liu et al. (2008) noted that some of the assumptions for these results are difficult to assess in practice and instead developed a Monte-Carlo method to generate the distribution for the score test statistic. They argued that results can be extended to multiple change points (although the changing coefficient β is still a scalar) and discussed one and two change points in a simulation study. In the next section we establish the weak convergence of the LRT to the supremum of a quadratic process with Gaussian components. Also we discuss how to compute critical values of the test. In Section 3, we illustrate our method using the Chronic Granulomatous Disease (CGD) data as in Fleming and Harrington (1991). We compare the performance of our asymptotic test to the use of Davies' approximation (Davies, 1987). In general, our test produces either comparable or better results. While Davies' approximation can be only applied to a χ 2 process, our results can be applicable to a wider class of quadratic processes where the χ 2 process is a special case. In Section 4, we explore the behavior of our tests under different alternatives. Here we note the potential for model selection by using multiple tests to indicate which parameter(s) may have changed. In Section 5 we briefly discuss our findings. 2. Main result We consider observing the time-to-event data under staggered entry and type I censoring. Staggered entry refers to situations when subjects enter a clinical trial at random times 0 ¼ τ0 o ⋯ o τN . We assume that the entry times τ1 o ⋯ o τn follow a Poisson process. Type I censoring occurs when event times are recorded if occurring within a time interval ½0; T R for some fixed T R 4 0, and censored otherwise. Such situations can occur not just in clinical trials, but also in reliability testing scenarios as well as in observational studies where the entry times cannot be arranged ahead of the study. Assume also that the ending time of the study TR is independent of the number of subjects or events. Conditioning on NðT R Þ ¼ n, the total number of subjects who entered the study before time TR, the censoring time of each subject has a uniform distribution, i.e., ðT R −τi Þ∼Uð0; T R Þ, where τi is the entering time for the ith subject. Suppose the underlying data t follow model (1), but we observe the pair fx; δg, where x ¼ minðt; T R −τÞ and δ is the indicator 1ft o T R −τg . The likelihood for an observation fx; δg is Lðθ; x; δÞ ¼ ðxk−1 expðzT β1 þ wT αÞÞδ exp½−xk expðzT β1 þ wT αÞ=k1fx o νg þ ðxk−1 expðzT β2 þ wT αÞÞδ exp ½−ðexpðzT β1 Þνk þexpðzT β2 Þðxk −νk ÞÞexpðwT αÞ=k1fx≥νg Under the null hypothesis and the uniform distribution of censoring times, f(x) the marginal density for x and f~ ðxÞ the joint density for fx; δ ¼ 1g (both conditioned on covariate vectors z; w) are the following: f ðxÞ ¼ expð−xk η=kÞð1 þ xk−1 ðT R −xÞη=kÞÞ=T R f~ ðxÞ ¼ expð−xk η=kÞðxk−1 ðT −xÞη=kÞÞ=T R
R
ð2Þ
where η ¼ expðz β þ w αÞ. a:s P d Let ⟶, ⟶, ⟶, and ⟹ represent almost sure convergence, convergence in probability, convergence in distribution, and weak convergence of a process, respectively. Let lν ðθÞi ¼ lν ðθ; xi ; δi ; zi ; wi Þ be the log-likelihood of the parameter vector θ ¼ ðβ1 ; β2 ; αÞT for model (1) with an observation fxi ; δi g, covariate vectors zi and wi and a given change point ν. Now define _l ν ðθ Þ and €l ν ðθ Þ as the first and second derivatives of the log-likelihood with respect to θ evaluated at θ ¼ θ . The vector θ 0 i 0 i 0 0 takes the null hypothesis values ðβ0 ; β0 ; α0 ÞT . Now consider fx; δg the paired vectors of observations fx1 ; δ1 g; …; fxn ; δn g with corresponding covariate matrices Z ¼ ½z′1 ; …; z′n ′ and W ¼ ½w′1 ; …; w′n ′. We reparameterize the log-likelihood lν ðθÞ ¼ ∑ni¼ 1 lν ðθÞi and approximate with a Taylor expansion about the null values θ0 : T
T
n pffiffiffi lν ðθÞ ¼ lν ðθ0 þ h= nÞ ¼ ∑ lν ðθ0 Þi þ W n ðθ0 ; hÞ i¼1
iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
3
where T
T
n h n h n € ∑ l ν ðθ0 Þi h þ ∑ Ri ðθ0 ; hÞ W n ðθ0 ; hÞ ¼ pffiffiffi ∑ _l ν ðθ0 Þi þ 2n i ¼ 1 ni¼1 i¼1
where the h vector is now the parameter of interest and Ri ðθ0 ; hÞ is the remainder term. The log-likelihood ratio test statistic (given ν) is now ! 2Λn ðνÞ ¼ −2
sup W n ðθ0 ; hÞ− sup W n ðθ0 ; hÞ
fh∈H0 g
fh∈Hg
where H0 is the restricted support of h under the null hypothesis and H is the full support. Define C ij ¼ Covð_l νi ðθ0 Þ; _l νj ðθ0 ÞÞ with 0 o a≤νi ≤νj ≤bo T R . The design parameters a and b specify the range of a change point. These are needed to avoid the boundary problem. For νi given, the Fisher information under the null hypothesis is −1 I θ0 ðνi Þ ¼ C ii . Define I β0 ðνi Þ as the subset of the Fisher information applying only to β1 and β2 . Now define Σ ij ¼ I −1 θ0 ðνi ÞC ij I θ0 ðνj Þ. Let X ¼ ðXβ1 ; Xβ2 ; Xα ÞT be a multivariate Gaussian process with mean zero and covariance Σ ij . Let Xβ ¼ ðXβ1 ; Xβ2 ÞT . Finally let I d;2;1 ¼ ½I d ; I d T where Id is the identity matrix of dimension d, the length of β0 . Theorem 1. For model (1) when H 0 : β1 ¼ β2 is true: for ν∈½a; b and 0 o a o bo T R , 2Λn ðνÞ converges weakly to the quadratic process XTβ QXβ ðνÞ, where Q ¼ I β0 −I β0 I d;2;1 ðI Td;2;1 I β0 I d;2;1 Þ−1 I Td;2;1 I β0 Proof. First we note that model (1) satisfies the conditions for local asymptotic normality (van der Vaart, 1998, p. 104). See Appendix A for details. Then for each fixed ν we have convergence in distribution (van der Vaart, 1998, p. 232). That is, d
2Λn ðνÞ⟶2ΛðνÞ where 2ΛðνÞ ¼ inf ðX−hÞT I θ0 ðX−hÞ− inf ðX−hÞT I θ0 ðX−hÞ fh∈H0 g
fh∈Hg
^ T I θ ðX−hÞ; ^ ¼ ðX−h^ 0 ÞT I θ0 ðX−h^ 0 Þ−ðX−hÞ 0
ð3Þ
pffiffiffi n _ d ^ nð∑i ¼ 1 l ν ðθ0 Þi =nÞ⟶I θ0 X. Details for solving h∈H and h^ 0 ∈H0 are provided in Appendix B. We then use and X∼Nð0; I −1 θ0 Þ, and empirical measure (Appendix A) to finally establish the weak convergence 2Λn ðνÞ⟹2ΛðνÞ for all ν∈½a; b. □ Corollary 2. Assume α in (1) is known. Then XTβ QXβ ðνÞ is a χ 2d process with d degrees of freedom. T 2 Proof. If α is known then X ¼ Xβ and I θ0 ¼ I β0 . Then Q I −1 θ0 is idempotent, which implies for each ν∈½a; b: Xβ QXβ is a χ d random variable with d degrees of freedom (Graybill, 1976).
Note that if α ¼ 0, then all covariate coefficients change after ν. Corollary 3. If β1 and β2 are scalars, then XTβ QXβ ðνÞ ¼ Z 2 ðνÞ where ZðνÞ is a mean zero Gaussian process with covariance ðfΣ ij g11 þ fΣ ij g22 Þ−ðfΣ ij g12 þ fΣ ij g21 Þ CovðZðνi Þ; Zðνj ÞÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 11 12 21 11 22 12 21 ððI i þ I 22 i Þ−ðI i þ I i ÞÞððI j þ I j Þ−ðI j þ I j ÞÞ
ð4Þ
−1 where I lk i ¼ fI β0 ðνi Þ glk for l; k ¼ 1; 2 where fglk refers to the elements with rows corresponding to βl and columns to βk .
Proof. The case of β1 and β2 as scalars is comparable to the results from Williams and Kim (2011). That is, the matrix I d;2;1 reduces to 1 ¼ ð1; 1ÞT . Q is a 2 2 matrix with 1T Q ¼ 0. Furthermore fQ g11 ¼ jI β0 j=1T I β0 1 ¼ 1=ððI 11 þ I 22 Þ−ðI 12 þ I 21 ÞÞ, where 2 −1 2 I lk i ¼ fI β0 ðνi Þ gl;k and jI β0 j is the determinant of the matrix I β0 . Then 2ΛðνÞ ¼ fQ g11 ðX β1 −X β2 Þ ¼ Z ðνÞ. □ Note that if β1 and β2 are scalars and α is known, then VarðZðνi Þ; Zðνj ÞÞ ¼ 1 and Z 2 ðνÞ is a χ 21 process. Corollary 4. For model (1) when H0 is true, sup f2Λn ðνÞ : a≤ν≤bg converges in distribution to supfXTβ QXβ ðνÞ : a≤ν≤bg where XTβ QXβ ðνÞ is as in Theorem 1. The significance level αn for a critical value c is then " # αn ¼ P sup XTβ QXβ ðνÞ 4 c : a≤ν≤b
Proof. It is clear by direct application of Slutsky's Theorem and the Continuous Mapping Theorem.
□
3. Illustration: chronic granulomatous disease data To illustrate our asymptotic test, we consider the chronic granulomatous disease (CGD) data from Fleming and Harrington (1991). CGD is a hereditary disease of the immune system which leads to serious, often recurring infections. We only consider the first infection after the treatment either by placebo (control group) or gamma interferon (treatment iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
4
group). In this clinical trial, patients entered the study at different start times (staggered entry). Responses are time to first infection after treatment or the right censoring time (rescaled to TR ¼1). There are 20 uncensored events out of 53 for the control, and only 7 of 55 for the treatment are uncensored. Therefore, heavy censoring occurred to both groups. Age and sex were included as covariates. From exploratory analysis (Williams and Kim, 2013), k¼2 was determined to be a reasonable assumption so it is assumed to be known. Additionally, the survival curves for the control group indicate a possible change point, while those for the treatment group do not. To test this formally, we assume (1) has the form ( hðt i ; zi ; wi Þ ¼
t i expðControli β1 þ ð1; Agei ; Sexi ÞT αÞ; t i expðControli β2 þ ð1; Agei ; Sexi ÞT αÞ;
ti o ν t i ≥ν
ð5Þ
where Controli ¼1 for patients receiving the placebo and Controli ¼ 0 for the treatment. This slight reversal of indicators (usually control is indicated by 0) allows for the baseline of the control curve to change without affecting the baseline treatment curve. We use the maximum likelihood estimate θ^ 0 and an interval ½a; b based on the 80% coverage of the control group data with a¼0.057 as the 10% quantile and b¼0.748 as the 90% quantile. To generate the critical values we create a grid of 1000 ν values spread evenly from a to b. We then generate realizations of the corresponding quadratic processes XTβ QXβ ðνÞ ¼ Z 2 ðνÞ (by Corollary 3) substituting θ^ 0 for θ0 (See Appendix B for details of constructing CovðXðνi Þ; Xðνj ÞÞÞ. We generate 10,000 realizations and find the maximum of each over the grid of ν. The upper 1% quantile of these maximum is then an estimate of the 1% critical value of the test. In this case, the observed test statistic of 37.79 exceeds the estimated 1% critical value of 13.8, so the null hypothesis is rejected at 1% significance level. The estimated change point for the control group is ν^ ¼ 0:072, or about 23 days after the placebo treatment. To assess type I error and power, data were simulated based on the maximum likelihood estimates under the null β^ 0 ¼ 1:54 and α^ 0 ¼ −0:15; −0:43; −0:45 and alternative hypotheses β^ 1 ¼ 4:56, β^ 2 ¼ 1:04, and α^ ¼ −0:27; −0:39; −0:28. Covariates z (control) and w (intercept, age, and sex) were resampled from the 108 observations in the CGD data. Likelihood ratio test statistics were calculated and compared to the critical values obtained from Corollary 3. This was repeated 1000 times for sample sizes of 100, 200, and 500. Type I error was calculated as the percentage of times that the test statistic exceeded the 5% and 1% critical values when the data was generated under the null hypothesis. Power was calculated as the percentage of times that the test statistic exceeded the 5% and 1% critical values when the data was generated under the alternative hypothesis. The rates obtained by using Corollary 3 are compared to rates from applying Davies' approximation (Davies, 1987), which allows us to approximate p-values for the supremum of χ 2 processes. While the result of Davies only strictly applies to χ 2 processes, it is possible that Z2 can be approximated by a χ 21 process if VarðZðνÞÞ≈1, as it is for the CGD example. The results for the Type I error simulations are summarized in Table 1. We see good convergence at a sample size of 100. Rates for our method (LAN) are comparable to those for Davies. Furthermore, our method is within a standard error of the nominal levels at each sample size. Interestingly, the test using Davies' approximation appears to be conservative for the Table 1 Type I error simulations using CGD data, comparing results from Corollary 3 (LAN) and Davies (1987) for 1000 replications with ½a; b ¼ ½0:057; 0:748. Test
Sample size
s.e. 5% s.e. 1% LAN 5% LAN 1% Davies 5% Davies 1%
100
200
500
0.022 0.010 0.047 0.012 0.053 0.019
0.015 0.007 0.035 0.007 0.038 0.010
0.010 0.004 0.043 0.006 0.020 0.003
Table 2 Power simulations using CGD data, comparing results from Corollary 3 (LAN) and Davies (1987) for 1000 replications with ½a; b ¼ ½0:057; 0:748. Test
LAN 5% LAN 1% Davies 5% Davies 1%
Sample size 100
200
500
0.950 0.895 0.961 0.911
0.998 0.998 0.998 0.998
1.000 1.000 1.000 1.000
iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
5
sample size n ¼500. Power simulations based on θ^ are displayed in Table 2. Power is high for both methods, and it is essentially the same. The close to nominal type I error combined with high power demonstrates that our test using both the LAN and Davies methods is appropriate and useful for this CGD data set. 4. Simulation study: exponential hazard regression We consider a constant hazard function with covariate information zi ∼Nð0; 1Þ: hðt i ; zi Þ ¼ expðβ0 þ β1 zi Þ
ð6Þ
We choose the simple exponential case ðk ¼ 1Þ, because CovðXðνi Þ; Xðνj ÞÞ is greatly simplified (Appendix B) and consistent estimators for parameters under some change point alternatives have already been studied (Dupuy, 2006). Based on our results from Section 2, we can construct tests for three different change point alternatives. 1. The intercept β0 changes at unknown ν : β01 ≠β02 2. The slope β1 changes at unknown ν : β11 ≠β12 3. Both the slope and intercept change at the same ν : β01 ≠β02 and β11 ≠β12 Letting β0 ¼ 1, β1 ¼ 1, and TR ¼1; we examine the behavior of these three tests. Under the null hypothesis of no change, we observe approximately 40% censoring. We choose the change point interval to be the expected 25% and 75% quantiles of the data: a¼0.06 and b¼0.37 respectively. To obtain critical values, we generate 10,000 realizations of each of the Z2 process for the single change tests and XTβ QXβ for the vector change. Each realization occurs over a grid of 1000 values of ν∈½a; b. The upper 5% and 1% quantiles from the supremum of each realization are then selected as the critical values. For example, we obtain 9.43, 8.77, and 11.33 as the 5% critical values and 13.20, 12.33, and 15.50 as the 1% for the test for changes in intercept, slope, and both, respectively. Even at sample size 100, the empirical type I error rates approach their nominal values. Interestingly results using the Davies (1987) method are quite conservative, resulting in empirical type I error rates less than 0.01 for a 5% level test. We next consider the behavior of the tests under the three alternatives. Table 3 displays the power for each of the three tests under all three alternatives. Only results for positive changes in β using 5% critical values are included. Results are similar for negative changes in β and for 1% critical values. Not surprisingly, a given test turns in the best result in terms of power when its specific alternative is true. On the other hand, the joint test retains good power when only one coefficient changes. When both parameters change, it appears that all of the three tests require much larger sample size. This is just an artifact of choosing a smaller effect size (1.5 instead of 2) so as not to have the power be close to 1.000 for both the 200 and 500 sample sizes. Together these results suggest that evaluating all three tests (perhaps with some correction for multiple testing) may be useful for model selection when determining not only if a change occurred but also which regression parameters are involved. We might interpret the joint test as a global test with reasonable power over a range of alternatives and consider the single test as follow-up with higher power for more specific alternatives. 5. Discussion In this paper we have proposed a likelihood-ratio test for an unknown change point in hazard regression models where one or more covariate coefficients change abruptly. The test allows for staggered entry and type I censoring so it is useful for Table 3 Power for single and joint change alternatives to model (6) using 5% critical values. Each scenario was simulated 1000 times for a standard error o 0:016. Test
Sample size 100
200
500
H A : β02 ¼ 2 Intercept change Covariate change Both change
0.904 0.060 0.875
0.999 0.099 0.994
1.000 0.225 1.000
H A : β12 ¼ 2 Intercept change Covariate change Both change
0.068 0.638 0.573
0.124 0.929 0.879
0.244 1.000 0.999
H A : β02 ¼ β12 ¼ 1:5 Intercept change Covariate change Both change
0.217 0.221 0.331
0.493 0.513 0.719
0.958 0.942 0.993
iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
6
situations where the entry times of time-to-event data are random, as shown in the study of CGD data in Section 3. We have derived the asymptotic distribution of the test statistic so time-consuming calculations can be avoided. Our procedure can test a change in a single parameter or multiple parameters in the covariate coefficient vector, it is possible to develop a battery of tests that are useful for model selection purposes. The LAN based approach combined with the weak convergence of the empirical measure enables us to approximate the profile log-likelihood ratio test statistic by a quadratic process, which greatly extends the use of tests when the limit is not a χ 2 process, as dictated in Davies' approximation. While we have focused on developing the likelihood-ratio test, our weak convergence results can also be used to build score and score-related tests. Our derivations and proofs use all the components needed for a score test (the score function and its covariance), so future use of score tests for model (1) can be facilitated by these results. Critical values for the test can be obtained in several ways, depending on the particular application. The most general method is to generate realizations of the multivariate Gaussian processes Xβ and then calculate the quadratic form XTβ QXβ ðνÞ from Theorem 1 and then use Corollary 4. If a test is desired to test only one coefficient in the covariate vector, we only need to generate the univariate Gaussian process ZðνÞ and use Corollary 3 and Corollary 4. Finally, if the covariate information is not used, the current nonstationary Gaussian process can be transformed to a stationary Ornstein–Uhlenbeck process (Kim et al., 2004; Williams and Kim, 2013) for which tail approximations for the supremum of the process are available (Leadbetter et al., 1983). For the example discussed in Section 4, we used four normal variates for each value of ν corresponding to two parameters before and two after the change. If more covariates are considered for β, only the dimension of the covariance matrix increases as a result; the test procedure remains identical. On the practical side, however, very large number of covariates in β may cause computational problems. Depending on the number of covariates and computing capability, it may require reduction of the number of grid points at which the covariance matrix is evaluated, or selection of a subset of covariates to be entered into β in the model. This work has focused on a change in the linear predictor of the hazard function when the Weibull shape parameter k is assumed to be known. We conducted computer simulations for the cases where the nuisance parameter k is unknown. We found that the maximum likelihood estimate for the unknown shape parameter is numerically unstable in some cases. Thus, further study is needed to get around the problem to obtain a stable estimate for the shape parameter. On the other hand, testing for a change from one known value of k to another, unknown value should be possible. For example, Williams and Kim (2011) consider the test of no change versus alternative hypothesis where k¼1 before the unknown change point, and k assumes some unknown value larger than 1 after the change. These results did not incorporate censoring or covariates, but they do consider tests under various forms of continuous hazard function. Focusing on abrupt change and using the techniques presented in this work should allow for the establishment of corresponding convergence results for likelihood-ratio tests.
Acknowledgments The first author was sponsored through the Institute for Critical Technology and Applied Science (ICTAS) Doctoral Scholars Program at Virginia Tech. The work of the second author was partially supported as part of the National Science Foundation Biocomplexity of Coupled Human and Natural Systems Program, Award No. BCS-0709671. Appendix A. Proof of weak convergence For the independent and identically distributed case (for example, see Williams and Kim, 2011), we would consider the likelihood as a density pθ ðx; δÞ ¼ Lðθ; x; δÞ. However, covariate information leads to nonidentical distributions. To circumvent this, we treat the covariates as random variables and consider the joint distribution (see the generalized linear model example in van der Vaart, 1998, p. 203): pθ ðx; δ; z; wÞ ¼ f ðx; δjz; wÞgðz; wÞ ¼ LðθÞL′ðγÞ where LðθÞ is the likelihood component for θ and L′ðγÞ is the likelihood component for the unrelated nuisance parameter vector γ. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi For our model pθ ðx; δ; z; wÞ is continuously differentiable with respect to θ for every ðx; δÞ, z, w, and ν. The Fisher Information matrix component for θ is I θ ¼ Ez;w ðEx;δjz;w ð€l ν ðθÞÞÞ: For model (1) under staggered entry, Ex;δjz;w ðl€ν ðθÞÞ is well-defined and continuous with respect to θ for any ν. Then for a compatible marginal distribution gðz; wÞ, I θ is also well-defined and continuous with respect to θ for any ν. We will assume z; w have such a compatiblep distribution ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiand will approximate it using the empirical distribution of the observed covariates. By the above properties of pθ ðx; δ; z; wÞ and I θ , our experiment is differentiable in quadratic mean (van der Vaart, 1998, p. 95). Therefore local asymptotic normality holds for each ν (van der Vaart, 1998, p. 104). For (1) we show that all elements of the vector-valued function _l ν ðθ0 Þ are in one Donsker class (van der Vaart, 1998, p. 269) indexed by ν. We construct an envelope function Mðx; δ; z; wÞ such that Mðx; δ; z; wÞ≥∥∂_l ν ðθ0 Þ=∂ν∥1 for ν∈½a; b and iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
7
EðjMðx; δ; z; wÞjÞ o ∞. The norm ∥ ∥1 is the sum of the absolute value of the elements of ∂l_ν ðθ0 Þ=∂ν. Note that Mðx; δ; z; wÞ is not a function of ν. The existence of this envelope function satisfies the conditions for a Donsker class (van der Vaart, 1998, p. 271). Based on the derivatives (See Appendix C) we choose Mðx; δ; z; wÞ ¼ jzexpðzT β0 þ wT α0 Þðνn Þk−1 j where νn maximizes νk−1 for fixed k 40 and ν∈½a; b. If 0 ok o1, then νn ¼ a. Otherwise νn ¼ b. We note that EðjMðx; δ; z; wÞjÞ ¼ Ez;w ðjMðx; δ; z; wÞjÞ since M is constant with respect to x and δ. If we can assume that the distribution gðz; wÞ is such that Ez;w ðjzexpðzT β0 þ wT α0 Þðνn Þk−1 jÞ o∞, then EðjMðx; δ; z; wÞjÞ o ∞. By our construction ! pffiffiffi 1 n _ ∑ l ν ðθ0 Þ ⟹I θ0 ðνÞXðνÞ; ν∈½a; b n ni¼1 where XðνÞ is a zero-mean Gaussian process with CovðXðνi Þ; Xðνj ÞÞ ¼ Σ ij .
For our model, we would like to show that each element in the matrix of functions €l ν ðθ0 Þ is Glivenko–Cantelli classes (van der Vaart, 1998, p. 269) indexed by ν∈½a; b. We construct an envelope function Fðx; δ; z; wÞ, such that Fðx; δ; z; wÞ≥∥€l ν ðθ Þ∥ 0
1
and EðjFðx; δ; z; wÞjÞ o ∞ for ν∈½a; b. This process is very similar to that of constructing Mðx; δ; z; wÞ. The existence of this envelope function satisfies the conditions for a Glivenko–Cantelli class (van der Vaart, 1998, p. 272). Alternatively, we can show membership of l€ν ðθ Þ in a Donsker class directly, which implies membership in a Glivenko–Cantelli class. For this we 0
consider ∥∂€l ν ðθ0 Þ=∂ν∥1 with envelope function M 2 ðx; δ; z; wÞ ¼ ð∥zzT ∥1 þ ∥2zwT ∥1 þ ∥wwT ∥1 ÞexpðzT β0 þ wT α0 Þðνn Þk−1 where νn maximizes νk−1 for fixed k 4 0 and ν∈½a; b. If 0 o k o 1, then νn ¼ a. Otherwise νn ¼ b. Again, M2 is a constant with respect to x and δ, so as long as the marginal density gðz; wÞ is compatible, we have established a Donsker class. Then P ∥n−1 ∑n €l ν ðθ Þ−I θ ∥ ⟶0. sup ν∈½a;b
i¼1
0
0
1
We next consider the remainder T Ri ðθ0 ; hÞ ¼ R1i ðθ0 ; hÞ−h €l ν ðθ0 Þh=2n
where pffiffiffi pffiffiffi T R1i ðθ0 ; hÞ ¼ lν ðθ0 þ h= nÞ−lν ðθ0 Þ−h _l ν ðθ0 Þ= n: pffiffiffi T T According to Taylor's Theorem: R1i ðθ0 ; hÞ ¼ h €l ν ðθ0 þ hϕ = nÞh=2n, where hϕ ¼ ðϕ1 h1 ; …; ϕd hd Þ for some ϕ1 ; …; ϕd ∈½0; 1 with d pffiffiffi the dimension of θ. For fixed θ0 and h and for n large enough, ½a; b ½θ0 ; θ0 þ h= n is well-defined. Furthermore, each pffiffiffi θ∈½θ0 ; θ0 þ h= n has a corresponding envelope M 2 ðθÞ which is an exponential function of a linear combination of θ. Thus pffiffiffi pffiffiffi there exists an envelope M2 ðθn Þ with θn ∈½θ0 ; θ0 þ h= n which dominates. Then €l ν ðθ0 þ hϕ = nÞ is Donsker (and Glivenko– pffiffiffi pffiffiffi Cantelli) since ½θ0 ; θ0 þ hϕ = n⊂½θ0 ; θ0 þ h= n. We can construct an appropriate envelope function for the difference €l ν ðθ0 þ pffiffiffi € P hϕ = nÞ−l ν ðθ0 Þ using the triangle inequality, so it is also Glivenko-Cantelli. Thus supν∈½a;b j∑ni¼ 1 Ri ðθ0 ; hÞj⟶0. We note here that our envelopes Mðx; δ; z; wÞ and M2 ðx; δ; z; wÞ are not functions of either x or δ. This suggests that any well-defined (local asymptotic normality holds for each ν) right-censoring regime for model (1) will possess the weak convergence properties demonstrated above. Of course the details of Cij will change, since f(x) and f~ ðxÞ from (2) depend on censoring structure, but we expect the results from Theorem 1 to hold for broader classes of censoring. One example is that of simple right censoring at time TR. It is a straight-forward exercise to show that the results of the previous sections apply to this simple case as well. Appendix B. Specifying the covariance structure To characterize the 2ΛðνÞ process, we need −1 CovðXðνi Þ; Xðνj ÞÞ ¼ I −1 θ0 ðνi ÞC ij I θ0 ðνj Þ;
where C ij ¼ Covð_l νi ðθ0 Þ; _l νj ðθ0 ÞÞ. From local asymptotic normality (LAN): Ex;δjz;w ð_l ν ðθ0 ÞÞ ¼ 0. Then C ij ¼ Ez;w ðEx;δjz;w ð_l νi ðθ0 Þ_l νj ðθ0 ÞÞÞ. The inner expectation (which we name Cijl) can be evaluated directly through integration using (2), the distributions of fxl ; δl g conditional on zl ; wl for observation l ¼ 1; …; n. The integrands change depending on the domain of x∈½0; νi Þ; ½νi ; νj Þ; ½νj ; T R . The outer expectation can be approximated by ð1=nÞ∑nl¼ 1 C ijl , the expectation with respect to empirical distribution of the covariates zl ; wl . The following are the covariance equations for Cijl conditioned on the covariate values zl ; wl : 2 3 C β1 ;β1 C β1 ;β2 C β1 ;α 6 C β ;β C β ;β C β ;α 7 C ijl ¼ 4 2 1 5 2 2 2 C α;β1 C α;β2 C α;α iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
8
with C β1 ;β1 ¼ zl zTl ½1 þ expð−νki ηl =kÞðνi =T R −1Þ þ ðΓ½1=k; νki ηl =k−Γ½1=kÞ=ððηl =kÞ1=k kT R Þ C β1 ;β2 ¼ 0 C β1 ;α ¼ zl wTl ½1 þ expð−νkj ηl =kÞðνj =T R −1Þ þ ðΓ½1=k; νkj ηl =k−Γ½1=kÞ=ððηl =kÞ1=k kT R Þ C β2 ;β1 ¼ zl zTl ½expð−νkj ηl =kÞðνj =T R −1Þ−expð−νki ηl =kÞðνi =T R −1Þ þ ðΓ½1=k; νkj ηl =k−Γ½1=k; νkj ηl =kÞ=ððηl =kÞ1=k kT R Þ C β2 ;β2 ¼ zl zTl ½−expð−νkj ηl =kÞðνj =T R −1Þ þ ðΓ½1=k; T kR ηl =k−Γ½1=k; νkj ηl =kÞ=ððηl =kÞ1=k kT R Þ C β2 ;α ¼ zl wTl ½−expð−νki ηl =kÞðνi =T R −1Þ þ ðΓ½1=k; T kR ηl =k−Γ½1=k; νki ηl =kÞ=ððηl =kÞ1=k kT R Þ C α;β1 ¼ wl zTl ½1 þ expð−νki ηl =kÞðνi =T R −1Þ þ ðΓ½1=k; νki ηl =k−Γ½1=kÞ=ððηl =kÞ1=k kT R Þ C α;β2 ¼ wl zTl ½−expð−νkj ηl =kÞðνj =T R −1Þ þ ðΓ½1=k; T kR ηl =k−Γ½1=k; νkj ηl =kÞ=ððηl =kÞ1=k kT R Þ C α;α ¼ wl zTl ½1−expð−T kR ηl =kÞ þ ðΓ½1=k; T kR ηl =k−Γ½1=kÞ=ððηl =kÞ1=k kT R Þ where ηl ¼ expðzTl β0 þ wTl α0 Þ, Γ½z ¼ function.
R∞ 0
t z−1 e−t dt is the gamma function, and Γ½z; a ¼
R∞ a
t z−1 e−t dt is the incomplete gamma
Now that we have Cij and therefore I θ0 ðνÞ, we solve for h^ 0 and h^ from (3). Under the alternative h is unrestricted, therefore h^ ¼ X making the second term in (3) equal to 0. Under the null hypothesis, h0 is restricted to be ðh00 ; h00 ; hα ÞT , Þ−1 I T I β Xβ . Then h^ ¼ ðh^ ; h^ ; where the third component is unrestricted. From linear algebra: h^ ¼ ðI T I β I 00
d;2;1
0
d;2;1
d;2;1
0
0
00
00
Xα Þ ¼ ðh^ 00 I d;2;1 ; Xα ÞT . Plugging in h^ 00 into the left term in (3) and factoring out Xβ , we get XTβ QXβ as in Theorem 1. It is T
clear that I Td;2;1 Q ¼ 0, making Q rank d. Appendix C. Derivatives of the log-likelihood We include the first and second derivatives evaluated at θ ¼ θ0 . From these, we can construct envelope functions Mðx; δ; z; wÞ, Fðx; δ; z; wÞ, and M2 ðx; δ; z; wÞ and calculate Cijl via integration. Let η ¼ expðzT β0 þ wT α0 Þ. _l ν ðθ Þj 0 β1 ¼ β2 ¼ β0 ;α ¼ α0 ∂lν ðθ0 Þ=∂β1 ¼ zðδ−ηxk =kÞ1fx o νg þ zð−ηνk =kÞ1fx≥νg ∂lν ðθ0 Þ=∂β2 ¼ zðδ−ηðxk −νk Þ=kÞ1fx≥νg ∂lν ðθ0 Þ=∂α ¼ wðδ−ηxk =kÞ ∂_l ν ðθ0 Þ=∂νjβ1 ¼ β2 ¼ β0 ;α ¼ α0 ∂2 lν ðθ0 Þ=∂ν∂β1 ¼ −zηνk−1 1fx≥νg ∂2 lν ðθ0 Þ=∂ν∂β2 ¼ zηνk−1 1fx≥νg ∂2 lν ðθ0 Þ=∂ν∂α ¼ 0 €l ν ðθ Þj 0 β1 ¼ β2 ¼ β0 ;α ¼ α0 ∂2 lν ðθ0 Þ=∂βT1 ∂β1 ¼ zzT ð−ηxk =kÞ1fx o νg þ zzT ð−ηνk =kÞ1fx≥νg ∂2 lν ðθ0 Þ=∂βT2 ∂β1 ¼ 0
∂2 lν ðθ0 Þ=∂αT ∂β1 ¼ zwT −ηxk =k 1fx o νg þ zwT ð−ηνk =kÞ1fx≥νg ∂2 lν ðθ0 Þ=∂βT2 ∂β2 ¼ zzT ð−ηðxk −νk Þ=kÞ1fx≥νg ∂2 lν ðθ0 Þ=∂αT ∂β2 ¼ wzT ð−ηðxk −νk Þ=kÞ1fx≥νg ∂2 lν ðθ0 Þ=∂αT ∂α ¼ −wwT ηxk =k ∂€l ν ðθ0 Þ=∂νjβ1 ¼ β2 ¼ β0 ;α ¼ α0 ∂3 lν ðθ0 Þ=∂ν∂βT1 ∂β1 ¼ zzT ð−ηνk−1 Þ1fx≥νg ∂3 lν ðθ0 Þ=∂ν∂βT2 ∂β1 ¼ 0 ∂3 lν ðθ0 Þ=∂ν∂αT ∂β1 ¼ zwT ð−ηνk−1 Þ1fx≥νg ∂3 lν ðθ0 Þ=∂ν∂βT2 ∂β2 ¼ zzT ðηνk−1 Þ1fx≥νg iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
9
∂3 lν ðθ0 Þ=∂ν∂αT ∂β2 ¼ wzT ðηνk−1 Þ1fx≥νg ∂3 lν ðθ0 Þ=∂ν∂αT ∂α ¼ 0
References Davies, R.B., 1987. Hypothesis testing when a nuisance parameter is present only under the alternatives. Biometrika 74, 33–43. Dupuy, J.-F., 2006. Estimation in a change-point hazard regression model. Statistics and Probability Letters 76, 182–190. Dupuy, J.-F., 2009. Detecting change in a hazard regression model with right-censoring. Journal of Statistical Planning and Inference 139, 1578–1586. Fleming, T.R., Harrington, D.P., 1991. Counting Processes and Survival Analysis. John Wiley Sons. Graybill, F.A., 1976. Theory and Application of the Linear Model. Duxbury Press. Kim, D., Woodroofe, M., Wu, Y., 2004. Testing for a change in the hazard rate with staggered entry. Communications in Statistics—Theory and Methods 33, 2041–2058. Invited paper for the special issue in honor of Professor Z. Govindarajulu on his 70th birthday. Leadbetter, M.R., Lindgren, G., Rootzén, H., 1983. Extremes and Related Properties of Random Sequences and Processes. Springer, New York. Liang, K.-Y., Self, S.G., Liu, X., 1990. The Cox proportional hazards model with change point: an epidemiologic application. Biometrics 46, 783–793. Liu, M., Lu, W., Shao, Y., 2008. A Monte Carlo approach for change-point detection in the cox proportional hazards model. Statistics in Medicine 27, 3894–3909. O'Quigley, J., Natarajan, L., 2004. Erosion of regression effect in a survival study. Biometrics 20, 344–351. van der Vaart, A.W., 1998. Asymptotic Statistics. Cambridge University Press. Williams, M.R., Kim, D., 2011. Likelihood ratio tests for continuous monotone hazards with unknown change point. Statistics and Probability Letters 81, 1599–1603. Williams, M.R., Kim, D., 2013. A test for an abrupt change in Weibull hazard functions with staggered entry and type I censoring. Communications in Statistics—Theory and Methods 42, 1922–1933.
iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
Contents lists available at SciVerse ScienceDirect
Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi
A test for abrupt change in hazard regression models with Weibull baselines Matthew R. Williams a,b,n, Dong-Yun Kim a a b
Department of Statistics, Virginia Tech, Blacksburg, VA, USA Research and Development Division, National Agricultural Statistics Service, United States Department of Agriculture, Fairfax, VA, USA
a r t i c l e i n f o
abstract
Article history: Received 5 June 2012 Received in revised form 18 January 2013 Accepted 18 April 2013
We develop a likelihood ratio test for an abrupt change point in Weibull hazard functions with covariates, including the two-piece constant hazard as a special case. We first define the log-likelihood ratio test statistic as the supremum of the profile log-likelihood ratio process over the interval which may contain an unknown change point. Using local asymptotic normality (LAN) and empirical measure, we show that the profile loglikelihood ratio process converges weakly to a quadratic form of Gaussian processes. We determine the critical values of the test and discuss how the test can be used for model selection. We also illustrate the method using the Chronic Granulomatous Disease (CGD) data. & 2013 Elsevier B.V. All rights reserved.
Keywords: Change point Hazard regression Likelihood ratio test Local asymptotic normality
1. Introduction Due to their flexibility, Weibull distributions are suitable for parametric modeling of time-to-event data and they have been widely used in scientific applications, e.g., in survival analysis, reliability and industrial engineering, and hydrology, among others. To incorporate covariate information into modeling, regression models including covariates with the Weibull baseline are often used for the estimation of the hazard function. There are cases where one or more covariate coefficients may change for some unknown reason, and one may wish to test whether such change exists. In this paper we propose a likelihood-ratio test (LRT) statistic that can be used to test whether there is an abrupt change at an unknown point in the covariate coefficient vector in the Weibull hazard function. Since a change in the covariate parameter vector can be tested using a single test statistic, it is much more convenient than other procedures that test one parameter at a time, and it is useful for model selection purposes. The test allows for staggered entry and type I censoring, as often found in large-scale clinical trials or experiments, and the censoring mechanism is explicitly incorporated into the test statistic. We derive the asymptotic distribution of the test statistic using local asymptotic normality (LAN) and the weak convergence of the empirical measure. An advantage of this approach is that it allows us to approximate the profile log-likelihood ratio test statistic by a quadratic process even when closed-form estimators for model parameters are not available, which is the case even when the change point is assumed to be known.
n Corresponding author at: Research and Development Division, National Agricultural Statistics Service, United States Department of Agriculture, Fairfax, VA, USA. Tel.: +1 703 877 8000x160. E-mail addresses: [email protected] (M.R. Williams), [email protected] (D.-Y. Kim).
0378-3758/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jspi.2013.04.006
iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
2
To be more specific, assume that the hazard function hðt i Þ is of the form ( k−1 t i expðzTi β1 þ wTi αÞ; t i oν hðt i ; zi ; wi Þ ¼ t k−1 expðzTi β2 þ wTi αÞ; t i ≥ν i
ð1Þ
where t 1 ; …; t n are a sequence of independent observations, and β1 , β2 , and α are unknown covariate coefficient vectors. The shape parameter k 4 0 is assumed to be known, but α may be unknown nuisance parameters. The null hypothesis H0: β1 ¼ β2 ¼ β0 leads to an irregular test, since the unknown change point ν vanishes under the null hypothesis, and in this case (1) becomes a Weibull model and the scale parameter is modeled by a regression model via a log link function. For the exponential case (k ¼1), Dupuy (2009) derived non-asymptotic bounds for significance and type II error rates for the LRT, assuming that all coefficients change and data are subject to random right censoring. If semi-parametric approaches using the partial likelihood are used, some computational complexity arising from the baseline hazard λ0 ðt i Þ ¼ t k−1 can be i avoided. In this framework, there are several results for the case when β is a scalar: Using Davies' approximation (Davies, 1987). O'Quigley and Natarajan (2004) obtained a score test for a case equivalent to when α ¼ 0. Liang et al. (1990) proved the weak convergence of the score test process to an Ornstein–Uhlenbeck process. Liu et al. (2008) noted that some of the assumptions for these results are difficult to assess in practice and instead developed a Monte-Carlo method to generate the distribution for the score test statistic. They argued that results can be extended to multiple change points (although the changing coefficient β is still a scalar) and discussed one and two change points in a simulation study. In the next section we establish the weak convergence of the LRT to the supremum of a quadratic process with Gaussian components. Also we discuss how to compute critical values of the test. In Section 3, we illustrate our method using the Chronic Granulomatous Disease (CGD) data as in Fleming and Harrington (1991). We compare the performance of our asymptotic test to the use of Davies' approximation (Davies, 1987). In general, our test produces either comparable or better results. While Davies' approximation can be only applied to a χ 2 process, our results can be applicable to a wider class of quadratic processes where the χ 2 process is a special case. In Section 4, we explore the behavior of our tests under different alternatives. Here we note the potential for model selection by using multiple tests to indicate which parameter(s) may have changed. In Section 5 we briefly discuss our findings. 2. Main result We consider observing the time-to-event data under staggered entry and type I censoring. Staggered entry refers to situations when subjects enter a clinical trial at random times 0 ¼ τ0 o ⋯ o τN . We assume that the entry times τ1 o ⋯ o τn follow a Poisson process. Type I censoring occurs when event times are recorded if occurring within a time interval ½0; T R for some fixed T R 4 0, and censored otherwise. Such situations can occur not just in clinical trials, but also in reliability testing scenarios as well as in observational studies where the entry times cannot be arranged ahead of the study. Assume also that the ending time of the study TR is independent of the number of subjects or events. Conditioning on NðT R Þ ¼ n, the total number of subjects who entered the study before time TR, the censoring time of each subject has a uniform distribution, i.e., ðT R −τi Þ∼Uð0; T R Þ, where τi is the entering time for the ith subject. Suppose the underlying data t follow model (1), but we observe the pair fx; δg, where x ¼ minðt; T R −τÞ and δ is the indicator 1ft o T R −τg . The likelihood for an observation fx; δg is Lðθ; x; δÞ ¼ ðxk−1 expðzT β1 þ wT αÞÞδ exp½−xk expðzT β1 þ wT αÞ=k1fx o νg þ ðxk−1 expðzT β2 þ wT αÞÞδ exp ½−ðexpðzT β1 Þνk þexpðzT β2 Þðxk −νk ÞÞexpðwT αÞ=k1fx≥νg Under the null hypothesis and the uniform distribution of censoring times, f(x) the marginal density for x and f~ ðxÞ the joint density for fx; δ ¼ 1g (both conditioned on covariate vectors z; w) are the following: f ðxÞ ¼ expð−xk η=kÞð1 þ xk−1 ðT R −xÞη=kÞÞ=T R f~ ðxÞ ¼ expð−xk η=kÞðxk−1 ðT −xÞη=kÞÞ=T R
R
ð2Þ
where η ¼ expðz β þ w αÞ. a:s P d Let ⟶, ⟶, ⟶, and ⟹ represent almost sure convergence, convergence in probability, convergence in distribution, and weak convergence of a process, respectively. Let lν ðθÞi ¼ lν ðθ; xi ; δi ; zi ; wi Þ be the log-likelihood of the parameter vector θ ¼ ðβ1 ; β2 ; αÞT for model (1) with an observation fxi ; δi g, covariate vectors zi and wi and a given change point ν. Now define _l ν ðθ Þ and €l ν ðθ Þ as the first and second derivatives of the log-likelihood with respect to θ evaluated at θ ¼ θ . The vector θ 0 i 0 i 0 0 takes the null hypothesis values ðβ0 ; β0 ; α0 ÞT . Now consider fx; δg the paired vectors of observations fx1 ; δ1 g; …; fxn ; δn g with corresponding covariate matrices Z ¼ ½z′1 ; …; z′n ′ and W ¼ ½w′1 ; …; w′n ′. We reparameterize the log-likelihood lν ðθÞ ¼ ∑ni¼ 1 lν ðθÞi and approximate with a Taylor expansion about the null values θ0 : T
T
n pffiffiffi lν ðθÞ ¼ lν ðθ0 þ h= nÞ ¼ ∑ lν ðθ0 Þi þ W n ðθ0 ; hÞ i¼1
iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
3
where T
T
n h n h n € ∑ l ν ðθ0 Þi h þ ∑ Ri ðθ0 ; hÞ W n ðθ0 ; hÞ ¼ pffiffiffi ∑ _l ν ðθ0 Þi þ 2n i ¼ 1 ni¼1 i¼1
where the h vector is now the parameter of interest and Ri ðθ0 ; hÞ is the remainder term. The log-likelihood ratio test statistic (given ν) is now ! 2Λn ðνÞ ¼ −2
sup W n ðθ0 ; hÞ− sup W n ðθ0 ; hÞ
fh∈H0 g
fh∈Hg
where H0 is the restricted support of h under the null hypothesis and H is the full support. Define C ij ¼ Covð_l νi ðθ0 Þ; _l νj ðθ0 ÞÞ with 0 o a≤νi ≤νj ≤bo T R . The design parameters a and b specify the range of a change point. These are needed to avoid the boundary problem. For νi given, the Fisher information under the null hypothesis is −1 I θ0 ðνi Þ ¼ C ii . Define I β0 ðνi Þ as the subset of the Fisher information applying only to β1 and β2 . Now define Σ ij ¼ I −1 θ0 ðνi ÞC ij I θ0 ðνj Þ. Let X ¼ ðXβ1 ; Xβ2 ; Xα ÞT be a multivariate Gaussian process with mean zero and covariance Σ ij . Let Xβ ¼ ðXβ1 ; Xβ2 ÞT . Finally let I d;2;1 ¼ ½I d ; I d T where Id is the identity matrix of dimension d, the length of β0 . Theorem 1. For model (1) when H 0 : β1 ¼ β2 is true: for ν∈½a; b and 0 o a o bo T R , 2Λn ðνÞ converges weakly to the quadratic process XTβ QXβ ðνÞ, where Q ¼ I β0 −I β0 I d;2;1 ðI Td;2;1 I β0 I d;2;1 Þ−1 I Td;2;1 I β0 Proof. First we note that model (1) satisfies the conditions for local asymptotic normality (van der Vaart, 1998, p. 104). See Appendix A for details. Then for each fixed ν we have convergence in distribution (van der Vaart, 1998, p. 232). That is, d
2Λn ðνÞ⟶2ΛðνÞ where 2ΛðνÞ ¼ inf ðX−hÞT I θ0 ðX−hÞ− inf ðX−hÞT I θ0 ðX−hÞ fh∈H0 g
fh∈Hg
^ T I θ ðX−hÞ; ^ ¼ ðX−h^ 0 ÞT I θ0 ðX−h^ 0 Þ−ðX−hÞ 0
ð3Þ
pffiffiffi n _ d ^ nð∑i ¼ 1 l ν ðθ0 Þi =nÞ⟶I θ0 X. Details for solving h∈H and h^ 0 ∈H0 are provided in Appendix B. We then use and X∼Nð0; I −1 θ0 Þ, and empirical measure (Appendix A) to finally establish the weak convergence 2Λn ðνÞ⟹2ΛðνÞ for all ν∈½a; b. □ Corollary 2. Assume α in (1) is known. Then XTβ QXβ ðνÞ is a χ 2d process with d degrees of freedom. T 2 Proof. If α is known then X ¼ Xβ and I θ0 ¼ I β0 . Then Q I −1 θ0 is idempotent, which implies for each ν∈½a; b: Xβ QXβ is a χ d random variable with d degrees of freedom (Graybill, 1976).
Note that if α ¼ 0, then all covariate coefficients change after ν. Corollary 3. If β1 and β2 are scalars, then XTβ QXβ ðνÞ ¼ Z 2 ðνÞ where ZðνÞ is a mean zero Gaussian process with covariance ðfΣ ij g11 þ fΣ ij g22 Þ−ðfΣ ij g12 þ fΣ ij g21 Þ CovðZðνi Þ; Zðνj ÞÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 11 12 21 11 22 12 21 ððI i þ I 22 i Þ−ðI i þ I i ÞÞððI j þ I j Þ−ðI j þ I j ÞÞ
ð4Þ
−1 where I lk i ¼ fI β0 ðνi Þ glk for l; k ¼ 1; 2 where fglk refers to the elements with rows corresponding to βl and columns to βk .
Proof. The case of β1 and β2 as scalars is comparable to the results from Williams and Kim (2011). That is, the matrix I d;2;1 reduces to 1 ¼ ð1; 1ÞT . Q is a 2 2 matrix with 1T Q ¼ 0. Furthermore fQ g11 ¼ jI β0 j=1T I β0 1 ¼ 1=ððI 11 þ I 22 Þ−ðI 12 þ I 21 ÞÞ, where 2 −1 2 I lk i ¼ fI β0 ðνi Þ gl;k and jI β0 j is the determinant of the matrix I β0 . Then 2ΛðνÞ ¼ fQ g11 ðX β1 −X β2 Þ ¼ Z ðνÞ. □ Note that if β1 and β2 are scalars and α is known, then VarðZðνi Þ; Zðνj ÞÞ ¼ 1 and Z 2 ðνÞ is a χ 21 process. Corollary 4. For model (1) when H0 is true, sup f2Λn ðνÞ : a≤ν≤bg converges in distribution to supfXTβ QXβ ðνÞ : a≤ν≤bg where XTβ QXβ ðνÞ is as in Theorem 1. The significance level αn for a critical value c is then " # αn ¼ P sup XTβ QXβ ðνÞ 4 c : a≤ν≤b
Proof. It is clear by direct application of Slutsky's Theorem and the Continuous Mapping Theorem.
□
3. Illustration: chronic granulomatous disease data To illustrate our asymptotic test, we consider the chronic granulomatous disease (CGD) data from Fleming and Harrington (1991). CGD is a hereditary disease of the immune system which leads to serious, often recurring infections. We only consider the first infection after the treatment either by placebo (control group) or gamma interferon (treatment iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
4
group). In this clinical trial, patients entered the study at different start times (staggered entry). Responses are time to first infection after treatment or the right censoring time (rescaled to TR ¼1). There are 20 uncensored events out of 53 for the control, and only 7 of 55 for the treatment are uncensored. Therefore, heavy censoring occurred to both groups. Age and sex were included as covariates. From exploratory analysis (Williams and Kim, 2013), k¼2 was determined to be a reasonable assumption so it is assumed to be known. Additionally, the survival curves for the control group indicate a possible change point, while those for the treatment group do not. To test this formally, we assume (1) has the form ( hðt i ; zi ; wi Þ ¼
t i expðControli β1 þ ð1; Agei ; Sexi ÞT αÞ; t i expðControli β2 þ ð1; Agei ; Sexi ÞT αÞ;
ti o ν t i ≥ν
ð5Þ
where Controli ¼1 for patients receiving the placebo and Controli ¼ 0 for the treatment. This slight reversal of indicators (usually control is indicated by 0) allows for the baseline of the control curve to change without affecting the baseline treatment curve. We use the maximum likelihood estimate θ^ 0 and an interval ½a; b based on the 80% coverage of the control group data with a¼0.057 as the 10% quantile and b¼0.748 as the 90% quantile. To generate the critical values we create a grid of 1000 ν values spread evenly from a to b. We then generate realizations of the corresponding quadratic processes XTβ QXβ ðνÞ ¼ Z 2 ðνÞ (by Corollary 3) substituting θ^ 0 for θ0 (See Appendix B for details of constructing CovðXðνi Þ; Xðνj ÞÞÞ. We generate 10,000 realizations and find the maximum of each over the grid of ν. The upper 1% quantile of these maximum is then an estimate of the 1% critical value of the test. In this case, the observed test statistic of 37.79 exceeds the estimated 1% critical value of 13.8, so the null hypothesis is rejected at 1% significance level. The estimated change point for the control group is ν^ ¼ 0:072, or about 23 days after the placebo treatment. To assess type I error and power, data were simulated based on the maximum likelihood estimates under the null β^ 0 ¼ 1:54 and α^ 0 ¼ −0:15; −0:43; −0:45 and alternative hypotheses β^ 1 ¼ 4:56, β^ 2 ¼ 1:04, and α^ ¼ −0:27; −0:39; −0:28. Covariates z (control) and w (intercept, age, and sex) were resampled from the 108 observations in the CGD data. Likelihood ratio test statistics were calculated and compared to the critical values obtained from Corollary 3. This was repeated 1000 times for sample sizes of 100, 200, and 500. Type I error was calculated as the percentage of times that the test statistic exceeded the 5% and 1% critical values when the data was generated under the null hypothesis. Power was calculated as the percentage of times that the test statistic exceeded the 5% and 1% critical values when the data was generated under the alternative hypothesis. The rates obtained by using Corollary 3 are compared to rates from applying Davies' approximation (Davies, 1987), which allows us to approximate p-values for the supremum of χ 2 processes. While the result of Davies only strictly applies to χ 2 processes, it is possible that Z2 can be approximated by a χ 21 process if VarðZðνÞÞ≈1, as it is for the CGD example. The results for the Type I error simulations are summarized in Table 1. We see good convergence at a sample size of 100. Rates for our method (LAN) are comparable to those for Davies. Furthermore, our method is within a standard error of the nominal levels at each sample size. Interestingly, the test using Davies' approximation appears to be conservative for the Table 1 Type I error simulations using CGD data, comparing results from Corollary 3 (LAN) and Davies (1987) for 1000 replications with ½a; b ¼ ½0:057; 0:748. Test
Sample size
s.e. 5% s.e. 1% LAN 5% LAN 1% Davies 5% Davies 1%
100
200
500
0.022 0.010 0.047 0.012 0.053 0.019
0.015 0.007 0.035 0.007 0.038 0.010
0.010 0.004 0.043 0.006 0.020 0.003
Table 2 Power simulations using CGD data, comparing results from Corollary 3 (LAN) and Davies (1987) for 1000 replications with ½a; b ¼ ½0:057; 0:748. Test
LAN 5% LAN 1% Davies 5% Davies 1%
Sample size 100
200
500
0.950 0.895 0.961 0.911
0.998 0.998 0.998 0.998
1.000 1.000 1.000 1.000
iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
5
sample size n ¼500. Power simulations based on θ^ are displayed in Table 2. Power is high for both methods, and it is essentially the same. The close to nominal type I error combined with high power demonstrates that our test using both the LAN and Davies methods is appropriate and useful for this CGD data set. 4. Simulation study: exponential hazard regression We consider a constant hazard function with covariate information zi ∼Nð0; 1Þ: hðt i ; zi Þ ¼ expðβ0 þ β1 zi Þ
ð6Þ
We choose the simple exponential case ðk ¼ 1Þ, because CovðXðνi Þ; Xðνj ÞÞ is greatly simplified (Appendix B) and consistent estimators for parameters under some change point alternatives have already been studied (Dupuy, 2006). Based on our results from Section 2, we can construct tests for three different change point alternatives. 1. The intercept β0 changes at unknown ν : β01 ≠β02 2. The slope β1 changes at unknown ν : β11 ≠β12 3. Both the slope and intercept change at the same ν : β01 ≠β02 and β11 ≠β12 Letting β0 ¼ 1, β1 ¼ 1, and TR ¼1; we examine the behavior of these three tests. Under the null hypothesis of no change, we observe approximately 40% censoring. We choose the change point interval to be the expected 25% and 75% quantiles of the data: a¼0.06 and b¼0.37 respectively. To obtain critical values, we generate 10,000 realizations of each of the Z2 process for the single change tests and XTβ QXβ for the vector change. Each realization occurs over a grid of 1000 values of ν∈½a; b. The upper 5% and 1% quantiles from the supremum of each realization are then selected as the critical values. For example, we obtain 9.43, 8.77, and 11.33 as the 5% critical values and 13.20, 12.33, and 15.50 as the 1% for the test for changes in intercept, slope, and both, respectively. Even at sample size 100, the empirical type I error rates approach their nominal values. Interestingly results using the Davies (1987) method are quite conservative, resulting in empirical type I error rates less than 0.01 for a 5% level test. We next consider the behavior of the tests under the three alternatives. Table 3 displays the power for each of the three tests under all three alternatives. Only results for positive changes in β using 5% critical values are included. Results are similar for negative changes in β and for 1% critical values. Not surprisingly, a given test turns in the best result in terms of power when its specific alternative is true. On the other hand, the joint test retains good power when only one coefficient changes. When both parameters change, it appears that all of the three tests require much larger sample size. This is just an artifact of choosing a smaller effect size (1.5 instead of 2) so as not to have the power be close to 1.000 for both the 200 and 500 sample sizes. Together these results suggest that evaluating all three tests (perhaps with some correction for multiple testing) may be useful for model selection when determining not only if a change occurred but also which regression parameters are involved. We might interpret the joint test as a global test with reasonable power over a range of alternatives and consider the single test as follow-up with higher power for more specific alternatives. 5. Discussion In this paper we have proposed a likelihood-ratio test for an unknown change point in hazard regression models where one or more covariate coefficients change abruptly. The test allows for staggered entry and type I censoring so it is useful for Table 3 Power for single and joint change alternatives to model (6) using 5% critical values. Each scenario was simulated 1000 times for a standard error o 0:016. Test
Sample size 100
200
500
H A : β02 ¼ 2 Intercept change Covariate change Both change
0.904 0.060 0.875
0.999 0.099 0.994
1.000 0.225 1.000
H A : β12 ¼ 2 Intercept change Covariate change Both change
0.068 0.638 0.573
0.124 0.929 0.879
0.244 1.000 0.999
H A : β02 ¼ β12 ¼ 1:5 Intercept change Covariate change Both change
0.217 0.221 0.331
0.493 0.513 0.719
0.958 0.942 0.993
iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
6
situations where the entry times of time-to-event data are random, as shown in the study of CGD data in Section 3. We have derived the asymptotic distribution of the test statistic so time-consuming calculations can be avoided. Our procedure can test a change in a single parameter or multiple parameters in the covariate coefficient vector, it is possible to develop a battery of tests that are useful for model selection purposes. The LAN based approach combined with the weak convergence of the empirical measure enables us to approximate the profile log-likelihood ratio test statistic by a quadratic process, which greatly extends the use of tests when the limit is not a χ 2 process, as dictated in Davies' approximation. While we have focused on developing the likelihood-ratio test, our weak convergence results can also be used to build score and score-related tests. Our derivations and proofs use all the components needed for a score test (the score function and its covariance), so future use of score tests for model (1) can be facilitated by these results. Critical values for the test can be obtained in several ways, depending on the particular application. The most general method is to generate realizations of the multivariate Gaussian processes Xβ and then calculate the quadratic form XTβ QXβ ðνÞ from Theorem 1 and then use Corollary 4. If a test is desired to test only one coefficient in the covariate vector, we only need to generate the univariate Gaussian process ZðνÞ and use Corollary 3 and Corollary 4. Finally, if the covariate information is not used, the current nonstationary Gaussian process can be transformed to a stationary Ornstein–Uhlenbeck process (Kim et al., 2004; Williams and Kim, 2013) for which tail approximations for the supremum of the process are available (Leadbetter et al., 1983). For the example discussed in Section 4, we used four normal variates for each value of ν corresponding to two parameters before and two after the change. If more covariates are considered for β, only the dimension of the covariance matrix increases as a result; the test procedure remains identical. On the practical side, however, very large number of covariates in β may cause computational problems. Depending on the number of covariates and computing capability, it may require reduction of the number of grid points at which the covariance matrix is evaluated, or selection of a subset of covariates to be entered into β in the model. This work has focused on a change in the linear predictor of the hazard function when the Weibull shape parameter k is assumed to be known. We conducted computer simulations for the cases where the nuisance parameter k is unknown. We found that the maximum likelihood estimate for the unknown shape parameter is numerically unstable in some cases. Thus, further study is needed to get around the problem to obtain a stable estimate for the shape parameter. On the other hand, testing for a change from one known value of k to another, unknown value should be possible. For example, Williams and Kim (2011) consider the test of no change versus alternative hypothesis where k¼1 before the unknown change point, and k assumes some unknown value larger than 1 after the change. These results did not incorporate censoring or covariates, but they do consider tests under various forms of continuous hazard function. Focusing on abrupt change and using the techniques presented in this work should allow for the establishment of corresponding convergence results for likelihood-ratio tests.
Acknowledgments The first author was sponsored through the Institute for Critical Technology and Applied Science (ICTAS) Doctoral Scholars Program at Virginia Tech. The work of the second author was partially supported as part of the National Science Foundation Biocomplexity of Coupled Human and Natural Systems Program, Award No. BCS-0709671. Appendix A. Proof of weak convergence For the independent and identically distributed case (for example, see Williams and Kim, 2011), we would consider the likelihood as a density pθ ðx; δÞ ¼ Lðθ; x; δÞ. However, covariate information leads to nonidentical distributions. To circumvent this, we treat the covariates as random variables and consider the joint distribution (see the generalized linear model example in van der Vaart, 1998, p. 203): pθ ðx; δ; z; wÞ ¼ f ðx; δjz; wÞgðz; wÞ ¼ LðθÞL′ðγÞ where LðθÞ is the likelihood component for θ and L′ðγÞ is the likelihood component for the unrelated nuisance parameter vector γ. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi For our model pθ ðx; δ; z; wÞ is continuously differentiable with respect to θ for every ðx; δÞ, z, w, and ν. The Fisher Information matrix component for θ is I θ ¼ Ez;w ðEx;δjz;w ð€l ν ðθÞÞÞ: For model (1) under staggered entry, Ex;δjz;w ðl€ν ðθÞÞ is well-defined and continuous with respect to θ for any ν. Then for a compatible marginal distribution gðz; wÞ, I θ is also well-defined and continuous with respect to θ for any ν. We will assume z; w have such a compatiblep distribution ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiand will approximate it using the empirical distribution of the observed covariates. By the above properties of pθ ðx; δ; z; wÞ and I θ , our experiment is differentiable in quadratic mean (van der Vaart, 1998, p. 95). Therefore local asymptotic normality holds for each ν (van der Vaart, 1998, p. 104). For (1) we show that all elements of the vector-valued function _l ν ðθ0 Þ are in one Donsker class (van der Vaart, 1998, p. 269) indexed by ν. We construct an envelope function Mðx; δ; z; wÞ such that Mðx; δ; z; wÞ≥∥∂_l ν ðθ0 Þ=∂ν∥1 for ν∈½a; b and iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
7
EðjMðx; δ; z; wÞjÞ o ∞. The norm ∥ ∥1 is the sum of the absolute value of the elements of ∂l_ν ðθ0 Þ=∂ν. Note that Mðx; δ; z; wÞ is not a function of ν. The existence of this envelope function satisfies the conditions for a Donsker class (van der Vaart, 1998, p. 271). Based on the derivatives (See Appendix C) we choose Mðx; δ; z; wÞ ¼ jzexpðzT β0 þ wT α0 Þðνn Þk−1 j where νn maximizes νk−1 for fixed k 40 and ν∈½a; b. If 0 ok o1, then νn ¼ a. Otherwise νn ¼ b. We note that EðjMðx; δ; z; wÞjÞ ¼ Ez;w ðjMðx; δ; z; wÞjÞ since M is constant with respect to x and δ. If we can assume that the distribution gðz; wÞ is such that Ez;w ðjzexpðzT β0 þ wT α0 Þðνn Þk−1 jÞ o∞, then EðjMðx; δ; z; wÞjÞ o ∞. By our construction ! pffiffiffi 1 n _ ∑ l ν ðθ0 Þ ⟹I θ0 ðνÞXðνÞ; ν∈½a; b n ni¼1 where XðνÞ is a zero-mean Gaussian process with CovðXðνi Þ; Xðνj ÞÞ ¼ Σ ij .
For our model, we would like to show that each element in the matrix of functions €l ν ðθ0 Þ is Glivenko–Cantelli classes (van der Vaart, 1998, p. 269) indexed by ν∈½a; b. We construct an envelope function Fðx; δ; z; wÞ, such that Fðx; δ; z; wÞ≥∥€l ν ðθ Þ∥ 0
1
and EðjFðx; δ; z; wÞjÞ o ∞ for ν∈½a; b. This process is very similar to that of constructing Mðx; δ; z; wÞ. The existence of this envelope function satisfies the conditions for a Glivenko–Cantelli class (van der Vaart, 1998, p. 272). Alternatively, we can show membership of l€ν ðθ Þ in a Donsker class directly, which implies membership in a Glivenko–Cantelli class. For this we 0
consider ∥∂€l ν ðθ0 Þ=∂ν∥1 with envelope function M 2 ðx; δ; z; wÞ ¼ ð∥zzT ∥1 þ ∥2zwT ∥1 þ ∥wwT ∥1 ÞexpðzT β0 þ wT α0 Þðνn Þk−1 where νn maximizes νk−1 for fixed k 4 0 and ν∈½a; b. If 0 o k o 1, then νn ¼ a. Otherwise νn ¼ b. Again, M2 is a constant with respect to x and δ, so as long as the marginal density gðz; wÞ is compatible, we have established a Donsker class. Then P ∥n−1 ∑n €l ν ðθ Þ−I θ ∥ ⟶0. sup ν∈½a;b
i¼1
0
0
1
We next consider the remainder T Ri ðθ0 ; hÞ ¼ R1i ðθ0 ; hÞ−h €l ν ðθ0 Þh=2n
where pffiffiffi pffiffiffi T R1i ðθ0 ; hÞ ¼ lν ðθ0 þ h= nÞ−lν ðθ0 Þ−h _l ν ðθ0 Þ= n: pffiffiffi T T According to Taylor's Theorem: R1i ðθ0 ; hÞ ¼ h €l ν ðθ0 þ hϕ = nÞh=2n, where hϕ ¼ ðϕ1 h1 ; …; ϕd hd Þ for some ϕ1 ; …; ϕd ∈½0; 1 with d pffiffiffi the dimension of θ. For fixed θ0 and h and for n large enough, ½a; b ½θ0 ; θ0 þ h= n is well-defined. Furthermore, each pffiffiffi θ∈½θ0 ; θ0 þ h= n has a corresponding envelope M 2 ðθÞ which is an exponential function of a linear combination of θ. Thus pffiffiffi pffiffiffi there exists an envelope M2 ðθn Þ with θn ∈½θ0 ; θ0 þ h= n which dominates. Then €l ν ðθ0 þ hϕ = nÞ is Donsker (and Glivenko– pffiffiffi pffiffiffi Cantelli) since ½θ0 ; θ0 þ hϕ = n⊂½θ0 ; θ0 þ h= n. We can construct an appropriate envelope function for the difference €l ν ðθ0 þ pffiffiffi € P hϕ = nÞ−l ν ðθ0 Þ using the triangle inequality, so it is also Glivenko-Cantelli. Thus supν∈½a;b j∑ni¼ 1 Ri ðθ0 ; hÞj⟶0. We note here that our envelopes Mðx; δ; z; wÞ and M2 ðx; δ; z; wÞ are not functions of either x or δ. This suggests that any well-defined (local asymptotic normality holds for each ν) right-censoring regime for model (1) will possess the weak convergence properties demonstrated above. Of course the details of Cij will change, since f(x) and f~ ðxÞ from (2) depend on censoring structure, but we expect the results from Theorem 1 to hold for broader classes of censoring. One example is that of simple right censoring at time TR. It is a straight-forward exercise to show that the results of the previous sections apply to this simple case as well. Appendix B. Specifying the covariance structure To characterize the 2ΛðνÞ process, we need −1 CovðXðνi Þ; Xðνj ÞÞ ¼ I −1 θ0 ðνi ÞC ij I θ0 ðνj Þ;
where C ij ¼ Covð_l νi ðθ0 Þ; _l νj ðθ0 ÞÞ. From local asymptotic normality (LAN): Ex;δjz;w ð_l ν ðθ0 ÞÞ ¼ 0. Then C ij ¼ Ez;w ðEx;δjz;w ð_l νi ðθ0 Þ_l νj ðθ0 ÞÞÞ. The inner expectation (which we name Cijl) can be evaluated directly through integration using (2), the distributions of fxl ; δl g conditional on zl ; wl for observation l ¼ 1; …; n. The integrands change depending on the domain of x∈½0; νi Þ; ½νi ; νj Þ; ½νj ; T R . The outer expectation can be approximated by ð1=nÞ∑nl¼ 1 C ijl , the expectation with respect to empirical distribution of the covariates zl ; wl . The following are the covariance equations for Cijl conditioned on the covariate values zl ; wl : 2 3 C β1 ;β1 C β1 ;β2 C β1 ;α 6 C β ;β C β ;β C β ;α 7 C ijl ¼ 4 2 1 5 2 2 2 C α;β1 C α;β2 C α;α iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
8
with C β1 ;β1 ¼ zl zTl ½1 þ expð−νki ηl =kÞðνi =T R −1Þ þ ðΓ½1=k; νki ηl =k−Γ½1=kÞ=ððηl =kÞ1=k kT R Þ C β1 ;β2 ¼ 0 C β1 ;α ¼ zl wTl ½1 þ expð−νkj ηl =kÞðνj =T R −1Þ þ ðΓ½1=k; νkj ηl =k−Γ½1=kÞ=ððηl =kÞ1=k kT R Þ C β2 ;β1 ¼ zl zTl ½expð−νkj ηl =kÞðνj =T R −1Þ−expð−νki ηl =kÞðνi =T R −1Þ þ ðΓ½1=k; νkj ηl =k−Γ½1=k; νkj ηl =kÞ=ððηl =kÞ1=k kT R Þ C β2 ;β2 ¼ zl zTl ½−expð−νkj ηl =kÞðνj =T R −1Þ þ ðΓ½1=k; T kR ηl =k−Γ½1=k; νkj ηl =kÞ=ððηl =kÞ1=k kT R Þ C β2 ;α ¼ zl wTl ½−expð−νki ηl =kÞðνi =T R −1Þ þ ðΓ½1=k; T kR ηl =k−Γ½1=k; νki ηl =kÞ=ððηl =kÞ1=k kT R Þ C α;β1 ¼ wl zTl ½1 þ expð−νki ηl =kÞðνi =T R −1Þ þ ðΓ½1=k; νki ηl =k−Γ½1=kÞ=ððηl =kÞ1=k kT R Þ C α;β2 ¼ wl zTl ½−expð−νkj ηl =kÞðνj =T R −1Þ þ ðΓ½1=k; T kR ηl =k−Γ½1=k; νkj ηl =kÞ=ððηl =kÞ1=k kT R Þ C α;α ¼ wl zTl ½1−expð−T kR ηl =kÞ þ ðΓ½1=k; T kR ηl =k−Γ½1=kÞ=ððηl =kÞ1=k kT R Þ where ηl ¼ expðzTl β0 þ wTl α0 Þ, Γ½z ¼ function.
R∞ 0
t z−1 e−t dt is the gamma function, and Γ½z; a ¼
R∞ a
t z−1 e−t dt is the incomplete gamma
Now that we have Cij and therefore I θ0 ðνÞ, we solve for h^ 0 and h^ from (3). Under the alternative h is unrestricted, therefore h^ ¼ X making the second term in (3) equal to 0. Under the null hypothesis, h0 is restricted to be ðh00 ; h00 ; hα ÞT , Þ−1 I T I β Xβ . Then h^ ¼ ðh^ ; h^ ; where the third component is unrestricted. From linear algebra: h^ ¼ ðI T I β I 00
d;2;1
0
d;2;1
d;2;1
0
0
00
00
Xα Þ ¼ ðh^ 00 I d;2;1 ; Xα ÞT . Plugging in h^ 00 into the left term in (3) and factoring out Xβ , we get XTβ QXβ as in Theorem 1. It is T
clear that I Td;2;1 Q ¼ 0, making Q rank d. Appendix C. Derivatives of the log-likelihood We include the first and second derivatives evaluated at θ ¼ θ0 . From these, we can construct envelope functions Mðx; δ; z; wÞ, Fðx; δ; z; wÞ, and M2 ðx; δ; z; wÞ and calculate Cijl via integration. Let η ¼ expðzT β0 þ wT α0 Þ. _l ν ðθ Þj 0 β1 ¼ β2 ¼ β0 ;α ¼ α0 ∂lν ðθ0 Þ=∂β1 ¼ zðδ−ηxk =kÞ1fx o νg þ zð−ηνk =kÞ1fx≥νg ∂lν ðθ0 Þ=∂β2 ¼ zðδ−ηðxk −νk Þ=kÞ1fx≥νg ∂lν ðθ0 Þ=∂α ¼ wðδ−ηxk =kÞ ∂_l ν ðθ0 Þ=∂νjβ1 ¼ β2 ¼ β0 ;α ¼ α0 ∂2 lν ðθ0 Þ=∂ν∂β1 ¼ −zηνk−1 1fx≥νg ∂2 lν ðθ0 Þ=∂ν∂β2 ¼ zηνk−1 1fx≥νg ∂2 lν ðθ0 Þ=∂ν∂α ¼ 0 €l ν ðθ Þj 0 β1 ¼ β2 ¼ β0 ;α ¼ α0 ∂2 lν ðθ0 Þ=∂βT1 ∂β1 ¼ zzT ð−ηxk =kÞ1fx o νg þ zzT ð−ηνk =kÞ1fx≥νg ∂2 lν ðθ0 Þ=∂βT2 ∂β1 ¼ 0
∂2 lν ðθ0 Þ=∂αT ∂β1 ¼ zwT −ηxk =k 1fx o νg þ zwT ð−ηνk =kÞ1fx≥νg ∂2 lν ðθ0 Þ=∂βT2 ∂β2 ¼ zzT ð−ηðxk −νk Þ=kÞ1fx≥νg ∂2 lν ðθ0 Þ=∂αT ∂β2 ¼ wzT ð−ηðxk −νk Þ=kÞ1fx≥νg ∂2 lν ðθ0 Þ=∂αT ∂α ¼ −wwT ηxk =k ∂€l ν ðθ0 Þ=∂νjβ1 ¼ β2 ¼ β0 ;α ¼ α0 ∂3 lν ðθ0 Þ=∂ν∂βT1 ∂β1 ¼ zzT ð−ηνk−1 Þ1fx≥νg ∂3 lν ðθ0 Þ=∂ν∂βT2 ∂β1 ¼ 0 ∂3 lν ðθ0 Þ=∂ν∂αT ∂β1 ¼ zwT ð−ηνk−1 Þ1fx≥νg ∂3 lν ðθ0 Þ=∂ν∂βT2 ∂β2 ¼ zzT ðηνk−1 Þ1fx≥νg iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006
M.R. Williams, D.-Y. Kim / Journal of Statistical Planning and Inference ] (]]]]) ]]]–]]]
9
∂3 lν ðθ0 Þ=∂ν∂αT ∂β2 ¼ wzT ðηνk−1 Þ1fx≥νg ∂3 lν ðθ0 Þ=∂ν∂αT ∂α ¼ 0
References Davies, R.B., 1987. Hypothesis testing when a nuisance parameter is present only under the alternatives. Biometrika 74, 33–43. Dupuy, J.-F., 2006. Estimation in a change-point hazard regression model. Statistics and Probability Letters 76, 182–190. Dupuy, J.-F., 2009. Detecting change in a hazard regression model with right-censoring. Journal of Statistical Planning and Inference 139, 1578–1586. Fleming, T.R., Harrington, D.P., 1991. Counting Processes and Survival Analysis. John Wiley Sons. Graybill, F.A., 1976. Theory and Application of the Linear Model. Duxbury Press. Kim, D., Woodroofe, M., Wu, Y., 2004. Testing for a change in the hazard rate with staggered entry. Communications in Statistics—Theory and Methods 33, 2041–2058. Invited paper for the special issue in honor of Professor Z. Govindarajulu on his 70th birthday. Leadbetter, M.R., Lindgren, G., Rootzén, H., 1983. Extremes and Related Properties of Random Sequences and Processes. Springer, New York. Liang, K.-Y., Self, S.G., Liu, X., 1990. The Cox proportional hazards model with change point: an epidemiologic application. Biometrics 46, 783–793. Liu, M., Lu, W., Shao, Y., 2008. A Monte Carlo approach for change-point detection in the cox proportional hazards model. Statistics in Medicine 27, 3894–3909. O'Quigley, J., Natarajan, L., 2004. Erosion of regression effect in a survival study. Biometrics 20, 344–351. van der Vaart, A.W., 1998. Asymptotic Statistics. Cambridge University Press. Williams, M.R., Kim, D., 2011. Likelihood ratio tests for continuous monotone hazards with unknown change point. Statistics and Probability Letters 81, 1599–1603. Williams, M.R., Kim, D., 2013. A test for an abrupt change in Weibull hazard functions with staggered entry and type I censoring. Communications in Statistics—Theory and Methods 42, 1922–1933.
iPlease cite this article as: Williams, M.R., Kim, D.-Y. A test for abrupt change in hazard regression models with Weibull baselines. Journal of Statistical Planning and Inference (2013), http://dx.doi.org/10.1016/j.jspi.2013.04.006