Compartmental models and competing risk

Compartmental Models and Competing Risk G. M. TALLIS Department of Statistics, The University of Adelaide, Adelaide, South Australia Received 30 September 1992; revised 10 June 1993

ABSTRACT General compartmental models are derived using competing risk arguments. When the risk variables are exponential, the results specialize to the standard stationary Markov compartmental model. Iterative methods of solving the fundamental integral equation are given, and the uniqueness of the solution is incidentally established. The analysis is extended to include fixed inputs, orderly and nonorderly stream infusions, and time dependency. The study is motivated by a biological system that evolves through various stages over time.

1.

INTRODUCTION

A compartmental model is a closed, evolutionary system that consists of a finite number of states Y1,pZ,. . . ,Pn. The system starts in one of the states, 3, say, at time t = 0, and the general properties of the system determine the probability Pjj(t) that it is in state j at time t. From this broad description it is clear that such models have a very wide application. However, this work emphasizes a particular biological system as illustration, although the results that are generated have quite general relevance. Specifically, the biological system involves an organism that goes through a number of stages during its development. The stages equate to states in a compartmental model, and q represents stage i. In this setting, z_ 1 comes chronologically before q for i = 1,2,. . . , n - 1, while Yn is the “stage” death. Thus the system is closed, and at each stage the organism can mature to the next stage or go directly to pn, that is, die, before it matures. In the case of pn_ i, of course, the only possible transition is to pn. To specify the process, the transition matrix P(t) = [Pij(t>] is needed, where Pii

MATHEMATICAL

BIOSCIENCES

=Pr(q,tIe,O).

121:111-122 (1994)

OElsevier Science Inc., 1994 655 Avenue of the Americas, New York, NY 10010

111 0025-5564/94/$7.00

112

G. M. TALLIS

Then, if Xi(t) is the number of organisms in q at t, and xi(O) = X,(O), the probability generating function for X,(t), . . . , X,(t) is

where gi(s, t) = C& , Pij(t)sj, under the assumption that each organism acts independently of all others according to dictates of the system. There has been a large amount of research done on the theory and application of compartmental models. For an introduction, see the book by Chiang 141; a more comprehensive review, with a substantial bibliography and examples, is given by Matis and Wherly [8]. The latter reference also considers problems of parameter estimation. Another introductory article on the subject emphasizing the semi-Markov approach to compartmental modeling is that of Mode and Pickens [ll], who give a useful algorithm for computing transition probabilities for discrete time processes. They mention that transitions out of any state can be viewed as competing risks, and note that Berman [3] observed this in 1963. There are many other, more technical papers that tend not to be so accessible or in as convenient a computational form for users. Agrafiotis [l] developed a semi-Markov approach for a restricted class of continuous-time compartmental models with infusion, and Weiner and Purdue [13] also used semi-Markov systems in their studies. Mehata and Selvam [lo], who provide a useful bibliography, studied a class of semi-Markov systems with fixed and Poisson inputs, emphasizing recurrence times and the numbers in each compartment. For many genuine and detailed applications, see [9, 7, 121. The purpose of this paper is to derive general competing risk techniques to specify the waiting time distributions for transitions to the various compartments and to use a renewal type of matrix integral equation to specify P(t). This equation can be solved numerically by an iterative method that converges exponentially quickly to the unique solution. Some problems of nonstationarity are then considered, following which generating function techniques are used to extend the analysis to orderly and nonorderly stream infusions. A conventional biological model motivates and illustrates the work. II.

NOTATION

AND

STANDARD

RESULTS

It is convenient to use the notation 3 + q to represent passage from 3 to 3. Then the specification of the standard, stationary,

COMPARTMENTAL

Markov L$ + q

MODELS

AND

compartmental model by t, is as follows.

COMPETING

leading

113

RISK

to the Pi,(t),

the probability

of

Pr(q,t+AIq,t)=&,A+o(A), Pr(q,t+AIq,t)=l-

c

&,A+o(A).

j#i

These

two equations 4ij

=

can be put into matrix form by writing

C

qii=-

Pji;

Q=[qijI;

Pij;

j#i

p’(t) = (PlwYd%w)~

Pi(l) =Pr(q,t); Then

where I is the IZx n identity matrix. Subtracting p(t) from both sides, dividing by A, and letting A + 0 gives (p’(t) = dp(t)/dt)

if(t) = Qdt)

(2)

with solution = eQ%(0) 9

P(l) where

e@ is the matrix

exponential

eQt=

defined

c

(3) by

(Qt)",

kc0

k!

It can be shown that

P(t) = eQ”, where Q’ is the transpose of Q. Thus, 9(~,t) can be found from the rows of eQ” in this simple case. If eij is the probability of x -+ 3, given that there is a transition, then, as shown later, eij

and the waiting time in 2 parameter - qii, 8( - qii).

=

Sji

/( -

has a negative

4ii)

exponential

distribution

with

114

G. M. TALLIS

These results are mathematically neat and convenient, but the underlying assumptions are restrictive. It is therefore important to see in what ways the structure can be modified to make it more useful. Nevertheless, as a first approximation the above model, or simple modifications of it, may be entirely appropriate. The ensuing developments will be in terms of risk variables and hazard functions. Let X be a nonnegative random variable with distribution function (df) and probability density function (pdf) F(x) and f(x), respectively. Now consider Pr(XE[x,x+dY)IX>x)=f(x)&/[l-F(x)]=h(x)&. A(x) is called the hazard function (hf) for X, and if X is the time to some event, h(x) is the rate at which it occurs, given that it has not occurred prior to x. The relationship

fW/~l-

F(x)1 = h(x)

can be integrated to give qx)

= l-

e-*(x);

f(x)

= h(x)d’“‘;

A= A’.

The function A(x) is called the cumulative hazard function (chf). III.

GENERAL

SETUP AND RESULTS

Let

{X12,...,Xlnj ,... ,{&,...,Xnn-J be 12 series of n - 1 independent, nonnegative random variables (rv’s). For Xij, the df and pdf are Fij(x) and hi(x), with hf and chf hii and Rij(x), respectively. Thus, E;.j( x) = l-

c-*&);

hj( x) = hij( x)epAJx).

Now put Xi = min {Xij} j#i

and note that

xi > x

iff all

Xii > x,

COMPARTMENTAL so

MODELS

AND

COMPETING

115

RISK

that Pr(Xi>x)=l-&(x)=exp

j#i ‘I]

- CR..(x) [

=e-*i(‘),

say. The set {Xi, ,..., Xii_i,Xii+i ,..., Xi,} is associated with the system starting in 3, and the Xii are the potential waiting times for the change S$ --f 3. Now, (i) the waiting time in e

before a move has df

Fi(X)

(ii) the event 2 -+ 3

e-*1@);

= l-

occurs in the interval [x,x + cix) if and only if

XijE[x,x+&j; By independence,

j’# j.

Xij! > x,

this event has probability exp [

j’#j 1 C Rij.(x)

Aij( x)e-*JX)du.

(iii) The probability that 2 + 3 tljj =

is

lij(

x)e-Al(x) dx.

/0

Moreover, Aij( x)e-

*qe,, =gij(

x)

is a pdf. It is the waiting distribution given that there is a move 3 -+ 3. Therefore, Yij(X)&

= Bijgij(x)dx

is the joint probability of 3 + q

=

h,j(x)e-*w?!x

and Xi E [x, x + dx).

It is now possible to write an integral equation for the matrix of transition probabilities p(t) based on the assumption that the system is stationary and has an embedded Markov system. Put

YCx)= [ Yijtx>]

y

y,i( x) = 0;

G. M. TALLIS

116

then

pij(r)

5

=

k=l

/'Y,k(Y)'kj(' 0

i#

-Y)dYy

j.

If F(t)=diag[l-F,(t),...,l-F,(t)], then, in matrix notation,

P(l) =j:Y(Y)P(I

-Y)dY

+flt>-

(4)

SPECL4 L CASE

As an illustration, consider the case where the Xii are distributed as 2% pjj), where the pij are as defined in Section II. Thus, (Xil,.. . , XJ are n - 1 independent, exponentially distributed random variables. Now make the identification yi,(x)

= qjieq~~*;

eij = qji_/meq~~x & 0

= *. II

Note also that if A = diag(q,, , . . . , q,,), then y(x)

= e*“[ Q’ - A]

in this case. Does P(t) = eo” satisfy (4)? That is, we wish to know whether eQ’f =

eQ’(t-Y)dy

= eA’

+ eAf

‘,-ncr-.“, [Q' - A]eQ'('-Y)dy /0

+ e”‘.

But x

eQ’q

= _ ep”“~

,Q’x

+ e-hx

Q’eQ”,

COMPARTMENTAL so

MODELS

AND

COMPETING

RISK

117

the integral is e*‘[ e - A.r eP’“]:, = ,Q'[ _ eht,

and the equality is verified. This demonstration establishes the waiting time properties of the elementary compartmental model. It also emphasizes that the process can be viewed as a competing risk model with exponentially distributed risk variables. IV.

THE SOLUTION

OF EQUATION

(4)

The above special case is a system that can be specified by risk variables that have constant hfs. Where the &;(x) are genuine functions of x, the solution of (4) must proceed iteratively. Put P@‘(t) =qt>, and define P’“‘(t)=@y)P’““(t-y)dy+F(t). Also, let

where for any it x

n matrix

A, llAll = C

laijl

i,j

(see Belman [2] for this and other results on matrix analysis) will be used as a matrix norm. Thus,

Gm(ll Pck)( y)- Pckl)(y)lldy, with

118

Hence,

G. M. TALLIS

if 1)Pck)( t) - Pckml)( t) 11G nmktk/k!,

which is certainly

true for k = 1, then

lip(k+l)( t)

- Pck)( t) II< nmkt’tk+‘/(

For all i, j, therefore, {P,:k’(t)} is a Cauchy Pi,(t), say. If P(t) = [Pl,(t)l, then

k + 1) !

sequence

and has a limit

lim Pck)( t) = P(t), k+m

and convergence

is exponentially

fast. Notice that P(t) satisfies

(4) since

lim Pck’(t)=)h+nm ~‘~(y)Pckpl)(t-y)dy+F(t),

k-rm

so that

p(t)

=p(Y)P(t - Y)dY +qt>.

To establish the uniqueness also satisfies it; then

of P(t) as a solution

to (41, suppose

G(t)

Now put m, = maxO d y 4 t II P(y) - G(y) II to show that IIP(t)-G(t)ll~m,mt; then an induction

argument

shows that

for all k, implying that P(t) = G(t). These results provide an iterative procedure for numerically calculating P(t) from an arbitrary set of hfs, hij(x). Such computational work can be readily programmed for modern high-speed computers (see Section V).

COMPARTMENTAL

V.

EXAMPLE

MODELS

AND

AND

COMPETING

RISK

119

EXTENSIONS

In the Introduction, the general compartmental system was specialized to model an organism that matured through it - 1 stages q, the nth stage, Yn, being death. For each 3, i = 1,2,. . . , n - 2, there are two hfs, Aii+ 1(y) and h,,(y). Of course, for Yn,- I there is the single hf h n_ r,,(y), and for Yn, none at all. These functions can be estimated from properly designed trials. A large number of organisms that survive to y in z must be observed for a short period, A, and the numbers dying or passing on to q+ 1 noted. These numbers divided by A and then expressed as a ratio of the total number at y estimate Ai, and hii+ r(y). Once plots of the hii are available, suitable functional forms can be fitted to them to assist with the numerics of solving (4). Care must be taken to ensure that the hii + 00 as y + 00 so that bona fide dfs emerge. This particular example emphasizes the convenience of the hf specification of the process. The maturation of organisms through various stages can be conveniently viewed in terms of competing risks. The individual makes it to the next stage provided the risk variable for death is not less than the waiting time variable for that stage. Additionally, the numerical specification of the fij(y) by means of the h,j(y) is especially appropriate as it allows an actuarial type of construction of the various distributions. To test the numerical properties of the iterative procedure of Section IV, a model of the above type involving four stages (n = 4) was used. Exponential holding times were specified so that P(t) could be determined from exp{Q’t} = Cz(Q’tjk/k! for Q suitably chosen. Then, with the corresponding y(y), the methods of Section IV were implemented. The matrices entering into the iteration were specified at points of a fine partition of [O,t], and the integration process was approximated by elementary finite numerical methods. Mainly, double precision was adequate, although quadruple was used where necessary, and convergence to P(t) occurred after five iterations for various values of t and the parameters. Accuracy and speed were remarkably good, the work being carried out on a SUNSPARC-server 2. The short program, which was written in Fortran 77, was uncomplicated, and there is reason to believe that bigger and more involved models can be adequately dealt with similarly. More sophisticated numerical approaches should be even more successful. NONSTATIONARITY

The assumption of stationarity can be lifted. In Section II, the fiij may be functions of time, Pij(t> sa.y, and hence Q will be also, Q(t). The

120

G. M. TALLIS

set of differential

equations

becomes p’(t)

which has the unique

solution At>

where

H(t)

is the unique H’(t)

H(t) can be calculated

= Q(t)At),

= H(t)p(O),

matrix

satisfying

H(0) = I.

=Q(t)ff(t), recursively

from

which converges exponentially quickly, commencing with H(‘)(t) = I. The proof of these statements is along the lines of Section IV. The general treatment of Section III can also be modified to handle nonstationarity in a cruder but effective way. Partition [O,t] into r subintervals [O,t,),[t,, t2), . . .,[t,_ Ir t>, so that on each interval it can be assumed that the process is essentially stationary. Each can then be treated in the same way as described in Sections III and IV, Equation (4) for the Ith interval, I[, assuming the form

where

Pl(tr - t,_ i) is the transition matrix pertinent to 1,, and yI(y), III but with the hazard functions, A’,f’(y), matrix for [O,tl is then the product

q(y) are defined as in Section pertinent to I!. The transition

P(t) =

fJqt,-t,_1),

l=l

t,=t.

The estimation of A’,f’(y) requires data observed specifically over I,, and the fact that individuals in q have been there for varying periods is ignored. This leads to estimates of “equilibrium” hfs and removes the necessity of making some mathematical adjustments of the type used in renewal theory (see, e.g., Cox [51). In summary, then, compartmental models of a very general character can be constructed by the methods outlined above. These should provide an excellent description of the underlying process as the hazard functions that specify it can be tailored to each practical application.

COMPARTMENTAL MIGRATION

MODELS

AND COMPETING

RISK

121

PROCESSES

As a final comment we consider modifications that are required to alter the closed compartmental system to allow for random inputs. Such systems are known as migration processes. Consider initially 3 at t = 0 and the associated pgf (Section I) 9$(s, t). Suppose new individuals enter 3 according to a “nonorderly stream” process (see Khintchine [6]) that is characterized by the infinitesimal generating function (l+

hi(tL)At[hi(9i(

s,t-tt,))-l]+o(At))

where [O,t) is partitioned into small subintervals of length At and t, belongs to the fth one. The pgf hi(s) = Cy= 1hiid specifies the probability of j new individuals starting in z given that at least one does, and the assumption of independent increments, which is required of streams, ensures that the ensuing pgf for the numbers in the various states as a consequence of “external migration” is

Note that results for “orderly streams” are obtained by putting hi(s) = s. Since the same argument holds for all i, the total pgf is

If the system started with xi(O) in q, i = 1,2,. ..,rz, then the pgf resulting from all migrations, and from the initial numbers, is

This emphasizes how very general processes can be constructed a basic compartmental model.

around

Z thank the two referees for providing references and suggestions that improved the original version of this paper. Computations involving equation (4) were carried out by Mr. P. Leppard, and his work is gratefully acknowledged. REFERENCES 1

G. K. Agrafiotis, On the stochastic theory of compartments: A semi-Markov approach for the analysis of the k-compartmental system, Bull. Math. Biol. 44309-817 (1982).

122 2 3

G. M. TALLIS

R. Belman, Introduction to Matrix Anatjks, McGraw-Hill, Sydney, (1970). S. M. Berman, A note on extreme values, competing risks, and semi-Markov processes, Ann. Math. Stat. 34:1104-1106 (1963). 4 C. L. Chiang, An Introduction to Stochastic Processes and Their Applications, Krieger, Huntington, N.Y., 1980. 5 D. R. Cox, Renewal Theory, Methuen, London, 1962. 6 A. Y. Khintchine, Mathematical Methods in the Theory of Queuing, Charles Griffin, London, 1969. A. H. Marcus, Semi-Markov compartmental models in ecology and environmental health, in Compartmental Analysis and Ecosystem Models, J. H. Matis, B. Patten, and G. C. White, Eds., International Co-operative, 1979. 8 J. H. Matis and T. E. Wherly, Stochastic models of compartmental systems, Biomehics 35:199-220 (1979). 9 J. H. Matis and T. E. Wherly, Generalized stochastic compartmental models with Erlang transit times, J. Pharmocokin. Biopharm. 18589-607 (1980). 10 K. M. Mehata and D. D. Selvam, A stochastic model for the n-compartment irreversible system, Bull. Math. Biol. 43:549-561 (1981). 11 C. J. Mode and G. T. Pickens, Computation methods for renewal theory and semi-Markov processes with illustrative examples, Am. Stat. 42:143-152 (1988). 12 H. D. Tolley, D. Burdick, K. G. Manton, and E. Stallard, A compartmental model approach to the estimation of tumor incidence and growth: Investigation of a model of cancer latency, Biometrics 34:377-389 (1978). 13 D. Weiner and P. Purdue, A semi-Markov approach to stochastic compartmental models, Commun. Stat. Theor. Math. 6:1231-1243 (1977).