Indirect inference for survival data

Statistics & Operations Research Transactions SORT 27 (1) January-June 2003, 79-94 Statistics & Operations Research Transactions Indirect inference for survival data Bruce W. Turnbull*and Wenxin Jiang Cornell University Abstract In this paper we describe the so-called “indirect” method of inference, originally developed from the econometric literature, and apply it to survival analyses of two data sets with repeated events. This method is often more convenient computationally than maximum likelihood estimation when handling such model complexities as random effects and measurement error, for example; and it can also serve as a basis for robust inference with less stringent assumptions on the data generating mechanism. The first data set concerns recurrence times of mammary tumors in rats and is modeled using a Poisson process model with covariates and frailties. The second data set involves times of recurrences of skin tumors in individual patients in a clinical trial. The methodology is applied in both parametric and semi-parametric regression analyses to accommodate random effects and covariate measurement error. MSC: 62N01, 62N02, 62P10 Keywords: Estimating equations, frailty, hazard rate regression, indirect inference, measurement error, naive estimators, overdispersion, quasi-likelihood, random effects, robustness 1 Introduction Methods of indirect inference (Gourieroux, Monfort and Renault, 1993) have been developed and used in the field of econometrics where they have proved valuable for parameter estimation in highly complex models. This paper recasts the basic technique in a likelihood-flavoured approach and illustrates some applications in biostatistics, in particular for survival and repeated events data. We begin by illustrating the steps involved in the indirect method in the following simple pedagogic example. ∗Address for correspondence: B. W. Turnbull. School of Operations Research and Department of Statistical Science, 227 Rhodes Hall, Cornell University, Ithaca, NY 14853-3801, USA. turnb[email protected] Received: October 2002 Accepted: January 2003 80 Indirect inference for survival data Example 1: exponential survival with censoring. Consider lifetimes {T1,...,Tn}, which are independent and identically distributed (i.i.d.) according to an exponential distribution with mean θ. The data are subject to Type I single censoring after fixed time c. Thus the observed data are {Y1,...,Yn}, where Yi=min(Ti,c),(i=1,...,n). We consider indirect inference based on the intermediate statistic ˆs=Y. This choice can be considered either as the basis for a method of moments estimator or as the MLE (maximum likelihood estimator) for a misspecified model M0in which the presence of censoring has been ignored. The naive estimator Yin fact consistently estimates not θ but the “naive” or “auxiliary” parameter s(θ) = θ[1−exp(−c/θ)],(1) the expectation of Y. The equation (1) is an example of what may be termed as a “bridge relation” (Jiang and Turnbull, 2001) or a “binding relation” (Gourieroux, et al., 1993). We can see the obvious effect of the misspecification, namely that ˆsunderestimates θ. However a consistent estimator ˆ θof θas n∞can be obtained by solving (1) for θwith s(θ)replaced by ˆs=Y. That is, ˆ θ=s−1(ˆs). (Note that s(·)is strictly increasing on ℜ+and thus invertible). We note also that ˆ θis not the MLE of θwhich is nY /[∑n i=1I(Yi<c)].) More generally, a consistent estimator can be constructed based on an intermediate statistic ˆsthat does not need to have the interpretation of a ‘naive’ estimator. For example, above we could have chosen perhaps ˆs=Y2=n−1∑n i=1Y2 iso that ˆ θ=s−1(ˆs)where s−1 is the inverse function of s(θ)≡E(Y2|θ). In fact,the dimension of ˆscan be greater than that of θ— e.g., we could take ˆs= (Y,Y2)Tin the above example. Now a consistent estimator ˆ θof θcan be found by using weighted least squares, ˆ θ=argmin θ{ˆs−s(θ)}TA{ˆs−s(θ)}, where an optimal choice of Ais the inverse of the estimated variance matrix of ˆs, as we will discuss later. This is the principal idea of indirect inference– statistical inference of θbased on an indirect data “summary” ˆs. The choice of ˆsis not unique, but in most applications there will natural one to use as we shall see. We will term ˆ θas the “indirect MLE”, since it can be viewed as the MLE using an approximate likelihood based on the indirect data summary ˆs. We will also see how to obtain the standard error for ˆ θ. 2 Indirect inference In general, the indirect MLE has properties similar to those of the usual MLE: consistency, asymptotic normality, and certain efficiency properties. In addition, chisquared goodness-of-fit tests can be based on the indirect likelihood. Bruce W. Turnbull and Wenxin Jiang 81 Advantages of the indirect method include ease of computation; robustness; and informativeness on the effect of model misspecification. We will summarize the framework of indirect inference below. 2.1 The basic approach Suppose we have a data set consisting of nindependent units. The essential ingredients of the indirect approach, when reformulated in a likelihood-flavoured treatment, are as follows. •There is a hypothesized true model M for data generation, with distribution P(θ) which depends on an unknown parameter θof interest which is of dimension p. •One first computes an intermediate or auxiliary statistic ˆsof dimension q≥p, which is asymptotically normal with mean s(θ), say, under model M. •An indirect likelihood L(θ|ˆs)is then constructed based on the normal approximation, so that, apart from an additive constant, −2logL(θ|ˆs) = {ˆs−s(θ)}Tv−1{ˆs−s(θ)}=H(θ),say,(2) where vis a consistent estimate of the asymptotic variance cvar(ˆs). A typical choice might be the ‘robust’ or ‘sandwich formula’, when ˆssolves an estimating equation (see e.g. Carroll, Ruppert and Stefanski 1995, Section A.3). •This indirect likelihood is then maximized to generate an indirect maximum likelihood estimate (indirect MLE) or adjusted estimate ˆ θ(ˆs)for θ. In the case when the dimension (q) of the intermediate statistic equals that (p) of the parameter θ and s(θ)is invertible, it can be seen from (2) that maximization of the indirect likelihood is equivalent to solving the “bridge” or “binding” equation s(θ) = ˆsfor θ, because then (2) can be made zero. In the “indirect” analysis of the pedagogic example of Section 1, M is the i.i.d. exponential model with censoring: Yi=min(Ti,c)and P(θ)(Ti≤t) = 1−e−t/θ, for t∈[0,∞),i=1,...,n. In the initial approach, the intermediate statistic was ˆs= n−1∑n i=1Yi, which is asymptotically normal as n−∞by the central limit theorem, with asymptotic mean s(θ) = θ[1−exp(−c/θ)]. The indirect likelihood L(θ|ˆs)is given by −2logL(θ|ˆs) = {n−1 n ∑ i=1 Yi−s(θ)}Tv−1{n−1 n ∑ i=1 Yi−s(θ)} where one can substitute the robust estimate ˆv={n(n−1)}−1∑n i=1(Yi−Y)2for the asymptotic variance v=var(ˆs). Finally the adjusted estimate (or indirect MLE) is ˆ θ=s−1(ˆs). 82 Indirect inference for survival data In summary, in this indirect approach, the data are first summarized by the intermediate statistic ˆs. Its asymptotic mean sis referred to as the auxiliary parameter. The auxiliary parameter is related to the original parameter by a relation s=s(θ), termed the bridge relation or binding function. The starting point is the choice of an intermediate statistic ˆs. This can be chosen as some set of sample moments, or the solution of some estimating equations, or the MLE based on some convenient model M0, say, termed the auxiliary (or naive) model. If the last, then the model M0is a simpler but misspecified or partially misspecified model. As stated previously, the choice of an intermediate statistic ˆsis not necessarily unique; however in any given situation there is often a natural one to use. 2.2 Intermediate statistics arising from estimating equations Most intermediate statistics can be defined implicitly as a solution, s=ˆs, of a (qdimensional) estimating equation of the form G(W,s) = 0, say. [Clearly this includes any statistic ˆs=ˆs(W)that has an explicit expression as a special case, by taking G=s−ˆs(W).] The estimating equation could be the normal equation from a leastsquares analysis, or the score equation based on some likelihood function. In such situations there is a parallel formulation of indirect inference in ‘implicit form’. For instance, one can state the ‘bridge relation’ s(θ)implicitly as F(θ,s) = 0 where F(θ,s)≡EW|θG(W,s), which is the limiting version of the estimating equation G(W,ˆs) = 0. Correspondingly, in the definition of the indirect likelihood L,Hcan be (asymptotically) equivalently defined by H(θ,ˆs) = F(θ,ˆs)Tv−1F(θ,ˆs). Here vis (a sample estimate of) the avar of F(θ,ˆs), which can be evaluated by the delta method (e.g. Bickel and Doksum (2001), Sec. 5.3.2), and found to be the same as var(G)evaluated at s=s(θ)(the auxiliary parameter). Then we define the adjusted estimator (or the indirect MLE)ˆ θto be the maximizer of L, or the minimizer of H. 2.3 Properties of indirect MLE In general, the indirect MLE has a set of properties analogous to those of the usual MLE. These include, under appropriate regularity conditions: (i) (Indirect Score Function). The asymptotic mean and variance of the indirect likelihood score function satisfy the usual relations E(∇θlogL) = 0 and var(∇θlogL)+ E(∇2 θlogL) = 0. (ii) (Asymptotic Normality). The adjusted estimator ˆ θis asymptotically normal (AN) with mean θ, and with asymptotic variance (avar) estimated by −(∇2 θlogL)−1or 2(∇2 θH)−1where consistent estimates are substituted for parameter values. Bruce W. Turnbull and Wenxin Jiang 83 (iii) (Tests). Likelihood-ratio statistics based on the indirect likelihood for testing simple and composite null hypotheses have the usual asymptotic χ2distributions. (iv) (Efficient use of indirect data). The adjusted estimator has smallest avar among all consistent AN estimators f(ˆs)of θ, which are constructed from the naive estimator ˆsby continuously differentiable mappings f. These results can be found in the references cited in Section 2.5, and are summarized in Jiang and Turnbull (2001, Proposition 1). When different intermediate statistics are used, the asymptotic efficiency can be different. In general the indirect MLE is not as efficient as the MLE based on the true model M; although there are situations that can be identified where the efficiency will be high, as in the example of Section 3.1. In a special case when the dimension of the intermediate statistic (q) equals that (p) of the parameter θ, and s(·)is a diffeomorphism on the parameter space Θof θ, maximization of Lis equivalent to the bias correction ˆ θ=s−1(ˆs)(from solving F(θ,ˆs) = 0), which is AN and consistent for θ. See, e.g., Kuk (1995), Turnbull et al. (1997) and Jiang et al. (1999) for biostatistical applications. When q<p, there are more unknown true parameters than ‘naive parameters’. In this case the bridge relation is many-to-one and does not in general permit the construction of adjusted estimates. It is mainly of interest for investigating the effects of misspecification when the naive estimators are constructed under misspecified models. However, in such situations it may be possible to construct consistent estimates for a subset of true parameters, which may be of interest. In other situations, some components of the higher-dimensional true parameter are known or can be estimated from other outside data sources. This enables the other components to be consistently estimated by inverting the bridge relation. Examples of this kind arising from errors-in-variables regression models are given in Sections 3.2 and 3.3. 2.4 Why consider the indirect method? This indirect approach offers the following advantages: 1. Ease of computation. The indirect method is typically computationally simpler and more convenient. For example, when ˆsis based on some simplified model M0, it can often be computed with available standard computer software. 2. Informativeness on the effect of model misspecification. When ˆsis a ‘naive estimate’ obtained from a naive model M0neglecting certain model complexities, the approach is very informative on the effect of model misspecification — the bridge relation s=s(θ)provides a dynamic correspondence between M0and M. For example, in errors-in-variable regression, such a relation is sometimes termed an 84 Indirect inference for survival data ‘attenuation relation’ (see e.g., Carroll, Ruppert and Stefanski 1995, Chapter 2), and tells how regression coefficients can be underestimated when neglecting the measurement error in a predictor. 3. Robustness. The validity of the inference based on an intermediate statistic essentially relies on the correct specification of its asymptotic mean. This is typically a less demanding assumption than the correct specification of a full probability model, which would be generally needed for a direct likelihood inference to be valid. Therefore inferences based on the adjusted estimate ˆ θcan remain valid despite some departure of the data generation mechanism from the hypothesized true model M. 2.5 Bibliography and notes The above very brief exposition of the indirect method of inference represents a summary of results that have appeared in the econometric and statistical literature in varying forms and generality and tailored for various applications. Examples include: the generalized method of moments (GMM: Hansen 1982); the method of linear forms and minimum χ2(Ferguson 1958); the regular best asymptotic normal estimates that are functions of sample averages (Chiang 1956, Theorem 3); simulated method of moments and indirect inference [McFadden (1989), Pakes and Pollard (1989), Gourieroux et al. (1993), Gallant and Tauchen (1996, 1999) Gallant and Long (1997)]. Newey and McFadden (1994, Chapters 6 and 8) discuss two-stage parametric and nonparametric estimation in the GMM context, where some ‘nuisance’ parameter, possibly infinite dimensional, is estimated from a preliminary consistent method. Applications of GMM in the settings of generalized estimating equations from biostatistics are discussed in Qu, Lindsay and Li (2000). McCullagh and Nelder (1989, p. 341), as referred to by Qin and Lawless (1994, p. 315), consider optimal linear combination of estimating equations, as is traditionally done in GMM literature. Qin and Lawless (1994) also provide an alternative but asymptotically equivalent way of combining estimating equations using empirical likelihood. The theory of estimators obtained from misspecified likelihoods goes back at least as far as Cox (1962), Berk (1966) and Huber (1967) and is summarized in the comprehensive monograph by White (1994). The use of ˆs(based on an auxiliary model M0) in indirect inference about θ(under model M) appears recently in the field of econometrics to treat complex time series and dynamic models, see, e.g., Gourieroux et al. (1993) and Gallant and Tauchen (1996, 1999); as well as in the field of biostatistics to treat regression models with random effects and measurement error, see e.g., Kuk (1995), Turnbull, Jiang and Clark (1997), and Jiang et al. (1999). This bibliography is far from exhaustive. A thorough review and synthesis of the methods of indirect inference are given in Jiang and Turnbull (2001). Bruce W. Turnbull and Wenxin Jiang 85 3 Three applications with survival data In this section we discuss three applications with recurrent event data which use models of increasing order of complexity: 3.1 A Poisson process regression model with random effects (“frailties” or “unexplained heterogeneity”); 3.2 A Poisson process regression model with random effects and covariate measurement error; 3.3 A semi-parametric intensity rate regression model with random effects and measurement error. The first example uses mammary tumor recurrence times from a rodent carcinogenicity experiment. The remaining two examples use data on skin cancer recurrences in the Nutritional Prevention of Cancer (NPC) trial — a long-term randomized clinical trial for cancer prevention (Clark et al. 1996). 3.1 Animal carcinogenicity data: multiple times to tumor Gail, Santner and Brown (1980, Table 1) present data on multiple mammary tumor incidence times from an experiment conducted by Thompson et al. (1978). Forty-eight female rats which remained tumor-free after sixty days of pre-treatment of a prevention drug (retinyl acetate) were randomized with equal probability into two groups. In Group 1 they continued to receive treatment (Z=1), in Group 2 they received placebo (Z=0). All rats were followed for an additional 122 days and the time of any newly diagnosed mammary tumor was recorded. The numbers of tumors diagnosed in individual rats ranged from 0 to 13. The objective of the study was to estimate the effect of the preventive treatment (Z) on tumor recurrence. Suppose we consider a model in which the tumors occur over time in a given subject (rat) according to a Poisson process with a constant intensity rate which depends on treatment Z, a fixed effect, and on subject, a random effect. If we define Yto be the number of tumors diagnosed in a particular rat during the 122 day followup time, the model M specifies that, given Zand ,Yis Poisson distributed with mean exp(α+Zβ+). Here the assigned treatment Zis observed, but represents an unobserved random effect modeled as normally distributed with zero mean and constant variance σ2, independent of Z. This random effect or “unexplained heterogeneity” could be considered to be caused by omitted covariates. We observe n=48 i.i.d. pairs Wi= (Yi,Zi),i=1,...,n. The likelihood for the observed data involves integration over and is difficult to compute. (However it is possible – see below.) Instead, we start by taking the indirect approach with an auxiliary statistic ˆs= (ˆa,ˆ b,ˆ t2)T, where 86 Indirect inference for survival data (ˆa,ˆ b)are the regression coefficient estimates maximizing a naive log-likelihood R= ∑n 1{Yi(a+Zib)−ea+Zib}, and ˆ t2=n−1∑n i=1Y2 iis the second sample moment. Here the auxiliary parameter is s=plim(ˆa,ˆ b,ˆ t2)T, whereas the true parameter to be estimated is θ= (α,β,σ2)T. The use of the naive log-likelihood Rcorresponds to a simplified model M0in which the presence of the random effect is neglected. The second sample moment is included in the intermediate statistic to provide information for estimation of the variance parameter. Therefore ˆsis solved from the estimating equation G(W,s) = 0, where (formally) G= (n−1∂aR,n−1∂bR,ˆ t2−t2)T, i.e. G=n−1n ∑ i=1 (Yi−ea+Zib,Zi(Yi−ea+Zib),Y2 i−t2)T=n−1n ∑ i=1 gi,say. The solution ˆs= ( ˆa,ˆ b,ˆ t2)Tcan be computed easily. For the rat carcinogenicity data we obtain the auxiliary estimates ˆa=1.7984; ˆ b=−0.8230; ˆ t2=31.875. The asymptotic variance var(ˆs)can be estimated by the sandwich formula (see e.g. Carroll, Ruppert and Stefanski 1995, Section A.3) v= (∇sG)−1cvar(G)(∇sG)−T|s=ˆs where cvar(G) = n−2∑n i=1gigT i|s=ˆs,∇sGis a 3×3 matrix with elements (∇sG)jk =∂skGj, j,k=1,2,3, and A−T= (A−1)Tfor a generic matrix A. The indirect likelihood L(θ|ˆs), up to an additive constant, satisfies −2logL(θ|ˆs) = {ˆs−s(θ)}Tv−1{ˆs−s(θ)}, where s(θ)is the asymptotic mean or large sample almost sure limit of ˆs. Since ˆssolves the estimating equation G=0, its limit is the solution of the limiting estimating equation F(θ,s) = EW|θG(W,s) = 0, which can be explicitly solved to obtain s=s(θ). This yields the bridge equation: s=plim(ˆa,ˆ b,ˆ t2)T =s(θ) =  α+σ2/2,β,1 2(1+eβ)eα+1 2σ2 +1 2(1+e2β)e2(α+σ2)  T . Because dim(s) = dim(θ) = 3 and s(θ)is a smooth invertible mapping, the indirect MLE ˆ θ=argmaxθL(θ|ˆs)can be obtained by solving ˆs=s(θ), which gives the adjusted estimates ˆ θ= (ˆα,ˆ β,ˆσ2) = s−1(ˆs). Thus ˆ β=ˆ b, and ˆα=ˆa−ˆσ2/2 where ˆσ2= log2ˆ t2−eˆa(1+eˆ b) e2 ˆa(1+e2ˆ b). For the rat data, this leads to adjusted estimates ˆα=1.6808(0.1589); ˆ β=−0.8230(0.1968); ˆσ=0.4850(0.1274). The estimated standard errors shown in parentheses are obtained using the delta method formula: cvar(ˆ θ) = (∇θs)−1v(∇θs)−T|θ=ˆ θ, and then taking the square roots of the 3 diagonal elements of this matrix. It is noted that this delta method expression is equivalent to deriving the variance by Bruce W. Turnbull and Wenxin Jiang 87 {−∇2 θlogL(θ|ˆs)}−1|θ=ˆ θ based on the ‘indirect likelihood Fisher information’, where ∇2 θrepresents the Hessian. This follows because the jkth element of the Hessian is, for j,k=1,2,3, {−∇2 θlogL(θ|ˆs)}jk|θ=ˆ θ=−∂θj∂θklogL(θ|ˆs)|θ=ˆ θ = (∂θjsT)v−1(∂θks)|θ=ˆ θ−(∂θj∂θksT)v−1(ˆs−s)|θ=ˆ θ = (∂θjsT)v−1(∂θks)|θ=ˆ θ+0 ={cvar(ˆ θ)−1}jk. If we wish to obtain the MLE of θ= (α,β,σ2)based on model M, then it can be found by a somewhat tedious iterative numerical maximization of the true likelihood which involves numerical integration over the distribution of . These estimates are: ˆαML =1.6717 (0.1560);ˆ βML =−0.8125 (0.2078); ˆσML =0.5034 (0.0859). For the MLEs, the estimated standard errors are based on the inverse of the Fisher information matrix, evaluated at the corresponding estimate values. The estimated standard errors suggest that the efficiency of indirect estimation of the treatment effect parameter βis high here in this example. Related results (Cox, 1983; Jiang et al., 1999) show that such high efficiency is achievable if the follow-up times are about the same across different subjects (which is true here), or if the overdispersion is small. Also it should be noted that the adjusted estimator ˆ βis robust, in the sense that it remains consistent, essentially as long as the mean function E(Y|Z,)is correctly specified and and Zare independent. (Its standard error estimate from the sandwich formula is also model-independent and robust.) In particular, the consistency property does not depend on the specification of a complete probability model, namely that Yis Poisson and is normal. Thus the indirect estimator enjoys a robustness advantage over the MLE. The indirect approach, although formulated from the different perspective of using naive model plus method of moments, is intimately related to the work of Breslow (1990) based on quasi-likelihood and method of moments. Breslow used a different linear combination of Yi’s based on quasi-likelihood (Wedderburn, 1974; McCullagh and Nelder, 1989), which enjoy general efficiency properties among linear estimating equations. However, (i) our approach can be interpreted as basing inference on the simple moments n−1∑Yi,n−1∑ZiYiand n−1∑Y2 i(which can be easily seen from the estimating equation G=0), and (ii) our approach shows clearly, by the use of bridge relations, the sensitivity and robustness of parameter estimates to the omission of over-dispersion in modeling. Also note that here we used a log-normal distribution to model the random effects and the variance parameter also enters the mean model (unconditional on ), whereas Breslow (1990) focused on the examples such as ones with gamma multiplicative random effects in which the mean model does not change.