W a, b {exp(βggi)}−1 a , kus a = 3,4 ja b = 134 on Weibulli jaotuse baasparameetrid, gi ∈ {0, 1, 2} on i-nda indiviidi 32 genotüüp ning βg = log(1,5) ≈0,405 on genotüübile vastav parameeter. Sel juhul saame Weibulli jaotuse kumulatiivse riskifunktsiooni kuju (4) põhjal, et f1(t) = t 134 3,4 . Liitumisvanused ti on valitud tagasipanekuta juhuvalikuga TÜ EGV geenidoonorite liitumis vanuste seast ja i-nda indiviidi haigestumise indikaatoriks on Yi = I(Ti 6 ti). Binaarsele haigestumise tunnusele on sobitatud täiend-log-log mudel, kus parameetrite väärtusteks hinnati βˆ0 = −6,28, βˆt = 0,0561 ja βˆ0 g = 0,404. Seega f2(t) = exp(−6,28 + 0,0561t). Näeme, et funktsioonide f1 ja f2 väärtused on genereeritud liitumisvanuste korral sarnased ja ka genotüübile vastava parameetri hinnang täiend-log-log mudelist on sarnane õigele parameetrile. Seega, praktikas võime öelda, et kui f1 ≈f2 vaatluse all olevate jälgimisaegade korral, siis ka βg ≈β0 g. Sedasama seost uurime aga ka simulatsiooniuuringus. 2.4 Meta-analüüs Meta-analüüsiks ehk analüüside analüüsiks nimetatakse enamasti erinevate uuringute tulemuste statistilist analüüsi eesmärgiga neid tulemusi koondada ja üldistada (Glass, 1976). Meta-analüüsi kasutatakse peamiselt selliste uuringute koondamiseks, kus on käsitletud samu teemasid, näiteks uuritud mingi kindla ravi

efektiivsust või geenivariantide mõju teatud haigusele. Üksikutel geneetilistel seoseuuringutel ei ole tihti piisavalt võimsust, et tuvastada väikese sage duse ja efektiga geenivariante. Tänu mitme uuringu kombineerimisele on võimalik saavutada suuremat statistilist võimsust ja tuvastada selliseid seoseid, mis üksikute analüüside puhul jääk sid leidmata. Selles töös kasutatakse meta-analüüsi, kombineerimaks binaarsete mudelite ja võrdeliste riskide 33 mudelite hinnanguid vastavalt jälgimiseelsete ja jälgimisaegsete juhtude analüüsist. Lihtsaim meta-analüüsi tüüp on fikseeritud mõjudega meta-analüüs, mille puhul eeldatakse, et igas kaasatavas uuringus on mingi uuritava töötluse (näiteks mõne ravimeetodi või geeni variandi) tegelik mõju β samasugune. Olgu meil vaatluse all n sõltumatut uuringut, kus on leitud hinnangud βˆi (i = 1, . . . , n) tege likule mõjule β. Olgu σi nende hinnangute standardvigade hinnangud ja tähistame wi = 1/σ2 i. Siis fikseeritud mõjudega mudeli korral avaldub kombineeritud hinnang tegelikule mõjule üksik hinnangute kaalutud keskmisena: Pn β ˆ = i=1 wiβˆ Pi i=1 wi, n kus i-nda (i = 1, . . . , n) uuringu kaaluks on wi/Pn j=1wj. Seega, mida täpsem on uuringus lei tud hinnang ehk mida väiksem on tema standardvea hinnang, seda suurema kaaluga ta kombi neeritud hinnangusse panustab. Selline kaalude valik tagab, et meta-analüüsi hinnang on nihketa (eeldusel, et esialgsed hinnangud βˆi (i = 1, . . . , n) on nihketa) ja vähima võimaliku dispersioo niga. Kombineeritud hinnangu β ˆ dispersioon on D(β ˆ ) = (Pn i=1wi)−1. Selles töös kasutamegi fikseeritud mõjudega meta-analüüsi ja selle läbiviimiseks

kasutame R-i paketti metafor (Viechtbauer, 2010). 34 3 Simulatsiooniuuring Erinevate analüüsimeetodite – Coxi ja martingaalijääkide mudeli ning logistilise ja täiend-log log mudeli – võrdlemiseks viiakse läbi simulatsiooniuuring. Et võrdeliste riskide mudel on kirja pandud riskifunktsioonide kaudu, kuid elukestusanalüüsi rakendamiseks on vaja teada konk reetseid elukestuseid, siis kõigepealt kirjeldatakse, kuidas genereerida haigestumisandmeid nii, et need vastaksid etteantud võrdeliste riskide mudelile. 3.1 Haigestumisandmete simuleerimine Simulatsiooniuuringu läbiviimiseks tuleb genereerida haigestumist kirjeldavad andmed. Haigestumise vanuse jaotusena kasutame Weibulli jaotust. Olgu n valimi suurus. Algoritm haigestumisandmete genereerimiseks nii, et efektialleel suurendaks või vähendaks haigestu mise riskimäära vastavalt kindlaksmääratud riskimäärade suhtele, on järgmine: iga indiviidi i = 1, . . . , n korral

1. genereeritakse genotüüp Gi võimalike väärtustega 0, 1, 2 (mis tähistavad efektialleelide arvu), kasutades binoomjaotust B(2, MAF), kus MAF on fikseeritud efektialleeli (harvema alleeli) sagedus;
2. genereeritakse elukestus Ti Weibulli jaotusest W(α, β);
3. genereeritakse haigestumise vanus Hi Weibulli jaotusest W(a, bg), kus bg on b0, b1 või b2 vastavalt sellele, kas inimese genotüüp on 0, 1 või 2;
4. valitakse vaatlusvanuseks Di minimaalne haigestumise vanusest ja elukestusest: Di = min{Hi, Ti};
5. märgitakse ta haigeks, kui tema vaatlusvanus on haigestumise vanus, seega haigestumist kirjeldab indikaator Yi = I(Hi 6 Ti) = I(Di = Hi). 35 Selliseks simuleerimiseks tuleb kõigepealt leida sobivad Weibulli jaotuse parameetrid α ja β elukestuse jaoks ning a, b0, b1 ja b2 haigestumise vanuse jaoks. Üks võimalus sobivate parameetrite leidmiseks on kasutada suurima tõepära meetodit ja hin nata need parameetrid näiteks TÜ EGV andmetelt. Selleks saab kasutada R-i paketi survival funktsiooni survreg. Selles töös kasutamegi elukestuse genereerimiseks geenivaramu andmetelt hinnatud parameetreid α = 9 ja β = 90. Joonisel 5 on näitena kujutatud elukestuse ja haigestumise vanuse tihedusfunktsioonid juhul, kui elukestuse parameetrid on α = 9 ja β = 90 ning haigestumise vanuse parameetrid on a = 7 ja b0 = 120. Selliste parameetrite korral on haigestumise mediaanvanus umbes 82 aastat ja haiguse levimus umbes 11%. s u de hi T 0,04 0,03

0,02 0,01 0,00 0 50 100 150 Vanus Elukestus Haigestumise vanus Joonis 5. Elukestuse ja haigestumise vanuse tihedusfunktsioonid, kui elukestus on Weibulli jaotusega W(9, 90) ja haigestumise vanus on Weibulli jaotusega W(7, 120) Haigestumise vanuse parameetrid on samuti võimalik hinnata huvipakkuva haiguse andmete põhjal. Et aga simulatsiooniuuringu käigus soovime näha ka seda, kuidas hinnangud käitu 36 vad haiguse erinevate levimuste korral, siis arvutame need parameetrid ise. Selleks fikseerime esmalt simuleeritava haiguse mediaanvanuse m haigestumise hetkel ja

soovitava haiguse levi muse p baasgrupis ehk nende inimeste seas, kellel Gi = 0. Kuna haigestumise vanuse jaotuse parameetrid sõltuvad inimese genotüübist, siis leitakse esmalt Weibulli jaotuse baasparameetrid a ja b0, millega genereeritakse haigestumise vanus neile, kellel ei ole ühtegi efektialleeli ehk kelle genotüüp on Gi = 0. Parameetrid b1 ja b2 saab seejärel ar vutada fikseeritud riskimäärade suhte kaudu. Olgu elukestust tähistav juhuslik suurus T ∼W(α, β) ja sellest sõltumatu haigestumise va nust baasgrupis tähistav juhuslik suurus H0 ∼ W(a, b0) jaotus- ning tihedusfunktsioonidega vastavalt FT , fT ja FH0, fH0. Kuna haigestumise vanus ja elukestus genereeritakse igale vaatlusele ning haigeks märgitakse need inimesed, kelle puhul on simuleeritud haigestumise vanus väiksem kui simuleeritud elu kestus, siis määrab baasgrupi levimuse p tõenäosus P(H0 < T). Avaldades selle võrduse elu kestuse ja haigestumise vanuse jaotus- ning tihedusfunktsioonide kaudu, leiame esimese seose a ja b0 vahel: p = P(H0 < T) = = P(H0 − T < 0) =

Z ∞ 0 Z ∞ Z t 0 fH0(h)fT (t) dh dt =

FH0(t)fT (t) dt =

0 Z ∞0 1 − exp − t b0 a α β t β α−1 exp − t β α dt. Haigestumise mediaanvanuse puhul võtame arvesse ainult neid inimesi, kes said haiguse ehk kelle haigestumise vanus on väiksem kui elukestus. Seega, kui soovime, et haigestumise mediaan vanus baasgrupis oleks m, peab kehtima P(H0 < m|H0 < T) = 0,5. Kasutades tingliku tõenäosuse valemit, saame nüüd avaldada ka teise seose

parameetrite a ja b0 37 vahel: 0,5 = P(H0 < m|H0 < T) = =P(H0 < m, H0 < T) P(H0 < T)=

1 pZ m 0 Z ∞h fT (t)fH0(h) dt dh =

1 pZ m (1 − FT (h))fH0(h) dh = 0

1 pZ m 0 exp − h β α a b0 h b0 a−1 exp − h b0 a dh. Oleme saanud kaks võrrandit a ja b0 määramiseks. Arvestades seda, et tegu on keeruliste integraalidega, millel üldjuhul lihtsat analüütilist kuju ei leidu, kasutame a ja b0 määramiseks numbrilist integreerimist. Selles töös on seda tehtud programmi MATLAB abil, kasutades mittelineaarsete võrrandisüsteemide numbriliseks lahendamiseks mõeldud fsolve funktsiooni (MATLAB Optimization Toolbox, 2021). Tabelis 1 on näitena toodud leitud lahendid haigestu mise mediaanvanuse m = 60 korral ja viie erineva baasgrupi levimuse p korral. Tabel 1. Haigestumise vanust kirjeldava Weibulli jaotuse kujuparameeter a ja skaalaparameeter b0 erinevate baaslevimuste korral, kui haigestumise mediaanvanus baasgrupis on 60 Baaslevimus p a b0 0,01 1,9254 936,7213 0,05 1,9562 391,9832 0,10 1,9974 264,8140 0,20 2,0897 175,7165 0,30 2,1989 136,7687

Selleks, et haigestumine sõltuks genotüübist – efektialleeliga isikute seas oleks haigestumise riskimäär kas suurem või väiksem kui baasgrupis – tuleb haigestumise vanust kirjeldava Weibulli jaotuse skaalaparameetrit muuta vastavalt efektialleelide arvule ja fikseeritud riski määrade suhtele. Selleks, et riskimäärade suhe ühe võrra suurema arvu efektialleelide suhtes oleks HR, peab 38 võrduse (5) põhjal kehtima b1 = b0 {exp(βg · 1)}−1a ja analoogiliselt b2 = b1 {exp(βg · 1)}−1a= b0 {exp(βg · 2)}−1a, kus βg on log(HR). Seega, üldjuhul genereeritakse i-ndale indiviidile haigestumise vanus Weibulli jaotusest W a, b0 {exp(βgGi)}−1 a , kus Gi ∈ {0, 1, 2} on inimese genotüüp. Kood selliste andmete simuleerimiseks R-is on toodud lisas 3. 3.2 Simulatsiooniplaan Simulatsiooniuuringus käsitleme kolme olukorda. Esimesel juhul eeldame, et meil on olemas sünnikohort, kus kõik inimesed on vaatluse all sünnist kuni surmani. Seega teame täpselt iga inimese haigestumise vanust ja elukestust, ning iga haiguse saanud inimene on jälgimisaegne juht. See on nii-öelda ideaalne olukord, mida päriselus tegelikult ette ei tule. Teisel juhul eel dame, et inimesed liituvad geenivaramuga oma elu jooksul ja osa juhtudest on seega jälgimis eelsed, osa jälgimisaegsed. Kolmanda

stsenaariumi puhul eeldame nagu teiseski stsenaariumis, et geenivaramuga liitutakse elu jooksul, kuid enam me jälgimiseelsetel ja jälgimisaegsetel juh tudel vahet ei tee. Simulatsioon viiakse läbi erinevate efektisuuruste HR ∈{0,9; 0,95; 1; 1,05; 1,1; 1,2; 1,5}, harvema alleeli sageduste MAF ∈{0,05; 0,1; 0,2; 0,3; 0,4; 0,5}, haigestumise mediaanvanuse m = 60 ja haiguse baasgrupi levimuste p ∈{0,05; 0,1; 0,2} korral. Valimimahuna kasutati N = 100 000, mis on ligikaudu võrdne TÜ EGV kohordi suurusega, kust on välja jäetud lä hisuguluses olevad geenidoonorid. Erinevate mudelite tulemuste võrdlemiseks hinnati iga para meetrite kombinatsiooni korral iga stsenaariumi ja iga mudeli jaoks nsim = 1000 simulatsiooni 39 keskmisena hinnangute empiiriline nihe: 1 nsim Xn sim k=1 (βˆk − log(HR )), ruutjuur keskmisest ruutveast (RMSE, root mean square error): vuu t 1nsimXn sim k=1 ja võimsus, kasutades olulisuse nivood 0,05: (βˆk −log(HR))2 1 nsim Xn sim k=1 I(Pk < 0,05), kus βˆk ja Pk on vastavalt k-nda (k = 1, . . . , nsim) valimi puhul leitud hinnang efektisuuruse logaritmile ja mudelist saadud statistilise testi p-väärtus. Juhul kui tegelik efektisuurus on HR = 1, saab viimase valemi järgi arvutada empiirilise I liiki vea tõenäosuse.

Kirjeldame nüüd täpsemalt, kuidas andmete genereerimine ja analüüs erinevate stsenaariumite puhul läbi viiakse. Olgu meil genereeritud andmed nii, nagu on kirjeldatud peatüki 3.1 alguses. Seega on meil iga indiviidi i = 1, . . . , n jaoks olemas genotüüp Gi, elukestus Tija genotüübist sõltuv haigestumise vanus Hi. Esimene stsenaarium Esimese stsenaariumi puhul hinnatakse Coxi võrdeliste riskide mudel, kasutades ajaskaalana inimese vanust haigestumise ajal ja ainsa argumenttunnusena kaasatakse genotüüp. Et inimesed liitusid sündides, siis inimese vanus haigestumise ajal on samal ajal ka aeg liitumisest haiges tumiseni. Lisaks lähendatakse sama mudelit martingaalijääkide abil. Selleks sobitatakse ilma kovariaatideta Coxi mudel ja hinnatakse sellelt skaleeritud martingaalijäägid, millele sobita takse seejärel lineaarne mudel, kus ainsaks seletavaks tunnuseks on genotüüp. 40 Teine stsenaarium Teise ja kolmanda stsenaariumi puhul genereeritakse igale inimesele lisaks ka sünniaasta Si ühtlasest jaotusest U(1920, 1980) ja liitumisaasta Li ühtlasest jaotusest U(2000, 2020). Sellise genereerimise puhul jäetakse andmestikku alles vaid need inimesed, kes elasid vähemalt liitu miseni ehk kelle elukestus on suurem kui liitumisvanus: Ti > Li−Si. Samuti valitakse nii-öelda katkestusaeg K, millest alates on kõik inimesed tsenseeritud. Geenivaramu andmete analüüsi misel on katkestusajaks tavaliselt kuupäev, mil saabusid viimased andmed meditsiiniandmete allikatest. Simulatsiooniuuringus on katkestusajaks valitud aasta 2050 – inimesed, kes jäid hai geks või surid pärast seda aastat, on tsenseeritud, ja nende vaatlusvanus on vanus aastal 2050. Seega on inimese vaatlusvanuseks minimaalne haigestumise vanusest, elukestusest ja vanusest katkestusajal: Di = min{Hi, Ti, K −Si}, ja haigeks märgitakse need, kelle vaatlusvanus on nende haigestumise vanus, seega haigestumist kirjeldab indikaator Yi = I(Di = Hi).

Teise stsenaariumi puhul jagatakse andmestik enne analüüsi haigestumis- ja liitumisaja järgi kaheks. Ühte andmestikku jäävad kõik jälgimiseelsed juhud ja juhuslikult valitud terved kont rollid nii, et iga juhu kohta on andmestikus neli kontrolli. Terveteks kontrollideks on need, kes ei olnud haigestunud enne ega pärast liitumist. Kontrollide ja juhtude suhe 4:1 on juht kontrolluuringutes laialt kasutusel, sest on näidatud, et sellest suurema suhte kasutamine testide statistilist võimsust oluliselt ei paranda (Hong ja Park, 2012). Teise andmestikku jäävad kõik jälgimisaegsed juhud ja kõik need terved kontrollid, keda ei valitud esimesse andmestikku. Esimesele andmestikule hinnatakse nii logistilise regressiooni mudel kui täiend-log-log mudel, kus funktsioontunnuseks on haiguse olemasolu kirjeldav binaarne tunnus ja seletavaks tunnu seks genotüüp. Kovariaadina lisatakse mudelisse ka inimese liitumisvanus. Liitumisvanus on see vanus, mille korral me juhtude puhul teame, et nad olid selleks ajaks haigeks jäänud, kont rollide puhul teame, et selles vanuses olid nad terved. Teisele ehk jälgimisaegsete juhtude andmestikule hinnatakse Coxi võrdeliste riskide mudel, kasutades ajaskaalana aega liitumisest: Di −(Li −Si). Seletavaks tunnuseks on genotüüp ja kovariaadina on arvesse võetud liitumisvanus. Seda mudelit lähendatakse ka martingaali 41 jääkide abil, sobitades esmalt Coxi mudeli, kuhu kovariaadina on lisatud vaid liitumisvanus. Selle mudeli skaleeritud martingaalijääkidele hinnatakse seejärel lineaarne mudel, kus seleta vaks tunnuseks on genotüüp. Kahe andmestiku analüüsi tulemused meta-analüüsitakse kokku, kasutades fikseeritud mõju dega meta-analüüsi mudelit. Nii Coxi mudeli kui martingaalijääkide mudeli hinnangud kom bineeritakse nii logistilise kui täiend-log-log mudeli hinnangutega, seega saadakse kokku neli meta-analüüsi hinnangut. Kolmas stsenaarium