Estimating medical costs from a transition model

IMS Collectio ns Beyond P arametr ics in Interdisciplinary Resear c h: F estsc hrift in Honor of Professor Pranab K. Sen V ol. 1 (2008) 350– 363 c  Institute of Mathematical Stat istics , 2008 DOI: 10.1214/ 193940307 000000266 Estimating medical cos ts from a transition mo del Joseph C. Gardiner 1 , Lin Liu 2 and Zheh ui Luo 3 Michigan State Unive rsi t y Abstract: Nonparametric estimators of the mean total cost hav e b een pro- posed in a v ar i et y o f se ttings. In clini cal trials it is generally impractical to follow up patien ts until all hav e resp onded, and therefore censoring of patien t outcomes and total cost will o ccur in practice. W e describ e a general longi- tudinal fr amew ork in which costs emanate f rom tw o streams, during so journ in health states and in transition fr om one health state to another. W e con- sider estimation of net present v alue f or expenditures i ncurred o ve r a finite time horizon f rom medical cost data that might b e incompletely ascertained in some patien ts. Because patien t s pecific demographic and clini cal c harac- teristics would influence total cost, we use a regression mo del to incorp orate co v ariates. W e discuss similar ities and differences betw een our net present v alue estimator and other widely used estimators of total medical costs. Our mo del can accommodate heterosceda sticity , ske wness and censoring in cost data and prov ides a flexible approac h to analyses of health care cost. 1. In tro duction Estimating cost from medical follo w-up studies has b een the focus of extensive metho dological research. Cost data in observ ationa l studies exhibit sev eral features such as heteroscedas ticit y , skewness and censor ing tha t must b e addres sed in sta- tistical modeling so that ensuing inference would b e v a lid. In clinical trials it is generally impra ctical to prolong a study until all patients hav e resp onded, and therefore inevitably censor ing of patient outcomes and total c o st will o ccur in prac- tice. Since costs are incurred over time, the cum ulativ e cost C ( t ) at time t is a nonnegative monotone function. Co st a c cum ulation ends at an ev en t time T , for example at dea th for lifetime cost, or at a sp ecified finite time ho rizon τ . Interest lies in estimating the mea n cost µ = E ( C ( T ∗ )) where T ∗ = min( T , τ ). Bec a use T could be pr ecluded from obser v ation by cens oring at time U , that is, when T > U , the cor resp onding cost would b e complete only if U ≥ T ∗ . Several nonpara metric estimators o f µ hav e b een prop osed in a v ariety of settings with regressio n mo d- els b eing the mainstay for assess ing the influence of patient-specific characteristics (eg, tr eatments, demo graphics, como rbidity) on co st (for example, Bang and Ts i- ∗ Supported by the Agency for Healthcare Research & Qualit y under gran t 1R01 HS14206. 1 Division of Bi ostatistics, Department of Epidemiology , B629 W est F ee Hall , Michigan Stat e Unive rsity , East Lansing, MI 48824, USA, e-mail : jgardine r@epi.msu .edu 2 Now at Eli Lil ly and Company , US Medical Division, Indianap olis, Indiana 46285, USA, e-mail: liu lin ll@lilly. com 3 Now at R TI In ternational, Beha vi oral Health Economics Program, 3040 Corn wa llis Rd., Re- searc h T r iangle Park, NC 27709, USA, e-mail: zluo@rti .org AMS 2000 subje ct classific ati ons: Pr i mary 62N01, 60J27 ; secondary 62G05 . Keywor ds and phr ases: censoring, Kaplan-Meier estimator, longitudinal data, M arko v mo del, inv erse-weigh ting, random-effect s. 350 Estimating me dic al c osts 351 atis [ 2 , 3 ], Ba ser et al. [ 4 , 5 ], Lin [ 1 3 , 14 ], Lin et al. [ 15 ], O’Hag an and Stevens [ 17 ], Strawderman [ 18 ] a nd Gardiner et al. [ 8 ]). This article adopts a broader view o f the cumulative co st { C ( t ) : 0 ≤ t ≤ τ } within the framework of a long itudina l mo del. Section 2 descr ib es a ll the substan- tive a sp ects of our models star ting with an underlying finite state sto chastic pro cess for the evolution of patient even ts as they o ccur over time. The sta tes ar e different health conditions that the patien t presents ov er the p erio d [0 , τ ]. Costs emanate from t w o streams , during so journ in health states and in transitio n from one health state to another. W e conside r e stimation of net present v alue (NPV) for e x pe ndi- tures incurred ov er [0 , τ ]. Regressio n mo dels for the event histor y pro ces s and for observed c osts are us ed to incorp ora te cov a riates. Section 3 outlines the metho d of estimation o f NP V from a pa tient sample of time-ce ns ored even t history data. W e then discuss similarities a nd differe nc e s b etw een our net present v alue estima- tor and other widely used estimator s of total medical co sts. Section 4 is a brief summary and co nclusion. 2. Sto c hastic mo del 2.1. T r ansi tion and sojour n c ost A sto chastic pro ces s X = { X ( t ) : t ∈ T } on the interv al T = [0 , τ ] where τ < ∞ , describ es the health states of a patient from the relev a nt p o pulation under study . The time τ is the maximum limit of observ ation for all co st and patient outcomes. The state spa ce o f X is finite and labeled E= { 0 , . . . , m } and co nsists of several transient states , such as “well”, “r ecov ery”, “ relapse” , and one or mo re absorbing states such a s “dead” or “ disabled”. A transient state is one whic h if visited w ill be exited after a finite so journ, wherea s a transition o ut of a n absorbing state is impo ssible. Costs ar e incurred while so jour ning in a transient health state and in transition b etw een states. If the pa tient is in sta te h at time t , that is, X ( t ) = h , the expe nditur e ra te is B ( t, h ). If a tra nsition o cc ur s fro m state h to state j a t time t , that is, X ( t − ) = h and X ( t ) = j, a cost C ( t, h, j ) is incurred. The nota tio n [ A ] denotes the indicator function of the even t A taking v alue 1 if A is true and 0 if A is false. F or example, to indicate the state of o ccupation just prior to time t w e write Y h ( t ) = [ X ( t − ) = h ]. The num b er of direct transitions h → j , h 6 = j in the time in terv al [0 , t ] is N hj ( t ) = # { s ≤ t : X ( s − ) = h, X ( s ) = j } . If r is the discount rate, the present v alue of exp enditure s asso cia ted with all h → j transitions in T is (1) C (1) hj = Z τ 0 e − r t C ( t, h, j ) dN hj ( t ) , and the pres ent v alue of exp enditures fo r all so journs in state h in T is (2) C (2) h = Z τ 0 e − r t B ( t, h ) Y h ( t ) dt. W e will in terpret all integrals as on the semi-op en in terv al (0 , τ ]. In practice we wan t to es timate the ex pe c ted v a lues (av erages ) of thes e tw o quantities. T o do so we imp ose a Mar ko v mo del on X to gov ern the tra nsitions b etw een states. 352 J. C. Gar diner, L. L iu and Z. Luo 2.2. Markov mo del W e ca ll X a n on-homo gene ous Markov pr o c ess if P [ X ( t ) = j | X ( s ) = h, X ( u ) : u < s ] = P [ X ( t ) = j | X ( s ) = h ] for all h, j ∈ E and a ll s ≤ t . The tr ansition pr ob abilities P hj ( s, t ), s ≤ t , of X a r e given by P hj ( s, t ) = P [ X ( t ) = j | X ( s ) = h ] and the tr ansition intensities α hj ( t ) b y α hj ( t ) = lim ∆ t ↓ 0 P [ X ( t + ∆ t ) = j | X ( t ) = h ] / ∆ t, j 6 = h with α hh = − P j 6 = h α hj . Throughout we a ssume that the α hj are integrable o n T = [0 , τ ]. The m × m matrices P = { P hj ( s, t ) } a nd α = { α hj } are related by the pro duct-integral formula P ( s, t ) = Q s 0 : X ( t ) = m } is the time to absor ption. The surviv al distr ibution, c o nditional on X (0) = i , is S mi ( t | z ) = P [ τ m > t | X (0) = i, z ] = 1 − P im (0 , t | z ) , and the unconditiona l s urviv a l distr ibution is S m ( t | z ) = 1 − X i 6 = m π i (0 | z ) P im (0 , t | z ) . In the sp ecial cas e o f one transient sta te 0 (“alive”) and one termina l state 1 (“dead”) w e g et the usual surviv al time T (= τ 1 ) and its surviv a l distribution S ( t | z ) = P [ T > t | z ] = P 00 (0 , t | z ). 3. Estimation Suppo se we observe the a forementioned pro cess es fo r ea ch o f n s ub jects in a lo n- gitudinal study . F or the i -th patient the basic cov aria te vector is z i ( t ), the initial state X i (0), the state indicator Y hi ( t ) = [ X i ( t − ) = h, U i ≥ t ] and the num b er of direct h → j transitions N hj i ( t ) = # { u ≤ t ∧ U i : X i ( u − ) = h, X i ( u ) = j } , h 6 = j. Conditionally on { z i (0) , X i (0) : 1 ≤ i ≤ n } a ssume pr o cesses { X i ( t ) : t ∈ T } are independent and that mo del ( 3 ) holds for ea ch individual with the same bas e line int ensities. F rom now on denote by N hj ( t ) and Y h ( t ), resp ectively , the aggr egated pro cesses P n i =1 N hj i ( t ) a nd P n i =1 Y hi ( t ). In this context estimation of the tr ansi- tion proba bilities P hj (0 , t | z ) and integrated in tensities A hj ( t | z ) at a fixed cov ariate profile z is well k nown (Andersen et al. [ 1 ]). Combining this with appropr iate esti- mation of costs would lead to estimator s of NPV. How ever, b efore we describ e an approach to estimation we first consider several ex amples. 3.1. Single tr ansition without c ovariates The only p ermissible transition 0 → 1 is asso ciated with a s ingle cost C ( T , 0 , 1) (denoted here by y ) where T deno tes the surviv al time. F rom ( 6 ) we hav e NPV = Z τ 0 e − r t c 01 ( t ) P 00 (0 , t − ) dA 01 ( t ) . T o estimate NPV we use the estimato rs ˆ P 00 (0 , t − ) = ˆ S ( t − ) and d ˆ A 01 ( t ) = { Y 0 ( t ) } − 1 dN 01 ( t ) , where ˆ S is the K aplan-Meier estimator of the sur viv al distribution o f the surviv al time T (= τ 1 ), Y 0 ( t ) = P n i =1 [ T i ∧ U i ≥ t ] and N 01 ( t ) = P n i =1 [ T i ≤ t ∧ U i ]. If ther e Estimating me dic al c osts 355 are no ties in the surviv al times T i , the na tural estimator of c 01 is ˆ c 01 ( T i ) = y i . Therefore our es timator o f NPV is (7) ∧ NPV = Z τ 0 e − r t ˆ c 01 ( t ) ˆ S ( t − ) Y 0 ( t ) dN 01 ( t ) = n X i =1 e − r T i y i ˆ S ( T i − ) Y 0 ( T i ) [ T i ≤ U i ∧ τ ] . If ˆ G denotes the Kaplan-Meier estimato r of the surviv al distribution o f U i , a nd using the fact that ˆ S ( t − ) ˆ G ( t − ) = n − 1 Y 0 ( t ) if there are no ties b et ween surviv al and censor ing times, then ( 7 ) can b e re written a s (8) ∧ NPV = n − 1 Z τ 0 e − r t ˆ c 01 ( t ) ˆ G ( t − ) dN 01 ( t ) = n − 1 n X i =1 e − r T i y i [ T i ≤ U i ∧ τ ] / ˆ G ( T i − ) . Using the consistency of the Ka pla n-Meier estimato r ˆ G w e see that N ˆ PV converges to E ( e − r T y [ T ≤ τ ]) provided G ( τ − ) > 0 . Ther efore in the abs e nce o f discounting N ˆ PV estimates the av erage co st restricted to τ . In this case ( 8 ) with r = 0 is the mean cost estimato r describ ed by Bang and Tsiatis [ 2 ] and Zhao a nd Tia n [ 2 1 ]. If there are ties in the sur viv al times a nd 0 < t ∗ 1 < . . . < t ∗ p ≤ τ are the distinct observed times, then ˆ c 01 ( t ∗ j ) = ¯ y ∗ j is the mean of the observed c osts at time t ∗ j and the right-hand side o f ( 8 ) is n − 1 X j : t ∗ j ≤ τ e − r t ∗ j d j ¯ y ∗ j / ˆ G ( t ∗ j − ) , where d j is the multiplicit y of t ∗ j . 3.2. Single soj ourn without c ovari ates A single so journ b egins in state 0 a nd ends with transition to sta te 1 at time T . So journ cost is incur red through time T ∗ = min( T , τ ). F rom ( 6 ) the NPV of in terest is (9) NPV = Z τ 0 e − r t S ( t − ) b 0 ( t ) dt = Z τ 0 S ( t − ) dm ( t ) where m ( t ) = R t 0 e − r u b 0 ( u ) du . Allowing for an initial cost at t = 0, integration-by- parts yields (10) NPV + m (0 ) = E ( m ( T ∗ )) = Z τ 0 m ( t )( − dS ( t )) + m ( τ ) S ( τ ) where m (0) is the exp ected initial cost. In the absence of discounting ( r = 0) and ignoring cov aria tes, Str awderman [ 18 ] consider s the nonparametric estimation of NPV based o n observ ations on (censored) surviv al times and accum ulating co sts in [0 , τ ]. F or the i -th sub ject the observed data are ( N i ( t ) , Y i ( t ) , V i ( t ) : t ≤ τ ), where V i ( t ) is the accumulated costs up to time t , N i ( t ) = [ T i ≤ t, T i ≤ U i ], and Y i ( t ) = [ T i ∧ U i ≥ t ]. Define R ( t, u ) = E ( V i ( t ) | T i ≥ u ) for t ≥ u and estimate m ( t ) = Z t 0 E ( R ( du, u ) | T ≥ u ) 356 J. C. Gar diner, L. L iu and Z. Luo by ˆ m ( t ) = n X i =1 Z t 0 { Y 0 ( u ) } − 1 Y i ( u ) dV i ( u ) . This leads to the es timator of NPV, (11) ∧ NPV = n X i =1 Z τ 0 ˆ S ( t − ) Y i ( t ) Y − 1 0 ( t ) dV i ( t ) where, as b efore ˆ S is the Ka pla n-Meier estimator of S and Y 0 ( t ) = P n i =1 Y i ( t ). T o include in ( 11 ) a c o st at t = 0 we could add the term n − 1 P n i =1 ˆ m i (0) as the estimator of m (0 ). Beca use ˆ S ( t ) − ˆ S ( t − ) = − ˆ S ( t − ) ∆ N 01 ( t ) Y 0 ( t ) the right hand s ide o f ( 10 ) would b e estimated b y n X i =1 ˆ m ( T i ) ˆ S ( T i − ) { Y 0 ( T i ) } − 1 [ T i ≤ U i ∧ τ ] + ˆ m ( τ ) ˆ S ( τ ) . Expressio n ( 11 ) is useful when the accumulating cost histor y is obser ved. 3.3. Single soj ourn without c ovari ates with r estricte d c ost history Suppo se the c o st accumulation pro cess V i ( t ) is observed a t fixed time p oints { a 0 , . . . , a G } where 0 = a 0 < a 1 < · · · < a G = τ . Let V ig = V i ( a g ) − V i ( a g − 1 ). If observ ation go es pas t a g then V ig is observed. If T i ∈ ( a g − 1 , a g ] then V i ( a g ) = V i ( T i ) and if T i ≤ a g − 1 , V ig = 0 . When censoring occur s in ( a g − 1 , a g ] the tr ue incremental cost in the interv al is not known. W e o nly obser ve ˜ V ig = V i ( U i ) − V i ( a g − 1 ). In all other c a ses we define ˜ V ig = V ig . Regarding dV i ( t ) in ( 11 ) as a dis crete mea sure w ith mass ˜ V ig at t = a g − 1 we obtain (12) ∧ NPV = G X g =1 ˆ S ( a g − 1 − ) Y − 1 0 ( a g − 1 ) n X i =1 Y i ( a g − 1 ) ˜ V ig . This estimator was in tro duced b y Lin et al. [ 15 ]. By the weak law of large n umbers and the indep endence o f U i with T i and V i ( t ) Y − 1 0 ( a g − 1 ) n X i =1 Y i ( a g − 1 ) V ig → E ( Y i ( a g − 1 ) V ig ) /E ( Y i ( a g − 1 )) = E ( V ig | T i ≥ a g − 1 ) . Because ˜ V ig differs from V ig when there is cens oring, ( 12 ) conv erges to G X g =1 S 1 ( a g − 1 − ) E ( V ig | T i ≥ a g − 1 ) − E ∗ = E ( V i ( τ )) − E ∗ where E ∗ = P G g =1 E { ( V i ( a g ) − V i ( U i ))[ U i ≤ T i ∧ a g ] | U i ≥ a g − 1 } . Hence there is down ward bias in estimating the mea n c ost E ( V i ( τ )). If censor ing do es o ccur clos e to the rig ht endpo int of the interv als this bias is likely to b e small. Estimating me dic al c osts 357 3.4. R e gr ession mo del-b ase d estimates of NPV F or the i -th sub ject let Y i = ( y i 1 , . . . , y in i ) ′ denote the costs for the so journs ending at chronologically or dered times t i = ( t i 1 , . . . , t in i ) ′ . If the las t so journ ha s not ended with transition to an absorbing state, w e will de fine t in i = τ so that y in i is the cost asso ciated with the p er io d [ t in i − 1 , τ ]. Censo r ing w ould also preclude observ ation o f some so journ co sts. Observ ation ends in one o f three ways: (1) censo ring o ccur s at U i befo re τ , (2 ) an absorbing s tate is reached b efore τ , or (3) observ atio n go es past τ . The cos t y ig asso ciated with the g -th s o journ interv al ( t ig − 1 , t ig ] is obser ved if s ig = 1 wher e s ig = [ U i ≥ t ig ∧ τ ]. Let s i denote the diagona l matrix of the { s ig , g = 1 , . . . , n i } and X i = ( x i 1 , . . . , x in i ) ′ be a n i × p ma trix of cov ariates asso c ia ted with Y i . The comp onents of x ig contain co v aria tes that are fixed ov e r time as well as cov aria tes tha t v a ry with time, but only thro ugh ( t i 1 , . . . , t ig ). In particula r x ig will contain functions o f t ig − 1 , t ig . The co nditio na l mean v ector and co v aria nce ma tr ix are denoted, resp ectively , by µ i = E ( Y i | X i ) , V i = E [( Y i − µ i )( Y i − µ i ) ′ | X i ]. W e imp ose strict exo geneity on the conditional means µ ig = E ( y ig | X i ) that requires µ ig to be a function of x ig only , that is, E ( y ig | x i 1 , . . . x in i ) = E ( y ig | x ig ) for all g = 1 , . . . , n i . Independence across s ub jects is assumed, in fa c t that { ( Y i , X i , s i ) : 1 ≤ i ≤ n } is a random sample. The total nu mber of records in the s ample is N = P n i =1 n i . Let h b e a link function such that h ( µ ig ) = x ′ ig β wher e β is a p × 1 v ector of unknown parameters . (The β her e is not the same as the regres sion parameter in the int ensity model ( 3 ) of section 2.3 .) The n i × p matrix D i of deriv ativ es ∂ µ i ∂ β ′ can b e expressed a s D i = D 0 i X i where D 0 i is the diagonal matrix with elements ( dh/dx ) − 1 ev aluated at x = µ ig . Assuming V i is p os itive definite we may wr ite V i = L i L ′ i where L i is the unique lo wer triangular matrix with p ositive diago nal elements. Make the transformations ˜ Y i = w 1 / 2 i L − 1 i Y i , ˜ µ i = w 1 / 2 i L − 1 i µ i where w i is the diagonal matrix with elements w ig = s ig /p ( t ig ∧ τ − , z i ) and p ( t, z i ) = P [ U i > t | z i ]. Here z i are fixed cov a riates that mo del the censor ing distribution. They may differ from the comp onents of X i . Given z i , a ssume U i is indep endent of ( Y i , X i , t i ). Then E ( s ig | Y i , X i , t i , z i ) = P [ U i ≥ t ig ∧ τ | z i ] and E ( w ig | Y i , X i , t i , z i ) = 1 under the assumption p ( τ − , z i ) > 0. An estimator of β is obtained by minimizing the sum of squares ˜ q ( Y i , w i , X i ) = P n i =1 ( ˜ Y i − ˜ µ i ) ′ ( ˜ Y i − ˜ µ i ) with resp ect to β which leads to the estimating equation (13) n X i =1 D ′ i ( L − 1 i ) ′ w i ( L − 1 i )( Y i − µ i ) = 0 . Because E [ D ′ i ( L − 1 i ) ′ w i ( L − 1 i )( Y i − µ i )] = E [ D ′ i ( L − 1 i ) ′ E ( w i | Y i , X i , t i , z i )( L − 1 i )( Y i − µ i )] = E [ D ′ i V − 1 i ( Y i − µ i )] = 0 , ( 13 ) pr ovides a consistent es timator ˆ β of β . The tra nsformation of Y i − µ i and D i by w 1 / 2 i L − 1 i preserves time o rder and effectiv ely uses only uncensored data in ( 13 ). In the absence o f censoring we would use the estima ting equation n X i =1 D ′ i V − 1 i ( Y i − µ i ) = 0 . 358 J. C. Gar diner, L. L iu and Z. Luo Hence ( 13 ) is the g eneralized estimating equatio ns (GEE) analo g for the selected sample { ( Y i , X i , s i ) : 1 ≤ i ≤ n } . F ollowing the standard GEE methodolo gy , n 1 / 2 ( ˆ β − β ) is asymptotica lly nor mal with zero mean and cov ariance matrix A − 1 BA − 1 where A = E  ∂ S i ( w i , Y i , X i , β ) ∂ β ′  , (14) B = E [ S i ( w i , Y i , X i , β ) S ′ i ( w i , Y i , X i , β )] and S i ( w i , Y i , X i , β ) = D ′ i ( L − 1 i ) ′ w i ( L − 1 i )( Y i − µ i ). Cons istent estimators of A and B a re o btained by r eplacing the ex p ecta tions in ( 14 ) by their sample averages and β by ˆ β . In addition, w e also need a consistent e s timator of V i = L i L ′ i and the censor ing distribution p ( t, z i ). Metho ds for their estimation are suggested in sp ecific contexts in Lin [ 13 , 1 4 ], Baser et al. [ 5 ] and Gardiner et a l. [ 8 ]. Another approa ch is to estimate a ra ndom-effects (RE) mo del for Y i (or a tra ns- formation of Y i ) given by (15) Y i = X i β + a i 1 i + u i where β is a n unknown p × 1 parameter, 1 i the n i × 1 vector with all elements equal to 1, a i an unobser ved patien t-sp ecific heterog eneity and u i is the n i × 1 vector of idiosyncratic erro r s. The c o mpo site error is v i = a i 1 i + u i . Assume Ω i = E ( v i v ′ i ) is p ositive definite a nd that the standar d RE a s sumptions (W o oldr idg e [ 20 ]) hold: (a) E ( u i | X i , a i ) = 0 , E ( a i | X i ) = 0, (b) r ank E ( X ′ i Ω − 1 i X i ) = p , (c) E ( u i u ′ i | X i , a i ) = σ 2 u I i , E ( a 2 i | X i ) = σ 2 a where σ 2 u and σ 2 a are constants and I i is the n i × n i ident ity matrix . Therefore E ( v i ) = 0 and Ω i = σ 2 u I i + σ 2 a J i where J i is the n i × n i matrix with all elements equal to 1. T o estimate β in ( 15 ) from censored observ a tions on costs we fir st tra nsform ( Y i , X i , v i ) to ( ˜ Y i , ˜ X i , ˜ v i ) wher e ˜ v i = w 1 / 2 i L − 1 i v i and ˜ Y i , ˜ X i are similar ly defined. Here L i is the unique lower triang ular matrix with p ositive diagonal elements such that Ω i = L i L ′ i . The ob jective function for estimating β is ˜ q ( Y i , w i , X i ) = { w 1 / 2 i ( L ′ i ) − 1 ( Y i − X i β ) } ′ { w 1 / 2 i L − 1 i ( Y i − X i β ) } . Spec ia lizing ( 13 ) leads to the generalized leas t-squares (GLS) weigh ted estimator ˆ β w given b y (16) ˆ β w = n X i =1 ˜ X ′ i ˜ X i ! − 1 n X i =1 ˜ X ′ i ˜ Y i ! . F rom ( 16 ) we get the c o nsistency of ˆ β w and (17) n 1 / 2 ( ˆ β w − β ) → N (0 , A − 1 BA − 1 ) where A = E ( X ′ i Ω − 1 i X i ) and B = E ( ˜ X ′ i ˜ v i ˜ v ′ i ˜ X i ). 3.5. Estimation of NPV F rom our mo del ( 15 ) for all tra nsition costs w e o btain estimates of c hj ( t | z ) for a cov a riate pr ofile z by sp ecifying the cov aria tes x 0 corres p o nding to column positions Estimating me dic al c osts 359 in X. The row v ector x ′ ij of X i in o ur mo del for y ij will con tain the fixed cov aria tes x i , dummies for transitio ns types, terms of mo deling the transition times such as t ij , t 2 ij and perha ps interactions b etw een these times and x i . Our s pe c ial x 0 will contain the des ir ed z, interactions b etw een z , t and t 2 , indicato r v a riables with v alue 1 for tra ns ition type h → j , and v a lue 0 for all other tr a nsition type s . Denoting this cov aria te pro file by x hj 0 ( t ) then c hj ( t | z )= x ′ hj 0 ( t ) β and fr om ( 16 ) we o btain the estimator (18) ˆ c hj ( t | z ) = x ′ hj 0 ( t ) ˆ β w . Although the consistency o f ˆ c hj ( t | z ) might seem immediate from ( 18 ) the final form o f the computable ˆ β w inv o lves the estimated Ω i and weights w i , the latter through the censor ing distribution G . A formal verification is not attempted her e, but see Ba s er et al. [ 5 ] for a simila r context. Now reca ll the ex pe c ted net present v alue E ( C (1) hj | X (0) = i, z ). Plugging in estimators for the ent ities on the r ig ht hand s ide in ( 4 ) leads to (19) ˆ E ( C (1) hj ( τ ) | X 0 = i, z ) = Z τ 0 e − r t ˆ c hj ( t | z ) ˆ P ih (0 , t − | z ) d ˆ A hj ( t | z ) . The estimation of E ( C (2) h | X (0) = i , z ) is entirely analogous except that one must deal with the quantit y b h ( t | z ) which is the exp ected mean rate of exp enditures at time t while so journing in state h . In practice it will not b e o bserv a ble unless discrete information is av ailable. Ins tead, we will know only the total c o st o f the so journ. F or example, consider hospital costs for patients undergoing coro nary artery b ypass surgery . Exp enditures are incurred in v arious care units such a s the intensiv e care unit, c a rdiac care unit and in recov ery . W e would know the en try a nd exit dates for each unit and the asso ciated cost of the leng th of stay in eac h unit, but not necessarily the cos t p er day . An applicatio n modeling tr eatment co st rates in cancer patients is discussed in Gardiner et al. [ 8 ] using a mo del for the lo g-transfo r med rate of cost a ccumulation y ij = y ∗ ij / ( t ij − t ij − 1 ) betw een consecutiv e tr ansition times t i 1 , t i 2 , . . . w he r e y ∗ ij the so jour n cost in [ t ij − 1 , t ij ). 3.6. Single tr ansition with c ovari ates Consider the same scenario discussed previous ly in 3.1 with all patien ts starting in state “0” and fo llow e d until they r each the terminal s tate “1” (dead). F or the i -th pa tient T i is the s urviv a l time a nd U i the censo ring time. O bserv a tio n cea ses at min ( T i , U i , τ ), that is, either at the failure time, or c ensoring time or the limit of o bs erv atio n. The only cost incurr ed is y i = y i ( T i ) at time T i which is observed if s i = 1 wher e s i = [ U i ∧ τ ≥ T i ]. Let x i denote a p - vector of fixed cov aria tes o f int erest and z i denote fixed cov a r iates used for modeling the cens o ring distribution. Assuming indep endent censor ing, that is, given z i , U i is indep endent of ( y i , x i , T i ) , we get P [ s i = 1 | z i , y i , x i , T i ] = P [ U i ≥ T i , T i ≤ τ | z i , T i ] = G ( T i − | z i )[ T i ≤ τ ] . Defining w i = s i /G ( T i −| z i ) we see tha t ( 13 ) reduces to minimizing with resp ect to β the ob jective function n − 1 P n i =1 q ( y i , w i , x i ) where q ( y i , w i , x i ) = σ − 2 u w i ( y i − x ′ i β ) 2 . This yields the estima to r ˆ β w in ( 16 ) which in this case is ˆ β w = n X i =1 w i x i x ′ i ! − 1 n X i =1 w i x i y i . 360 J. C. Gar diner, L. L iu and Z. Luo This is the same e stimator describe d b y Lin [ 13 ] e xcept for a sligh t difference in the weigh ts. Because Lin [ 13 ] uses a model in which costs are incurre d through time T i ∧ τ his censoring indicator is s ∗ i = [ U i ≥ T i ∧ τ ] and w eight w ∗ i = s ∗ i /G ( T i ∧ τ − | z i ). In our case the cost is realized at time T i if and only if s i = 1. Let z 0 denote a fixed cov ar iate at which NP V( z 0 ) is to b e estimated. Since “ 0” is the initial s tate a nd the o nly transitio n is 0 → 1 , ˆ P 00 (0 , t − | z 0 ) = ˆ S ( t − | z 0 ) and ˆ S ( t | z 0 ) = exp( − ˆ A 01 ( t | z 0 )). Her e ˆ S is the es tima to r of the surviv al distribution S of T , the time of transition. F rom ( 18 ) and ( 19 ) our estimator of NPV( z 0 ) is ∧ NPV( z 0 ) = ˆ β ′ w Z τ 0 e − r t x 0 ( t ) {− d ˆ S ( t | z 0 ) } , where x 0 ( t ) is the cov ariate vector der ived from z 0 and ter ms us e d to mo del time (such as t, t 2 ) in the co st equation y i = x ′ i β + u i . Here a sing le co st y i = y i ( T i ) is incurred at T i , if obser ved b y time τ . Then E [ y i ( t ) | z 0 , T i = t ] = x ′ 0 ( t ) β and NPV( z 0 ) = Z τ 0 e − r t E [ y i ( t ) | z 0 , T i = t ] {− dS ( t | z 0 ) } simplifies to E ( e − r T i y i ( T i )[ T i ≤ τ ] | z 0 ). Since ˆ β w → β w in probability , and uniformly o n [0, τ ], ˆ S ( ·| z 0 ) → S ( ·| z 0 ) in probability , if S ( τ | z 0 ) > 0 , w e obtain the c o nsistency of ∧ NPV( z 0 ) provided R τ 0 e − r t x 0 ( t ) dS ( t | z 0 ) is finite. Also in estimating β w we require P [ T i ≤ τ | z 0 ] = 1 − S ( t | z 0 ) > 0 , bec ause other wise the c o st equatio n will b e v acuous since with pro bability 1 no transition takes plac e in [0, τ ]. 3.7. Single soj ourn with c ovariates Suppo se the interv al [0 , τ ] is par titioned by the fixed p o ints a j , j = 0 , . . . , K with 0 = a 0 < a 1 < . . . < a K = τ . If the exp ected rate of cost accumulation is constant in the interv als ( a j − 1 , a j ) with v a lue s b j we hav e NPV( z 0 ) = K X j =1 b j Z a j a j − 1 e − r t S ( t | z 0 ) dt. The integral is the increment ov er ( a j − 1 , a j ) in discounted life ex pec ta ncy , LE ( z 0 , t ) = Z t 0 e − r u S ( t | z 0 ) du. F ollowing Baser et al. [ 5 ] we could use a RE mo del for cost y i = ( y i 1 , . . . , y iK ) ′ incurred by the i -th pa tient . Here y ij is the co st incur red in in terv al ( a j − 1 , a j ) whic h is observed pro vided s ij = 1 where s ij = [ T i ∧ U i ≥ a j ] + [ a j − 1 < T i < U i ∧ a j ]. This reflects the tw o ca ses: (1) the pa tient neither died nor was censor ed b efore a j , or (2) death was obser ved in ( a j − 1 , a j ). Under the assumed indep endence of censoring P [ s ij = 1 | z i , T i ] = G ( T ∗ ij − | z i )[ T i ≥ a j − 1 ] wher e T ∗ ij = min( T i , a j ). The r egress ion mo del for y ij will include int erv a l-sp ecific elaps ed time T ∗ ij − a j − 1 which will yield an estimator of b j that may dep end on z 0 . Estimating me dic al c osts 361 4. Discussion and sum mary The estimation o f medical costs has received co nsiderable attent ion b ecause of its impo rtance in ass essing c o st-effectiveness of medical interven tions and tr e a tment s. F acing constra ined healthcar e budgets gov ernment planners and po licy makers are forced to cons ider the costs o f co mpe ting int erven tions in addition to claims of their clinical efficacy . The difference in the expected co st of tw o co mpe ting in terven tio ns is the n umerator o f the cos t-effectiveness ratio, the denominator being the inc r emental health b enefit as mea sured by life ex p ecta ncy or b y quality-adjusted life y ears (Chen and Sen [ 6 , 7 ] and Gardiner et a l. [ 9 ]). The cost-effectiveness ratio can b e used to compare comp eting interven tions with resp ect to b o th their health b enefits as well as their cos t. Obtaining reliable and v alid estimates of costs is imp e r ative. In this ar ticle we adopted a longitudinal framew ork in w hich patient co sts are manifested dynami- cally ov e r time. An under ly ing finite-state s to chastic pro ces s describ e s the e volving patient history with costs incurr e d at transitio n times b etw een states and during so journ in states. In this framework we sho wed how net pres ent v a lues ar e defined, following the basic notions o f a ctuarial v alues used extensively in the insurance a nd finance literature (Norb erg [ 16 ]). F or example, in the cla ssic disa bilit y mo del there are “ a ble” p erio ds and “disabled” p e r io ds. The individua l holding a disability in- surance p olicy would receive a fixe d payment str eam ov er the p er io d of his or her disability . In able p erio ds the individual would pay the fixed premium in accordance with the p olicy . There ar e three p olicy states –“a ble”, “disa bled”, and “dead” . A fundamen tal difference in our context is tha t cos ts a re not fixed but r andom. One may re g ard the total cost over a s pe cified p erio d as the sum o f all tra nsition costs and so jo urn cos ts . Several metho ds hav e b een propo sed to estimate medica l co st from follow up data. The prima ry fo cus has be e n on a single cost measure tha t might b e incom- pletely ascertained due to time censoring (Ba ng and Tsiatis [ 2 ], B a ser et al. [ 4 ], L in et al. [ 15 ], Strawderman [ 1 8 ] and O’Haga n a nd Stevens [ 17 ]). Regress io n a nalyses allow for assessing the influence of explana tory v ar iables on some measure of the cost distribution, s uch a s the mean or median (Ba ng a nd Tsiatis [ 3 ], Baser et al. [ 5 ], Lin [ 1 3 , 14 ] and Gardiner et al. [ 8 ]). Apart from address ing the incomplete- ness of co st data, the a bilit y to o bs erve costs ov er finer time p erio ds can serve to strengthen ensuing ana lyses. F or instance , cons ider the cost of a treatment which is a ssumed to last at most for one year and costs a re monitor ed monthly . If either the endp oint is reac hed b efor e the end o f the year, or observ atio n lasts one y ear, the total cos t is obser ved. If there is censoring of the endpoint b efor e the end of the year, we could use the monthly costs, except for the la st mo nth of obser v ation, to improv e our estimate of the av erage cos t of treatment. The metho ds disc us sed here for analyses of medical costs may be adapted to estimate other summar y mea sures used in cos t- effectiveness analyses (Gardiner et al. [ 9 ]). F or example, quality-adjusted surviv a l is defined by using a q uality weigh t q ( h, t ) whic h r e pr esents the utilit y , relative to the state of p erfect health, of ea ch unit of time spent in state h = X ( t ) at time t . Perfect health has a quality w eight 1, while death or states judged equiv a lent to death g et a qualit y weigh t of 0. The total quality adjuste d time in [0 , τ ] is P h ∈ E R τ 0 e − r t q ( h, t ) Y h ( t ) dt . Hence conditional on X (0 ) = i we define the exp e cte d quality adjuste d life ye ars , QALY i ( z ) = P h ∈ E R τ 0 e − r t q ( h, t ) P ih (0 , t − | z ) dt . The unconditional v ersion is given by QALY ( z ) = P i ∈ E π i (0 | z ) QALY i ( z ). This expression is similar to the second term of the NPV in ( 6 ). 362 J. C. Gar diner, L. L iu and Z. Luo The transition mo del a dopted in this a r ticle extends the s impler t wo-state sur- viv al mo del with a s ing le transition and so journ. The underlying analys is of surviv a l times is no w r eplaced b y the ana lysis of multiple even t times which is facilitated by using a no n-homogeneo us Markov mo del to govern the movemen t b etw een states , and a multiplicative in tensity mo del to inco rp orate cov aria te effects. F o r the ana l- ysis of longitudinal c o st da ta, techniques such as inv erse probability weigh ting to account for censo ring can be applied (Willan et al. [ 19 ]) but a more careful con- sideration is requir ed to combine the tw o pa rts of the mo del, the transition mo del for the even t times and a r egressio n mo del for costs. Methods for joint mo deling of lo ngitudinal obser v ations and ev ent times could b e adapted for this purp ose (Henderson et al. [ 10 ] a nd Hogan and Lair d [ 1 1 , 12 ]). Ac knowledgmen ts. W e thank the referees for their car eful reading of the man uscript. Their comments have improv ed the presentation of the pap er. References [1] Andersen, P . K., Borgan, O. , Gil l, R. D. a nd Keiding, N . (19 93). Sta- tistic al Mo dels Base d on Counting Pr o c esses . Springer, New Y ork. MR11988 84 [2] Bang, H. and Tsia tis, A. A. (2000). Estimating medica l costs with censor ed data. Biometrika 87 3 29–34 3. MR1782 482 [3] Bang, H. and Tsia tis, A. A. (2 0 02). Median r egressio n with censored cos t data. Biometrics 58 643–6 49. MR19261 17 [4] Baser, O. , Gardiner, J. C., Bradl ey, C. J. and Given, C. W. (2004 ). Estimation from censored medical cost data. Biometric al Journal 46 35 1–363 . MR20620 00 [5] Baser, O. , Gardiner, J. C., Bradley, C. J., Yuce, H. and Given, C. (2006). Longitudinal ana lysis of censo r ed medical cost data. He alth Ec onomics 15 513– 5 25. [6] Chen, P. L. an d Sen, P. K. (200 1). Quality-adjusted surv iv al estima tio n with p erio dic o bserv a tio ns. Biometrics 57 868 –874. MR18634 51 [7] Chen, P. L. an d Sen, P. K. (200 4). Quality-adjusted surv iv al estima tio n with p erio dic obser v ations: A multistate surviv a l analysis approach. Comm. Statist. The ory Metho ds 33 1327 –133 9. MR20695 71 [8] Gardiner, J. C., Luo, Z., Bradley, C. J. , Sirbu, C. M. and Given, C. W. (2006a ). A dynamic mo del for estimating changes in health status a nd costs. Statistics in Me dicine 25 36 48–3 667. MR22524 17 [9] Gardiner, J. C., Luo, Z., Liu, L. and Bradley, C. J. (2006b). A s to c has- tic framework for estimation of summary mea sures in cost-e ffectiveness analy- ses. Exp ert R eview of Pharmac o e c onomics and Outc omes R ese ar ch 6 347–35 8. [10] Henderson, R ., Diggle, P. and Dobson, A. (2000). Joint mo delling of longitudinal measurements a nd event time data. Biostatistics 1 46 5–480 . [11] Hogan, J. W. and Laird, N. M. ( 1997 a). Mixture models for the join t distribution o f re p ea ted meas ur es and even t times. Statistics in Me dicine 16 239–2 57. [12] Hogan, J. W. and Laird, N. M. (19 97b). Mo del-based appr oaches to analysing inco mplete longitudinal a nd failure time data. Statist ics in Me dici ne 16 259– 2 72. [13] Lin, D. Y. (2000). Linear r egressio n of censored medical costs. Biostatistics 1 35–47 . Estimating me dic al c osts 363 [14] Lin, D. Y. (2 003). Regression analysis of incomplete medical co st data. Statis- tics in Me dicine 22 11 81–1 2 00. [15] Lin, D. Y. , Feuer, E. J., Etz ioni, R. and W ax, Y. (1997). Estimating medical costs from inco mplete fo llow-up data. Biometrics 53 419– 434. [16] Norberg, R. (1995). Differential-equations fo r moments of present v alues in life-insurance. Insur anc e Mathematics and Ec onomics 17 171 –180 . MR13636 47 [17] O’Hagan, A. and Stevens, J. W. (2004 ). On e stimators o f medical costs with censored data. J ournal of He alth Ec onomics 23 6 15–6 25. [18] Stra wderman, R. L. (2 000). Estimating the mean of an incr e asing sto chas- tic pro cess at a censored stopping time. Journal of the A meric an Statistic al Asso ciation 95 1192 –120 8. MR18042 43 [19] Willan, A . R., Lin, D. Y., Cook, R. J. a nd Chen, E. B. (2 002). Using inv er se-weigh ting in cost- effectiveness ana lysis with censo red data. Statistic al Metho ds in Me dic al R ese ar ch 11 539 –551 . [20] Wooldridge, J. M. (200 2). Ec onometric Analysis of Cr oss Se ction and Panel Data . MIT Pr ess, Ca mbridge, MA. [21] Zhao, H. W. and Tian, L . L. (2001). On e s timating medical cost and incr e- men tal co s t-effectiveness ra tios with censor ed data. Biometrics 57 1002– 1 008. MR19504 17