Kui on leitud selle mudeli parameetrite hinnangud βˆj (j = 1, . . . , p) ja Nelson-Aaleni hinnang kumulatiivsele baasriskifunktsioonile, siis saab arvutada vastavad martingaalijäägid: mˆ i = δi − Hˆ0(ti) exp(βˆ1z1i + . . . + βˆpzpi). Therneau et al. (1990) on näidanud, et selliselt arvutatud martingaalijäägid on seotud mudelist välja jäänud argumenttunnusega (antud juhul on selleks genotüüp g). Kui genotüüp g tuleks sarnaselt mudelis olevate tunnustega lisada mudeli lineaarsesse osasse ilma teisenduseta, siis on ka tema seos martingaalijäägiga lineaarne, ehk mˆ i c= β0 + β˜ggi + εi, kus c = d/n on sündmuste arvu ja kõigi vaatluste arvu suhe, β0 on lineaarse mudeli vabaliige, β˜g ≈βg ja εi ∼N(0, σ2) on sõltumatud vead. Seega on selle asemel, et iga huvipakkuva SNP-ga hinnata Coxi mudel, mis on ülegenoomsete seoseuuringute puhul väga aja- ja arvutus mahukas, võimalik hinnata üks Coxi nullmudel, kus ühegi SNP mõju ei ole arvesse võetud, ja siis hinnata iga SNP-ga lineaarne mudel, kus uuritavaks tunnuseks on nullmudeli põhjal arvu tatud skaleeritud martingaalijääk. Teoreetiliselt on teada, et selline lähendamine töötab ainult kindlates piirides. Selles töös uuri me, millal ja kui hästi sobib martingaalijääkide meetodi rakendamine Coxi mudelite hindamise asemel. 22 2.2.2 Weibulli jaotusega võrdeliste riskide mudel

Coxi võrdeliste riskide mudel on poolparameetriline mudel, sest baasriskifunktsiooni kuju ei määrata. Seetõttu ei saa Coxi mudelit kasutada elukestusandmete simuleerimiseks. Seega, kuigi elukestusandmete analüüsimiseks kasutatakse enamasti Coxi mudelit, siis elukestusandmete simuleerimiseks kasutatakse tavaliselt parameetrilisi mudeleid. Parameetriliste võrdeliste riskide mudeli puhul, kus baasriskifunktsiooni kuju on määratud, kasutatakse kõige enam Weibulli jaotust (Kleinbaum ja Klein, 2012: 304). Kirjeldame siin Weibulli jaotuse parametrisatsiooni, mida kasutatakse R-i funktsioonides dweibull, pweibull, rweibull ja qweibull. Weibulli jaotusel on kaks parameetrit: kujuparameeter a > 0 ja skaala parameeter b > 0. Olgu juhuslik suurus T Weibulli jaotusega T ∼ W(a, b). Sel juhul on T tihedusfunktsiooni ja jaotusfunktsiooni kujud t > 0 korral vastavalt f(t) = a b t b a−1 exp − t b a ja F(t) = 1 −exp Weibulli jaotuse riskifunktsioo n on kujul − t b a . 1 − F(t)= a b t b a−1= ab−ata−1, t > 0. h(t) = f(t)

Seega on a > 1 korral tegemist kasvava riskiga, a < 1 korral kahaneva riskiga ja a = 1 korral 23 on risk konstantne. Weibulli jaotuse kumulatiivne riskifunktsioon on Z t H(t) = = 0 Z t 0 h(u)du = ab−aua−1d u = t b a , t > 0. (4) Joonisel 2 on näited Weibulli jaotuse riskifunktsiooni kujude kohta erinevate kujuparameetrite a korral. ) t ( h n o o i s t k n u f i k s i R a=5 2,0 a=2 1,5 1,0

a=1,5 0,5 a=1 a=0,5 0,0 0 1 2 3 t Joonis 2. Weibulli jaotuse W(a, b = 2) riskifunktsioon h(t) erinevate kujuparameetrite a korral Tänu lihtsale riski- ja üleelamisfunktsioonile ning erinevatele võimalikele riskifunktsiooni kuju dele on Weibulli jaotus parameetrilises elukestusanalüüsis laialdaselt kasutusel. Vaatame nüüd võrdeliste riskide mudelit juhul, kui elukestused T on Weibulli jaotusega. Kui baaselukestus on Weibulli jaotusega W(a, b), siis i-nda indiviidi, kelle argumenttunnuste vektor

24 on xi, riskifunktsioon avaldub kujul hi(t) = h0(t) exp(βTxi) = = ab−ata−1exp(βTxi) = b exp(βTxi) − 1 a −ata−1 = ab−a ita−1, = a kus bi = b exp(βTxi) − 1a. (5) Seega on i-nda indiviidi elukestus Weibulli jaotusega W(a, bi). Weibulli jaotusega võrdeliste riskide mudelite puhul eeldataksegi, et kujuparameeter a on igal indiviidil samasugune, aga skaalaparameeter bi sõltub baasparameetritest a ja b, parameetrite vektorist β ning argument tunnuste vektorist xi. 2.3 Binaarse tunnuse modelleerimine Jälgimiseelsete haigusjuhtude uurimisel vaadeldakse haigusseisundit binaarse tunnusena, mille väärtuseks on kas 1 või 0 vastavalt diagnoosi esinemisele või puudumisele. Kirjeldame siin kahte binaarse tunnuse modelleerimise meetodit (nimetame neid edaspidi lihtsuse mõttes binaarseteks mudeliteks) – logistilist regressiooni ja täiend-log-log regressiooni – ning nende mudelite efektisuuruste seost võrdeliste riskide mudelite efektisuurusega ehk riskimäärade suhtega. Olgu meil valimis n vaatlust ja p seletavat tunnust. Olgu Yi (i = 1, . . . , n) binaarne juhuslik suurus, mis kirjeldab meile huvipakkuva sündmuse toimumist ehk diagnoosi

esinemist i-ndal indiviidil, ja olgu pi = P(Yi = 1) selle sündmuse toimumise tõenäosus. 25 2.3.1 Logistiline mudel Logistilise regressiooni korral kasutatakse uuritava sündmuse tõenäosuse kirjeldamiseks sõltu matute argumenttunnuste lineaarkombinatsiooni kaudu logit-seosefunktsiooni: logit(pi) = log pi 1 − pi = β0 + βTxi, kus β0 ja β = (β1, . . . , βp)T on mudeli tundmatud parameetrid, xi = (x1i, . . . , xpi)T on i-nda indiviidi argumenttunnuste vektor ja pi/(1 −pi) on sündmuse toimumise šanss. Sündmuse toi mumise šanss näitab, mitu korda on uuritava sündmuse toimumine tõenäolisem kui sündmuse mittetoimumine. Tundmatud parameetrid βj (j = 0, . . . , p) hinnatakse tavaliselt suurima tõepära meetodil ja hinnangute arvutamiseks kasutatakse iteratiivseid meetodeid, näiteks R-i funktsioonis glm on kasutusel Fisheri skoorimeetod. Logistilise regressioonimudeli parameetreid interpreteeritakse tavaliselt šansside suhte (OR, odds ratio) kaudu. Erinegu indiviidid a ja b vaid j-nda (j = 1, . . . , p) argumenttunnuse väär tuse poolest: xja = xjb + 1. Siis on nende indiviidide šansside suhte ja j-nda argumenttunnuse parameetri βj seos järgmine: pa/(1 − pa) pb/(1 − pb)=exp(β0 + β1x1a + . . . + βj (xjb + 1) + . . . + βpxpa) exp(β0 + β1x1a + . . . + βjxjb + . . . + βpxpa)= exp(βj ). See tähendab, et kui j-nda argumenttunnuse väärtus muutub ühe ühiku võrra (ja

teised väärtused jäävad samaks), siis šansid muutuvad exp(βj ) korda. Šansside suhe OR ja riskimäärade suhe HR Eespool nägime, et võrdeliste riskide mudeli puhul saab mudeli parameetreid interpreteerida kui logaritme riskimäärade suhetest, logistilise mudeli puhul aga kui logaritme šansside suhe test. Riskimäärade suhe ja šansside suhe on kaks sagedasti raporteeritavat suurust, kuid nende 26 interpreteerimisel ollakse sageli hooletu, kusjuures mõlemaid kasutatakse tihti hoopis kolmanda näitaja – suhtelise riski (relative risk, risk ratio) – iseloomustamiseks. Teatud juhtudel – näiteks, kui huvipakkuva sündmuse toimumise tõenäosus on väike või jälgimisaeg lühike – ongi need näitajad sarnased. Samuti on kõigi nende kolme näitaja suund alati sama: kui teame, et üks neist on suurem kui 1, on seda ka teised kaks, ja vastupidi. (Davies et al., 1998; Bangdiwala, 2010; Sutradhar ja Austin, 2018) Vaatame lähemalt, kuidas on omavahel seotud šansside suhe (OR) ja riskimäärade suhe (HR). Olgu meil lihtsuse mõttes võrdluses kaks gruppi: kontrollgrupp, kus argumenttunnus x = 0 ja riskigrupp, kus x = 1, riskifunktsioonidega vastavalt h0 ja h1. Olgu β tunnusele x vastav Coxi mudeli regressioonikordaja ja olgu jälgimisaja pikkuseks t > 0. Võrdeliste riskide mudeli kehtimise korral avaldub riskigrupi riskifunktsioon h1 kontrollgrupi riskifunktsiooni h0 kaudu: h1(t) = h0(t) exp(βx) = h0(t) exp(β) ja kahe grupi riskimäärade suhe on h0(t)=h0(t) exp(β) HR = h1(t) h0(t)= exp(β).

Tähistame nüüd huvipakkuva sündmuse toimumise tõenäosuse ajavahemikus algmomendist hetkeni tkontrollgrupisp0 ja riskigrupisp1 . Kumulatiivse baasriskifunktsiooni H0 kaudu saame need tõenäosusedkirja panna järgmiselt: p0 = 1− S0 (t) = 1−e −H0 (t) , p1 = 1− S1 (t) = 1−e −H1 (t) = 1−e −H0 (t) exp(β) . Kahe grupi šansside suhte saame nende tõenäosuste kaudu arvutada järgmiselt: 27 OR = p1 /(1− p1 ) p0 /(1− p0 )= p1 (1− p0 ) p0 (1− p1 )

= (1−e −H0 (t) exp(β) ) e −H0 (t) (1−e −H0 (t) ) e −H0 (t) exp(β)

= 1−e −H0 (t) exp(β) (e H0 (t) −1) e −H0 (t) exp(β)

= e H0 (t) exp(β) −1 e H0 (t) −1

= (1− p0 ) −exp(β) −1 (1− p0 ) −1 −1= (1− p0 ) −HR −1 (1− p0 ) −1 −1 , seega avaldub šansside suhe baasgrupi levimusep0 ja riskimäärade suhte HR =exp(β) kaudu ning kumulatiivset baasriskifunktsiooni ei ole vaja teada (Green jaSymons, 1983). Seda seost on illustreeritud joonisel 3. Sellel on näidatud, kuidas sõltub šansside suhe OR kontrollgrupi levimusest p0 erinevate riskimäärade suhete HR korral. Paneme tähele, et  OR > HR, kui HR > 1 ehk β > 0, OR < HR, kui HR < 1 ehk β < 0,  OR = HR, kui HR = 1 ehk β = 0. Samuti näeme, et HR ja OR väärtused on sarnased, kui HR on ühe lähedal ning väärtused on seda erinevamad, mida suurem on p0. Et exp(x) ≈1 + x, kui x on absoluutväärtuselt väike, siis saame šansside suhet ka lähendada: OR = eH 0 (t)HR −1 eH 0 (t) −1≈ ≈1 + H0(t)HR −1 1 + H0(t) −1= HR, juhul kui vaatlusaeg t on lühike ja seega kumulatiivne riskifunktsioon H0(t) väike, ning sünd muse esinemise tõenäosus vaatlusaja jooksul ega riskifaktori mõju ei ole suured (Green ja Symons, 1983). 28 R O 5 4 3 2 1 0

0,00 0,25 0,50 0,75 1,00 p0 HR 0,5 0,9 1 1,05 1,5 2 3 Joonis 3. Šansside suhe OR erinevate riskimäärade suhete HR ja baasgrupi levimuste p0 korral 2.3.2 Täiend-log-log mudel Analoogiliselt logistilise regressiooniga kasutatakse ka täiend-log-log mudelit, kirjeldamaks uuritava sündmuse toimumise tõenäosust ja selle muutumist sõltuvalt argumenttunnuste väär tuste muutumisest. Selle mudeli puhul rakendatakse täiend-log-log (cloglog, complementary log-log) seosefunktsiooni: cloglog(pi) = log(−log(1 − pi)) = β0 + βTxi, kus β0 ja β = (β1, . . . , βp)T on mudeli tundmatud parameetrid, xi = (x1i, . . . , xpi)T on i-nda indiviidi argumenttunnuste vektor, ja millest sündmuse toimumise tõenäosus pi avaldub kujul pi = 1 −exp(−exp(β0 + βTxi)). 29

Täiend-log-log mudeli saab R-is hinnata funktsiooniga glm, kus argumendina tuleb kasutada family = binomial(link = "cloglog"). Täiend-log-log mudeli parameetrid hinnatakse enamasti suurima tõepära meetodil, funktsioonis glm kasutatakse selleks Fisheri skoorimeetodit. Vaatame, kuidas on seotud Coxi võrdeliste riskide elukestusmudel täiend-log-log mudeliga. Olgu meil vaatluse all n indiviidi ja vastaku haiguse diagnoosimise vanused Ti (i = 1, . . . , n) Coxi võrdeliste riskide mudelile ehk olgu haigestumise vanusele Ti vastav kumulatiivne riski funktsioon Hi(t) = H0(t) exp(βTxi), kus β on võrdeliste riskide mudeli parameetrite vektor, xi on i-nda indiviidi argumenttunnuste vektor ja H0 on kumulatiivne baasriskifunktsioon, mille kuju jäetakse Coxi mudeli puhul mää ramata. Olgu jälgimisaja pikkuseks i-ndal indiviidil ti > 0 (see võib olla näiteks liitumisvanus). Hai gestumise indikaator on siis Yi = I(Ti 6 ti). Kasutades seost (1), saame tõenäosuse, et i-ndal indiviidil on haigus diagnoositud ajaks ti, kirjutada kujul pi = P(Yi = 1| xi, ti) = P(Ti 6 ti|xi) = = 1 − Si(ti) = = 1 −exp(−Hi(ti)) = = 1 −exp(−H0(ti) exp(βTxi)) = = 1 −exp(−exp(log(H0(ti)) + βTxi)). Juhul kui H0(ti) = exp(β0 + βtti), saame pi = 1 −exp(−exp(β0 + βtti + βTxi)). See vastab täpselt täiend-log-log mudelile vabaliikmega β0 ja parameetrite vektoriga

(βt; β) 30 (jälgimisaeg ti on lisatud argumenttunnusena). Seega saab ka täiend-log-log mudeli para meetreid interpreteerida kui riskimäärade suhteid. See kehtib aga vaid siis, kui on täidetud eel dus, et kumulatiivne risk avaldub kui H0(ti) = exp(β0 + βtti). Viimane eeldus ei pruugi aga alati täidetud olla. Meid huvitaval juhul, kus haigestumine sõltub genotüübist, avaldub i-nda indiviidi haigestumise tõenäosus võrdeliste riskide mudeli kehtimise korral järgmiselt: P(Ti 6 ti| gi) = 1 − Si(ti) = = 1 −exp(−H0(ti) exp(βggi)), kus gi ∈ {0, 1, 2} on i-nda indiviidi genotüüp, H0 on kumulatiivne baasriskifunktsioon ja βg on genotüübile vastav parameeter võrdeliste riskide mudelis. Tähistades haigestumise indikaatori Yi = I(Ti 6 ti), avaldub haigestumise tõenäosus täiend log-log mudeli kehtimise korral järgmiselt: P(Yi = 1|ti, gi) = 1 −exp(−exp(β0 + β0 ggi + βtti)), kus β0, β0 gja βt on vastavalt täiend-log-log mudeli vabaliige, genotüübile vastav parameeter ja liitumisvanusele vastav parameeter. Näeme, et mõlema mudeli korral saame üleelamistõenäosuse kirjutada liitumisvanusest sõltuva funktsiooni ja genotüübist sõltuva funktsiooni kaudu: võrdeliste riskide mudeli puhul saame selle kirjutada kujul P(Ti 6 ti| gi) = 1 −exp(−f1(ti) exp(βggi)), kus f1(ti) = H0(ti), ja täiend-log-log mudeli puhul kujul P(Yi = 1|ti, gi) = 1 −exp(−exp(β0 + βtti) exp(β0 ggi)) =

= 1 −exp(−f2(ti) exp(β0 ggi)), 31 kus f2(ti) = exp(β0 +βtti). Kuna mõlema mudeli puhul saame hinnata ühtesama haigestumise tõenäosust, siis nüüd juhul, kui f1(ti) = f2(ti) valimi liitumisvanuste ti korral, siis ka βg = β0 g. Et f1 on kumulatiivne baasriskifunktsioon ja võib olla kuitahes keerulise kujuga, kuid f2 on alati e aste lineaarsest liitumisvanuse funktsioonist, siis üldjuhul need funktsioonid võrdsed ei ole. Praktikas võivad need funktsioonid siiski sarnased olla. Näitena on joonisel 4 kujutatud f1 ja f2 juhul, kui haigestumisandmed on simuleeritud vastavalt Weibulli jaotusega võrdeliste riskide mudelile. ) t (f 0,25 0,20 0,15 0,10 0,05 0,00 30 50 70 90 Liitumisvanus t f1(t) = (t/b)a f2(t) = exp(β^0 + β^tt) Joonis 4. Liitumisvanusest t sõltuvad funktsioonid f1 ja f2 vastavalt Weibulli jaotusega võrdeliste riskide mudeli korral ja täiend-log-log mudeli korral. Parameetrite väärtusteks on a = 3,4, b = 134, βˆ0 = −6,28, βˆt = 0,0561. Selle näite jaoks on i-nda (i = 1, . . . , 1000) indiviidi haigestumise vanus Ti genereeritud Weibulli jaotusest