行政院國家科學委員會專題研究計畫 成果報告
非正比例涉險模式變動係數之估計與診斷
計畫類別： 個別型計畫
計畫編號： NSC93-2118-M-039-001-
執行期間： 93 年 08 月 01 日至 94 年 07 月 31 日
執行單位： 中國醫藥大學公共衛生學系
計畫主持人： 吳宏達
報告類型： 精簡報告
處理方式： 本計畫可公開查詢
中 華 民 國 94 年 9 月 29 日
Report of Grant NSC 93-2118-M-039-001
Hong-Dar Isaac Wu,
School of Public Health, China Medical University, 91 Hsueh-Shih Rd., Taichung 404, Taiwan.
E-mail: [email protected] September 28, 2005
Title: A Combined Modeling of Cure-Fraction and Multiple Cross-Effect
through A Varying-Coefficient Hazards Regression
Abstract
We consider a piecewise-constant varying-coefficient model for the analysis of survival data to deal with, simultaneously, multiple crossing-hazards phenomenon and the existence of cure prob-ability. Artificial examples for illustrating the use of the model and analysis of actual data are provided. Moreover, we construct a class of tests for a nested model structure. Under this structure, we can test the need of imposing varying coefficients for specific covariates (a special case: test for varying heteroscedasticity), the test for varying-coefficient PH setting, and the test for equality between the distributions of subpopulations with or without a cure fraction. Estimating procedures and simulation study are also reported.
Key Words: time-varying effect, heteroscedasticity, multiple crossings, proportional hazards, non-proportional hazards, cure fraction.
1
Introduction
In event-history data analysis where the effect of a specific variable is of the main interest, modeling and estimation of time-varying effect gets more important in the recent years. In contrast with the Cox proportional hazards model (Cox, 1972), many authors have devoted to the study of varying-coefficient Cox (PHvc) model. For example, see Murphy and Sen (1991), Murphy (1993), Martinussen and Scheike (2001), and Tien, Zucker, and Wei (2005) among oth-ers. Without regard to the space of ’time’, the PHvc model basically still estimate homogeneity effect over the space of covariate(s). By this, homogeneity means there is a common ’effect’ (be-tween two covariate-specific subpopulations) measured at different values of the covariate or of the ’configurations’ of several covariates (Wu and Hsieh, 2005). On the contrary, heterogeneity states that the effect is different and diverse over the covariate space for a fixed time point. As was studied by Hsieh (2001), the heteroscedastic hazards regression (HHR) model is a model suitable for capturing heterogeneity effect. Note that Hsieh’s HHR model can be derived from the transformation model with heteroscedastic errors. In addition to modeling heterogeneity, the HHR model also has the merit in modeling time-varying effect. The HHR model proposed in Hsieh (2001) has the following hazard function:
λ(t; z, x) = λ0(t)eγ Tx(t)
{Λ0(t)}e
γT x(t)−1
eβTz(t), (1)
where z(t) and x(t) are two sets of predictable time-dependent covariates, and Λ0(t) = Rt
0λ0(u)du
is the baseline cumulative hazard. In view of the intrinsic time-varying property of the hazard ra-tio (which may be a measure of effect of interest) implied by (1), it is possible to extend the HHR model to incorporate varying-coefficient settings at both of the ’heteroscedasticity component’ (exp(γTx(t))) and the ’risk function’ (exp(βTz(t))). Hereafter we denote the varying-coefficient HHR model as a HHRvc model. Motivation of this extension can be interpreted by the following two examples. First, the inter-relation among groups, in terms of survivor or cumulative hazard functions, may be ’diverse’ in time. A more apparent phenomenon is the multiple cross-effect (MCE) also studied by Bagdonaviˇcius and Nikulin (2005). Second, cure-fraction (CF) appeared in many clinical and oncological studies in which survivals of cancer patients receiving surgery followed by chemo- and/or radio-therapies are of concern. Though, the definition of ’cure’ still needs to be clarified.
The HHRvc model can deal with survival data with time-diversity (e.g., MCE) and cure-fraction simultaneously, at a loose definition. The purpose of this paper is to study the applica-tions of the HHRvc model, in particular in the aspects of model validity. Here we should note that with a more-extended model, some ’focused’ tests can be studied; including the test of
time-varying effect of specific covariate(s), the validity of PHvc model, the equality of distributions with or without underlying cure fraction, and so on.
The use of a more-extended model (say the current HHRvc model) than a smaller class (say the PH or PHvc models) should be strongly motivated, at least in several aspects. First, when the data cannot be suitably described by the simpler model; on the other hand, the extended model do fit the data well enough. Second, the extended model is easy to interpret; at least, it models the hazards ratio directly although the parameters β may not have parallel meaning with that of conventional PH model. Third, diagnostic plot(s) and model validity problem can be easily implemented with the model estimated.
By an inspection of the application of the HHRvc model and taking PH model as a com-parison basis, the most important parts are certainly (i) to make feasible the incorporation of parameter γ; (ii) to persuade the time-varying characteristic of β and γ; and (iii) to further investigate the goodness-of-fit problems. With respect to the previous explanations, we seek in this study either to test for a sub-model of HHRvc against a general alternative (i.e., with an omnibus flavour), or, simply to test the validity of the embedded sub-model (null hypothesis) by assuming the HHRvc model as the alternative hypothesis. For further elucidation, consider a class of models Pθ indexed by θ∈ Θ ⊂ Rp and a sub-model Pθ˜, ˜θ∈Θe ⊂ Θ. Further denotes
θk= θ\ ˜θ. We define in this paper two types of test: the type-O test, ’O’ stands for ’omnibus’,
and the type-N test, where N means that the null and alternative hypotheses actually form a nested structure. For the first type, the statistical hypotheses are: H0: the underlying model
is Pθ˜, against Ha: the underlying model is not Pθ˜. However, construction of the test statistic
relies on the property that Pθ˜⊂ Pθ (see Section 4).
For the second type, we assume the true model belongs to the class of Pθand H0∗: θ = ˜θtθk0
versus Ha∗ : θ = ˜θt θk, where θk0 is some prespecified value of the unknown θk. The type-N
hypotheses can be written even more simpler as H0∗ : θk = θk0 and Ha∗ : θk 6= θk0. Finally we
claim that although the varying-coefficient model we investigate is new, the employed tests are commonly seen but were not often discussed together under a complex class.
2
Illustration of The Piecewise Constant Model
Piecewise-constant setting
For a simpler exposition, we adopt notations only with ’univariate’ case. Model (1) can be extended to allow for varying coefficients:
λz,x(t) = λ0(t){Λ0(t)}e γ(t)z−1
eβ(t)z+γ(t)x. (2)
Because the partial likelihood does not eliminate out the baseline hazard, there are three time-dependent parameters to be estimated simulataneously, that is, Λ0(t) (or λ0(t)), β(t), and γ(t).
For the conventional varying-coefficient Cox model, there are at least two approaches to es-timation of the varying-coefficient(s): the piecewise-constant approximation method (Mur-phy and Sen, 1991; Mur(Mur-phy, 1993; Marzec and Marzec, 1997) and the smoothing estimate (Martinussen and Scheike, 2001; Tien at al., 2005). We use the former setting because it is convenient to bypass the work of Hsieh (2001), based on his OEE approach. Let [0, τ ] be the observational period, and τ0 < τ1 < . . . < τm be a set of cutoff points with τ0 = 0 and τm = τ .
We use the following piecewise-constant approximation to the three functions: Λ0(t) = Z t 0 m X 0 αj1(τj−1<u≤τj)du, β(t) = m X 0 βj1(τj−1<t≤τj), γ(t) = m X 0 γj1(τj−1<t≤τj). (3)
So in this paper we consider the piecewise-constant HHR model, called the HHRvc model, which has the following hazard and cumulative hazard:
λ(t; z) = αj{Λ0(t)}σj−1σjµj, τj−1 < t≤ τj,
Λ(t; z) = Λ(τj−1; z) + [{Λ0(t)}σj− {Λ0(τj−1)}σj]µj, τj−1 < t≤ τj, (4)
where λ(·) denotes the approximation of λ(·), σj = eγ T jz, µ j = eβ T jz, and Λ(τ 0; z) = Λ0(τ0) = 0.
Formula (4) is very useful to the understanding of HHRvc and to the forthcoming random number generation according to HHRvc, because we can simply use the relation S(·) = exp{Λ(·)} and equate it to a Uniform(0,1) random number.
Of course the main reason is due to the varying-coefficient heteroscedasticity component. (I) HHR without varying coefficient can only result in a one-time crossing. (II) The ordinary PH model with varying-coefficient risk function can make multiple crossings, but the inter-subpopulation effect is still homogeneous for any fixed time point. (III) HHR with varying-coeffcient risk function and heteroscedasticity component is a natural extension. It can offer a multiple-crossing modeling as well as including heterogeneity.
Example 2.1: The case of MCE modeling
Let’s consider the Weibull regression as an example for multiple crossings. With this parametric regression model, we take the covariates z and x to be totally equal (z = x). If we take the baseline cumulative hazard as Λz(t) = a(t)tb(t), with a(t) = eβ(t)z and b(t) = eγ(t)z. For
simplicity, moreover, we choose z = 0 and 1 to discuss the two-sample case. For z = 0 and 1, we have Λ0(t) = t and Λ1(t) = eβ(t)texp{γ(t)}, respectively. Now if we take β(t) = log 2,∀t, and
γ(t) = log 2, 0 < t≤ 1 2, = 0,1
2 < t≤ 1, = − log 2, 1 < t,
then the two survival (or cumulative-hazard) curves have two crossings (see Fig.1). With the same spirit, a case with more crossings also can be designed. It is obvious that multiple crossings can also come from β(t), when it fluctuates around the horizontal line β(t) = 0. Nevertheless, if the inter-relations among several (≥ 3) groups (obtained from suitably grouping a continuous z-variable) are diverse over time, the multiple cross-effect is more adequately captured by a time-varying γ(t) than by β(t). So this example illustrate a commonly encountered situation that the two survival curves cross at an early stage, and then separate each other, and cross (or tend to cross) again at a later stage. A more realistic example with actual data will be provided in Section 5.
[Figure 1 about here.]
3
Estimation under the Piecewise-constant Setting
Suppose there are n ordered, randomly right-censored observations (survival times or censored times) T1,T2, . . ., Tn. Let λi(t; Z, X) be the intensity process and Ni(t) be the counting process
for the ith individual and Yi(t) be the associated at-risk indicator at time t. Further denote
SJ(t) = (1/n) n X i=1
Yi(t)J (t)eβ(t)zi+γ(t)xi{Λ0(t)}eγ(t)xi−1, (5)
with time-dependent covariates J (t) = zi(t), xi(t), or vi(t)≡ xi(t){1 + eγ(t)xilog Λ0(t)}.
As was discussed in Wu and Hsieh (2005), we have the following Breslow-type estimating equation for the baseline cumulative hazard and estimating equations for βjs and γjs:
Λ0(t) = XZ t 0 dNi(u) nS1(u) , (6) M2j ≡ 1 √n j n X i=1 Z τj τj−1 {zi− SZ S1}dNi (u) = 0, j = 1, 2, . . . , m, and M3j ≡ 1 √n j n X i=1 Z τj τj−1 {vi− SV S1}dNi (u) = 0, j = 1, 2, . . . , m, (7)
Proposition 3.1.
Asymptotically (n↑ ∞, nj = O(n 2
3)), estimates of the piecewise parameters θj ≡ (βj, γj)0 solved
from the above equations have the following property: √n j(θbj− θj)−→ N(0, A−1j ), where Aj = A22,j A23,j A32,j A33,j ! , with A22,j = lim 1 nj XZ τj τj−1 {Zi(u)− SZ(u) S1(u)} 2dN i(u), A23,j = A32,j = lim 1 nj XZ τj τj−1 {Zi(u)− SZ(u) S1(u)}{Vi (u)−SV(u) S1(u)}dNi (u), A33,j = lim 1 nj XZ τj τj−1 {Vi(u)− SV(u) S1(u)} 2dN i(u). Points to consider
• How to choose the cut-off points. And, more importantly, how to detect (sequentially) the crossing points, in practice.
• For a regression set-up with multiple regressors, of course not all variables have varying effect; and not all the varying coefficients have the same ’crossing point(s)’.
• The Weibull model can also have a possibility to modeling cure fraction, by adjusting the ’scale parameter(s)’ (Carroll, 2002?, Controlled CLinical Trial). However, it is only a varying-coefficient ’homogeneous’ model, not being capable of capturing all patterns of crossings, in particular for those stemed from time-varying heteroscedasticity.
4
The Tests
4.1 Starting from A Degenerate Statistic
In this section we study the varying-coefficient HHRvc model, starting from considering a test proposed in Hsieh (2001) for the case of HHR model without varying coefficients. Preliminarily we calculate the following statistic:
Tdegen= m X j=1 {(Mc2j,Mc3j)Abj−1(Mc2j,Mc3j)0} (βbj,bγj,bΛ0)
with the parameters of interest being evaluated piecewisely at (βj, γj) = (βbj,γbj),∀j, whereβbj
The statistic is called degenerate because all degrees of freedom were consumed out at each segments. However, it offers an important clue to the construction of several tests for model validity. The tests studied in Wu et al. (2002) can be treated as a special case because the HHR model is a submodel of the HHRvc model. By this perspective, the Tdegen-statistic can
be amended to augment the degrees of freedom to 2m− 2 for the purpose of testing the validity of the HHR model (being nested within HHRvc), simply by repalceing all βjs and γjs with the
overall estimate β andb γ respectively.b
In this section, we propose several type-O and type-N tests with the possibility of viewing the tested submodels as being nested in the HHRvc class.
4.2 Several tests
As stated in Section 1, there are two types of test constructed under two levels. The type-O test has the flavour of an omnibus test. Although it is constructed by employing the HHRvc setting, the alternative hypothesis has no regard to the HHRvc model. So this kind of test is obtained only by suitably amending the Tdegen statistic to fit in for the following statistical hypotheses:
H0 : The true model is a sub-model of HHRvc model, against
Ha: The true model is NOT a sub-HHRvc model.
The type-N test, on the other hand, is constructed by assuming that the HHRvc model is true, and then test for a subset of the parameters at a specified value. That is, we will have the following hypotheses:
H0∗ : The true model is a sub-model of HHRvc model with some θk= θk0, against
Ha∗ : The true model is the HHRvc model, and θk6= θk0.
We can make the statement more specifically. Consider a sub-HHRvc model with the hazard: λHHR(t; z) = λ0(t){Λ0(t)}e
γ(t)z−1
e{β(t)z+γ(t)x}. (8)
We express θ(t) = (β(t)T, γ(t)T)T, and β(t)Tz+γ(t)Tx = θ(t)Tw, so that, for testing a θk= θk0,
the hypotheses according to the above two types of test can be re-written as follows: (1) Without regard to the HHRvc-model.
That means, we don’t assume the HHRvc class or subclass as a generally valid model for H0SHa,
but only for H0; and moreover, θk= θk0:
Ha: The true model is not HHRvc(θk=θk0).
(2) With regard to the HHRvc-model.
That is, we do assume the HHRvc class, and, under which, we are going to test for H0∗: θk= θk0, versus
Ha∗: θk6= θk0.
For the type-O test, we the proposed statistic T{·} is defined as
Tk= m X j=1 {McTjAb−1j Mcj} b θ(k),θk0,
for Mj = (M2j, M3j)T. For the type-N test, we should further consider a projection between
the spaces of information partitioned by the parameters θ(k) and θk. So the corresponding test
statistic is: Wk= m X j=1 {McTjAb∗ −1 j Mcj} b θ(k),θk0, for which A∗−1 j = (−1){Mj,kk− Mj,k(k)M−1j,(k)(k)Mj,(k)k}−1.
The application of W{·}-test needs further explanations. If an assumed model is correct, it
should then be noted that −E{∂ 2log(lik θ) ∂θk∂θk0 } = E{( ∂ log(likθ) ∂θk )(∂ log(likθ) ∂θk0 )}. (9)
Both the full likelihood and the partial likelihood (with Johansen’s (1983) meaning) satisfy the above equation. It is studied by Lin (1991, Statist. Sinica) that an omnibus test for the PH model can be constructed on the basis of measuring discrepancy of the second-order observed information matrix and the covariation matrix. For the type-N test, because the HHRvc model is assumed for H0 and Hatogether, we are indeed constructing the test statistic under H0SHa
within which the identity (9) holds. This test is said to be used under the meaning of with regard to the HHRvc-model. It is natural to consider an alternative -type-N test when A∗j−1 is
substituted by
A◦j−1 ={Aj,kk− Aj,k(k)A−1j,(k)(k)Aj,(k)k}−1.
Another useful test to compare with the type-N test is the (full-) likelihood ratio test. A similar concern can be referred to the likelihood ratio (LR) test of piecewisely varying-coefficient PH model proposed in Murphy (1993). Although the LR test does not follow a χ2 law in some nonstandard situations (for example in those the nuisance parameters become nonidentifiable in the Ha-space). Fortunately, our model-setting can have parallel results with those of Murphy
as long as the technical conditions of Section 3 are assumed. We will study and compare theW and LR tests through simulations reported later.
Example 4.1: Test for varying effect of a specific covariate
If the HHRvc is the underlying model and piecewise-constant approximation is utilized, then the Tdegen-statistic can be adopted, with a small amendment, to test for ’varying effect’
with respect to a specific covariate. As an illustrative example, if we want to test for constant heteroscedasticity (that is, γ(t) = γ0,∀t, for some constant γ0), the testing statistic can be
constructed as Thetvc = m X j=1 {McTjAb−1j Mcj} (βbj,γ0,bΛ0).
In practice, γ0 is substituted by an overall estimate bγ. In this case, we are exactly testing for
the piecewise-effect HHRvc model without varying heteroscedasticity (H0). If rejected, we only
know the constant-heteroscedasticity HHRvc model is rejected. That is the reason it is claimed having the flavor of an ’omnibus’ test. As for the second type test, we have
Whetvc = m X j=1 {McTjAb∗ −1 j Mcj} (βbj,γ0,bΛ0).
That is, assuming the HHRvc model, our hypotheses are H0∗ : γ1 = . . . = γm = γ0 versus
Ha∗ : γj’s are not all equal. If we set γ0 = γ, the statistic Tb hetvc will be a χ 2
m−1-variate
approximately.
Example 4.2: Test for the varying-coefficient PH model
There are tests and diagnostic plots proposed to check for varying effects under the varying-coefficients PH (PHvc-) model (Murphy, 1993; Martinussen and Scheike, 2001? 2002?; Valsecchi et al., 1996.) We are here at the stage of studying tests for a PHvc model. As before, we investigate two types of test statistics,T{·} andW{·}. For the first type,T{·}, we can compare it
with the performance of several tests proposed in Marzec and Marzec (1997); for the second one, the power ofW{·} is interestingly to be compared with that of the LR test . For both types, our
test statistics have the same form except for being evaluated at (cβj, 0,Λb0) for the j-th segment.
It is also worth of noting that the second-type test should not be compared to those of Murphy (1993) because, at the level of this example, we are actually inserting the HHRvc-model as an outer model of the PHvc-model. Under the hypotheses H0∗, γj = 0,∀j versus Ha∗, at least one of
the γjs is not equal to 0, theTphvc-statistic is distributed as χ2mfor large n. We can see that the
Thetvc-statistic is amended into Tphvc to further test whether the common heteroscedasticity is
zero (γ0 = 0). Because now there is no need to replace all γjs by an ’overall’ estimate γ, web
Likewise, the proposed test statistics can be modified slightly to test for equality among groups represented by different covariate-values. For example, for the second-type test, the null hypothesis and alternative hypothesis are: H0∗, βj = γj = 0,∀j; and Ha∗, at least one of βjs
and γjs is not equal to 0; and the test statistics are both distributed as χ22m asymptotically.
Furthermore, it is also redundant to introduce a test for PH-assumption under the HHRvc-model. The purpose is easily achieved by calculating the above two statistics at (βb0, 0,Λb0) when
b
β0 is a global estimate obtained under the hypothesis H0∗ : βj = β0,∀j, for some constant β0.
A commonly used test for equality is the logrank test in K-sample problem. It is a special case of score test under the PH regression setting. It is then also passed into the situation when there is cure-fraction (CF). When the CF probabilities are large for distinct groups, a genuine difference among groups could be masked (or ignored) by these large probabilities of cure-fraction. So, in this paper, the test Tequal can have greater power in testing for difference.
So it is also interesting to compare, in this case, the performance of the present Tequal with the
’modified’ score test studied by Bagdonavicius and Nikulin (2005) under their multiple cross-effect (MCE) model.
• A further point to consider: How to construct an omnibus test for goodness-of-fit of HHRvc model? Ref. Lin, 1991, JASA; Marzec and Marzec, 1997, Biometrika; Lin, Wei, and Ying, 1993, Biometrika.)
