- Email: [email protected]
Nonparametric hazard versus nonparametric frailty distribution in modelling recurrence of breast cancer
ELSEVIER
Journal of Statistical Planning and Inference 47 (1995) 111 127
joumal of statistical planning and inference
Nonparametric hazard versus nonparametric frailty distribution in modelling recurrence of breast cancer Dirley M. dos Santos, Richard B. Davies*, Brian Francis Centre.for Applied Statistics, Fylde College, Lancaster University, Lancaster LA l 4 YF, U K
Abstract
In extending survival models to include frailty effects, the relative merits of parametric and nonparametric formulations are unclear. The epidemiological emphasis (Clayton, Statist. in Med. 7 (1988) 819 841) has been upon nonparametric specification of conditional effects while, within econometrics (Heckman and Singer, Longitudinal Analysis of Labor Market Data (Cambridge University Press, Cambridge, 1985) 39 110) there has been some concern with the nonparametric specification of residual heterogeneity. This paper illustrates and examines some of the issues that arise in applying frailty models to the recurrence of breast cancer. The results suggest marginal advantage for the 'econometric' approach provided that there is no interest per se in the frailty distribution. However, it is important to adopt a flexible parametric specification for the conditional hazard and to allow different shaped hazard functions for successive durations. Key words: Frailty; Survival model; Piecewise exponential; Profile log-likelihood; Nonparametric mixing distribution, Recurrent events
1. Introduction
It has long been recognised that failure to control for unmeasured differences between individuals (usually termed 'frailty' in the epidemiological context and 'residual heterogeneity' in the social sciences) can result in misleading inference about temporal dependence (e.g. Bates and Neyman, 1952). The basic problem is that recent outcomes can be similar to earlier outcomes either because of a causal link (temporal dependence) or because some individuals are more prone to specific outcomes (spurious temporal dependence due to population heterogeneity). Early studies of the practical implications of this phenomenon include Massy et al. (1970) in the field of consumer behaviour; Ginsberg (1973) in sociology; Heckman and Willis (1977) and Chamberlain (1979) in econometrics; and Clayton (1978) in epidemiology, The
* Corresponding author. 0378-3758/95/$09.50 © 1995--Elsevier Science B.V. All rights reserved SSDI 0 3 7 8 - 3 7 5 8 ( 9 4 ) 0 0 1 2 5 - 1
112
D.M. dos Santos et al./ Journal of Statistical Planning and Inference 47 (1995) 111-127
importance of controlling for omitted variables in micro-level longitudinal analysis has been reinforced by both theoretical and empirical results which demonstrate that uncontrolled heterogeneity can also bias the estimated effects of the included explanatory variables (e.g. Lancaster 1979, 1985; Lancaster and Nickell, 1980; Heckman and Singer, 1985; Davies and Pickles, 1985; Ridder, 1987). In cancer studies there is every reason to expect frailty variation because the tendency for malignant cell transformation differs between individuals due to different exposures to carcinogens, different genetic dispositions, and so on. Even for cancers where almost all patients succumb eventually, the speed with which the disease develops may vary considerably. Appropriate epidemiological models for survival analysis with heterogeneous populations have been proposed in recent years (e.g. Clayton, 1978, 1988; Hougaard, 1984, 1986; Clayton and Cuzick, 1985a, b; Zeger and Karim, 1991) but have not been applied to many real datasets. There is therefore little experience of the practical problems which arise in their use. In this paper we compare several approaches for modelling the distributions of the hazard and frailty in order to identify the effects of time dependence and covariates on breast cancer recurrence at an individual level. We are therefore focusing upon methods for modelling the epidemiological processes in breast cancer recurrence and our primary objective is to investigate empirically the main parametric and nonparametric approaches which are feasible for this purpose. We emphasise that, in contrast to most comparable social science application areas, the unconditional hazard (i.e. the hazard obtained without allowing for frailty) is also important in medical research. The practitioner, in particular, does not know the individual's frailty and must make judgements based upon unconditional measures. A supplementary objective is therefore to investigate the differences between results obtained with and without allowance for heterogeneity. The most general model considered has the individual level (i.e. conditional) proportional hazards formulation
hij(tjlvi) = hoj(tjlvi) exp(fl' xij) vi, where hlj(. I') is the conditional baseline hazard for the jth recurrence, duration tj is measured from primary treatment (j = 1) or when the disease responded to treatment after the first or second recurrence (j = 2 or 3), fl'xij is the linear predictor, and vl is an individual-specific frailty parameter. However, the empirical investigations we report begin with the familiar Cox (1972) proportional hazards model and additional complexities are introduced progressively.
2. The data
The data used in this research cover the women referred to the Christie Hospital, Manchester, UK, with breast cancer between 1980 and 1985, and their subsequent monitoring until July 1991.
D.M. dos Santos et al./ Journal of Statistical Planning and ln['erence 47 (1995) 111 127
113
The data were checked for values out of range, incorrect sequence of events or dates~ and logical inconsistencies such as discrepancy between date of death and date of follow-up. Patients with bilateral disease (left and right breast with cancer) were excluded from the study. The youngest patient in the study was aged 21 yr, the oldest 88; mean age was 56.3 yr (standard deviation, 13.1 yr). The analyses reported in this paper use two explanatory variables which have proved to be improved prognostic indicators for breast cancer: Age and Stage of the disease at initial treatment. In preliminary work we investigated a three level menopause variable (premenopausal, menopausal, and postmenopausal) as an alternative to Age but it proved to be a less effective explanatory variable. No information on grading of the tumour was given in the dataset and nodal status was often missing but the data included a simple and consistent tumour classification: Stage 1 turnout size ~< 2 cm and no nodes involved (646 cases); Stage 2 - - tumour size > 2 cm to ~< 5 cm and no nodes, or tumour size ~< 5 cm with nodes (230 cases); Stage 3 tum o u r size > 5 cm with or without nodes (41 cases); and Stage 4 any tumour size with the presence of distant metastasis, i.e. disease elsewhere other than breast or local nodes. Classification o f t u m o u r stage was missing for 27 cases. In principle, these cases could have been included in the analysis by defining an additional level for Stage. This was not attempted because of the small numbers involved. As recurrence is not relevant for w o m a n with Stage 4 classification at initial presentation, these cases were also excluded from the analysis, leaving a total sample size of 917. The outcome of interest is duration to recurrence after treatment. Death from other causes was treated as right censoring and therefore assumed to be noninformative. This pragmatic assumption allows us to focus upon the methodological issues involved in modelling recurrence. It could be relaxed by adopting a competing risks formulation. The definition of duration to recurrence is illustrated in Diagram 1. Here A is the date of primary treatment; B, D, F are dates of first, second and third recurrences, respectively; C and E indicate dates when the disease is reported to have responded to treatment, after first and second recurrences; and G is date of death. The total survival time survival time since 2nd recurrence (
State 1
State 2
State 3 A
I
I
L
I
L--
I
I t
./.-
survival time since 1st recurrence A
B
C
D
E
F
Diagram 1. Illustration of a complete case.
G
114
D.M. dos Santos et aL / Journal of Statistical Planning and lnference 47 (1995) 111-127
duration to first recurrence is from A to B; second recurrence from C to D; and third recurrence from E to F. The duration to first recurrence is present for all women in the dataset, although it may be right censored. For a duration to a further recurrence to be included in the disease history, progression must have been followed by the recording of 'no disease' or 'complete response' (for overall response, soft tissue, bone, liver, pleura, lung and nodes). Durations to second and third recurrences are present for 403 and 116 women, respectively. In modelling repeated durations, we include the endogenous variable State (factor with 3 levels) to enable the intercept in the linear predictor to vary between durations.
3. Conventional proportional hazards models In this section we consider the conventional proportional hazards model without frailty and for single duration data given by (1)
hi(t) = h o ( t ) e x p ( f l ' x i ) ,
where hi(. ) is the unconditional hazard for case i, h0(. ) is the unconditional baseline hazard, fl is the parameter vector, xi is the covariate vector and t is the duration. This model is straightfoward to interpret because the explanatory variables have a multiplicative effect upon the baseline hazard. By using the well-known Cox (1975) partial likelihood approach, it is possible to estimate the ~ parameters without making any parametric assumptions about the baseline hazard. The results in Table 1 for duration to first recurrence use an equivalent model fitting method (the piecewise exponential; Breslow, 1974). Duration was measured from date of primary treatment. The results were obtained using GLIM4; many other software packages would be equally effective. From Table 1, we conclude that the unconditional hazard decreases with Age and that Stage is highly significant. It seems that younger women tends to have a faster rate of recurrence than older women. This is consistent with studies of breast cancer in women under 40 yr which have shown that this group of patients has a poorer prognosis than the older group (e.g. Ribeiro and Swindell, 1981). Also, the hazard for women for Stage 2 and Stage 3 severity is 2.3 and 4.8 times greater, respectively, than the hazard for women with Stage 1 severity. Table 1 Modelling the length of time to first recurrence Variable Age x 10-2 Stage 2 Stage 3
Estimate - 0.61 0.82 1.57
S.E.
Multiplier
0.38 0.11 0.18
2.3 4.8
115
D.M. dos Santos et al./ Journal of Statistical Planning and Inlerence 47 (1995) 111 127
-7.25
-7.50
"...Weibull
-7.75
-8.00 N -8.25 ¢-.
'~.~eibulI-Normal
-8.50
-8.75
l ~iecewise
-9.00
Exp.
-9.25
-9.50 2.00
I
I
I
I
I
I
I
3.00
4.00
5.00
6.00
7.00
8.00
9.00
log t i m e to first r e c u r r e n c e
Fig. 1. Unconditional baseline log-hazards for time to first recurrence.
Fig. 1 shows the baseline log-hazard function for first recurrence. To smooth this plot, we took each tenth distinct recurrence time for the piecewise representation. It is evident that the unconditional log-hazard declines after an initial decrease. Superimposed on Fig. 1 is the estimated hazard of first recurrence from a conventional Weibull model fitted to the single duration data.
4. Modelling the frailty and hazard for repeated durations In this section we examine extensions of the conventional proportional hazards model, allowing for repeated durations and frailty. However, the baseline hazard is assumed to be identical for each duration.
116
D.M. dos Santos et al. / Journal of Stat&tical Planning and Inference 47 (1995) 111-127
Most frailty models considered in the literature have a proportional hazards specification with a scalar and time-invariant frailty variable. Lancaster (1979) was one of the first to generalise the proportional hazards model by introducing a positive multiplicative disturbance to represent unobserved heterogeneity. The corresponding model for repeated durations is given by the conditional hazard hq(tjlvi) = ho~(t~lvi) exp(fl'xu) vi,
(2)
where vi is the individual-specific frailty variable, independent of the exogenous explanatory variables. This formulation is consistent with the notion of frailty being due to omitted variables and therefore having an additive effect on the linear predictor: hu(tjlvi) = hoj(tj]vi) exp(fl'xij + el),
(3)
where vi = exp(e/). These formulations provide the basis of a 'mixture' model approach in which the nuisance parameter vi or e~ is integrated out of the likelihood for each case. The absence of appropriate computer software has prevented the widespread application of mixture models. Researchers have used G L I M macros, NAG subroutines and have written their own software, but these approaches are difficult to operationalise and tend to be remarkably slow. Most researchers have used a fully parametric specification (parametric distributions for the hazard and frailty) for model (2), to minimise the model fitting problems. In general, the gamma distributions is used for the frailty (e.g. Lancaster and Nickell, 1980; Manton et al., 1986; Follmann and Goldberg, 1988) and the Weibull for the hazard. These are conjugate distributions which lead to an analytically tractable model specification. Some researchers have reported the parameters of the frailty distribution to be very sensitive to the choice of the hazard function. We use the parsimonious and popular Weibull distribution initially for the conditional hazard as a base against which to compare more complex models. We report the results in this section. Some authors have used a nonparametric approach for the conditional hazard to avoid this problem; we consider this further in Section 5. 4.1. The W e i b u l l - g a m m a frailty m o d e l
Assume that recurrence times for the ith woman are drawn from a Weibull distribution with scale parameter vi and shape parameter ~. Assume further that the scale parameter vl represents a random draw from a gamma mixing distribution, g(vi) = (O~/F(7)) v~ 1 e x p ( - Ovi), where F(-) denotes the gamma function and that the recurrence times are assumed independent (conditional on v~). Then, the total
D.M. dos Santos et al./ Journal of Statistical Planning and In[erence 47 (1995) 111 127
117
likelihood equals the product of the contributions of all n cases, and is L(fl)= fif i=1
~2I vi~t'[k-aexp(flkx~k)exp( - ~'itik~ exp(flkXik) ) ] d'k k=l
0 7
× [ e x p ( - Viti~kexp(flkXik))] 1 -d,~
V~" I e x p ( - Ovi)dvi, F(7) '
(4)
where dik is the censoring indicator (equals 0 for censoring, 1 otherwise) and x~k denotes the vector of explanatory variables corresponding to recurrence time t~k, and K~ is the number of durations for case i. The intercept on the log linear scale corresponds to ln(#), where # = 7/0 is the mean of the mixing distribution. Eq. (4) simplifies to i=1
F(7)
~ •
exp •
\k
= 1
dikflkXik ~Jl 1~ t(,~ 1)d,k, (5) -
where di = }~'-ldik' and ~i, t (~ = ~K, (~ k = 1 tik exp(flkXik). If~ = 1, Eq. (4) is the likelihood for a gamma mixture of exponentials, whereas if the gamma mixing distribution is degenerate, (4) reduces to a Weibull distribution. Thus, Eq. (4) allows for the two extreme causes of apparent decreasing hazard rates: due solely to heterogeneity and due solely to temporal dependence. If K~ = l, (4) reduces to the model for single durations considered by Lancaster (1979) and Lancaster and Nickell (1980) amongst others. This, in turn, is a generalisation of the two parameter Burr (1942) distribution.
4.2!. The Weibull-normal frailty model In the Weibull-gamma model, the {vl} in Eq. (2) are considered to be gamma distributed. The Weibull-normal model is similar except that the {ai} in Eq. (3) are assumed to be normal distributed and, as this does not lead to an analytically tractable expression, quadrature methods are used to approximate the integral. The unconditional density function for single duration may be written as
.~(t) = Ifi(tl~i) deb(ei/w) mo
~ {et a 1exp(flo + fl'xi + W~k)exp[ -- t'exp(flo + fl'xl + W~k)]} pk,
(6)
k=l
where 4~(. ) is the standard normal distribution function, the {~, p} are the quadrature point locations and masses, respectively, mo is the number of quadrature points (taken to be 8 in the analysis we report later), rio is the location parameter and w the scale parameter for the normal mixing distribution.
118
D.M. dos Santos et al. / Journal of Statistical Planning and Inference 47 (1995) 111-127
Similarly, the unconditional survival function, Si, is given by trl O
Si(t) = ~ exp[--Hexp(flo + fl'xi + W~k)]Pk,
(7)
k=l
thus, the unconditional hazard is calculated as
jS(t) hi(t) = S,(t--)"
(8)
It must be noted that this unconditional hazard includes the sample selection bias effect whereby the less fail are progressively over represented with duration. These equations are readily extended to include multiple durations for each case. 4.3. The Weibull-nonparametric frailty model Evidence that results are not robust to choice of parametric form for the frailty distribution, at least for single duration data (see, in particular, Heckman and Singer, 1982), has stimulated interest in nonparametric representation of frailty (e.g. Laird, 1978; Heckman and Singer, 1984; Davies and Crouchley, 1984). The main approach uses the theoretical result, dating back to Kiefer and Wolfowitz (1956), that the nonparametric characterisation of the frailty distribution within maximum likelihood estimation is confined to a finite number of mass points. It also draws upon the empirical experience that the number of mass points tends to be small (Davies, 1993). For single duration data, the Weibull-nonparametric model may be written as mO
f ( t ) = ~ {~t " - t exp(fl'xi + ~ k ) e x p [ - - t ' e x p ( f l % + ~k)]} Pk. k=t
(9)
Note that this is algebraically similar to the Weibull-normal model with quadrature (Eq. (6)) but the m0 and {~,p} are no longer specified externally and have to be estimated. The location and scale parameters have been absorbed into the nonparametric representation. Apart from results on consistency, large sample theory has yet to be developed for the nonparametric frailty approach. However, limited simulation evidence (Heckman and Singer, 1982; Davies, 1987) suggests standard large sample properties and, in particular, provides some justification for the use of likelihood ratio tests and standard errors derived from the estimated information matrix. 4.4. Results on breast cancer recurrence The results of fitting the models discussed above to the breast cancer data are shown in Table 2. The software MIXTURE (Ezzet and Davies, 1988) was used to fit these models although a modified version was required for the Weibull-gamma. The Weibull-nonparametric model required just four mass points. Of course, the nonparametric characterisation of the frailty distribution may require further small mass
D.M. dos Santos et al. / Journal o f Statistical Planning and Inference 47 (1995) 111 127
119
Table 2 Conventional and mixture models Variable
Weibull
Weibull-gamma
Estimate S . E . Age x 10 2 Stage 2 Stage 3 State 2 State 3 [/o w 0 7 log-likelihood
- 0.76 0.66 1.07 1.53 1.81 0.75 -- 6.17
0,25 0,06 0.12 0.07 0.11 0.04 0.25
- 6786.91
Estimate S.E. - 1.30 1.20 1.84 1.25 1.47 1.05
0.40 0.11 0.14 0.08 0.10 0.04
1002.82 0.26 0.47 0.10 - 6729.79
Weibull-normal
Estimate S . E . 1.00 1.26 2.01 1.44 1.74 1.05 -- 8.74 1.39
0.50 0.15 0.33 0.09 0.14 0.04 0.46 0.12
- 6735.71
Weibullnonparametric Estimate S.E. - 1.40 1.17 2.39 1.29 1.47 1.07
0.50 0.16 0.37 0.1 I 0.15 0.04
6726.47
p o i n t s b u t these w o u l d n o t alter the l o g - l i k e l i h o o d a n d structural p a r a m e t e r estimates within the convergence criteria used (primarily, a r e d u c t i o n of less t h a n 10 _4 in the l o g - l i k e l i h o o d after a N e w t o n - R a p h s o n i t e r a t i o n with step length adjustment); the four mass p o i n t s p r o v e an a d e q u a t e if not c o m p l e t e n o n p a r a m e t r i c c h a r a c t e r i s a t i o n of the mixing distribution. T h e r e is an a p p r e c i a b l e increase in the l o g - l i k e l i h o o d when frailty is represented explicitly in the m o d e l (Table 2) and, a l t h o u g h formal l i k e l i h o o d ratio tests are not a p p r o p r i a t e , the n o n p a r a m e t r i c r e p r e s e n t a t i o n of frailty results in a n o t a b l e i m p r o v e m e n t u p o n the fully p a r a m e t r i c models. Theoretically, i g n o r i n g h e t e r o g e n e i t y leads to u n d e r e s t i m a t i o n of the effects of e x o g e n o u s e x p l a n a t o r y variables a n d their s t a n d a r d errors a n d also to s p u r i o u s r e d u c t i o n of the W e i b u l l p a r a m e t e r ( L a n c a s t e r a n d Nickell, 1980). This is consistent with the results for Age a n d Stage in T a b l e 2 a n d it a p p e a r s that the W e i b u l l m o d e l w i t h o u t heterogeneity p r o v i d e s a quite misleading picture of the h a z a r d declining with time; the h e t e r o g e n e o u s m o d e l s suggest that the c o n d i t i o n a l h a z a r d might actually increase, a l t h o u g h the results are n o t inconsistent with a c o n s t a n t hazard. The results for State are m o r e difficult to i n t e r p r e t because this is an e n d o g e n o u s r a t h e r than e x o g e n o u s e x p l a n a t o r y variable. T h e m a i n c o n s e q u e n c e of allowing for frailty is a r e d u c t i o n in the e s t i m a t e d h a z a r d following first a n d second recurrences. This is a s a m p l e selection effect; those who have experienced recurrence will tend to be those w h o are m o r e frail. W i t h o u t c o n t r o l for frailty, the e s t i m a t e d h a z a r d of further recurrence o r d e a t h will be b i a s e d u p w a r d . S u p e r i m p o s e d on Fig. 1 is the u n c o n d i t i o n a l e s t i m a t e d rate of first recurrence using the W e i b u l l - n o r m a l m o d e l (see Eq. (8)) with e x p l a n a t o r y variables set equal to zero. It is therefore a p a r a m e t r i c equivalent of the piecewise e x p o n e n t i a l in Fig. 1 a l t h o u g h it
120
D.M. dos Santos et al./ Journal of Statistical Planning and Inference 47 (1995) 111-127
has been obtained from the estimated conditional distribution rather than estimated directly from the data. It is noted that, whereas the conditional hazard for each woman increases over time in Weibull-normal model, the unconditional rate of recurrence shows the typical decline after an initial increase. This demonstrates the statistical effect whereby the less frail progressively dominate the survivors. Conversely, these results provide no evidence that the declines in the recurrence rate shown in Fig. 1 are attributable to any regenerative process.
5. Different baseline hazards for each recurrence type
In this section we allow the conditional baseline hazard to be different for each duration. Both parametric and nonparametric characterisation of the hazards are considered. A development due to Clayton (1988) enables estimation of models with nonparametric conditional hazards and multiple durations for each case. Frailty is assumed to have a gamma distribution. The approach is based upon a piecewise exponential formulation. As for the mass-point semiparametric approach, formal large sample theory has not been developed. But simulation evidence is encouraging with respect to consistency, likelihood ratio tests, and the normal distribution of parameter estimates (Nielsen et al., 1992) and Clayton and Cuzick (1985a) provide a heuristic justification for partial likelihood inference, using the results of Voelkel and Crowley (1984), if the model has a semi-Markov structure with ordered states such that transition from state k to state d is only possible if d > k. This condition is not overly restrictive for processes, such as many diseases, with an explicit progression. In particular, it enables the breast cancer data to be analysed with a different nonparametric hazard for each of the three ordered durations corresponding to State 1 (from primary treatment), State 2 (from first recurrence), and State 3 (from second recurrence). More specifically, the conditional hazard function for State j is given by hii(t~lvi) = hoj(tj[vi) exp(fl'xi~) vi,
(1o)
where tj is measured from initial treatment (j = 1) or when the disease last responded to treatment (j = 2 or 3), hot('l') is the conditional baseline hazard, and, in the Clayton method, {vj} are independent gamma variates with unit mean and variance tr 2. The problem is to estimate /3, tr 2, and the three unknown functions h0~('l-), j = 1,2,3. Epochs for the piecewise exponential formulation are constructed separately for each State. For example, the durations to first recurrence are divided into epochs on the basis of all the durations from initial treatment to first recurrence (or censoring) in the sample. Each baseline hazard h 0 / . ].) is represented in the partial likelihood by
D.M. dos Santo.; et al./ Journal of Statistical Planning and Inference 47 (t995) l l l 127
a step function
{hjk}
121
where k indexes epochs. The partial likelihood may be written
as
L(fl) =
lq [(Aijk)a'J"exp(--aijkUjk)] i=1
lexp(--Tvi)dvi ,
(11)
=
where Kij is the number of epochs for individual i in State j, Aijk = vihjkexp(fl'xO, U;k is the epoch duration, dljk is an outcome indicator variable (equals 1 if case i has a recurrence in epoch k of state j, and 0 otherwise), and 7 = 1/°-2, The GLIM macros given by Clayton (1988) use an EM algorithm to estimate the parameters in Eq. (11), estimating the v~ at the E step and the fl parameters and the hjk at the M step. To obtain a maximum likelihood solution it is necessary to repeat the algorithm for different values of cr2. Clayton emphasises that this profile log-likelihood approach is computationally simpler than seeking to estimate the variance simultaneously with the other parameters. Nielsen et al. (1992) propose a similar approach. The model fitting results for the Clayton approach including Age and Stage as explanatory variables are shown in Table 3. The piecewise conditional baseline log-hazards are plotted in Fig. 2(a), smoothed by taking every tenth distinct recurrence time for each State. It is evident that, at least for State 1 the hazard increases Table 3 Results of the semiparametric models Variable
Weibull non-parametric frailty
Clayton Estimate Agex Stage Stage State State
10 2 2 3 2 3
O-2
1.12 0.99 1.55
S.E.
Estimate
S.E.
0.28 0.07 0.13
-- 1.00 1.27 1.94 2.55 2.71
0.50 0.17 0.38 0.52 0.63
1.22 1.00 1.00
0.07 0.05 0.08
Estimate -- 1.00 1.14 1.80 3.47 4.05
S.E. 0.40 0.16 0.26 0.69 0.88
0.95
~2 ~3 71 × 103 ~'2 × 103 ]~3 × 103
-- 6.89 -- 8.08 -- 9.71
~2
G Pl ,o2 P3 P4
log-likelihood
Weibull Gompertz non-parametric frailty
6396.08
-- 13.28 0.07 0.26 0.33 0.34 -- 6723.34
1.54 1.22 1.14 1.02 1.19 0.88 8.69 1O.06 11.58
0.10 0.07 0.13 0.16 0.29 0.60
0.06 0.33 0.55 0.06 6692.24
122
D.M. dos Santos et al./ Journal of Statistical Planning and Inference 47 (1995) 111-127 (a)
Clayton approach
(b) Weii~ulFgompertz
-55
-55 ~=~..-.-! ....
i
i!;~"~ ~ i •
-60
lii - , I. ~ ~ i
- -
-
-.
STATE3
-6S
- -
STATE I
......
STATE 3
STARE 2 .........................
STATE t STATE 2
-'~
-60
~5
i.]
'"... -70
-7O
-75
-75
-~ O
~o
-aS
-8 5
-~0
-e 0
3OO
4OO
5OO
6OO
7OO
8OO
3o0
4o0
log time
6oo
600
700
a.o0
log lima
Fig. 2. Conditional baseline log-hazards. and then decreases. Moreover, the visual evidence suggests that the three hazards vary in shape. To fit a Weibull-nonparametric model with different shape parameters for each State we used N A G subroutines. The Weibull results, also shown in Table 3, are for four mass points because attempts at adding a fifth mass point did not improve the log-likelihood and the models showed exactly the results of the previous solution. The results indicate, in contrast to those from the Clayton method, increasing or constant hazards for the three states. The parameter estimates for Stage are about 25% larger than those for the Clayton model. The standard error estimates are also much larger for the Weibull model but we know that, because of the iterative algorithm, the Clayton model standard errors will tend to be underestimated (see, for example, Aitkin and Clayton, 1980). There is an improvement of only 3.13 in the log-likelihood with two additional parameters in comparison with the Weibull-nonparameteric model with identical baseline hazards for each recurrence (see Table 2).
6. Resolving the discrepancies The problem in reconciling the different conditional hazard and parameter estimate results for the Weibull-nonparametric and the Clayton nonparametric hazard-gamma models is that either or both may be misspecified. In particular, the Weibull cannot represent a nonmonotonic hazard and we have no assurance that the gamma is an appropriate distribution for the frailty (see Hougaard, 1986). Ideally, we would fit a model with nonparametric hazard and nonparametric frailty to provide a nested structure for the comparison but such models would require unrealistic computing resources and there are doubts about their identifiability. We therefore proceed in this
D.M. dos Santos et al./ Journal o f Statistical Planning and In/erence 47 (1995) 111 127
123
section to investigate the discrepancies between the Weibull-nonparametric and Clayton models by allowing more flexible parametric modelling of the hazards. The three parameter Weibull-Gompertz distribution (Lee, 1980) provides appreciable additional flexibility. The conditional baseline hazard in Eq. (10) is given by
hoj(tjlvi)
= c~jt~' l exp(~flj),
t~ >7 0,
(12)
with ~ j > 0 and - G o < ~ < ~ . The model reduces to the Weibull when 7 j = 0 . When ~j > 1, 7j > 0, ho~(tjlvi) is monotonic increasing; when ~j > 1, ~ < 0, hoj(tjlvi) increases and then decreases; when ~j < l, 7i > 0, hoj(tj[vi) decreases from infinity to a minimum and then increases to infinity; and when ~j < 1, ),j< 0, h0j(t~lv3 is monotonic decreasing. N A G optimisation subroutines were used to fit the Weibull-Gompertz model with nonparametric frailty. As the density and survival functions of the Weibull-Gompertz do not have explicit forms, they were approximated by Taylor series expansions with sufficient terms to ensure that the log-likelihood is accurate to at least the fourth decimal place. The survival function (dropping subscripts for simplicity) is
S(t)=exp
-~t ~
+
7t + + + --. c~ + 1 2!(~ + 2) 3!(c~+ 3)
.
(13)
The results are shown in Table 3. With ~ > 1 and 7 < 0 for all three hazards, these results are consistent with the pattern of increasing and then decreasing hazards suggested by the Clayton model. Also, the 31.10 increase in the log-likelihood with_ just three parameters compared with the Weibull-nonparametric (Table 3) model confirms that the Weibull hazard specification is untenable. Constraining the Weibull-Gompertz hazards to the same shape for each State reduces the log-likeli-. hood by 10.23 with four fewer parameters. It is now evident that the hazard functions do differ significantly. The estimated mass points give frailty variances of 2.33 and 1.74 for the Weibull and Weibull-Gompertz model, respectively, after scaling to a unit mean to ensure comparability with the Clayton approach. Unsurprisingly, the estimated frailty variance is reduced by improving the flexibility of the hazard specification. The baseline hazards for the Weibull-Gompertz model with nonparametric frailty are plotted in Fig. 2(b). To ensure comparability with the piecewise hazards in Fig. 2(a), they have been evaluated at the mean of the transformed mass point locations, exp(¢k). Thus the baseline hazard for State j is given by 3
ho~(tjlv,)
=
~ ~j ^ tj-~j 1 exp(~3jti + /~j)exp(~k),bk, k=l
where fli is the estimate of the intercept in the linear predictor for State j. Although they have the same 'sickle' shape, it is emphasised that these conditional hazards still differ markedly from the unconditional hazards. This is illustrated in Fig. 3 which compares the conditional hazard for first recurrence with the unconditional hazard estimated from the Weibull-Gompertz results (cf. Section 4.4 for the Weibull-normal
124
D.M. dos Santos et al./Journal of Statistical Planning and Inference 47 (1995) 111-127
-6.5
-7.0
. . .,. . .'"'"'°"""
..°.. .............. ...o... "'"'""-.,
.,°
-7.5 "O N m t-
-8.0
.go -8.5
-9.0
-9.5 I
I
I
I
I
3.00
4.00
5.00
6.00
7.00
8.00
log time Fig. 3. Estimated baseline hazards for first recurrence, Weibu]]-Gompertz.
model). As expected, the unconditional hazard is progressively lower than the conditional hazard. After 4 yr, the conditional hazard is more than twice the unconditional Moreover, the unconditional hazard begins to decline earlier (after less than 11 months compared with about 18 months for the conditional hazard); any regenerative process is operating more slowly than suggested by the unconditional h a z a r d
7. Concluding remarks The results in this paper have confirmed the importance of recognising frailty effects in any modelling of breast cancer recurrence and of distinguishing between conditional and unconditional hazards; the hazards may have different shape
D.M. dos Santos et al./ Journal ~?f Statistical Planning and ln[erence 47 (1995) 111 127
125
characteristics and be influenced differently by explanatory variables. The results also demonstrate that there is no technical impediment to the effective application of' semiparametric frailty models to epidemiological data. Flexible specification of hazard functions is a distinctive feature of the epidemiological tradition in survival analysis with, in particular, a strong emphasis upon nonparametric functions. This distinguishes it from the econometric tradition in which theoretical arguments are often adduced to justify the adoption of monotonic, parametric hazard functions. One constraint on including frailty effects within the epidemiological tradition is the empirical evidence of identifiability problems when a nonparametric specification is adopted for both the hazard and the frailty distribution; model fitting tends to converge to a degenerate, zero frailty solution even for simulated data with known frailty characteristics (Meyer, 1990). More seriously, Clayton (1988) has noted a similar tendency for nonparametric hazards with parametric frailty. This contrasts with the frailty models favoured by econometricians; empirical experience is encouraging and identifiability has been established under comparatively weak conditions for fully parametric models and models with parametric hazards and nonparametric frailty distributions, even for single duration data (Elbers and Ridder, 1982; Heckman and Singer, 1984). Results reported elsewhere using subsets of the breast cancer data (dos Santos, 1994) provide further evidence of the convergence to zero frailty variance with fully nonparametric specifications. However, no such problems were encountered in fitting the nonparametric hazard-parametric frailty models discussed in this paper. We would speculate that, at least for parametric frailty, this effect arises through model misspecification and/or inadequate sample sizes. Nevertheless, a more satisfactory compromise in many situations may be to use a parametric hazard and a nonparametric frailty distribution. The advantage of this approach is that it confines the nonparametric specification to the nuisance parameters; the frailty distribution is rarely of interest in its own right and, in any case, is difficult to estimate with precision by parametric (Nielsen eta., 1992) or nonparametric (Heckman and Singer, 1984) characterisations. Conversely, this approach enables conventional inferential methods to be applied to the characteristics of the hazards, to compare shapes, identify turning points, and so on. Concerns about the flexibility of the hazard specification can be addressed by using parametric forms which are more complex than those often recommended for survival analysis. The three parameter WeibuU Gompertz distribution used in this paper appears to have seen sufficiently flexible to represent the main features of the hazard for breast cancer recurrence.
Acknowledgements This research was supported, in part, by a grant from CAPES-Brazil, proc. 189/90-1 and the ESRC (Research project R000233850). The authors are also grateful to an anonymous referee for helpful comments.
126
D.M. dos Santos et al./ Journal of Statistical Planning and Inference 47 (1995) 111-127
References Aitkin, M. and D. Clayton (1980). The fitting of exponential, Weibull and extreme value distributions to complex censored survival data using GLIM. Appl. Statist. 29, 156-163. Bates, G. and J. Neyman (1952). Contributions to the theory of accident proneness II: true or false contagion. Univ. California Publ. Statist. 1, 255-275. Breslow, N. (1974). Covariance analysis of censored survival data. Biometric 30, 89-99. Burr, I.W. (1942). Cumulative frequency functions. Ann. of Math. Statist. 13, 215-232. Chamberlain, G. (1979). Heterogeneity, omitted variable bias and duration dependence. Harvard Institute of Economic Research, DP691. Cambridge, MA. Clayton, D.G. (1978). A model for association in bivariate life-tables and its application in epidemiological studies of familial tendency in chronic disease incidence. Biometrika 65, 141-151. Clayton, D.G. (1988). The analysis of event history data: a review of progress and outstanding problems. Statist. in Med. 7, 819-841. Clayton, D.G. and J. Cuzick (1985a). Multivariate generalizations of the proportional hazards model (with discussion). J. Roy. Statist. Soc. Set. A 148(2), 82-117. Clayton, D.G. and J. Cuzick (1985b). The semi-parametric pareto model for regression analysis of survival times. Bull. Inst. Statist. Inst. 51(23.3), 1 18. Cox, D.R. (1972). Regression models and life tables (with discussion). J. Roy. Statist. Soc. Ser. B 34, 187 220. Cox, D.R. (1975). Partial likelihood. Biometrika 62(2), 269-276. Davies, R.B. (1987). Mass point methods for dealing with nuisance parameters. In: R. Crouchley ed., Longitudinal Data Analysis, Aldershot, Avebury. Davies, R.B. (1993). Nonparametric control for residual heterogeneity in modelling recurrent behaviour. Comput. Statist. Data Anal. 16, 143-160. Davies, R.B. and R. Crouchley (1984). Calibrating longitudinal models of residential mobility and migration. Regional. Sci. Urban Econom. 14, 231-247. Davies, R.B. and A.R. Pickles (1985). Longitudinal versus cross-sectional methods for behavioural research: a first round knockout. Environ. Plann. A 17, 1315-1329. dos Santos, D.M. (1994). Distinguishing the effects of time dependence from the effects of frailty in multiple spell breast cancer data. Unpublished Ph.D. thesis, Lancaster University, Lancaster, UK. Elbers, C. and G. Ridder (1982). True and spurious duration dependence: the identifiability of the proportional hazard model. Rev. Econom. Stud. 49, 403-409. Ezzet, F.L. and R.B. Davies (1988). A manual for MIXTURE. Centre for Applied Statistics, University of Lancaster, Lancaster, England. Follmann, D.A. and M,S. Goldberg (1988). Distinguishing heterogeneity from decreasing hazard rates. Technometrics 30(4), 389-396. Francis, B., M. Green and C. Payne (1992). The GLIM System: Release 4 Manual. Oxford University Press, Oxford. Ginsberg, D.B. (1973). Stochastic models of residential geographic mobility for heterogeneous populations. Environ. Plann. 5, 113-124. Heckman, J.J. and B. Singer (1982). The identification problem in econometric models for duration data. In: W. Hildenbrand, Ed., Advances in Econometrics, Proc. World Meetings of the Econometric Soc 1980. Cambridge University Press, Cambridge. Heckman, J.J and B. Singer (1984). A method for minimizing the impact of distributional assumptions in econometric models for duration data. Econometriea 52(2), 271-320. Heckman, J.J. and B. Singer (1985). Social science duration analysis. In: J.J. Heckman and B. Singer, Eds., Longitudinal Analysis of Labor Market Data. Cambridge University Press, Cambridge, 39-110. Heckman, J.J. and R.J. Willis (1977). A beta-logistic model for the analysis of sequential labor force participation by married women. J. Political Economy 85, 27-58. Hougaard, P. (1984). Life table methods for heterogeneous populations: distributions describing the heterogeneity. Biometrika 71, 75 83. Hougaard, P. (1986). Survival models for heterogeneous populations derived from stable distributions. Biometrika 73(2), 387-396. Kiefer, J. and J. Wolfowitz (1956). Consistency of the maximum likelihood estimator in the presence of infinitely many incidental parameters. Ann. of Math. Statist. 27, 887-906.
D.M. dos Santos et al./ Journal of Statistical Planning and ln[erence 47 (1995) 11l 127
127
Laird, N. (1978). Nonparametric maximum likelihood estimation of a mixing distribution. J. Amer. Statist. Assoc. 73, 805-811. Lancaster, T. (1979). Econometric methods for the duration of unemployment. Econometrica 47, 939 956. Lancaster, T. (1985). Generalised residuals and heterogeneous duration models: with applications to the Weibull model. J. Econometrics 28(1), 155--169. Lancaster, T. and S. Nickell (1980). The analysis of re-employment probabilities for the unemployed. J. Roy Statist. Ser. Soc. A 143, 141 152. Lee, L. (1980). Testing adequacy of the Weibull and log linear rate models for a Poisson process Technometrics 22(2), 195 199. Manton, K.G., E. Stallard and J.W. Vaupel (1986). Alternative models for the heterogeneity of mortality risks among the aged. J. Amer. Statist. Assoc. 81,395, 635-644. Massy, W.F., D.M. Montgomery and D.G. Morrison (1970). Stochastic Models of Buying Behaviour. MIT Press, Cambridge, MA. Meyer, B.D. (1990). Unemployment insurance and unemployment spells. Econometrica 58(4), 757 782. NAG (1983). Numerical Algorithms Group Library Manual. Oxford, England. Nielsen, G.G., R.D. Gill, P.K. Andersen and T.I.A. Sorensen (1992). A counting process approach to maximum likelihood estimation in frailty models. Scand. J. Statist. 19, 25-43. Ribeiro, G.G. and R. Swindell (1981). The prognosis of breast carcinoma in women aged less than 40 years. Clin Radiol. 32, 231 236. Ridder, G. (1987). The sensitivity of duration models to misspecified unobserved heterogeneity and duration dependence. Working Paper, University of Amsterdam. Voelkel, J.G. and J. Crowley (1984). Non-parametric inference for a class of semi-markov models with censored observations. Ann. of Statist. 12, 142-160. Zeger, S.L. and M.R. Karim (1991). Generalized linear models with random effects, a Gibbs sampling approach. J. Amer. Statist. Assoc. 86, 413; Theory and methods, 79-86.
Journal of Statistical Planning and Inference 47 (1995) 111 127
joumal of statistical planning and inference
Nonparametric hazard versus nonparametric frailty distribution in modelling recurrence of breast cancer Dirley M. dos Santos, Richard B. Davies*, Brian Francis Centre.for Applied Statistics, Fylde College, Lancaster University, Lancaster LA l 4 YF, U K
Abstract
In extending survival models to include frailty effects, the relative merits of parametric and nonparametric formulations are unclear. The epidemiological emphasis (Clayton, Statist. in Med. 7 (1988) 819 841) has been upon nonparametric specification of conditional effects while, within econometrics (Heckman and Singer, Longitudinal Analysis of Labor Market Data (Cambridge University Press, Cambridge, 1985) 39 110) there has been some concern with the nonparametric specification of residual heterogeneity. This paper illustrates and examines some of the issues that arise in applying frailty models to the recurrence of breast cancer. The results suggest marginal advantage for the 'econometric' approach provided that there is no interest per se in the frailty distribution. However, it is important to adopt a flexible parametric specification for the conditional hazard and to allow different shaped hazard functions for successive durations. Key words: Frailty; Survival model; Piecewise exponential; Profile log-likelihood; Nonparametric mixing distribution, Recurrent events
1. Introduction
It has long been recognised that failure to control for unmeasured differences between individuals (usually termed 'frailty' in the epidemiological context and 'residual heterogeneity' in the social sciences) can result in misleading inference about temporal dependence (e.g. Bates and Neyman, 1952). The basic problem is that recent outcomes can be similar to earlier outcomes either because of a causal link (temporal dependence) or because some individuals are more prone to specific outcomes (spurious temporal dependence due to population heterogeneity). Early studies of the practical implications of this phenomenon include Massy et al. (1970) in the field of consumer behaviour; Ginsberg (1973) in sociology; Heckman and Willis (1977) and Chamberlain (1979) in econometrics; and Clayton (1978) in epidemiology, The
* Corresponding author. 0378-3758/95/$09.50 © 1995--Elsevier Science B.V. All rights reserved SSDI 0 3 7 8 - 3 7 5 8 ( 9 4 ) 0 0 1 2 5 - 1
112
D.M. dos Santos et al./ Journal of Statistical Planning and Inference 47 (1995) 111-127
importance of controlling for omitted variables in micro-level longitudinal analysis has been reinforced by both theoretical and empirical results which demonstrate that uncontrolled heterogeneity can also bias the estimated effects of the included explanatory variables (e.g. Lancaster 1979, 1985; Lancaster and Nickell, 1980; Heckman and Singer, 1985; Davies and Pickles, 1985; Ridder, 1987). In cancer studies there is every reason to expect frailty variation because the tendency for malignant cell transformation differs between individuals due to different exposures to carcinogens, different genetic dispositions, and so on. Even for cancers where almost all patients succumb eventually, the speed with which the disease develops may vary considerably. Appropriate epidemiological models for survival analysis with heterogeneous populations have been proposed in recent years (e.g. Clayton, 1978, 1988; Hougaard, 1984, 1986; Clayton and Cuzick, 1985a, b; Zeger and Karim, 1991) but have not been applied to many real datasets. There is therefore little experience of the practical problems which arise in their use. In this paper we compare several approaches for modelling the distributions of the hazard and frailty in order to identify the effects of time dependence and covariates on breast cancer recurrence at an individual level. We are therefore focusing upon methods for modelling the epidemiological processes in breast cancer recurrence and our primary objective is to investigate empirically the main parametric and nonparametric approaches which are feasible for this purpose. We emphasise that, in contrast to most comparable social science application areas, the unconditional hazard (i.e. the hazard obtained without allowing for frailty) is also important in medical research. The practitioner, in particular, does not know the individual's frailty and must make judgements based upon unconditional measures. A supplementary objective is therefore to investigate the differences between results obtained with and without allowance for heterogeneity. The most general model considered has the individual level (i.e. conditional) proportional hazards formulation
hij(tjlvi) = hoj(tjlvi) exp(fl' xij) vi, where hlj(. I') is the conditional baseline hazard for the jth recurrence, duration tj is measured from primary treatment (j = 1) or when the disease responded to treatment after the first or second recurrence (j = 2 or 3), fl'xij is the linear predictor, and vl is an individual-specific frailty parameter. However, the empirical investigations we report begin with the familiar Cox (1972) proportional hazards model and additional complexities are introduced progressively.
2. The data
The data used in this research cover the women referred to the Christie Hospital, Manchester, UK, with breast cancer between 1980 and 1985, and their subsequent monitoring until July 1991.
D.M. dos Santos et al./ Journal of Statistical Planning and ln['erence 47 (1995) 111 127
113
The data were checked for values out of range, incorrect sequence of events or dates~ and logical inconsistencies such as discrepancy between date of death and date of follow-up. Patients with bilateral disease (left and right breast with cancer) were excluded from the study. The youngest patient in the study was aged 21 yr, the oldest 88; mean age was 56.3 yr (standard deviation, 13.1 yr). The analyses reported in this paper use two explanatory variables which have proved to be improved prognostic indicators for breast cancer: Age and Stage of the disease at initial treatment. In preliminary work we investigated a three level menopause variable (premenopausal, menopausal, and postmenopausal) as an alternative to Age but it proved to be a less effective explanatory variable. No information on grading of the tumour was given in the dataset and nodal status was often missing but the data included a simple and consistent tumour classification: Stage 1 turnout size ~< 2 cm and no nodes involved (646 cases); Stage 2 - - tumour size > 2 cm to ~< 5 cm and no nodes, or tumour size ~< 5 cm with nodes (230 cases); Stage 3 tum o u r size > 5 cm with or without nodes (41 cases); and Stage 4 any tumour size with the presence of distant metastasis, i.e. disease elsewhere other than breast or local nodes. Classification o f t u m o u r stage was missing for 27 cases. In principle, these cases could have been included in the analysis by defining an additional level for Stage. This was not attempted because of the small numbers involved. As recurrence is not relevant for w o m a n with Stage 4 classification at initial presentation, these cases were also excluded from the analysis, leaving a total sample size of 917. The outcome of interest is duration to recurrence after treatment. Death from other causes was treated as right censoring and therefore assumed to be noninformative. This pragmatic assumption allows us to focus upon the methodological issues involved in modelling recurrence. It could be relaxed by adopting a competing risks formulation. The definition of duration to recurrence is illustrated in Diagram 1. Here A is the date of primary treatment; B, D, F are dates of first, second and third recurrences, respectively; C and E indicate dates when the disease is reported to have responded to treatment, after first and second recurrences; and G is date of death. The total survival time survival time since 2nd recurrence (
State 1
State 2
State 3 A
I
I
L
I
L--
I
I t
./.-
survival time since 1st recurrence A
B
C
D
E
F
Diagram 1. Illustration of a complete case.
G
114
D.M. dos Santos et aL / Journal of Statistical Planning and lnference 47 (1995) 111-127
duration to first recurrence is from A to B; second recurrence from C to D; and third recurrence from E to F. The duration to first recurrence is present for all women in the dataset, although it may be right censored. For a duration to a further recurrence to be included in the disease history, progression must have been followed by the recording of 'no disease' or 'complete response' (for overall response, soft tissue, bone, liver, pleura, lung and nodes). Durations to second and third recurrences are present for 403 and 116 women, respectively. In modelling repeated durations, we include the endogenous variable State (factor with 3 levels) to enable the intercept in the linear predictor to vary between durations.
3. Conventional proportional hazards models In this section we consider the conventional proportional hazards model without frailty and for single duration data given by (1)
hi(t) = h o ( t ) e x p ( f l ' x i ) ,
where hi(. ) is the unconditional hazard for case i, h0(. ) is the unconditional baseline hazard, fl is the parameter vector, xi is the covariate vector and t is the duration. This model is straightfoward to interpret because the explanatory variables have a multiplicative effect upon the baseline hazard. By using the well-known Cox (1975) partial likelihood approach, it is possible to estimate the ~ parameters without making any parametric assumptions about the baseline hazard. The results in Table 1 for duration to first recurrence use an equivalent model fitting method (the piecewise exponential; Breslow, 1974). Duration was measured from date of primary treatment. The results were obtained using GLIM4; many other software packages would be equally effective. From Table 1, we conclude that the unconditional hazard decreases with Age and that Stage is highly significant. It seems that younger women tends to have a faster rate of recurrence than older women. This is consistent with studies of breast cancer in women under 40 yr which have shown that this group of patients has a poorer prognosis than the older group (e.g. Ribeiro and Swindell, 1981). Also, the hazard for women for Stage 2 and Stage 3 severity is 2.3 and 4.8 times greater, respectively, than the hazard for women with Stage 1 severity. Table 1 Modelling the length of time to first recurrence Variable Age x 10-2 Stage 2 Stage 3
Estimate - 0.61 0.82 1.57
S.E.
Multiplier
0.38 0.11 0.18
2.3 4.8
115
D.M. dos Santos et al./ Journal of Statistical Planning and Inlerence 47 (1995) 111 127
-7.25
-7.50
"...Weibull
-7.75
-8.00 N -8.25 ¢-.
'~.~eibulI-Normal
-8.50
-8.75
l ~iecewise
-9.00
Exp.
-9.25
-9.50 2.00
I
I
I
I
I
I
I
3.00
4.00
5.00
6.00
7.00
8.00
9.00
log t i m e to first r e c u r r e n c e
Fig. 1. Unconditional baseline log-hazards for time to first recurrence.
Fig. 1 shows the baseline log-hazard function for first recurrence. To smooth this plot, we took each tenth distinct recurrence time for the piecewise representation. It is evident that the unconditional log-hazard declines after an initial decrease. Superimposed on Fig. 1 is the estimated hazard of first recurrence from a conventional Weibull model fitted to the single duration data.
4. Modelling the frailty and hazard for repeated durations In this section we examine extensions of the conventional proportional hazards model, allowing for repeated durations and frailty. However, the baseline hazard is assumed to be identical for each duration.
116
D.M. dos Santos et al. / Journal of Stat&tical Planning and Inference 47 (1995) 111-127
Most frailty models considered in the literature have a proportional hazards specification with a scalar and time-invariant frailty variable. Lancaster (1979) was one of the first to generalise the proportional hazards model by introducing a positive multiplicative disturbance to represent unobserved heterogeneity. The corresponding model for repeated durations is given by the conditional hazard hq(tjlvi) = ho~(t~lvi) exp(fl'xu) vi,
(2)
where vi is the individual-specific frailty variable, independent of the exogenous explanatory variables. This formulation is consistent with the notion of frailty being due to omitted variables and therefore having an additive effect on the linear predictor: hu(tjlvi) = hoj(tj]vi) exp(fl'xij + el),
(3)
where vi = exp(e/). These formulations provide the basis of a 'mixture' model approach in which the nuisance parameter vi or e~ is integrated out of the likelihood for each case. The absence of appropriate computer software has prevented the widespread application of mixture models. Researchers have used G L I M macros, NAG subroutines and have written their own software, but these approaches are difficult to operationalise and tend to be remarkably slow. Most researchers have used a fully parametric specification (parametric distributions for the hazard and frailty) for model (2), to minimise the model fitting problems. In general, the gamma distributions is used for the frailty (e.g. Lancaster and Nickell, 1980; Manton et al., 1986; Follmann and Goldberg, 1988) and the Weibull for the hazard. These are conjugate distributions which lead to an analytically tractable model specification. Some researchers have reported the parameters of the frailty distribution to be very sensitive to the choice of the hazard function. We use the parsimonious and popular Weibull distribution initially for the conditional hazard as a base against which to compare more complex models. We report the results in this section. Some authors have used a nonparametric approach for the conditional hazard to avoid this problem; we consider this further in Section 5. 4.1. The W e i b u l l - g a m m a frailty m o d e l
Assume that recurrence times for the ith woman are drawn from a Weibull distribution with scale parameter vi and shape parameter ~. Assume further that the scale parameter vl represents a random draw from a gamma mixing distribution, g(vi) = (O~/F(7)) v~ 1 e x p ( - Ovi), where F(-) denotes the gamma function and that the recurrence times are assumed independent (conditional on v~). Then, the total
D.M. dos Santos et al./ Journal of Statistical Planning and In[erence 47 (1995) 111 127
117
likelihood equals the product of the contributions of all n cases, and is L(fl)= fif i=1
~2I vi~t'[k-aexp(flkx~k)exp( - ~'itik~ exp(flkXik) ) ] d'k k=l
0 7
× [ e x p ( - Viti~kexp(flkXik))] 1 -d,~
V~" I e x p ( - Ovi)dvi, F(7) '
(4)
where dik is the censoring indicator (equals 0 for censoring, 1 otherwise) and x~k denotes the vector of explanatory variables corresponding to recurrence time t~k, and K~ is the number of durations for case i. The intercept on the log linear scale corresponds to ln(#), where # = 7/0 is the mean of the mixing distribution. Eq. (4) simplifies to i=1
F(7)
~ •
exp •
\k
= 1
dikflkXik ~Jl 1~ t(,~ 1)d,k, (5) -
where di = }~'-ldik' and ~i, t (~ = ~K, (~ k = 1 tik exp(flkXik). If~ = 1, Eq. (4) is the likelihood for a gamma mixture of exponentials, whereas if the gamma mixing distribution is degenerate, (4) reduces to a Weibull distribution. Thus, Eq. (4) allows for the two extreme causes of apparent decreasing hazard rates: due solely to heterogeneity and due solely to temporal dependence. If K~ = l, (4) reduces to the model for single durations considered by Lancaster (1979) and Lancaster and Nickell (1980) amongst others. This, in turn, is a generalisation of the two parameter Burr (1942) distribution.
4.2!. The Weibull-normal frailty model In the Weibull-gamma model, the {vl} in Eq. (2) are considered to be gamma distributed. The Weibull-normal model is similar except that the {ai} in Eq. (3) are assumed to be normal distributed and, as this does not lead to an analytically tractable expression, quadrature methods are used to approximate the integral. The unconditional density function for single duration may be written as
.~(t) = Ifi(tl~i) deb(ei/w) mo
~ {et a 1exp(flo + fl'xi + W~k)exp[ -- t'exp(flo + fl'xl + W~k)]} pk,
(6)
k=l
where 4~(. ) is the standard normal distribution function, the {~, p} are the quadrature point locations and masses, respectively, mo is the number of quadrature points (taken to be 8 in the analysis we report later), rio is the location parameter and w the scale parameter for the normal mixing distribution.
118
D.M. dos Santos et al. / Journal of Statistical Planning and Inference 47 (1995) 111-127
Similarly, the unconditional survival function, Si, is given by trl O
Si(t) = ~ exp[--Hexp(flo + fl'xi + W~k)]Pk,
(7)
k=l
thus, the unconditional hazard is calculated as
jS(t) hi(t) = S,(t--)"
(8)
It must be noted that this unconditional hazard includes the sample selection bias effect whereby the less fail are progressively over represented with duration. These equations are readily extended to include multiple durations for each case. 4.3. The Weibull-nonparametric frailty model Evidence that results are not robust to choice of parametric form for the frailty distribution, at least for single duration data (see, in particular, Heckman and Singer, 1982), has stimulated interest in nonparametric representation of frailty (e.g. Laird, 1978; Heckman and Singer, 1984; Davies and Crouchley, 1984). The main approach uses the theoretical result, dating back to Kiefer and Wolfowitz (1956), that the nonparametric characterisation of the frailty distribution within maximum likelihood estimation is confined to a finite number of mass points. It also draws upon the empirical experience that the number of mass points tends to be small (Davies, 1993). For single duration data, the Weibull-nonparametric model may be written as mO
f ( t ) = ~ {~t " - t exp(fl'xi + ~ k ) e x p [ - - t ' e x p ( f l % + ~k)]} Pk. k=t
(9)
Note that this is algebraically similar to the Weibull-normal model with quadrature (Eq. (6)) but the m0 and {~,p} are no longer specified externally and have to be estimated. The location and scale parameters have been absorbed into the nonparametric representation. Apart from results on consistency, large sample theory has yet to be developed for the nonparametric frailty approach. However, limited simulation evidence (Heckman and Singer, 1982; Davies, 1987) suggests standard large sample properties and, in particular, provides some justification for the use of likelihood ratio tests and standard errors derived from the estimated information matrix. 4.4. Results on breast cancer recurrence The results of fitting the models discussed above to the breast cancer data are shown in Table 2. The software MIXTURE (Ezzet and Davies, 1988) was used to fit these models although a modified version was required for the Weibull-gamma. The Weibull-nonparametric model required just four mass points. Of course, the nonparametric characterisation of the frailty distribution may require further small mass
D.M. dos Santos et al. / Journal o f Statistical Planning and Inference 47 (1995) 111 127
119
Table 2 Conventional and mixture models Variable
Weibull
Weibull-gamma
Estimate S . E . Age x 10 2 Stage 2 Stage 3 State 2 State 3 [/o w 0 7 log-likelihood
- 0.76 0.66 1.07 1.53 1.81 0.75 -- 6.17
0,25 0,06 0.12 0.07 0.11 0.04 0.25
- 6786.91
Estimate S.E. - 1.30 1.20 1.84 1.25 1.47 1.05
0.40 0.11 0.14 0.08 0.10 0.04
1002.82 0.26 0.47 0.10 - 6729.79
Weibull-normal
Estimate S . E . 1.00 1.26 2.01 1.44 1.74 1.05 -- 8.74 1.39
0.50 0.15 0.33 0.09 0.14 0.04 0.46 0.12
- 6735.71
Weibullnonparametric Estimate S.E. - 1.40 1.17 2.39 1.29 1.47 1.07
0.50 0.16 0.37 0.1 I 0.15 0.04
6726.47
p o i n t s b u t these w o u l d n o t alter the l o g - l i k e l i h o o d a n d structural p a r a m e t e r estimates within the convergence criteria used (primarily, a r e d u c t i o n of less t h a n 10 _4 in the l o g - l i k e l i h o o d after a N e w t o n - R a p h s o n i t e r a t i o n with step length adjustment); the four mass p o i n t s p r o v e an a d e q u a t e if not c o m p l e t e n o n p a r a m e t r i c c h a r a c t e r i s a t i o n of the mixing distribution. T h e r e is an a p p r e c i a b l e increase in the l o g - l i k e l i h o o d when frailty is represented explicitly in the m o d e l (Table 2) and, a l t h o u g h formal l i k e l i h o o d ratio tests are not a p p r o p r i a t e , the n o n p a r a m e t r i c r e p r e s e n t a t i o n of frailty results in a n o t a b l e i m p r o v e m e n t u p o n the fully p a r a m e t r i c models. Theoretically, i g n o r i n g h e t e r o g e n e i t y leads to u n d e r e s t i m a t i o n of the effects of e x o g e n o u s e x p l a n a t o r y variables a n d their s t a n d a r d errors a n d also to s p u r i o u s r e d u c t i o n of the W e i b u l l p a r a m e t e r ( L a n c a s t e r a n d Nickell, 1980). This is consistent with the results for Age a n d Stage in T a b l e 2 a n d it a p p e a r s that the W e i b u l l m o d e l w i t h o u t heterogeneity p r o v i d e s a quite misleading picture of the h a z a r d declining with time; the h e t e r o g e n e o u s m o d e l s suggest that the c o n d i t i o n a l h a z a r d might actually increase, a l t h o u g h the results are n o t inconsistent with a c o n s t a n t hazard. The results for State are m o r e difficult to i n t e r p r e t because this is an e n d o g e n o u s r a t h e r than e x o g e n o u s e x p l a n a t o r y variable. T h e m a i n c o n s e q u e n c e of allowing for frailty is a r e d u c t i o n in the e s t i m a t e d h a z a r d following first a n d second recurrences. This is a s a m p l e selection effect; those who have experienced recurrence will tend to be those w h o are m o r e frail. W i t h o u t c o n t r o l for frailty, the e s t i m a t e d h a z a r d of further recurrence o r d e a t h will be b i a s e d u p w a r d . S u p e r i m p o s e d on Fig. 1 is the u n c o n d i t i o n a l e s t i m a t e d rate of first recurrence using the W e i b u l l - n o r m a l m o d e l (see Eq. (8)) with e x p l a n a t o r y variables set equal to zero. It is therefore a p a r a m e t r i c equivalent of the piecewise e x p o n e n t i a l in Fig. 1 a l t h o u g h it
120
D.M. dos Santos et al./ Journal of Statistical Planning and Inference 47 (1995) 111-127
has been obtained from the estimated conditional distribution rather than estimated directly from the data. It is noted that, whereas the conditional hazard for each woman increases over time in Weibull-normal model, the unconditional rate of recurrence shows the typical decline after an initial increase. This demonstrates the statistical effect whereby the less frail progressively dominate the survivors. Conversely, these results provide no evidence that the declines in the recurrence rate shown in Fig. 1 are attributable to any regenerative process.
5. Different baseline hazards for each recurrence type
In this section we allow the conditional baseline hazard to be different for each duration. Both parametric and nonparametric characterisation of the hazards are considered. A development due to Clayton (1988) enables estimation of models with nonparametric conditional hazards and multiple durations for each case. Frailty is assumed to have a gamma distribution. The approach is based upon a piecewise exponential formulation. As for the mass-point semiparametric approach, formal large sample theory has not been developed. But simulation evidence is encouraging with respect to consistency, likelihood ratio tests, and the normal distribution of parameter estimates (Nielsen et al., 1992) and Clayton and Cuzick (1985a) provide a heuristic justification for partial likelihood inference, using the results of Voelkel and Crowley (1984), if the model has a semi-Markov structure with ordered states such that transition from state k to state d is only possible if d > k. This condition is not overly restrictive for processes, such as many diseases, with an explicit progression. In particular, it enables the breast cancer data to be analysed with a different nonparametric hazard for each of the three ordered durations corresponding to State 1 (from primary treatment), State 2 (from first recurrence), and State 3 (from second recurrence). More specifically, the conditional hazard function for State j is given by hii(t~lvi) = hoj(tj[vi) exp(fl'xi~) vi,
(1o)
where tj is measured from initial treatment (j = 1) or when the disease last responded to treatment (j = 2 or 3), hot('l') is the conditional baseline hazard, and, in the Clayton method, {vj} are independent gamma variates with unit mean and variance tr 2. The problem is to estimate /3, tr 2, and the three unknown functions h0~('l-), j = 1,2,3. Epochs for the piecewise exponential formulation are constructed separately for each State. For example, the durations to first recurrence are divided into epochs on the basis of all the durations from initial treatment to first recurrence (or censoring) in the sample. Each baseline hazard h 0 / . ].) is represented in the partial likelihood by
D.M. dos Santo.; et al./ Journal of Statistical Planning and Inference 47 (t995) l l l 127
a step function
{hjk}
121
where k indexes epochs. The partial likelihood may be written
as
L(fl) =
lq [(Aijk)a'J"exp(--aijkUjk)] i=1
lexp(--Tvi)dvi ,
(11)
=
where Kij is the number of epochs for individual i in State j, Aijk = vihjkexp(fl'xO, U;k is the epoch duration, dljk is an outcome indicator variable (equals 1 if case i has a recurrence in epoch k of state j, and 0 otherwise), and 7 = 1/°-2, The GLIM macros given by Clayton (1988) use an EM algorithm to estimate the parameters in Eq. (11), estimating the v~ at the E step and the fl parameters and the hjk at the M step. To obtain a maximum likelihood solution it is necessary to repeat the algorithm for different values of cr2. Clayton emphasises that this profile log-likelihood approach is computationally simpler than seeking to estimate the variance simultaneously with the other parameters. Nielsen et al. (1992) propose a similar approach. The model fitting results for the Clayton approach including Age and Stage as explanatory variables are shown in Table 3. The piecewise conditional baseline log-hazards are plotted in Fig. 2(a), smoothed by taking every tenth distinct recurrence time for each State. It is evident that, at least for State 1 the hazard increases Table 3 Results of the semiparametric models Variable
Weibull non-parametric frailty
Clayton Estimate Agex Stage Stage State State
10 2 2 3 2 3
O-2
1.12 0.99 1.55
S.E.
Estimate
S.E.
0.28 0.07 0.13
-- 1.00 1.27 1.94 2.55 2.71
0.50 0.17 0.38 0.52 0.63
1.22 1.00 1.00
0.07 0.05 0.08
Estimate -- 1.00 1.14 1.80 3.47 4.05
S.E. 0.40 0.16 0.26 0.69 0.88
0.95
~2 ~3 71 × 103 ~'2 × 103 ]~3 × 103
-- 6.89 -- 8.08 -- 9.71
~2
G Pl ,o2 P3 P4
log-likelihood
Weibull Gompertz non-parametric frailty
6396.08
-- 13.28 0.07 0.26 0.33 0.34 -- 6723.34
1.54 1.22 1.14 1.02 1.19 0.88 8.69 1O.06 11.58
0.10 0.07 0.13 0.16 0.29 0.60
0.06 0.33 0.55 0.06 6692.24
122
D.M. dos Santos et al./ Journal of Statistical Planning and Inference 47 (1995) 111-127 (a)
Clayton approach
(b) Weii~ulFgompertz
-55
-55 ~=~..-.-! ....
i
i!;~"~ ~ i •
-60
lii - , I. ~ ~ i
- -
-
-.
STATE3
-6S
- -
STATE I
......
STATE 3
STARE 2 .........................
STATE t STATE 2
-'~
-60
~5
i.]
'"... -70
-7O
-75
-75
-~ O
~o
-aS
-8 5
-~0
-e 0
3OO
4OO
5OO
6OO
7OO
8OO
3o0
4o0
log time
6oo
600
700
a.o0
log lima
Fig. 2. Conditional baseline log-hazards. and then decreases. Moreover, the visual evidence suggests that the three hazards vary in shape. To fit a Weibull-nonparametric model with different shape parameters for each State we used N A G subroutines. The Weibull results, also shown in Table 3, are for four mass points because attempts at adding a fifth mass point did not improve the log-likelihood and the models showed exactly the results of the previous solution. The results indicate, in contrast to those from the Clayton method, increasing or constant hazards for the three states. The parameter estimates for Stage are about 25% larger than those for the Clayton model. The standard error estimates are also much larger for the Weibull model but we know that, because of the iterative algorithm, the Clayton model standard errors will tend to be underestimated (see, for example, Aitkin and Clayton, 1980). There is an improvement of only 3.13 in the log-likelihood with two additional parameters in comparison with the Weibull-nonparameteric model with identical baseline hazards for each recurrence (see Table 2).
6. Resolving the discrepancies The problem in reconciling the different conditional hazard and parameter estimate results for the Weibull-nonparametric and the Clayton nonparametric hazard-gamma models is that either or both may be misspecified. In particular, the Weibull cannot represent a nonmonotonic hazard and we have no assurance that the gamma is an appropriate distribution for the frailty (see Hougaard, 1986). Ideally, we would fit a model with nonparametric hazard and nonparametric frailty to provide a nested structure for the comparison but such models would require unrealistic computing resources and there are doubts about their identifiability. We therefore proceed in this
D.M. dos Santos et al./ Journal o f Statistical Planning and In/erence 47 (1995) 111 127
123
section to investigate the discrepancies between the Weibull-nonparametric and Clayton models by allowing more flexible parametric modelling of the hazards. The three parameter Weibull-Gompertz distribution (Lee, 1980) provides appreciable additional flexibility. The conditional baseline hazard in Eq. (10) is given by
hoj(tjlvi)
= c~jt~' l exp(~flj),
t~ >7 0,
(12)
with ~ j > 0 and - G o < ~ < ~ . The model reduces to the Weibull when 7 j = 0 . When ~j > 1, 7j > 0, ho~(tjlvi) is monotonic increasing; when ~j > 1, ~ < 0, hoj(tjlvi) increases and then decreases; when ~j < l, 7i > 0, hoj(tj[vi) decreases from infinity to a minimum and then increases to infinity; and when ~j < 1, ),j< 0, h0j(t~lv3 is monotonic decreasing. N A G optimisation subroutines were used to fit the Weibull-Gompertz model with nonparametric frailty. As the density and survival functions of the Weibull-Gompertz do not have explicit forms, they were approximated by Taylor series expansions with sufficient terms to ensure that the log-likelihood is accurate to at least the fourth decimal place. The survival function (dropping subscripts for simplicity) is
S(t)=exp
-~t ~
+
7t + + + --. c~ + 1 2!(~ + 2) 3!(c~+ 3)
.
(13)
The results are shown in Table 3. With ~ > 1 and 7 < 0 for all three hazards, these results are consistent with the pattern of increasing and then decreasing hazards suggested by the Clayton model. Also, the 31.10 increase in the log-likelihood with_ just three parameters compared with the Weibull-nonparametric (Table 3) model confirms that the Weibull hazard specification is untenable. Constraining the Weibull-Gompertz hazards to the same shape for each State reduces the log-likeli-. hood by 10.23 with four fewer parameters. It is now evident that the hazard functions do differ significantly. The estimated mass points give frailty variances of 2.33 and 1.74 for the Weibull and Weibull-Gompertz model, respectively, after scaling to a unit mean to ensure comparability with the Clayton approach. Unsurprisingly, the estimated frailty variance is reduced by improving the flexibility of the hazard specification. The baseline hazards for the Weibull-Gompertz model with nonparametric frailty are plotted in Fig. 2(b). To ensure comparability with the piecewise hazards in Fig. 2(a), they have been evaluated at the mean of the transformed mass point locations, exp(¢k). Thus the baseline hazard for State j is given by 3
ho~(tjlv,)
=
~ ~j ^ tj-~j 1 exp(~3jti + /~j)exp(~k),bk, k=l
where fli is the estimate of the intercept in the linear predictor for State j. Although they have the same 'sickle' shape, it is emphasised that these conditional hazards still differ markedly from the unconditional hazards. This is illustrated in Fig. 3 which compares the conditional hazard for first recurrence with the unconditional hazard estimated from the Weibull-Gompertz results (cf. Section 4.4 for the Weibull-normal
124
D.M. dos Santos et al./Journal of Statistical Planning and Inference 47 (1995) 111-127
-6.5
-7.0
. . .,. . .'"'"'°"""
..°.. .............. ...o... "'"'""-.,
.,°
-7.5 "O N m t-
-8.0
.go -8.5
-9.0
-9.5 I
I
I
I
I
3.00
4.00
5.00
6.00
7.00
8.00
log time Fig. 3. Estimated baseline hazards for first recurrence, Weibu]]-Gompertz.
model). As expected, the unconditional hazard is progressively lower than the conditional hazard. After 4 yr, the conditional hazard is more than twice the unconditional Moreover, the unconditional hazard begins to decline earlier (after less than 11 months compared with about 18 months for the conditional hazard); any regenerative process is operating more slowly than suggested by the unconditional h a z a r d
7. Concluding remarks The results in this paper have confirmed the importance of recognising frailty effects in any modelling of breast cancer recurrence and of distinguishing between conditional and unconditional hazards; the hazards may have different shape
D.M. dos Santos et al./ Journal ~?f Statistical Planning and ln[erence 47 (1995) 111 127
125
characteristics and be influenced differently by explanatory variables. The results also demonstrate that there is no technical impediment to the effective application of' semiparametric frailty models to epidemiological data. Flexible specification of hazard functions is a distinctive feature of the epidemiological tradition in survival analysis with, in particular, a strong emphasis upon nonparametric functions. This distinguishes it from the econometric tradition in which theoretical arguments are often adduced to justify the adoption of monotonic, parametric hazard functions. One constraint on including frailty effects within the epidemiological tradition is the empirical evidence of identifiability problems when a nonparametric specification is adopted for both the hazard and the frailty distribution; model fitting tends to converge to a degenerate, zero frailty solution even for simulated data with known frailty characteristics (Meyer, 1990). More seriously, Clayton (1988) has noted a similar tendency for nonparametric hazards with parametric frailty. This contrasts with the frailty models favoured by econometricians; empirical experience is encouraging and identifiability has been established under comparatively weak conditions for fully parametric models and models with parametric hazards and nonparametric frailty distributions, even for single duration data (Elbers and Ridder, 1982; Heckman and Singer, 1984). Results reported elsewhere using subsets of the breast cancer data (dos Santos, 1994) provide further evidence of the convergence to zero frailty variance with fully nonparametric specifications. However, no such problems were encountered in fitting the nonparametric hazard-parametric frailty models discussed in this paper. We would speculate that, at least for parametric frailty, this effect arises through model misspecification and/or inadequate sample sizes. Nevertheless, a more satisfactory compromise in many situations may be to use a parametric hazard and a nonparametric frailty distribution. The advantage of this approach is that it confines the nonparametric specification to the nuisance parameters; the frailty distribution is rarely of interest in its own right and, in any case, is difficult to estimate with precision by parametric (Nielsen eta., 1992) or nonparametric (Heckman and Singer, 1984) characterisations. Conversely, this approach enables conventional inferential methods to be applied to the characteristics of the hazards, to compare shapes, identify turning points, and so on. Concerns about the flexibility of the hazard specification can be addressed by using parametric forms which are more complex than those often recommended for survival analysis. The three parameter WeibuU Gompertz distribution used in this paper appears to have seen sufficiently flexible to represent the main features of the hazard for breast cancer recurrence.
Acknowledgements This research was supported, in part, by a grant from CAPES-Brazil, proc. 189/90-1 and the ESRC (Research project R000233850). The authors are also grateful to an anonymous referee for helpful comments.
126
D.M. dos Santos et al./ Journal of Statistical Planning and Inference 47 (1995) 111-127
References Aitkin, M. and D. Clayton (1980). The fitting of exponential, Weibull and extreme value distributions to complex censored survival data using GLIM. Appl. Statist. 29, 156-163. Bates, G. and J. Neyman (1952). Contributions to the theory of accident proneness II: true or false contagion. Univ. California Publ. Statist. 1, 255-275. Breslow, N. (1974). Covariance analysis of censored survival data. Biometric 30, 89-99. Burr, I.W. (1942). Cumulative frequency functions. Ann. of Math. Statist. 13, 215-232. Chamberlain, G. (1979). Heterogeneity, omitted variable bias and duration dependence. Harvard Institute of Economic Research, DP691. Cambridge, MA. Clayton, D.G. (1978). A model for association in bivariate life-tables and its application in epidemiological studies of familial tendency in chronic disease incidence. Biometrika 65, 141-151. Clayton, D.G. (1988). The analysis of event history data: a review of progress and outstanding problems. Statist. in Med. 7, 819-841. Clayton, D.G. and J. Cuzick (1985a). Multivariate generalizations of the proportional hazards model (with discussion). J. Roy. Statist. Soc. Set. A 148(2), 82-117. Clayton, D.G. and J. Cuzick (1985b). The semi-parametric pareto model for regression analysis of survival times. Bull. Inst. Statist. Inst. 51(23.3), 1 18. Cox, D.R. (1972). Regression models and life tables (with discussion). J. Roy. Statist. Soc. Ser. B 34, 187 220. Cox, D.R. (1975). Partial likelihood. Biometrika 62(2), 269-276. Davies, R.B. (1987). Mass point methods for dealing with nuisance parameters. In: R. Crouchley ed., Longitudinal Data Analysis, Aldershot, Avebury. Davies, R.B. (1993). Nonparametric control for residual heterogeneity in modelling recurrent behaviour. Comput. Statist. Data Anal. 16, 143-160. Davies, R.B. and R. Crouchley (1984). Calibrating longitudinal models of residential mobility and migration. Regional. Sci. Urban Econom. 14, 231-247. Davies, R.B. and A.R. Pickles (1985). Longitudinal versus cross-sectional methods for behavioural research: a first round knockout. Environ. Plann. A 17, 1315-1329. dos Santos, D.M. (1994). Distinguishing the effects of time dependence from the effects of frailty in multiple spell breast cancer data. Unpublished Ph.D. thesis, Lancaster University, Lancaster, UK. Elbers, C. and G. Ridder (1982). True and spurious duration dependence: the identifiability of the proportional hazard model. Rev. Econom. Stud. 49, 403-409. Ezzet, F.L. and R.B. Davies (1988). A manual for MIXTURE. Centre for Applied Statistics, University of Lancaster, Lancaster, England. Follmann, D.A. and M,S. Goldberg (1988). Distinguishing heterogeneity from decreasing hazard rates. Technometrics 30(4), 389-396. Francis, B., M. Green and C. Payne (1992). The GLIM System: Release 4 Manual. Oxford University Press, Oxford. Ginsberg, D.B. (1973). Stochastic models of residential geographic mobility for heterogeneous populations. Environ. Plann. 5, 113-124. Heckman, J.J. and B. Singer (1982). The identification problem in econometric models for duration data. In: W. Hildenbrand, Ed., Advances in Econometrics, Proc. World Meetings of the Econometric Soc 1980. Cambridge University Press, Cambridge. Heckman, J.J and B. Singer (1984). A method for minimizing the impact of distributional assumptions in econometric models for duration data. Econometriea 52(2), 271-320. Heckman, J.J. and B. Singer (1985). Social science duration analysis. In: J.J. Heckman and B. Singer, Eds., Longitudinal Analysis of Labor Market Data. Cambridge University Press, Cambridge, 39-110. Heckman, J.J. and R.J. Willis (1977). A beta-logistic model for the analysis of sequential labor force participation by married women. J. Political Economy 85, 27-58. Hougaard, P. (1984). Life table methods for heterogeneous populations: distributions describing the heterogeneity. Biometrika 71, 75 83. Hougaard, P. (1986). Survival models for heterogeneous populations derived from stable distributions. Biometrika 73(2), 387-396. Kiefer, J. and J. Wolfowitz (1956). Consistency of the maximum likelihood estimator in the presence of infinitely many incidental parameters. Ann. of Math. Statist. 27, 887-906.
D.M. dos Santos et al./ Journal of Statistical Planning and ln[erence 47 (1995) 11l 127
127
Laird, N. (1978). Nonparametric maximum likelihood estimation of a mixing distribution. J. Amer. Statist. Assoc. 73, 805-811. Lancaster, T. (1979). Econometric methods for the duration of unemployment. Econometrica 47, 939 956. Lancaster, T. (1985). Generalised residuals and heterogeneous duration models: with applications to the Weibull model. J. Econometrics 28(1), 155--169. Lancaster, T. and S. Nickell (1980). The analysis of re-employment probabilities for the unemployed. J. Roy Statist. Ser. Soc. A 143, 141 152. Lee, L. (1980). Testing adequacy of the Weibull and log linear rate models for a Poisson process Technometrics 22(2), 195 199. Manton, K.G., E. Stallard and J.W. Vaupel (1986). Alternative models for the heterogeneity of mortality risks among the aged. J. Amer. Statist. Assoc. 81,395, 635-644. Massy, W.F., D.M. Montgomery and D.G. Morrison (1970). Stochastic Models of Buying Behaviour. MIT Press, Cambridge, MA. Meyer, B.D. (1990). Unemployment insurance and unemployment spells. Econometrica 58(4), 757 782. NAG (1983). Numerical Algorithms Group Library Manual. Oxford, England. Nielsen, G.G., R.D. Gill, P.K. Andersen and T.I.A. Sorensen (1992). A counting process approach to maximum likelihood estimation in frailty models. Scand. J. Statist. 19, 25-43. Ribeiro, G.G. and R. Swindell (1981). The prognosis of breast carcinoma in women aged less than 40 years. Clin Radiol. 32, 231 236. Ridder, G. (1987). The sensitivity of duration models to misspecified unobserved heterogeneity and duration dependence. Working Paper, University of Amsterdam. Voelkel, J.G. and J. Crowley (1984). Non-parametric inference for a class of semi-markov models with censored observations. Ann. of Statist. 12, 142-160. Zeger, S.L. and M.R. Karim (1991). Generalized linear models with random effects, a Gibbs sampling approach. J. Amer. Statist. Assoc. 86, 413; Theory and methods, 79-86.