Univerzita Karlova v Praze Matematicko-fyzik´aln´ı fakulta
´ RSK ˇ ´ PRACE ´ BAKALA A
ˇ Ivana Zohov´ a Aplikace Markovov´ ych proces˚ u pˇ ri modelov´ an´ı pr˚ ubˇ ehu choroby HIV Katedra pravdˇepodobnosti a matematick´e statistiky
Vedouc´ı bakal´aˇrsk´e pr´ace: Mgr. Michal Kulich, Ph.D. Studijn´ı program: matematika, obecn´a matematika
2006
Na tomto m´ıstˇe bych chtˇela podˇekovat vedouc´ımu bakal´aˇrsk´e pr´ace Mgr. Michalu Kulichovi, Ph.D. za jeho pomoc pˇri vypracov´av´an´ı t´eto pr´ace.
Prohlaˇsuji, ˇze jsem svou bakal´aˇrskou pr´aci napsala samostatnˇe a v´ yhradnˇe s pouˇzit´ım citovan´ ych pramen˚ u. Souhlas´ım se zap˚ ujˇcov´an´ım pr´ace a jej´ım zveˇrejˇ nov´an´ım. ˇ Ivana Zohov´ a
V Praze dne 5. 8. 2006
2
Obsah Seznam znaˇ cen´ı
5
´ Uvod
6
1 Choroba zp˚ usoben´ a virem HIV
7
2 Modelov´ an´ı inkubaˇ cn´ıch a infekˇ cn´ıch dob 2.1 Stochastick´e modely . . . . . . . . . . . . . . . . . . . . . . 2.2 Vlastnosti model˚ u. . . . . . . . . . . . . . . . . . . . . . . .
8 8 10
3 Markovovy procesy a jejich aplikace 3.1 Teorie Markovov´ ych proces˚ u . . . . 3.2 Modelov´an´ı pr˚ ubˇehu choroby HIV . 3.3 Vhodnost Markovov´ ych model˚ u . .
14 14 18 27
Literatura
pˇ ri modelov´ an´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
3
N´ azev pr´ ace: Aplikace Markovov´ ych proces˚ u pˇri modelov´an´ı pr˚ ubˇehu choroby HIV ˇ Autor: Ivana Zohov´ a Katedra: Katedra pravdˇepodobnosti a matematick´e statistiky Vedouc´ı bakal´ aˇ rsk´ e pr´ ace: Mgr. Michal Kulich, Ph.D. e-mail vedouc´ıho:
[email protected] Abstrakt: V pˇredloˇzen´e pr´aci se zab´ yv´ame stochastick´ ym modelov´an´ım pr˚ ubˇehu choroby HIV, pˇredevˇs´ım pak inkubaˇcn´ı doby AIDS. Hlavn´ı pozornost je vˇenov´ana modelov´an´ı pomoc´ı homogenn´ıch Markovov´ ych proces˚ u se spojit´ ym ˇcasem a koneˇcnou mnoˇzinou stav˚ u. V pr´aci nejprve studujeme nˇekter´a spojit´a rozdˇelen´ı, a to pˇredevˇs´ım z hlediska jejich vhodnosti k modelov´an´ı inkubaˇcn´ı doby AIDS. Dalˇs´ı ˇc´ast je vˇenov´ana konstrukci Markovov´ ych model˚ u. Pr˚ ubˇeh choroby HIV je nejprve modelov´an pomoc´ı dvou konkr´etn´ıch model˚ u a pot´e pomoc´ı zobecnˇen´eho modelu. V modelech jsou odvozena rozdˇelen´ı inkubaˇcn´ı doby AIDS. Kl´ıˇ cov´ a slova: choroba HIV, inkubaˇcn´ı doba AIDS, homogenn´ı Markov˚ uv proces, Markov˚ uv model
Title: Applications of Markov processes to modeling HIV disease progression ˇ Author: Ivana Zohov´ a Department: Department of Probability and Mathematical Statistics Supervisor: Mgr. Michal Kulich, Ph.D. Supervisor’s e-mail address:
[email protected] Abstract: In the present work we study stochastic modeling of HIV disease progression, and especially of the AIDS incubation period. The main attention is paid to models based on continuous time homogeneous Markov processes with a finite number of states. Firstly, in this work, some continuous distributions are considered and their suitability to model the AIDS incubation period is investigated. In the next section, we attend to the construction of Markov models. HIV disease progression is at first patterned through two concrete models, then a general model is cosidered. In models the AIDS incubation distributions are derived. Keywords: HIV disease, AIDS incubation period, homogeneous Markov process, Markov model
4
Seznam znaˇ cen´ı N0 R p pT ei B BT I Ik δij f ∗g
mnoˇzina cel´ ych nez´aporn´ ych ˇc´ısel mnoˇzina re´aln´ ych ˇc´ısel sloupcov´ y vektor vektor transponovan´ y k vektoru p jednotkov´ y sloupcov´ y vektor s ˇc´ıslem 1 na i-t´em m´ıstˇe matice matice transponovan´a k matici B jednotkov´a matice jednotkov´a matice ˇr´adu k Kroneckerovo delta konvoluce funkc´ı f a g
5
´ Uvod Onemocnˇen´ı AIDS (syndrom z´ıskan´eho imunodeficitu) je posledn´ı f´az´ı choroby zp˚ usoben´e virem HIV (d´ale jen choroby HIV). Prvn´ı zdokumentovan´e pˇr´ıpady AIDS se objevily na poˇc´atku 80. let 20. stolet´ı mezi homosexu´aly ve velk´ ych americk´ ych mˇestech jako New York, Los Angeles a San Francisco. Choroba HIV a zejm´ena jej´ı posledn´ı f´aze AIDS se brzy staly celosvˇetovou hrozbou takov´ ych rozmˇer˚ u, ˇze hovoˇr´ıme o pandemii HIV/AIDS. Podle u ´daj˚ u Svˇetov´e zdravotnick´e organizace (WHO) bylo doposud virem HIV infikov´ano kolem 50 milion˚ u lid´ı, z toho jiˇz 16 milion˚ u zemˇrelo na AIDS. V Evropˇe ˇzije v souˇcasn´e dobˇe kolem 2 milion˚ u infikovan´ ych jedinc˚ u. Z´avaˇznost pandemie vedla mimo jin´e k nutnosti zab´ yvat se ot´azkou modelov´an´ı pr˚ ubˇehu a ˇs´ıˇren´ı choroby HIV. Pˇr´ıpadn´e modely mohou b´ yt konstruov´any na z´akladˇe deterministick´eho nebo stochastick´eho pˇr´ıstupu. Protoˇze vˇsak jevy, se kter´ ymi se v tˇechto modelech pracuje, jsou pevnˇe spjaty s pojmy jako variabilita, neurˇcitost, n´ahodnost apod., je pˇri formulov´an´ı model˚ u ˇcasto vhodnˇejˇs´ı stochastick´ y pˇr´ıstup. Hlavn´ım c´ılem t´eto pr´ace je popsat vyuˇzit´ı homogenn´ıch Markovov´ ych proces˚ u se spojit´ ym ˇcasem pˇri modelov´an´ı pr˚ ubˇehu choroby HIV. Pr´ace je rozdˇelena do tˇr´ı kapitol. V prvn´ı kapitole jsou uvedeny klinick´e projevy choroby HIV a zp˚ usoby jej´ıho pˇrenosu. Druh´a a tˇret´ı kapitola jsou vˇenov´any zejm´ena problematice stochastick´eho modelov´an´ı inkubaˇcn´ı doby AIDS (d´ale jen inkubaˇcn´ı doby). Ve druh´e kapitole jsou pops´ana nˇekter´a statistick´a rozdˇelen´ı, zamˇeˇrujeme se pˇredevˇs´ım na posouzen´ı jejich vhodnosti jako model˚ u inkubaˇcn´ı doby. Tˇret´ı kapitola se zab´ yv´a konstrukc´ı Markovov´ ych model˚ u. V u ´vodn´ı ˇc´asti je shrnuta z´akladn´ı teorie potˇrebn´a pˇri modelov´an´ı. Pr˚ ubˇeh choroby HIV je modelov´an pomoc´ı dvou konkr´etn´ıch model˚ u, pozornost je vˇenov´ana zejm´ena odvozen´ı hustoty inkubaˇcn´ı doby vˇcetnˇe jej´ıho odhadu a odvozen´ı pravdˇepodobnost´ı pˇrechodu. D´ale je v pr´aci pops´an zobecnˇen´ y model a v tomto modelu odvozen explicitn´ı tvar hustoty inkubaˇcn´ı doby. V d˚ ukazu tohoto odvozen´ı lze nal´ezt postup, jak k tomuto vyj´adˇren´ı doj´ıt. Nicm´enˇe, aplikace tˇechto v´ ysledk˚ u v praxi je numericky pomˇernˇe n´aroˇcn´a.
6
Kapitola 1 Choroba zp˚ usoben´ a virem HIV Choroba HIV zaˇc´ın´a obvykle jako lehk´e chˇripkov´e onemocnˇen´ı (tzv. akutn´ı infekce). Po 2 aˇz 3 t´ ydnech nast´av´a latentn´ı f´aze, kter´a trv´a nˇekolik let. V tomto obdob´ı pacienti nem´ıvaj´ı ˇz´adn´e vnˇejˇs´ı pˇr´ıznaky, doch´az´ı vˇsak u nich k postupn´emu oslabov´an´ı imunitn´ıho syst´emu. D˚ usledkem toho se objevuj´ı nˇekter´e obt´ıˇze (zduˇren´e lymfatick´e uzliny, u ´nava, ztr´ata tˇelesn´e hmotnosti, pr˚ ujem), kter´e jsou pˇr´ıznakem dalˇs´ı, tzv. ARC f´aze (AIDS-Related Complex). Zhorˇsuj´ıc´ı se obranyschopnost organizmu vede k selh´av´an´ı imunity a rozv´ıj´ı se onemocnˇen´ı AIDS (Acquired Immunodeficiency Syndrome). Tato f´aze je charakterizov´ana rozvojem jedn´e nebo v´ıce chorob v d˚ usledku tˇeˇzk´eho poˇskozen´ı imunitn´ıho syst´emu (napˇr. pneumonie, tuberkul´oza, nˇekter´e nervov´e poruchy, nˇekter´e typy n´ador˚ u). Pacienti s projevy rozvinut´eho AIDS pˇreˇz´ıvaj´ı obvykle 1 aˇz 2 roky. P˚ uvodcem nemoci je virus HIV (Human Immunodeficiency Virus) patˇr´ıc´ı mezi retroviry. V napaden´em organizmu tento virus pˇrednostnˇe infikuje CD4+ T lymfocyty (d´ale jen T4 buˇ nky). V nˇekolika prvn´ıch mˇes´ıc´ıch po n´akaze se poˇcet T4 bunˇek v´ yraznˇe sn´ıˇz´ı z norm´aln´ıho stavu kolem 1125 T4 bunˇek v µl krve na asi 800 bunˇek v µl krve. V dalˇs´ım pr˚ ubˇehu se u ´bytek tˇechto bunˇek zpomal´ı, nicm´enˇe pokraˇcuje d´ale. Mnoˇzstv´ı T4 bunˇek v krvi proto slouˇz´ı jako v´ yznamn´ y ukazatel pr˚ ubˇehu choroby HIV. Dosud jsou zn´amy dvˇe formy viru HIV: HIV-1 a HIV-2. Virus HIV-1 byl identifikov´an jako prvn´ı a je pˇr´ıˇcinou vˇetˇsiny onemocnˇen´ı HIV v USA a dalˇs´ıch rozvinut´ ych zem´ıch. Virus HIV-2 byl poprv´e izolov´an u pacient˚ u v z´apadn´ı Africe. V porovn´an´ı s HIV-1 je m´enˇe patogenick´ y a pr˚ ubˇeh nemoci je m´ırnˇejˇs´ı. HIV-2 je bl´ızce pˇr´ıbuzn´ y jin´emu viru oznaˇcovan´emu jako SIV (Simian Immunodeficiency Virus). Tento virus byl objeven v krvi nˇekter´ ych druh˚ u africk´ ych opic. Podle vˇseobecnˇe pˇrij´ıman´e domnˇenky maj´ı SIV a viry HIV spoleˇcn´eho pˇredch˚ udce. Choroba HIV se ˇs´ıˇr´ı nechr´anˇen´ ym pohlavn´ım stykem, u narkoman˚ u spoleˇcn´ ym uˇz´ıv´an´ım injekˇcn´ıch jehel, moˇzn´ y je i pˇrenos z infikovan´e ˇzeny na jej´ı plod. V´ıce informac´ı o viru HIV, chorobˇe HIV a nemoci AIDS lze nal´ezt na str´ank´ach 3. l´ekaˇrsk´e fakulty Univerzity Karlovy [10] nebo v [7, str. 1–20]. 7
Kapitola 2 Modelov´ an´ı inkubaˇ cn´ıch a infekˇ cn´ıch dob 2.1
Stochastick´ e modely
Ve druh´e kapitole se budeme vˇenovat problematice stochastick´eho modelov´an´ı inkubaˇcn´ı a infekˇcn´ı doby choroby HIV a inkubaˇcn´ı doby AIDS. Jelikoˇz se definice tˇechto dob mohou v odborn´e literatuˇre liˇsit, upˇresn´ıme nyn´ı pojmy, se kter´ ymi budeme d´ale pracovat. Inkubaˇcn´ı dobou choroby HIV oznaˇc´ıme obdob´ı mezi infikac´ı organizmu virem HIV a prvn´ımi pˇr´ıznaky tzv. akutn´ı infekce. Toto obdob´ı trv´a obvykle 1 aˇz 3 mˇes´ıce. Infekˇcn´ı dobou choroby HIV budeme naz´ yvat obdob´ı, bˇehem nˇehoˇz m˚ uˇze infikovan´ y jedinec pˇren´aˇset chorobu na ostatn´ı. Toto obdob´ı zaˇc´ın´a t´emˇeˇr okamˇzitˇe po vniknut´ı viru do organizmu a trv´a aˇz do konce pacientova ˇzivota. Inkubaˇcn´ı dobou AIDS oznaˇc´ıme dobu mezi n´akazou organizmu virem HIV a rozvojem klinick´eho st´adia AIDS. (Pr´avˇe toto obdob´ı je nˇekdy oznaˇcov´ano jako inkubaˇcn´ı doba choroby HIV.) D´elky v´ yˇse definovan´ ych dob se mohou u jednotliv´ ych jedinc˚ u liˇsit. C´ılem t´eto podkapitoly je zkoumat, jak tuto variabilitu v d´elce dob popsat pomoc´ı nˇejak´eho vhodn´eho statistick´eho rozdˇelen´ı. Protoˇze postup, kter´ y zde pop´ıˇseme, lze aplikovat na vˇsechny tˇri v´ yˇse definovan´e doby (a obecnˇe pˇri studiu jak´ekoli infekˇcn´ı nemoci), budeme se vˇenovat pouze modelov´an´ı inkubaˇcn´ı doby AIDS (d´ale jen inkubaˇcn´ı doby). Inkubaˇcn´ı dobu budeme reprezentovat spojitou n´ahodnou veliˇcinou T s oborem hodnot RT = [0, ∞) pˇredstavuj´ıc´ı ˇcas. Definujme jej´ı distribuˇcn´ı funkci F (t) = P (T ≤ t), t ∈ R a pˇredpokl´adejme, ˇze plat´ı F (0) = 0. Hustotu veliˇciny T oznaˇc´ıme f (t) a funkci pˇreˇzit´ı S(t). Pˇri hled´an´ı vhodn´eho rozdˇelen´ı veliˇciny T je tˇreba m´ıt na pamˇeti, ˇze inkubaˇcn´ı doba je relativnˇe dlouh´a (pr˚ umˇernˇe kolem 10 let) a znaˇcnˇe variabiln´ı. Je tak´e pˇrirozen´e oˇcek´avat, ˇze ˇc´ım d´ele inkubaˇcn´ı doba trv´a, t´ım vyˇsˇs´ı je podm´ınˇen´a pravdˇepodobnost, ˇze se v dalˇs´ım okamˇziku rozvine klinick´e st´adium AIDS. 8
Nyn´ı uvedeme nˇekter´a statistick´a rozdˇelen´ı a posoud´ıme, zda jsou jako modely inkubaˇcn´ı doby vhodn´a. U nˇekter´ ych rozdˇelen´ı zm´ın´ıme tvar rizikov´e funkce, tedy funkce, pro kterou plat´ı f (t) F (t + h) − F (t) P (t < T ≤ t + h | T > t) = lim = h→0+ h→0+ h hS(t) S(t)
θ(t) = lim
pro vˇsechna t ∈ R, pro kter´a limita existuje. Tento vztah lze ekvivalentnˇe vyj´adˇrit diferenci´aln´ı rovnic´ı d ln S(t) dt s poˇc´ateˇcn´ı podm´ınkou S(0) = 1. Protoˇze v naˇsem pˇr´ıpadˇe je f (t) = 0 pro vˇsechna t ≤ 0, budeme hustotou nad´ale rozumˇet hustotu na intervalu (0, ∞). Obdobnˇe pro rizikovou funkci. θ(t) = −
Weibullovo a exponenci´ aln´ı rozdˇ elen´ı Necht’ α > 0, β > 0. Weibullovo rozdˇelen´ı W(α, β) m´a hustotu danou f (t) = αβtα−1 exp (−βtα ) , t > 0. Rizikov´a funkce m´a potom tvar θ(t) = αβtα−1 , t > 0. Vlastnosti rizikov´e funkce urˇcuje pˇredevˇs´ım hodnota parametru α. Je-li α > 1, je rizikov´a funkce rostouc´ı. T´ımto pˇr´ıpadem se budeme zab´ yvat v dalˇs´ı podkapitole. Naopak pro 0 < α < 1 je rizikov´a funkce klesaj´ıc´ı, takov´e rozdˇelen´ı proto nem˚ uˇze b´ yt pro modelov´an´ı inkubaˇcn´ı doby vhodn´e. Je-li α = 1, jedn´a se o exponenci´ aln´ı rozdˇelen´ı Exp(β) s konstantn´ı rizikovou funkc´ı. Pˇrestoˇze se exponenci´aln´ı rozdˇelen´ı ˇcasto pouˇz´ıv´a k modelov´an´ı n´ahodn´e doby ˇcek´an´ı na nˇejakou ud´alost, jako model inkubaˇcn´ı doby nen´ı pˇr´ıliˇs vhodn´e. Exponenci´aln´ı rozdˇelen´ı je totiˇz rozdˇelen´ım bez pamˇeti: M´a-li veliˇcina T exponenci´aln´ı rozdˇelen´ı, plat´ı P (T > t + h | T > t) = P (T > h) pro kaˇzd´e t > 0, h > 0. Gama rozdˇ elen´ı Necht’ α > 0, β > 0. Gama rozdˇelen´ı Ga(α, β) je urˇceno hustotou f (t) =
β α α−1 t exp(−βt), t > 0, Γ(α)
kde Γ(α) je gama funkce definovan´a pro α > 0. Je-li α = 1, z´ısk´ame stejnˇe jako u Weibullova rozdˇelen´ı exponenci´aln´ı rozdˇelen´ı Exp(β). Pro 0 < α < 1 je rizikov´a funkce gama rozdˇelen´ı klesaj´ıc´ı a pro α > 1 je rostouc´ı. V obou pˇr´ıpadech rizikov´a funkce konverguje k asymptotˇe β pro t → ∞. Pro velk´a t m´a tedy tato funkce podobn´e vlastnosti jako rizikov´a funkce rozdˇelen´ı Exp(β). 9
Log-logistick´ e rozdˇ elen´ı Necht’ α ∈ R, β > 0. Logistick´e rozdˇelen´ı L(α, β) m´a hustotu tvaru g(x) =
β exp[−(α + βx)] , x ∈ R. {1 + exp[−(α + βx)]}2
Takov´e rozdˇelen´ı nen´ı k modelov´an´ı inkubaˇcn´ı doby vhodn´e, protoˇze jeho hustota je kladn´a vˇsude na R, zat´ımco veliˇcina reprezentuj´ıc´ı inkubaˇcn´ı dobu je nez´aporn´a. M˚ uˇzeme vˇsak uvaˇzovat log-logistick´e rozdˇelen´ı. Definujme veliˇcinu T vztahem T = exp(X), kde X ∼ L(α, β). Potom m´a veliˇcina T log-logistick´e rozdˇelen´ı, tedy ln T ∼ L(α, β), a pro hustotu T plat´ı f (t) =
dF (t) dP (T ≤ t) dP (ln T ≤ ln t) dG(ln t) 1 = = = = g(ln t), t > 0, dt dt dt dt t
kde g je hustota veliˇciny X a G, F jsou distribuˇcn´ı funkce veliˇcin X,T v tomto poˇrad´ı. Log-norm´ aln´ı rozdˇ elen´ı Necht’ µ ∈ R, σ > 0. Uvaˇzujme norm´aln´ı rozdˇelen´ı N(µ, σ 2 ) se stˇredn´ı hodnotou µ a rozptylem σ 2 . Stejnˇe jako logistick´e rozdˇelen´ı m´a norm´aln´ı rozdˇelen´ı hustotu kladnou na cel´em R. Necht’ veliˇcina X m´a rozdˇelen´ı N(µ, σ 2 ) s hustotou g. Potom n´ahodn´a veliˇcina T = exp(X) m´a log-norm´aln´ı rozdˇelen´ı. Jej´ı hustotu odvod´ıme analogicky jako u log-logistick´eho rozdˇelen´ı, tedy je 1 1 (ln t − µ)2 , t > 0. f (t) = g(ln t) = √ exp − t 2σ 2 2π σt Log-Cauchyovo rozdˇ elen´ı Necht’ α ∈ R, β > 0. Cauchyovo rozdˇelen´ı C(α, β) m´a hustotu g(x) =
β 1 , x ∈ R. 2 π β + (x − α)2
Pro toto rozdˇelen´ı neexistuje koneˇcn´a stˇredn´ı hodnota. M´a-li veliˇcina X Cauchyovo rozdˇelen´ı, obdobnˇe jako v pˇredchoz´ıch dvou pˇr´ıpadech definujeme nez´apornou veliˇcinu T = exp(X) s log-Cauchyov´ym rozdˇelen´ım.
2.2
Vlastnosti model˚ u
Pˇri posuzov´an´ı, jak jsou jednotliv´a rozdˇelen´ı vhodn´a k modelov´an´ı inkubaˇcn´ı doby, je tˇreba vz´ıt v u ´vahu nˇekolik skuteˇcnost´ı. Pˇredevˇs´ım je nutn´e si uvˇedomit, ˇze obvykle m´ame k dispozici pouze mal´e soubory dat a ˇze pozorov´an´ı, ze kter´ ych tato data poch´az´ı, prob´ıhala relativnˇe kr´atkou dobu v porovn´an´ı s d´elkou inkubaˇcn´ı doby. D˚ usledkem toho mohli b´ yt z tˇechto pozorov´an´ı vylouˇceni jedinci s pˇrirozenˇe dlouhou inkubaˇcn´ı dobou. Proto 10
m˚ uˇze b´ yt zjiˇstˇen´a pr˚ umˇern´a d´elka inkubaˇcn´ı doby podhodnocena. Je tedy dobr´e vˇedˇet, zda a v jak´e m´ıˇre jednotliv´a rozdˇelen´ı pˇripouˇstˇej´ı odlehl´a pozorov´an´ı. D´ale je tˇreba m´ıt na pamˇeti, ˇze pr˚ ubˇeh nemoci ovlivˇ nuje cel´a ˇrada faktor˚ u. Zmiˇ nme napˇr´ıklad rozd´ıln´ y pr˚ ubˇeh chorob HIV-1 a HIV-2. Pr˚ ubˇeh choroby zˇrejmˇe z´avis´ı i na vˇeku pacienta a pravdˇepodobnˇe tak´e na pohlav´ı a zp˚ usobu n´akazy. Nav´ıc novˇe vyv´ıjen´e l´eky dok´aˇz´ı podstatnˇe prodlouˇzit trv´an´ı inkubaˇcn´ı doby. Jednotliv´a rozdˇelen´ı m˚ uˇzeme porovn´avat pomoc´ı kvantil˚ u. Parametry porovn´avan´ ych rozdˇelen´ı lze zvolit r˚ uzn´ ymi zp˚ usoby. U nˇekter´ ych rozdˇelen´ı m˚ uˇzeme postupovat tak, ˇze vyj´adˇr´ıme jejich stˇredn´ı hodnoty (pˇr´ıpadnˇe i rozptyly) jako funkce nezn´am´ ych parametr˚ u. Hodnoty tˇechto parametr˚ u zvol´ıme tak, aby stˇredn´ı hodnoty byly pro vˇsechna porovn´avan´a rozdˇelen´ı stejn´e. Pro nˇekter´a rozdˇelen´ı vˇsak nelze stˇredn´ı hodnotu element´arnˇe vyj´adˇrit a u nˇekter´ ych rozdˇelen´ı nemus´ı koneˇcn´a stˇredn´ı hodnota v˚ ubec existovat. V takov´em pˇr´ıpadˇe m˚ uˇzeme hodnoty parametr˚ u urˇcit tak, aby se pro porovn´avan´a rozdˇelen´ı shodovaly vybran´e kvanily. Nyn´ı porovn´ame nˇekter´a rozdˇelen´ı (pro zvolen´e parametry) z hlediska odlehl´ ych pozorov´an´ı. Pˇ r´ıklad 2.1. V tomto pˇr´ıkladu porovn´ame Weibullovo rozdˇelen´ı s gama rozdˇelen´ım, rozdˇelen´ım exponenci´aln´ım, log-norm´aln´ım, log-logistick´ ym a log-Cauchyov´ ym. Parametry tˇechto rozdˇelen´ı zvol´ıme tak, aby vˇsechna rozdˇelen´ı mˇela stejn´ y medi´an a pˇr´ıpadnˇe tak´e stejn´ y tˇret´ı kvartil (0,75kvantil). U rozdˇelen´ı s takto zvolen´ ymi parametry urˇc´ıme nˇekter´e dalˇs´ı kvantily a pomoc´ı nich rozdˇelen´ı vz´ajemnˇe porovn´ame. Zd˚ uraznˇeme, ˇze takto urˇcen´e parametry zde slouˇz´ı pouze pro ilustrativn´ı porovn´an´ı. Hodnoty budeme v pˇr´ıkladu zaokrouhlovat na 3 desetinn´a m´ısta Vyjdeme z v´ yzkumu, kter´ y provedli Lui a kol. Ti sledovali d´elku inkubaˇcn´ı doby AIDS u skupiny homosexu´aln´ıch muˇz˚ u. Inkubaˇcn´ı dobu modelovali pomoc´ı Weibullova rozdˇelen´ı W(α, β) a jeho parametry odhadli z dat metodou maxim´aln´ı vˇerohodnosti: α b = 2, 571000 a βb = 0, 003807. b medi´an tohoto rozdˇelen´ı Uvaˇzujme tedy Weibullovo rozdˇelen´ı W(b α, β), oznaˇc´ıme w0,50 a tˇret´ı kvartil w0,75 . U exponenci´aln´ıho rozdˇelen´ı Exp(β) zvol´ıme parametr β tak, aby byl pˇr´ısluˇsn´ y medi´an roven w0,50 , tedy β = = ln 2/w0,50 . U gama rozdˇelen´ı Ga(α, β) budeme vzhledem k poˇcetn´ı n´aroˇcnosti poˇzadovat pouze shodu jeho medi´anu s w0,50 . Z teorie matematick´e statistiky (transformace n´ahodn´ ych veliˇcin) je zn´amo, ˇze, m´a-li veliˇcina X rozdˇelen´ı Ga(α, β), potom veliˇcina Y = kX, kde k > 0, m´a rozdˇelen´ı Ga(α, β/k). Oznaˇcme F , G distribuˇcn´ı funkce veliˇcin X, Y v tomto poˇrad´ı, xq q-kvantil veliˇciny X a yq q-kvantil veliˇciny Y , q ∈ (0, 1). Jelikoˇz plat´ı y yq q G(yq ) = P (Y ≤ yq ) = P X ≤ =F = F (xq ) , k k 11
Tabulka 2.1: Vybran´e q-kvantily v roc´ıch pro Weibullovo, gama, exponenci´aln´ı, log-norm´aln´ı, log-logistick´e a log-Cauchyovo rozdˇelen´ı se stejn´ ymi medi´any q 0,25 0,50 0,75 0,95 0,99 0,999
Weibull. 5,378 7,571 9,913 13,378 15,813 18,514
gama 5,634 7,571 9,911 14,037 17,503 21,971
exponen. 3,142 7,571 15,141 32,720 50,298 75,447
log-norm. 5,782 7,571 9,913 14,610 19,185 26,036
log-logis. 5,782 7,571 9,913 15,593 23,381 41,232
log-Cauchy. 5.782 7,571 9,913 41,532 4, 026 · 104 1, 409 · 1038
jsou kvantily yq urˇceny vztahem yq = kxq . Chceme-li pro veliˇcinu Y urˇcit parametry tak, aby platilo y0,50 = w0,50 , m˚ uˇzeme postupovat n´asleduj´ıc´ım ˜ ˜ zp˚ usobem: Zvol´ıme nˇejak´e α ˜ > 0, β > 0 a pro veliˇcinu X ∼ Ga(˜ α, β) urˇc´ıme kvantil x0,50 (napˇr. pomoc´ı programu R). Veliˇcinu Y definujeme potom vztahem Y = kX, kde konstanta k splˇ nuje rovnici kx0,50 = w0,50 . My ˜ jsme zvolili α ˜ = 6 a β = 1, tˇret´ı kvartil y0,75 m´a v tomto pˇr´ıpadˇe hodnotu bl´ızkou hodnotˇe w0,75 . Poznamenejme jeˇstˇe, ˇze jsme z´amˇernˇe volili α ˜ > 1. Pro α ˜ = 1 bychom totiˇz dostali exponenci´an´ı rozdˇelen´ı, kter´ ym jsme se jiˇz zab´ yvali v´ yˇse. Pro α ˜ < 1 bychom dostali gama rozdˇelen´ı s klesaj´ıc´ı rizikovou funkc´ı, takov´e rozdˇelen´ı nen´ı k modelov´an´ı inkubaˇcn´ı doby AIDS vhodn´e. U log-norm´aln´ıho, log-logistick´eho a log-Cauchyova rozdˇelen´ı budeme poˇzadovat, aby se shodovaly jejich medi´any s w0,50 a tˇret´ı kvartily s w0,75 . Pro norm´aln´ı, logistick´e i Cauchyovo rozdˇelen´ı plat´ı, ˇze line´arn´ı transformac´ı se typ rozdˇelen´ı nezmˇen´ı. Toho nyn´ı vyuˇzijeme. Urˇc´ıme pouze hledan´e log-norm´aln´ı rozdˇelen´ı, pˇri odvozov´an´ı log-logistick´eho a log-Cauchyova rozdˇelen´ı lze postupovat analogicky. Uvaˇzujme veliˇcinu Z ∼ N (0, 1), jej´ı kvantily oznaˇc´ıme zq . Pro libovoln´e α, β ∈ R m´a veliˇcina α + βZ norm´aln´ı rozdˇelen´ı a nez´aporn´a veliˇcina T = exp(α + βZ) odpov´ıdaj´ıc´ı log-norm´aln´ı rozdˇelen´ı. Kvantily tq veliˇciny T splˇ nuj´ı vztah tq = exp(α + βzq ), nebot’ pro tq > 0 plat´ı
H(tq ) = P (T ≤ tq ) = P
ln tq − α Z≤ β
=F
ln tq − α β
= F (zq ),
kde F , H znaˇc´ı distribuˇcn´ı funkce veliˇcin Z a T v dan´em poˇrad´ı. Kvantily z0,50 a z0,75 urˇc´ıme pomoc´ı programu R. Konstanty α, β mus´ı splˇ novat rovnice exp(α + βz0,50 ) = w0,50 , exp(α + βz0,75 ) = w0,75 , odkud jiˇz snadno α a β vypoˇcteme. ´ Uloha je vyˇreˇsena pomoc´ı softwaru R. Poznamenejme, ˇze program R f (α, γ) s parametry α, γ = β − α1 , pracuje s Weibullov´ ym rozdˇelen´ım W kde α, β znaˇc´ı parametry n´ami uvaˇzovan´eho rozdˇelen´ı W(α, β). V´ ysledky 12
u ´lohy jsou zaps´any v tabulce 2.1. Porovn´an´ım druh´eho a tˇret´ıho sloupce tabulky je vidˇet, ˇze pro vhodnˇe zvolen´e parametry maj´ı Weibullovo a gama rozdˇelen´ı podobn´e kvantily. Porovn´ame-li s tˇemito rozdˇelen´ımi exponenci´aln´ı rozdˇelen´ı ve ˇctvrt´em sloupci, vid´ıme, ˇze pro dan´ y parametr jsou jeho qkvantily pro q ≥ 0.95 nˇekolikan´asobnˇe vyˇsˇs´ı neˇz odpov´ıdaj´ıc´ı q-kvantily Weibullova a gama rozdˇelen´ı. Ve vzorku generovan´em z takto urˇcen´eho rozdˇelen´ı by pˇribliˇznˇe 1% dat pˇrekroˇcilo hodnotu 50 rok˚ u. Z v´ ysledk˚ u d´ale plyne, ˇze log-norm´aln´ı a pˇredevˇs´ım log-logistick´e rozdˇelen´ı obsahuj´ı pro dan´e parametry v´ıce odlehl´ ych pozorov´an´ı neˇz uvaˇzovan´e Weibullovo rozdˇelen´ı. Tato b l´epe zohledˇ rozdˇelen´ı tedy v porovn´an´ı s rozdˇelen´ım W(b α, β) nuj´ı jedince s velmi dlouhou inkubaˇcn´ı dobou. Log-Cauchyovo rozdˇelen´ı v posledn´ım sloupci k modelov´an´ı inkubaˇcn´ı doby AIDS vhodn´e nen´ı, nebot’ obsahuje v´ yraznˇe odlehl´a pozorov´an´ı. Ve vzorku generovan´em z tohoto rozdˇelen´ı by bylo pˇribliˇznˇe 5% dat vˇetˇs´ıch neˇz 41,5 rok˚ u a pˇribliˇznˇe 1% dat by dokonce pˇres´ahlo hodnotu 40260 rok˚ u, tedy dobu, kter´a v´ yraznˇe pˇrekraˇcuje d´elku lidsk´eho ˇzivota. ♦ Dalˇs´ı moˇznosti, jak porovn´avat statistick´a rozdˇelen´ı, pˇredevˇs´ım z hlediska odlehl´ ych pozorov´an´ı, lze nal´ezt v [7, str. 37–42].
13
Kapitola 3 Markovovy procesy a jejich aplikace pˇ ri modelov´ an´ı 3.1
Teorie Markovov´ ych proces˚ u
V´ yznamnou roli mezi stochastick´ ymi modely pouˇz´ıvan´ ymi v epidemiologii hraj´ı n´ahodn´e procesy, mezi nˇe patˇr´ı i tzv. Markovovy procesy. Definice 3.1. Necht’ (Ω, A, P ) je pravdˇepodobnostn´ı prostor a T ⊂ R. Syst´em n´ahodn´ ych veliˇcin {X(t), t ∈ T } definovan´ ych na (Ω, A, P ) se naz´ yv´a n´ ahodn´y proces. Parametr t m´a v´ yznam ˇcasu. Je-li T mnoˇzina diskr´etn´ıch hodnot, hovoˇr´ıme o procesu s diskr´etn´ım ˇcasem. Je-li T interval, hovoˇr´ıme o procesu se spojit´ym ˇcasem. Hodnoty n´ahodn´ ych veliˇcin X(t), t ∈ T se naz´ yvaj´ı stavy, jejich mnoˇzinu budeme znaˇcit S. Nab´ yvaj´ı-li veliˇciny X(t) pouze diskr´etn´ıch hodnot, ˇr´ık´ame, ˇze {X(t), t ∈ T } je proces s diskr´etn´ımi stavy. Je-li S interval, hovoˇr´ıme o procesu se spojit´ymi stavy. V dalˇs´ı ˇc´asti se budeme vˇenovat Markovov´ ym proces˚ um se spojit´ ym ˇcasem a diskr´etn´ımi stavy. Pˇripomeneme pouze z´akladn´ı terminologii a nˇekter´e d˚ uleˇzit´e vlastnosti, kter´e vyuˇzijeme v druh´e podkapitole pˇri modelov´an´ı. Ucelenou teorii o tˇechto procesech lze nal´ezt napˇr. v [9]. Definice 3.2. Necht’ je d´an pravdˇepodobnostn´ı prostor (Ω, A, P ). Syst´em celoˇc´ıseln´ ych n´ahodn´ ych veliˇcin {X(t), t ≥ 0} definovan´ ych na (Ω, A, P ) se naz´ yv´a Markov˚ uv proces se spojit´ ym ˇcasem a mnoˇzinou stav˚ u S, jestliˇze P (X(t) = j | X(s) = i, X(tn ) = in , . . . , X(t0 ) = i0 ) = = P (X(t) = j | X(s) = i)
(3.1)
pro libovoln´e n ∈ N0 , pro vˇsechna i0 , . . . , in , i, j ∈ S a pro vˇsechna 0 ≤ t0 < < t1 < · · · < tn < s < t, pro kter´a je P (X(s) = i, X(tn ) = in , . . . , X(t0 ) = = i0 ) > 0. 14
Vztah (3.1) vyjadˇruje tzv. markovskou vlastnost, kter´a je z´akladn´ı charakteristikou Markovov´ ych proces˚ u. Podm´ınˇen´e pravdˇepodobnosti P (X(t) = = j | X(s) = i) = pij (s, t) (jsou-li definov´any) budeme naz´ yvat pravdˇepodobnosti pˇrechodu ze stavu i v ˇcase s do stavu j v ˇcase t. Pokud pro kaˇzd´e i, j ∈ S plat´ı pij (s, t) = pij (t − s), ∀s, t : 0 ≤ s < t, ˇr´ık´ame, ˇze Markov˚ uv proces je homogenn´ı. Nad´ale se budeme vˇenovat pouze teorii homogenn´ıch proces˚ u. Homogenn´ı Markov˚ uv proces se spojit´ ym ˇcasem budeme v t´eto podkapitole d´ale oznaˇcovat jako (CTHMP, continuous time homogeneous Markov process). Pravdˇepodobnosti pˇrechodu m˚ uˇzeme pro kaˇzd´e t ≥ 0 sestavit do tzv. matice pravdˇepodobnost´ı pˇrechodu P(t) = {pij (t), i, j ∈ S}. Zˇrejmˇe pro kaˇzd´e i, j ∈ S a kaˇzd´e t ≥ 0 plat´ı X pij (t) ≥ 0; pij (t) = 1; pij (0) = δij , j∈S
kde δij znaˇc´ı Kroneckerovo delta. V d˚ usledku markovsk´e vlastnosti splˇ nuj´ı pravdˇepodobnosti pˇrechodu tzv. Chapmanovu-Kolmogorovovu rovnici.
Vˇ eta 3.1 (Chapmanova-Kolmogorovova rovnice). Necht’ je d´an (CTHMP) s mnoˇzinou diskr´etn´ıch stav˚ u S. Potom pro pravdˇepodobnosti pˇrechodu plat´ı X pij (s + t) = pik (s)pkj (t), ∀i, j ∈ S, ∀s > 0, ∀t > 0. k∈S
D˚ ukaz. Viz [3, str. 115]. Protoˇze pˇri modelov´an´ı budeme pracovat pouze s procesy s koneˇcnou mnoˇzinou stav˚ u, v dalˇs´ı teorii se omez´ıme na tento pˇr´ıpad. D˚ uleˇzitou roli v teorii Markovov´ ych proces˚ u se spojit´ ym ˇcasem hraj´ı Kolmogorovy diferenci´ aln´ı rovnice. Ty lze odvodit z Chapmanovy-Kolmogorovy rovnice. K tomu je vˇsak nejprve nutn´e definovat tzv. intenzity pˇrechodu. Vˇ eta 3.2. Necht’ je d´an (CTHMP) s koneˇcnou mnoˇzinou stav˚ u S. Potom pro kaˇzd´e i, j ∈ S existuj´ı limity 1 − pii (∆) = −qii , − ∞ < qii ≤ 0, ∆→0+ ∆ pij (∆) = qij , 0 ≤ qij < ∞, kde i 6= j, lim ∆→0+ ∆ lim
15
a pro takto definovan´e konstanty qij plat´ı gii = −
P
j6=i qij .
D˚ ukaz. D˚ ukaz existence konstant P qii a qij , j 6= i je uveden v [4, str. 131–133], vˇety 4 a 5. D˚ ukaz vztahu gii = − j6=i qij lze nal´ezt v [9, str. 74]. Definice 3.3. Konstanty qij , i, j ∈ S definovan´e v pˇredchoz´ı vˇetˇe se naz´ yvaj´ı intenzity pˇrechodu ze stavu i do stavu j. Konstanty qi , i ∈ S definovan´e pˇredpisem qi = −qii se naz´ yvaj´ı celkov´e intenzity. Matici Q = {qij , i, j ∈ S} ˇr´ık´ame matice intenzit pˇrechodu.
Vˇ eta 3.3 (Kolmogorovy diferenci´ aln´ı rovnice). Necht’ je d´an (CTHMP) s koneˇcnou mnoˇzinou stav˚ u S a intenzitami pˇrechodu qij , i, j ∈ S. Potom pravdˇepodobnosti pˇrechodu pij (t) splˇ nuj´ı pro vˇsechna i, j ∈ S a t > 0 syst´emy diferenci´ aln´ıch rovnic p′ij (t) =
X
pik (t)qkj (prospektivn´ı rovnice),
(3.2)
qik pkj (t) (retrospektivn´ı rovnice)
(3.3)
k∈S
p′ij (t)
=
X k∈S
s poˇc´ateˇcn´ı podm´ınkou pij (0) = δij . D˚ ukaz. Viz vˇeta 3.9 v [9, str. 82]. Pˇrirozenˇe n´as zaj´ım´a, zda a jak´e maj´ı soustavy (3.2) a (3.3) ˇreˇsen´ı. Nejprve vˇsak pˇripomeˇ nme pojem exponenci´ aln´ı maticov´e funkce. Tato funkce je definov´ana pro libovolnou ˇctvercovou matici B vztahem exp(Bt) =
∞ X Bk tk k=0
k!
,
kde B0 je jednotkov´a matice odpov´ıdaj´ıc´ıho ˇr´adu. Nyn´ı jiˇz m˚ uˇzeme uv´est vˇetu t´ ykaj´ıc´ı se ˇreˇsen´ı soustav (3.2) a (3.3). Vˇ eta 3.4. Necht’ je d´ana ˇctvercov´ a matice P koneˇcn´eho ˇr´adu Q = {qij }, pro jej´ıˇz prvky plat´ı qij ≥ 0, i 6= j, qii = − j6=i qij . Pak existuje jedin´e ˇreˇsen´ı, spoleˇcn´e pro obˇe soustavy (3.2) a (3.3), s poˇc´ateˇcn´ı podm´ınkou P(0) = I, kde I znaˇc´ı jednotkovou matici pˇr´ısluˇsn´eho ˇr´adu. Toto ˇreˇsen´ı pˇredstavuje syst´em pravdˇepodobnost´ı pˇrechodu (CTHMP) s koneˇcnou mnoˇzinou stav˚ u. Maticovˇe lze ˇreˇsen´ı zapsat ve tvaru P(t) = exp(Qt).
16
D˚ ukaz. Viz vˇeta 3.10 v [9, str. 85]. Na z´akladˇe intenzit pˇrechodu m˚ uˇzeme klasifikovat stavy procesu. Zadefinujme nejprve ˇcas prvn´ıho n´ avratu do stavu i ∈ S jako τi = inf {t ≥ J : X(t) = i}, kde J znaˇc´ı ˇcasov´ y okamˇzik, kdy doˇslo k prvn´ımu pˇrechodu, tj. J = inf {t > 0 : X(t) 6= X(0)}. Definice 3.4. Necht’ je d´an (CTHMP) s koneˇcnou mnoˇzinou stav˚ u S. Stav i ∈ S se naz´ yv´a trval´y, jestliˇze bud’ qi = 0 (absorpˇcn´ı stav), nebo qi > 0 a z´aroveˇ n P (τi < ∞ | X(0) = i) = 1. Stav i ∈ S se naz´ yv´a pˇrechodn´y, jestliˇze qi > 0 a z´aroveˇ n P (τi < ∞ | X(0) = i) < 1. Vˇ eta 3.5. Necht’ je d´an (CTHMP) s koneˇcnou mnoˇzinou stav˚ u S. Je-li qi > 0, m´a doba, po kterou proces setrv´ a ve stavu i ∈ S, exponenci´ aln´ı rozdˇelen´ı Exp(qi ). Je-li qi = 0, potom pii (t) = 1 pro vˇsechna t ≥ 0. D˚ ukaz. Viz vˇeta 3.5 v [9, str. 77]. Nyn´ı se jeˇstˇe vr´at´ıme ke Kolmogorovov´ ym diferenci´aln´ım rovnic´ım. C. L. Chiang odvodil za urˇcit´ ych pˇredpoklad˚ u jejich explicitn´ı ˇreˇsen´ı. Lemma 3.1. Uvaˇzujme (CTHMP) s koneˇcnou mnoˇzinou stav˚ u S a intenzitami qij , i, j ∈ S. Necht’ plat´ı S = U ∪ T , kde U je nepr´ azdn´ a mnoˇzina pˇrechodn´ych stav˚ u a T je nepr´ azdn´ a mnoˇzina absorpˇcn´ıch stav˚ u. Bez u ´jmy na obecnosti pˇredpokl´adejme, ˇze je U = {1, . . . , n}. Definujme matici V = = {qij , i, j ∈ U }. Necht’ tato matice m´a re´aln´ a, navz´ajem r˚ uzn´a vlastn´ı ˇc´ısla ρ1 , . . . , ρn . Pak pravdˇepodobnosti pˇrechodu maj´ı pro t > 0 tvar
pij (t) =
n X l=1
pij (t) =
ATij (l) exp(ρl t), i, j ∈ U, m=1 (ρl − ρm )
Qn
m6=l n n XX l=1 k=1
ρl
ATik (l) [exp(ρl t) − 1] qkj , i ∈ U, j ∈ T, m=1 (ρl − ρm )
Qn
(3.4) (3.5)
m6=l
ATij (l) znaˇc´ı algebraick´y doplnˇek prvku aTij (l) v matici AT (l) = ρl I − VT . D˚ ukaz. Viz [3, str. 153–158] s odkazem na [3, str. 135–137].
17
3.2
Modelov´ an´ı pr˚ ubˇ ehu choroby HIV
Jedinec infikovan´ y virem HIV proch´az´ı r˚ uzn´ ymi st´adii choroby. Je proto pˇrirozen´e modelovat pr˚ ubˇeh nemoci pomoc´ı stavov´ ych model˚ u, jejichˇz stavy ztotoˇzn´ıme se st´adii nemoci. V t´eto podkapitole budeme pr˚ ubˇeh choroby HIV modelovat pomoc´ı homogenn´ıch Markovov´ ych proces˚ u se spojit´ ym ˇcasem a koneˇcnou mnoˇzinou stav˚ u. Moˇznost´ı, jak definovat syst´em st´adi´ı choroby HIV, je samozˇrejmˇe mnoho. Zde uvedeme pˇr´ıklady dvou syst´em˚ u, kter´e d´ale vyuˇzijeme pˇri konstrukci model˚ u. Prvn´ı syst´em vyuˇz´ıv´a f´aze popsan´e v prvn´ı kapitole, uveden je v tabulce 3.1. Druh´ y syst´em je pops´an v tabulce 3.2. St´adia choroby jsou v nˇem definov´ana na z´akladˇe poˇctu T4 bunˇek v krvi. Jak bylo uvedeno v prvn´ı kapitole, tyto buˇ nky virus HIV infikuje ihned po n´akaze a jejich poˇcet kles´a po celou dobu onemocnˇen´ı. Tabulka 3.1: Symptomatick´a st´adia choroby HIV Stadium E1 E2 E3 E4
Symptomy infikov´an virem HIV, ale s´eronegativn´ı s´eropozitivn´ı, ale bez zˇrejm´ ych pˇr´ıznak˚ u f´aze ARC klinick´e st´adium AIDS
Tabulka 3.2: St´adia choroby HIV zaloˇzen´a na poˇctu T4 bunˇek v krvi Stadium E1 E2 E3 E4 E5 E6
Poˇcet T4 bunˇek / mm3 v´ıce neˇz 899 700-899 500-699 350-499 200-349 0-199
Nyn´ı uvedeme pˇr´ıklady konkr´etn´ıch model˚ u pr˚ ubˇehu choroby HIV. V modelech budeme pro zjednoduˇsen´ı pˇredpokl´adat, ˇze se u kaˇzd´eho jedince infikovan´eho virem HIV ˇcasem rozvine AIDS. Pˇ r´ıklad 3.1. Model popsan´ y v tomto pˇr´ıkladu je pˇrevzat z [6]. Uvaˇzujme homogenn´ı Markov˚ uv proces {X(t), t ≥ 0} se spojit´ ym ˇcasem a mnoˇzinou stav˚ u S = {E1 , . . . , E5 }, kde E1 , . . . , E4 jsou stavy definovan´e v tabulce 3.1 a stav E5 reprezentuje u ´mrt´ı v d˚ usledku AIDS. Jin´e pˇr´ıˇciny smrti zde 18
nebudeme uvaˇzovat. Pˇredpokl´ad´ame, ˇze stavy E1 , . . . , E4 jsou pˇrechodn´e a stav E5 absorpˇcn´ı. Matice intenzit Q = {qij } necht’ m´a tvar
−α1 α1 0 0 0 0 −α2 α2 0 0 0 −α3 α3 0 Q= 0 , kde αi > 0 pro kaˇzd´e i = 1, . . . , 4. 0 0 0 −α4 α4 0 0 0 0 0
Pˇrechody mezi stavy jsou tedy urˇceny n´asleduj´ıc´ım diagramem: α
α
α
α
1 2 3 4 E1 −−−− −→ E2 −−−− −→ E3 −−−− −→ E4 − −−− − → E5 .
N´ahodnou dobu, kterou proces str´av´ı ve stavu Ei , oznaˇc´ıme Ti , kde i = = 1, . . . , 4. Jelikoˇz uvaˇzujeme homogenn´ı Markov˚ uv proces se spojit´ ym ˇcasem, jsou Ti nez´avisl´e n´ahodn´e veliˇciny, pro kter´e plat´ı Ti ∼ Exp(αi ). Veliˇcina Ti m´a tedy hustotu fi (t) = αi exp(−αi t), t > 0. Pˇr´ısluˇsn´e distribuˇcn´ı funkce oznaˇc´ıme Fi (t) a funkce pˇreˇzit´ı Si (t). Necht’ je proces v ˇcase 0 ve stavu Ei , i = 1, . . . , 4. N´ahodnou veliˇcinu, kter´a znaˇc´ı dobu ˇcek´an´ı na vstup do stavu Ej , j = i + 1, . . . , 5, oznaˇc´ıme Tij . Jelikoˇz v modelu pˇripouˇst´ıme pouze jednosmˇern´e pˇrechody, m˚ uˇzeme ps´at Tij = Ti + · · · + Tj−1 , kde Tj−1 = 0 pro j = i + 1. Vzhledem k nez´avislosti n´ahodn´ ych dob lze hustotu veliˇciny Tij zapsat pomoc´ı konvoluce gij (t) = fi ∗ · · · ∗ fj−1 (t), t > 0.
(3.6)
S odkazem na druhou kapitolu lze inkubaˇcn´ı dobu AIDS reprezentovat n´ahodnou veliˇcinou T14 a infekˇcn´ı dobu choroby HIV veliˇcinou T15 . D´ale se budeme zab´ yvat pouze rozdˇelen´ım inkubaˇcn´ı doby AIDS. Hustotu veliˇciny T14 vypoˇcteme aplikac´ı vztahu (3.6): g14 (t) = α1 α2 α3
3 X l=1
exp(−αl t) , t > 0. Q3 m=1 (αl − αm )
(3.7)
m6=l
Pˇr´ısluˇsnou distribuˇcn´ı funkci oznaˇc´ıme G14 (t). Plat´ı
G14 (t) =
Z
0
t
g14 (τ )dτ = α1 α2 α3
3 X l=1
1 − exp(−αl t) , t > 0. Q αl 3m=1 (αl − αm )
(3.8)
m6=l
Jelikoˇz je T14 = T1 + T2 + T3 a T1 , T2 , T3 jsou nez´avisl´e, plat´ı pro stˇredn´ı hodnotu a rozptyl vztah
19
E(T14 ) = var(T14 ) =
3 X
i=1 3 X
E(Ti ) =
3 X
αi−1 ,
i=1 3 X
var(Ti ) =
(3.9)
αi−2 .
i=1
i=1
Nyn´ı urˇc´ıme pravdˇepodobnosti pˇrechodu. V souladu se zaveden´ ym oznaˇcen´ım definujme pij (t) = P (X(t) = Ej |X(0) = Ei ), kde i, j ∈ {1, . . . , 5}. Zˇrejmˇe plat´ı pij (0) = δij . Budeme se tedy nad´ale zab´ yvat pˇr´ıpadem t > 0. Z vlastnost´ı modelu ihned plyne, ˇze pij (t) = 0 pro i > j. D´ale plat´ı pii (t) = = P (Ti > t) = Si (t) pro kaˇzd´e i = 1, . . . , 4 a p55 (t) = 1. Nyn´ı jeˇstˇe zb´ yv´a vyˇreˇsit pˇr´ıpad i < j. M˚ uˇzeme vyuˇz´ıt napˇr. vzorce (3.4) a (3.5). Jsou-li α1 , . . . , α4 navz´ajem r˚ uzn´a ˇc´ısla, potom dle tˇechto vzorc˚ u dostaneme
j−i
pij (t) = (−1)
αi · · · αj−1
pi5 (t) = (−1)4−i αi · · · α4
j X l=i
4 X l=i
exp(−αl t) , j 6= 5, Qj m=i (αl − αm )
(3.10)
m6=l
1 − exp(−αl t) . Q αl 4m=i (αl − αm )
(3.11)
m6=l
M˚ uˇzeme ale postupovat i jin´ ym zp˚ usobem, kter´ y nevyˇzaduje pˇredpoklad r˚ uznosti intenzit αi . Uvaˇzujme, ˇze proces v ˇcase 0 vstoupil do stavu Ei a v ˇcase t > 0 se nach´az´ı ve stavu Ej . Potom proces vstoupil do stavu Ej v nˇejak´em ˇcase u ∈ (0, t] s pravdˇepodobnost´ı gij (u)du a z˚ ustal v nˇem po dobu t − u ˇcasov´ ych jednotek s pravdˇepodobnost´ı Sj (t − u) pro j 6= 5 a s pravdˇepodobnost´ı 1 pro j = 5. Plat´ı tedy pij (t) =
Z
t
gij (u)Sj (t − u)du, j 6= 5,
0
pi5 (t) =
Z
t
gi5 (u)du.
0
Nezn´am´e parametry αi , i = 1, . . . , 4 m˚ uˇzeme odhadnout z dat metodou maxim´aln´ı vˇerohodnosti. Uvaˇzujme, ˇze m´ame k dispozici data poch´azej´ıc´ı od n sledovan´ ych jedinc˚ u. Necht’ u kaˇzd´eho jedince k, k = 1, . . . , n zn´ame stavy ykm ∈ S, ve kter´ ych se nach´azel v ˇcasech τkm , m = 0, 1, . . . , nk , 0 = τk0 < τk1 · · · < τknk . Oznaˇcme α = (α1 , α2 , α3 , α4 ). Jelikoˇz pracujeme s homogenn´ım Markovov´ ym procesem a osudy sledovan´ ych jedinc˚ u povaˇzujeme za nez´avisl´e, m´a vˇerohodnostn´ı funkce tvar
L(α) =
n Y
k=1
Lk (α), kde Lk (α) =
nk Y
m=1
20
pyk(m−1) ykm (τkm − τk(m−1) ),
1
0.8
0.6
0.4
0.2
0
50
100
150
200
250
ˆ 14 (t). Na vodorovn´e ose je vynesen ˇcas v mˇes´ıc´ıch. Obr´azek 3.1: Graf funkce G
kde Lk (α) znaˇc´ı vˇerohodnostn´ı funkci pˇr´ısluˇsej´ıc´ı jedinci k a pij (t) jsou d´any vztahy (3.10) a (3.11). Pro u ´plnost uved’me hodnoty maxim´alnˇe vˇerohodn´ ych odhad˚ u αˆi intenzit αi , kter´e spoˇcetli Longini a kol. (v mˇes´ıc−1 ): αˆ1 = 0, 4571, αˆ2 = 0, 0190, αˆ3 = 0, 0159, αˆ4 = 0, 0424. Dosad´ıme-li α ˆ i za αi v (3.9), dostaneme odhad stˇredn´ı d´elky inkubaˇcn´ı doby AIDS v mˇes´ıc´ıch:
\ E(T 14 ) =
3 X i=1
\ E(T i) =
3 X
α ˆ i = 2, 1877 + 52, 6316 + 62, 8931 = 117, 7124,
i=1
coˇz je pˇribliˇznˇe 9, 8 roku. Odhad hustoty gˆ14 (t) a odhad distribuˇcn´ı funkce ˆ 14 (t) inkubaˇcn´ı doby AIDS z´ısk´ame dosazen´ım α G ˆ i za αi v (3.7) a (3.8). Na ˆ obr´azku 3.1 je zn´azornˇen graf funkce G14 (t). Graf byl sestrojen pomoc´ı programu MAPLE. ♦ St´adia pouˇzit´a v pˇredeˇsl´em modelu n´am umoˇzn ˇovala bez jak´ ychkoli omezen´ı pˇredpokl´adat jednosmˇern´e pˇrechody mezi stavy. V n´asleduj´ıc´ım modelu to uˇz platit nebude. Pˇ r´ıklad 3.2. Stejnˇe jako v pˇredchoz´ım pˇr´ıkladu modelujeme pr˚ ubˇeh choroby HIV homogenn´ım Markovov´ ym procesem se spojit´ ym ˇcasem a koneˇcnou mnoˇzinou stav˚ u. Pro mnoˇzinu stav˚ u S necht’ plat´ı S = {E1 , . . . , E7 }, 21
kde E1 , . . . , E6 jsou pˇrechodn´e stavy definovan´e v tabulce 3.2 a E7 absorpˇcn´ı stav reprezentuj´ıc´ı klinick´e st´adium AIDS. Pˇri pouˇzit´ı takto definovan´ ych stav˚ u je tˇreba m´ıt na pamˇeti, ˇze zjiˇstˇen´ y poˇcet T4 bunˇek m˚ uˇze u jedince znaˇcnˇe kol´ısat, a to jak v d˚ usledku chyb mˇeˇren´ı, tak v d˚ usledku bˇeˇzn´ ych fyziologick´ ych zmˇen organizmu. V modelu proto pˇripust´ıme pˇrechody v obou smˇerech. Nicm´enˇe, d´elka interval˚ u definovan´ ych v tabulce 3.2 je dostateˇcnˇe velk´a, aby se vliv takov´eto fluktuace eliminoval. Pˇredpokl´adejme tedy, ˇze pˇrechody mezi pˇrechodn´ ymi stavy se ˇr´ıd´ı n´asleduj´ıc´ım sch´ematem: Ei → Ej , kde j = i − 1, i + 1 pro i = 2, . . . , 5 a E1 → E2 , E6 → E5 . D´ale je tˇreba vz´ıt v u ´vahu, ˇze seznam chorob definuj´ıc´ıch AIDS byl v pr˚ ubˇehu ˇcasu ofici´alnˇe zmˇenˇen. Proto pˇripust´ıme tak´e pˇr´ım´e pˇrechody typu Ei → E7 , i > 1. Vzhledem k d´elce inkubaˇcn´ı doby AIDS je vˇsak nepravdˇepodobn´e, ˇze by se tato nemoc rozvinula pˇr´ımo z ˇcasn´eho st´adia. V modelu m˚ uˇzeme napˇr´ıklad pˇredpokl´adat, ˇze pˇrechody Ei → E7 nast´avaj´ı pouze pro i = 1, 2, 3. Intenzity pˇrechodu ze stavu Ei do stavu Ej oznaˇc´ıme qij , i, j ∈ {1, . . . , 7}. Matice intenzit Q = {qij } necht’ m´a tvar −q1 α1 0 0 0 0 0 β2 −q2 α2 0 0 0 0 0 β3 −q3 α3 0 0 0 , 0 0 β −q α 0 γ Q= 4 4 4 4 0 0 0 β5 −q5 α5 γ5 0 0 0 0 β6 −q6 γ6 0 0 0 0 0 0 0 P7 kde α6 = γ6 , αi > 0, i = 1, . . . , 6 a qi = j=1 qij . j6=i
Opˇet n´as zaj´ım´a rozdˇelen´ı inkubaˇcn´ı doby AIDS. Je-li jedinec infikov´an virem HIV v ˇcase 0, lze inkubaˇcn´ı dobu AIDS reprezentovat veliˇcinou T17 . Pˇri urˇcov´an´ı jej´ıho rozdˇelen´ı jiˇz nelze postupovat jako v pˇredchoz´ım pˇr´ıkladu. Tomuto probl´emu je vˇenov´ana n´asleduj´ıc´ı teorie. ♦
Model z pˇr´ıkladu 3.2 zobecn´ıme. Budeme uvaˇzovat homogenn´ı Markov˚ uv proces {X(t), t ≥ 0} se spojit´ ym ˇcasem a koneˇcnou mnoˇzinou stav˚ u S = = {E1 , E2 , . . . , Ek , Ek+1 }, k ∈ N. Stavy E1 , . . . , Ek jsou pˇrechodn´e, stav E1 reprezentuje poˇc´atek infekce a stavy E2 , . . . , Ek znaˇc´ı nˇejak´a infekˇcn´ı st´adia. Stav Ek+1 je absorpˇcn´ı stav reprezentuj´ıc´ı klinick´e st´adium AIDS. V modelu budeme pˇredpokl´adat pˇrechody Ei → Ei+1 , i = 1, . . . , k s intenzitami αi > 0 a obecnˇe pˇripust´ıme i pˇrechody Ei → Ei−1 , i = 2, . . . , k s intenzitami βi ≥ 0 a pˇrechody Ei → Ek+1 , i = 1, . . . , k s intenzitami γi , kde γk = αk . Model m˚ uˇzeme shrnout sch´ematem na obr´azku 3.2. Pˇredpokl´adejme, ˇze proces je v ˇcase 0 ve stavu Ei0 , i0 = 1, . . . , k a oznaˇcme Ti0 (k+1) n´ahodn´ y ˇcas, kdy je stav Ei0 absorbov´an stavem Ek+1 . Hustotu veliˇciny Ti0 (k+1) oznaˇc´ıme gi0 (t) a distribuˇcn´ı funkci Gi0 (t). Inkubaˇcn´ı doba AIDS je tedy reprezentov´ana rozdˇelen´ım s hustotou g1 (t). Rozdˇelen´ı veliˇciny Ti0 (k+1) urˇcuje n´asleduj´ıc´ı lemma. 22
α
αk−2
αk−1
→ −−−−→ −−− −−− −−− −−− −−−− −−− −−− −− −− −− −− −− −−−→ −−−−−−→
−−−−−−−→
α
−−−−−−−→
1 2 −−−−− −−→ −−−−− −−→ · · · · · · −−−−−−−→ Ek−1 −−−−−−−→ Ek E1 − −−− β2 E2 −− β3 βk−1 βk −−− − γ2 −− γ γ1 −−− γk k−1 −−− −− −−− −− −−− −→ −→ Ek+1
Obr´azek 3.2: Sch´ema zobecnˇen´eho modelu
Lemma 3.2. Necht’ plat´ı v´yˇse zaveden´e znaˇcen´ı. Potom pro kaˇzd´e i0 = = 1, . . . , k je gi0 (t) =
k k−1 X X
ciij0 tj exp(−ρi t), t > 0,
j=0 i=1
kde ciij0 , ρi jsou konstanty.
Neˇz nap´ıˇseme d˚ ukaz, uvedeme nˇekolik pomocn´ ych tvrzen´ı. Pro n´asleduj´ıc´ı text definujme β1 = 0. Lemma 3.3. Uvaˇzujme ˇctvercovou matici Ak ˇr´adu k ve tvaru λ1 −α1 0 ... ... ... 0 −β2 λ2 −α2 0 ... ... 0 0 −β3 λ3 −α3 0 . . . 0 Ak = . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . , 0 ... ... 0 −βk−1 λk−1 −αk−1 0 ... ... ... 0 −βk λk
αi > 0, βi ≥ 0, λi ≥ αi + βi pro i = 1, . . . , k. Potom det(Ak ) > 0. D˚ ukaz. Viz vˇeta 2.2 v [8, str. 16].
Lemma 3.4. (viz [5, str. 301]) Necht’ jsou d´any ˇctvercov´e matice A, B ˇr´adu k. Potom plat´ı exp(A + B) = exp(A) exp(B), je-li AB = BA, (3.12) d exp(At) = A exp(At), (3.13) dt V−1 exp(A)V = exp(V−1 AV), kde V je libovoln´a regul´ arn´ı matice ˇr´adu k. (3.14) 23
Uvaˇzujme matici A = Ak , kde Ak je definov´ana v lemmatu 3.3 a pro kterou plat´ı λi = αi + βi + γi , i = 1, . . . , k − 1 a λk = αk + βk . Konstanty αi , βi , γi necht’ maj´ı v´ yznam intenzit v´ yˇse definovan´eho procesu {X(t), t ≥ 0} (viz obr´azek 3.2). Oznaˇcme 0k nulov´ y sloupcov´ y vektor o k prvc´ıch a γ = T = (γ1 , . . . , γk ) . Potom matice −A γ Q= 0T 0 k
je matice intenzit procesu {X(t), t ≥ 0}.
Lemma 3.5. Uvaˇzujme matice A, Q definovan´e v´yˇse. Je-li A invertibiln´ı, potom plat´ı g(t) = exp(−At)γ, t > 0,
(3.15)
kde g(t) = (g1 (t), . . . , gk (t))T . D˚ ukaz. Zˇrejmˇe plat´ı j
j
Q = (−1) Dle vˇety 3.4 je
Aj −Aj−1 γ , j = 1, 2, . . . . 0T 0 k
∞ ∞ X X 1 j j (−1)j j Aj −Aj−1 γ P(t) = exp(Qt) = t Q = Ik+1 + t = 0T 0 j! j! k j=0 j=1 P∞ (−1)j j j ! P (−1)j j j −1 t A −A Ik + ∞ j=1 j! t A γ j=1 j! = = T 0k 1 exp(−At) A−1 [Ik − exp(−At)]γ = , t > 0, (3.16) 0T 1 k kde Ik+1 znaˇc´ı jednotkovou matici ˇr´adu k + 1, obdobnˇe pro Ik . Oznaˇcme pk+1 (t) = (p1(k+1) (t), . . . , pk(k+1) (t))T . Pro kaˇzd´e i = 1, . . . , k plat´ı pi(k+1) (t) = P (X(t) = k + 1, Ti(k+1) > t | X(0) = i) + P (X(t) = k+ + 1, Ti(k+1) ≤ t | X(0) = i) = P (X(t) = k + 1, Ti(k+1) ≤ t | X(0) = i) = = P (Ti(k+1) ≤ t) = Gi (t). Je tedy pk+1 (t) = (G1 (t), . . . , Gk (t))T , t > 0. S vyuˇzit´ım (3.16) a (3.13) m˚ uˇzeme ps´at g(t) =
d d pk+1 (t) = {A−1 [Ik − exp(−At)]γ} = exp(−At)γ, t > 0. dt dt
24
Jeˇstˇe pˇripomeneme pojem Jordanovy matice. Problematice Jordanov´ ych matic vˇcetnˇe metod jejich nalezen´ı se vˇenuje napˇr. [1, str. 235-264]. Definice 3.5. Necht’ jsou d´any ˇctvercov´e matice A,B stejn´eho ˇr´adu defiˇ novan´e na t´emˇze tˇelese F . Rekneme, ˇze matice A,B jsou podobn´e, jestliˇze existuje regul´arn´ı matice V nad tˇelesem F takov´a, ˇze A = V−1 BV.
ˇ Definice 3.6. Ctvercov´ a matice tvaru a 1 0 ... ... ... 0 0 a 1 0 . . . . . . 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . J(a) = 0 . . . . . . . . . 0 a 1 0 ... ... ... ... 0 a
se naz´ yv´a Jordanova buˇ nka pˇr´ısluˇsn´a prvku a. Diagon´aln´ı blokov´a matice J, jej´ıˇz bloky na diagon´ale jsou Jordanovy buˇ nky, se naz´ yv´a Jordanova matice. Jordanovu matici m˚ uˇzeme definovat ekvivalentn´ım zp˚ usobem jako matici a1 d1 0 . . . . . . . . . 0 0 a2 d2 0 . . . . . . 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . , (3.17) 0 . . . . . . . . . 0 ak−1 dk−1 0 ... ... ... ... 0 ak
kde pro kaˇzd´e i = 1, . . . , k − 1 je di = 0 nebo di = 1 a, je-li di = 1, potom ai = ai+1 .
Vˇ eta 3.6. Necht’ A je ˇctvercov´ a matice ˇr´adu k s vlastn´ımi ˇc´ısly ρ1 , . . . , ρk . Potom matice A je podobn´ a Jordanovˇe matici. D˚ ukaz. Viz vˇeta 17.8 v [2, str. 116-118]. Vlastn´ı ˇc´ısla ρ1 , . . . , ρk matice A nemus´ı b´ yt nutnˇe r˚ uzn´a. Necht’ ρil , l = 1, . . . , n jsou vˇsechna navz´ajem r˚ uzn´a vlastn´ı ˇc´ısla, jejich n´asobnosti oznaˇc´ıme kil . Potom m´a Jordanova matice Jordanovy buˇ nky J(ρil ) ˇr´adu kil , kter´e budeme d´ale znaˇcit Jl , l = 1, . . . , n . Jordanova matice nen´ı ve vˇetˇe 3.6 d´ana jednoznaˇcnˇe. Plat´ı totiˇz, ˇze dvˇe Jordanovy matice liˇs´ıc´ı se pouze poˇrad´ım Jordanov´ ych bunˇek jsou podobn´e (viz D˚ usledky 18.12 (i) v [1, str. 244]).
25
D˚ ukaz. (Lemma 3.2) Necht’ matice A je definov´ana jako v lemmatu 3.5. Dle lemmatu 3.3 je matice A regul´arn´ı a dle lemmatu 3.5 tedy plat´ı (3.15). Vlastn´ı ˇc´ısla matice A oznaˇcme ρ1 , . . . , ρk . Dle vˇety 3.6 existuje Jordanova matice J a regul´arn´ı matice V tak, ˇze A = V−1 JV. Necht’ plat´ı znaˇcen´ı zaveden´e pˇred t´ımto d˚ ukazem. Uvaˇzujme matici J, kter´a se skl´ad´a z Jordanov´ ych bunˇek J1 , . . . , Jn (v tomto poˇrad´ı). Pro kaˇzd´e l = 1, . . . , n oznaˇcme Ml diagon´aln´ı matici ˇr´adu kil s diagon´aln´ım prvkem ρil a definujme matici Nl = Jl −Ml . Odpov´ıdaj´ıc´ı diagon´aln´ı blokovou matici s diagon´aln´ımi bloky M1 , . . . , Mn (v tomto poˇrad´ı) oznaˇc´ıme M, analogicky z blok˚ u N1 , . . . , Nn zkonstruujeme matici N. Potom plat´ı J = M + N. Jelikoˇz pro kaˇzd´e l je Ml Nl = Nl Ml , je zˇrejmˇe tak´e MN = NM. Oznaˇc´ıme-li M = V−1 MV a N = V−1 NV, plat´ı A = V−1 JV = V−1 (M + N)V = M + N, MN = V
−1
MNV = V
−1
(3.18)
NMV = NM.
(3.19)
(−1)
Oznaˇc´ıme V = {vij } a V−1 = {vij }. Pro s = 1, 2, . . . plat´ı Ms = (s) (s) (s) = {mij }, kde mii = ρsi a mij = 0 jinak, a tedy plat´ı exp(−Mt) = {µij }, (s) kde µii = exp(−ρi t) a µij = 0 jinak. Pro Ns = {nij } je Qi+s−1 aroveˇ n j ≤ k, (s) m=i dm je-li j = i + s a z´ nij = 0 jinak,
kde dm jsou prvky matice J ve smyslu (3.17). Plat´ı Ns = 0, s ≥ k a pro exp(−Nt) = {νij } je (−1)p p (p) n j = i + p, p! t nij je-li j > i a z´aroveˇ νij = 1 je-li j = i, 0 je-li j < i.
Oznaˇcme pro kaˇzd´e i0 = 1, . . . , k jednotkov´ y sloupcov´ y vektor s ˇc´ıslem 1 na i0 -t´em m´ıstˇe jako ei0 . Nyn´ı vyuˇzijeme (3.15) a, jelikoˇz plat´ı (3.18) a (3.19), tak´e (3.12) a (3.14). Pro t > 0 plat´ı
T T −1 gi0 (t) = eT exp(−Mt) exp(−Nt)V(γ1 , . . . , γk )T i0 g(t) = ei0 exp(−At)γ = ei0 V !T k k X X (−1) (−1) = vi0 1 exp(−ρ1 t), . . . , vi0 k exp(−ρk t) exp(−Nt) vki γi v1i γi , . . . , i=1
i=1
2 X
k X
2−r
(−1)k−r k−r (−1) (−1) (2−r) t2−r exp(−ρr t)vi0 r nr2 , . . . , t · (2 − r)! (k − r)! r=1 r=1 !T k k k k−1 X X XX (−1) (k−r) vki γi ciij0 tj exp(−ρi t). · exp(−ρr t)vi0 r nrk ) v1i γi , . . . , = (−1)
= (exp(−ρ1 t)vi0 1 ,
i=1
i=1
26
j=0 i=1
Pro konstanty ciij0 , kde i, j, i0 = 1, . . . , k plat´ı
ciij0 = kde
Qj
(
m=i
P (−1)j (−1) Qi+j−1 k je-li i + j ≤ k, v p=1 γp v(i+j)p m=i dm j! i0 i 0 jinak,
dm = 1 pro j < i.
D˚ ukaz n´am d´av´a n´avod, jak spoˇc´ıst pˇr´ısluˇsn´e konstanty z lemmatu 3.2. Nicm´enˇe, uvaˇzujeme-li model s obecnˇe kladn´ ymi intenzitami αi , βi a γi , je i pro malou mnoˇzinu stav˚ u S takov´ yto v´ ypoˇcet numericky n´aroˇcn´ y.
3.3
Vhodnost Markovov´ ych model˚ u
Pro pouˇzit´ı Markovov´ ych model˚ u pˇri modelov´an´ı pr˚ ubˇehu choroby HIV hovoˇr´ı nˇekolik skuteˇcnost´ı. Pˇri modelov´an´ı je pˇrirozen´e definovat jednotliv´a st´adia choroby, ta ztotoˇznit se stavy nˇejak´eho stavov´eho modelu a probl´em redukovat na modelov´an´ı pˇrechod˚ u mezi stavy tohoto modelu. D´ale je nutno si uvˇedomit, ˇze data, se kter´ ymi pracujeme, jsou cenzorov´ana, a to jak intervalovˇe, tak i zleva a zprava. V´ıme totiˇz pouze, v jak´ ych st´adi´ıch se pacienti nach´azeli ve vybran´ ych ˇcasov´ ych okamˇzic´ıch, ale nezn´ame pˇresn´e ˇcasy pˇrechod˚ u mezi st´adii choroby (interval censoring). U pacient˚ u obvykle nezn´ame ani pˇresn´ y ˇcas infikace (left censoring). D´ale mus´ıme poˇc´ıtat s t´ım, ˇze vzhledem k dlouh´e inkubaˇcn´ı dobˇe se u nˇekter´ ych pacient˚ u bˇehem pozorov´an´ı nerozvine klinick´e st´adium AIDS (right censoring). Z tˇechto d˚ uvod˚ u nelze pˇri anal´ yze dat pouˇz´ıt bˇeˇzn´e statistick´e metody. Pˇrirozen´ y charakter Markovov´ ych proces˚ u vˇsak umoˇzn ˇuje dobˇre odhadovat nezn´am´e parametry i z dat, kter´a jsou silnˇe zat´ıˇzena cenzorov´an´ım. Nicm´enˇe, jak je zn´amo z teorie, n´ahodn´e doby mezi pˇrechody maj´ı v Markovov´ ych modelech exponenci´aln´ı rozdˇelen´ı, tedy rozdˇelen´ı bez pamˇeti. Tento fakt je zˇrejmˇe pˇri konstrukci model˚ u omezuj´ıc´ı. Probl´em lze ˇreˇsit pouˇzit´ım semi-Markovov´ ych proces˚ u, jejichˇz doby mezi pˇrechody jsou tak´e nez´avisl´e, ale mohou m´ıt libovoln´e spojit´e rozdˇelen´ı. Teorii semi-Markovov´ ych proces˚ u se vˇsak v t´eto pr´aci vˇenovat nebudeme. Pˇrehledn´ yu ´vod do t´eto problematiky lze nal´ezt napˇr. v [7, str. 86–95].
27
Literatura [1] Beˇcv´aˇr J.: Line´arn´ı algebra, MATFYZPRESS, Praha, 2002. [2] Bican L.: Line´arn´ı algebra a geometrie, Academia, Praha, 2002. [3] Chiang C. L.: Introduction to Stochastic Processes in Biostatistics, John Wiley & Sons, Inc., New York, 1968. [4] Chung K. L.: Markov Chains with Stationary Transition Probabilities, Springer-Verlag, New York, 1967. [5] Curtis C. W.: Linear Algebra: An Introductory Approach, SpringerVerlag, New York, 1993. [6] Longini I. M. a kol.: Statistical Analysis of the Stages of HIV Infection Using a Markov Model, Statistics in Medicine 8 (1989) 831–843. [7] Mode C. J. , Sleeman C. K.: Stochastic Processes in Epidemiology, World Scientific, Singapore, 2000. [8] Neammanee K. a kol.: On the AIDS Incubation Distribution, Statistics & Probability Letters 73 (2005) 13–23. [9] Pr´aˇskov´a Z., Lachout P.: Z´ aklady n´ ahodn´ych proces˚ u, Nakladatelstv´ı Karolinum, Praha, 2005. [10] HIV infekce na internetov´e adrese http://www.lf3.cuni.cz/studium/materialy/infekce/HIV infekce.doc
28
