UNIVERZITA PALACKÉHO V OLOMOUCI P°írodov¥decká fakulta Katedra matematické analýzy a aplikací matematiky
DIPLOMOVÁ PRÁCE
Pouºití logistické regrese pro diagnostiku výskytu rakoviny prostaty
Vedoucí diplomové práce:
Vypracovala:
Mgr. Ond°ej Vencálek, Ph.D.
Bc. Kamila Fa£evicová
Rok odevzdání: 2012
AME, 2. ro£ník
Prohlá²ení Prohla²uji, ºe jsem diplomovou práci zpracovala samostatn¥ pod vedením pana Ond°eje Vencálka a výhradn¥ s pouºitím uvedené literatury. V Olomouci dne 25. b°ezna 2012
Pod¥kování Na tomto míst¥ bych ráda pod¥kovala vedoucímu práce Ond°eji Vencálkovi za jeho trp¥livost, rady i cenné p°ipomínky. Zárove¬ chci pod¥kovat rodin¥ a p°átel·m za jejich podporu v pr·b¥hu celého studia. Touto cestou d¥kuji i doktoru Michalu Greplovi, ºe mi prost°ednictvím pana Tomá²e Fürsta dal k dispozici rozsáhlá data, která byla hlavní inspirací pro vznik této práce.
Obsah
1 Úvod
5
2 Rakovina prostaty
6
2.1
Incidence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.2
Mortalita (úmrtnost) . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.3
Lifetime Risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.4
Rizikové faktory . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.4.1
V¥k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.4.2
Rasa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.4.3
Rodinná anamnéza . . . . . . . . . . . . . . . . . . . . . .
12
2.4.4
Obezita . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.4.5
Sexuální aktivita . . . . . . . . . . . . . . . . . . . . . . .
13
Diagnostika . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.5.1
Prostatický speci cký antigen PSA . . . . . . . . . . . . .
13
2.5.2
Digitální rektální vy²et°ení DRE
. . . . . . . . . . . . . .
14
2.5.3
Transrektální sonogra e TRUS . . . . . . . . . . . . . . .
15
2.5.4
Biopsie . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
Lifetime risk - odhad . . . . . . . . . . . . . . . . . . . . . . . . .
15
2.5
2.6
3 Základní poznatky o kategoriálních datech
18
3.1
Kategoriální prom¥nná . . . . . . . . . . . . . . . . . . . . . . . .
18
3.2
Binomické rozd¥lení . . . . . . . . . . . . . . . . . . . . . . . . . .
19
3.2.1
Odhad parametru π
. . . . . . . . . . . . . . . . . . . . .
20
3.2.2
Testy odhadu π ˆ . . . . . . . . . . . . . . . . . . . . . . . .
21
3.2.3
Intervaly spolehlivosti
. . . . . . . . . . . . . . . . . . . .
23
Multinomické rozd¥lení . . . . . . . . . . . . . . . . . . . . . . . .
25
3.3.1
Odhad parametru πj . . . . . . . . . . . . . . . . . . . . .
26
3.3.2
Testy multinomického parametru . . . . . . . . . . . . . .
27
3.3
3
4 Kontingen£ní tabulky 4.1
31
Vztahy mezi prom¥nnými v tabulce 2x2 . . . . . . . . . . . . . . .
35
4.1.1
Senzitivita a speci cita testu . . . . . . . . . . . . . . . . .
36
4.1.2
Rozdíl proporcí, relative risk . . . . . . . . . . . . . . . . .
39
4.1.3
Pom¥r ²ancí . . . . . . . . . . . . . . . . . . . . . . . . . .
42
5 Logistická regrese
46
5.1
Interpretace parametr· . . . . . . . . . . . . . . . . . . . . . . . .
47
5.2
Odhad parametru β
. . . . . . . . . . . . . . . . . . . . . . . . .
49
5.2.1
Newton-Raphsonova metoda . . . . . . . . . . . . . . . . .
50
5.2.2
Varian£ní matice odhadu parametr· . . . . . . . . . . . . .
53
5.2.3
Intervaly spolehlivosti a testování parametru βj . . . . . .
54
Hodnocení kvality modelu . . . . . . . . . . . . . . . . . . . . . .
54
5.3.1
Deviance . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
5.3.2
Rezidua . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
5.3.3
Modi kace koe cientu determinace . . . . . . . . . . . . .
58
5.3.4
Akaikeho informa£ní kritérium . . . . . . . . . . . . . . . .
58
5.3
6 P°íklad pouºití logistické regrese - diagnostika rakoviny prostaty 60 6.1
Hledání modelu . . . . . . . . . . . . . . . . . . . . . . . . . . . .
60
6.2
Odhad parametr· . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
6.2.1
Po£áte£ní odhady . . . . . . . . . . . . . . . . . . . . . . .
63
6.2.2
Iterace Newton-Raphsonovy metody . . . . . . . . . . . .
65
6.2.3
Intervaly spolehlivosti a testování významnosti parametr·
67
6.3
Ov¥°ení kvality modelu . . . . . . . . . . . . . . . . . . . . . . . .
69
6.4
Interpretace výsledk· . . . . . . . . . . . . . . . . . . . . . . . . .
70
6.5
Srovnání s modelem pro hyperplazii prostaty . . . . . . . . . . . .
73
7 Záv¥r
75
Pouºitá literatura a webové odkazy
77
4
1
Úvod
I v dne²ní dob¥, kdy je medicína jiº na velmi vysoké úrovni, existuje mnoho chorob, jejichº vznik, chování a lé£ba nejsou je²t¥ dokonale pochopeny a popsány a jsou tak stále p°edm¥tem v¥deckého bádání. Typickým p°íkladem takové choroby je rakovina, o které se sice jiº mnoho ví, ale která svojí závaºností stále vyºaduje pozornost mnoha odborník·. Jedním ze základních p°edpoklad· úsp¥²né lé£by rakoviny je její v£asné odhalení a toto tvrzení platí dvojnásobn¥ v p°ípad¥ zhoubného novotvaru prostaty, který je v rané fázi velmi dob°e lé£itelný, jak ukazují £etné studie. Cílem této práce je s vyuºítím záznam· z vy²et°ení, která prob¥hla na Urologické klinice Fakultní nemocnice Olomouc, kvanti kovat, jaký vliv mají jednotlivé p°i vy²et°ení zji²´ované ukazatele na výsledek biopsie resp. na pravd¥podobnost, ºe biopsie odhalí karcinom prostaty a pokusit se tak p°isp¥t k v£asné detekci tohoto novotvaru. Úkolem druhé kapitoly je seznámení £tená°e s problematikou rakoviny prostaty. Kapitola tedy obsahuje informace o po£tu nov¥ diagnostikovaných p°ípad· a také po£tu úmrtí v d·sledku této choroby. Tyto údaje jsou srovnány s údaji sv¥tovými a také s obdobnými ukazateli vztahujícími se k jiným typ·m karcinom·. Druhá kapitola navíc obsahuje informace o rizikových faktorech a moºnostech diagnostiky, která zahrnuje vyhodnocení ukazatel·, které budou pouºity v praktické £asti práce. Kapitolu t°etí a £tvrtou lze nazvat p°ípravnými kapitolami pro kapitolu pátou, v¥nující se logistické regresi. V této £asti práce jsou tedy popsány základní vlastnosti kategoriálních prom¥nných a také metody práce s kontingen£ními tabulkami. Kaºdý nový pojem je objasn¥n na jednoduchém p°íkladu. Pátá kapitola je v¥novaná logistické regresi, která bude v následující kapitole pouºita ke zpracování vstupních dat. Kapitola vyuºívá pojm· zavedených ve t°etí a £tvrté kapitole a v¥nuje se jak tvorb¥ vhodného modelu a ov¥°ení jeho kvality, tak i metod¥ odhadu jeho parametr· a jejich interpretaci. V této kapitole jiº p°íklady chybí, problematika pouºití logistické regrese je totiº objasn¥na v praktické £ásti práce. Ve v²ech t°ech teoretických kapitolách jsou pouºita v²eobecn¥ známá tvrzení, která jsou zde uvád¥na bez d·kazu a která lze nalézt nap°íklad v [15] nebo v [17]. Poslední ²está kapitola je jiº ryze praktická a zabývá se zejména hledáním modelu vhodného pro popis vztahu mezi veli£inami m¥°enými p°i vy²et°ení a pravd¥podobností pozitivní biopsie. Po nalezení takového modelu a odhadu jeho parametr· následuje jejich interpretace.
5
2
Rakovina prostaty
Rakovina prostaty je jedním z nej£ast¥j²ích nádorových onemocn¥ní muº· jak u nás, tak i ve sv¥t¥, jejíº po£et nových p°ípad· se neustále zvy²uje. Vzhledem k této vysoké incidenci má také pom¥rn¥ vysokou mortalitu (úmrtnost), která má ale naopak klesající trend vývoje. Nádor odhalený v po£áte£ním stadiu je pom¥rn¥ dob°e lé£itelný a proto je d·leºité dbát na jeho v£asnou diagnostiku, zaloºenou a´ uº na hodnotách prostatického speci ckého antigenu a z n¥j odvozených charakteristik, tak i na analýze alternativních indikátor·.
2.1
Incidence
Vzhledem k dob¥ pot°ebné ke sb¥ru a zpracování dat, je zatím posledním zve°ejn¥ným údajem o po£tu p°ípad· v eské republice údaj z roku 2008. V tomto roce bylo diagnostikováno 5332 nových p°ípad· rakoviny prostaty, tedy asi 104 p°ípad· na 100 000 muº· ºijících v R [1]. Jak je vid¥t z Tabulky 1, jedná se o zatím nejvy²²í po£et nových onemocn¥ní. Jiº na tomto míst¥ je ale d·leºité zmínit, ºe tento nár·st nemusí být nutn¥ zp·soben jen ur£itou v¥t²í agresivitou tohoto druhu rakoviny, ale m·ºe za ním být i v¥t²í informovanost muº· o chorob¥, d·sledn¥j²í diagnostika spojená se zavedením sledování hladiny PSA do b¥ºné praxe a také prodluºování délky ºivota. Rok Incidence
1977
1978
1979
1980
1981
1982
1983
celkem
1115
1251
1147
1245
1268
1466
1381
1352
na 100 tisíc
22,51
25,12
22,92
24,82
25,37
29,29
27,57
26,96 1992
Rok Incidence
1985
1986
1987
1988
1989
1990
1991
celkem
1356
1487
1533
1561
1634
1607
1785
1952
na 100 tisíc
27,02
29,61
30,49
31,01
32,45
31,9
35,66
38,94 2000
Rok Incidence
1993
1994
1995
1996
1997
1998
1999
celkem
1986
2174
2287
2659
2707
2895
2916
2869
na 100 tisíc
39,57
43,3
45,59
53,05
54,05
57,84
58,31
57,39 2008
Rok Incidence
1984
2001
2002
2003
2004
2005
2006
2007
celkem
3194
3499
3853
4397
4939
4902
5188
5332
na 100 tisíc
64,29
70,45
77,45
88,28
98,73
97,53
102,07
103,81
Tabulka 1: Vývoj incidence rakoviny prostaty v R v letech 1977 aº 2008 Ve srovnání s jinými typy nádorových onemocn¥ní tvo°ila rakovina prostaty asi 19 % v²ech u nás hlá²ených muºských novotvar· a v roce 2008 tak byla nej£ast¥j²ím nádorovým onemocn¥ním muº· p°ed rakovinami plic a tlustého st°eva, kterých bylo zji²t¥no 4686 a 2696, tedy 16 % a 9 % [1]. Z tohoto srovnání jsou ov²em vy°azeny koºní nádory, které jsou sice nejv¥t²í, ale také velmi speci ckou skupinou. Na Obrázku 1. je zobrazeno srovnání vývoje incidencí rakovin prostaty, plic a tlustého st°eva. 6
Obrázek 1: Srovnání vývoje incidencí rakovin prostaty, plic a tlustého st°eva v R v letech 1977 aº 2008, p°epo£tených na 100 000 muº· Vývoj ve sv¥t¥ je podobný jako u nás. Ve Spojených státech amerických se odhaduje, ºe v roce 2011 bude podíl rakoviny prostaty na v²ech onemocn¥ních £init 29 % oproti 14 % a 9 % pro rakoviny plic resp. tlustého st°eva, viz [3]. Podle údaj· z roku 2002 se po p°epo£tu incidence karcinomu prostaty na 100 000 osob, °adí eská republika na 18. místo v Evrop¥. Nejvy²²í incidence je pak moºná trochu p°ekvapiv¥ v severských zemích, zejména védsku, viz [1]. Podrobn¥j²í srovnání evropských zemí je na Obrázku 2 [1].
Obrázek 2: Srovnání incidence karcinomu prostaty v Evrop¥ v roce 2002 p°epo£tené na 100 000 osob
7
I v Olomouckém kraji lze pozorovat výrazný rostoucí trend vývoje, zatímco v roce 1998 byla rakovina prostaty diagnostikována 161 muº·m, o deset let pozd¥ji uº jich bylo 337, coº znamená 51,5 resp. 107,42 p°ípad· na 100 000 muº· [2]. Podrobn¥j²í údaje o incidenci v kraji jsou uvedené v Tabulce 2. Rok Incidence
celkem na 100 tisíc
Rok Incidence
celkem na 100 tisíc
Rok Incidence
celkem na 100 tisíc
Rok Incidence
celkem na 100 tisíc
1977
1978
1979
1980
1981
1982
1983
91
103
79
105
86
121
92
1984 114
29,58
33,26
25,38
33,62
27,63
38,82
29,51
36,56
1985
1986
1987
1988
1989
1990
1991
1992
93
100
113
117
117
109
119
145
29,79
32,06
36,24
37,51
37,52
34,9
38,36
46,62
1993
1994
1995
1996
1997
1998
1999
2000
131
156
170
227
191
161
175
180
42,07
50,02
54,56
72,4
61
51,5
56,06
57,69
2001
2002
2003
2004
2005
2006
2007
2008
189
224
264
273
272
247
305
337
60,88
72,37
85,3
88,37
87,42
79,2
97,36
107,42
Tabulka 2: Vývoj incidence rakoviny prostaty v Olomouckém kraji v letech 1977 aº 2008
2.2
Mortalita (úmrtnost)
Vzhledem k vysoké incidenci je i po£et zem°elých v d·sledku rakoviny prostaty pom¥rn¥ vysoký. V roce 2008 zem°elo v eské republice na následky této choroby 1521 muº· [1]. Vývoj po£tu úmrtí zp·sobených karcinomem prostaty je uveden v Tabulce 3. Vy²²í mortalitu m¥la jen rakovina plic, které v roce 2008 podlehlo 4017 muº· a zatímco karcinomy prostaty zap°í£inily asi desetinu v²ech 15333 úmrtí muº· na následky rakoviny, u karcinomu plic to byla více neº £tvrtina [2]. Rok Mortalita
1977
1981
1982
1983
1984
260
477
543
671
677
733
840
841
5,25
9,58
10,85
13,38
13,55
14,65
16,77
16,77
1985
1986
1987
1988
1989
1990
1991
1992
827
897
990
952
1047
987
1044
1058
16,48
17,86
19,69
18,91
20,79
19,6
20,85
21,1 2000
celkem
Rok
1993
1994
1995
1996
1997
1998
1999
celkem
1109
1192
1262
1278
1284
1372
1369
1359
na 100 tisíc
22,09
23,74
25,16
25,5
25,64
27,41
27,37
27,18 2008
Rok Mortalita
1980
celkem
na 100 tisíc
Mortalita
1979
na 100 tisíc Rok Mortalita
1978
2001
2002
2003
2004
2005
2006
2007
celkem
1413
1338
1357
1454
1443
1475
1467
1521
na 100 tisíc
28,44
26,94
27,28
29,19
28,84
29,35
28,86
29,61
Tabulka 3: Vývoj úmrtnosti na rakovinu prostaty v R od roku 1977 do 2008
8
Pro srovnání je²t¥ uve¤me, ºe na následky t°etí nej£ast¥j²í rakoviny tlustého st°eva zem°elo v roce 2008 1352 muº· [1]. Vývoj mortality v²ech t°í chorob od roku 1977 je zachycen na Obrázku 3. P°i srovnání t¥chto údaj· je pot°eba vzít v úvahu, ºe rakovina prostaty postihuje zejména star²í muºe, jejichº úmrtí m·ºe být zp·sobeno mnoºstvím jiných chorob a u nichº tedy pokro£ilé stadium karcinomu nemusí být v·bec odhaleno.
Obrázek 3: Srovnání vývoje úmrtností na rakoviny prostaty, plic a tlustého st°eva v R, p°epo£tených na 100 000 muº· v letech 1977 aº 2008 Sv¥tový vývoj je op¥t podobný. V USA je karcinom prostaty s 11 % druhou nej£ast¥j²í p°í£inou smrti onkologických pacient·. Tou nej£ast¥j²í je op¥t rakovina plic, která se na po£tu zem°elých podílí z 28 %. T°etí nej£ast¥j²í p°í£inou úmrtí u muº· jsou pak op¥t nádory v okolí tlustého st°eva, viz [3]. Na Obrázku 4 je zachyceno srovnání evropských zemí na základ¥ údaj· o úmrtnosti na rakovinu prostaty z roku 2002, které jsou p°epo£teny na 100 000 osob. I v tomto srovnání se eská republika °adí na 18. pozici. Nelichotivé první místo op¥t pat°í védsku, tedy zemi s nejvy²²í incidencí rakoviny prostaty na 100 000. Za pov²imnutí také stojí Finsko, které má sice vysokou incidenci, ale z hlediska mortality je aº 12. nebo naopak Dánsko, které se na Obrázku 2 v·bec nevyskytovalo, ale nyní se nachází jiº na 4. míst¥. Zatímco incidence se v Olomouckém kraji rychle zvy²uje, jak jiº bylo uvedeno v p°edchozí kapitole, úmrtnost roste jen málo. V roce 2008 tak na následky nádoru prostaty zem°elo 95 muº·, tj. 30,28 na 100 000 muº· v kraji. Pro srovnání je²t¥ uve¤me, ºe v roce 1998, kdy byl po£et nov¥ zji²t¥ných p°ípad· polovi£ní, jich bylo 77 tedy 24,63 na 100 000 [1]. P°ehled v²ech hodnot je uveden v Tabulce 4.
9
Obrázek 4: Srovnání úmrtnosti na rakovinu prostaty v Evrop¥ v roce 2002 p°epo£tené na 100 000 muº· Rok Mortalita
celkem na 100 tisíc
Rok Mortalita
celkem na 100 tisíc
Rok Mortalita
celkem na 100 tisíc
Rok Mortalita
celkem na 100 tisíc
1977
1978
1979
1980
1981
1982
1983
21
34
37
32
49
39
67
1984 58
6,83
10,98
11,88
10,25
15,74
12,51
21,49
18,6
1985
1986
1987
1988
1989
1990
1991
1992
53
56
56
62
74
61
60
68
16,98
17,95
17,96
19,88
23,73
19,53
19,34
21,87
1993
1994
1995
1996
1997
1998
1999
2000
71
74
71
81
70
77
70
78
22,8
23,73
22,79
25,84
22,36
24,63
22,42
25
2001
2002
2003
2004
2005
2006
2007
2008
84
88
97
88
88
93
104
95
27,06
28,43
31,34
28,49
28,28
29,82
33,2
30,28
Tabulka 4: Vývoj úmrtnosti na rakovinu prostaty v Olomouckém kraji v letech 1977 aº 2008
2.3
Lifetime Risk
Pojem lifetime risk lze interpretovat jako pravd¥podobnost, se kterou jedinec v pr·b¥hu celého svého ºivota zaºije n¥jakou sledovanou událost. Touto událostí bude jednak to, ºe muº onemocní rakovinou a pak také, ºe v d·sledku tohoto nádorového onemocn¥ní zem°e. Americká studie (viz [4]) uvádí, ºe pravd¥podobnost, ºe bude muºi v pr·b¥hu jeho ºivota alespo¬ jednou diagnostikována rakovina libovolného typu je 44,29 %. Pokud jsou jednotlivé typy onemocn¥ní brány odd¥len¥, je riziko nejvy²²í práv¥ pro karcinom prostaty, pro n¥jº je rovno 16,22 %. Pro srovnání riziko onemocn¥ní rakovinou plic je 7,67 % a tlustého st°eva 5,3 %. Riziko, ºe muº zem°e na rako10
vinu je odhadnuto na 23,2 %. Pravd¥podobnost, ºe bude smrt zp·sobena nádorem prostaty je 2,79 %. Z tohoto hlediska se jako nejrizikov¥j²í ukazují nádory plic, k nimº náleºí pravd¥podobnost 6,95 %. Pravd¥podobnost úmrtí v d·sledku rakoviny tlustého st°eva je uvedena jako 2,17 %. V zemích Evropské unie je pravd¥podobnost onemocn¥ní rakovinou prostaty 6 %, pravd¥podobnost, ºe na ni muº zem°e je pak 1,1 %, viz [5]. I v eské republice je hodnota rizika vzniku rakoviny nejvy²²í pro karcinom prostaty, jeho odhad £iní 9,2 %. Pro nádory plic a tlustého st°eva je to pak 7,7 % a 4,7 %. Pravd¥podobnost, ºe u nás muº onemocní rakovinou libovolného typu, je 49,4 %. Pravd¥podobnosti úmrtí jsou pro prostatu 2,6 %, pro plíce 6,5 % a pro tlusté st°evo 2,3 %. A£koliv je lifetime risk £asto pouºívanou veli£inou, jedná se pouze o informativní ukazatel, jehoº pouºívání má n¥kolik úskalí. Ta budou spole£n¥ s metodou jeho odhadu uvedena v kapitole 2.6.
2.4
Rizikové faktory
Faktor· o kterých se p°edpokládá, ºe mají vliv na výskyt rakoviny prostaty je celá °ada. Zatímco u n¥kterých je jejich vliv dob°e znám a popsán, u jiných tento vliv stále je²t¥ není zcela jasný a pátrání po n¥m dokonce ob£as vede k protich·dným záv¥r·m. Hlavními faktory jsou v¥k, rodinná anamnéza a rasa. Mezi ty mén¥ významné a zárove¬ h·°e popsané pak pat°í obezita, t¥lesná a sexualní aktivita, p°íjem ºivo£i²ných tuk· a alkoholu, podstoupení vazektomie, kou°ení £i ºivotní prost°edí.
2.4.1 V¥k Jedním z nejvýznamn¥j²ích rizikových faktor· je v¥k. Zatímco u mladých muº· se toto onemocn¥ní tém¥° nevyskytuje (v roce 2008 byl diagnostikován jen jeden p°ípad u muºe mlad²ího 40 let a 29 p°ípad· u muº· v rozmezí 40 aº 49 let), po p°ekro£ení hranice padesáti let po£et p°ípad· rychle nar·stá, jak je vid¥t i z Tabulky 5, která zachycuje po£ty nádor· prostaty diagnostikovaných v roce 2008, rozd¥lené podle v¥ku pacienta [2]. Nár·st po£tu p°ípad· zachycený v Tabulce 1 tedy m·ºeme £áste£n¥ vysv¥tlit prodlouºením délky ºivota a i do budoucna lze o£ekávat, ºe pokud tato délka dále poroste, poroste i incidence tohoto typu nádorového onemocn¥ní. Vzhledem k tomu jak významný rizikový faktor v¥k je, doporu£uje se, aby byli muºi star²í 50-ti let pravideln¥ kontrolováni, aby byl p°ípadný novotvar v£as odhalen. 11
V¥k Incidence
0-39 1
40-44 6
45-49 23
50-54 174
55-59 505
60-64 994
65-69 1015
70-74 1019
75-79 902
80-84 452
85-89 241
Tabulka 5: Incidence rakoviny prostaty v R v roce 2008 rozd¥lená podle v¥ku
2.4.2 Rasa Rozdíl v incidenci v závislosti na rase je pom¥rn¥ dob°e patrný na populaci v USA. Srovnání incidence podle rasy je uvedeno v [3]. Zatímco b¥loch·m náleºí hodnota 143,8 diagnóz na 100 000 muº·, u afroameri£an· jich je 230. Nejniº²í incidence byla zpozorována u ameri£an· asijského p·vodu, pro které má hodnotu 81 p°ípad· na 100 000. Podobn¥ vyznívá i srovnání úmrtnosti, ta je op¥t nejvy²²í u afroameri£an· (54,2/100 tis.) a nejniº²í u asiat· (10,6/100 tis.), mortalita b¥loch· na nádory prostaty je 22,8 úmrtí na 100 000 muº·. Vliv rasy na rakovinu je tedy z°ejmý, jeho p°í£ina ov²em ne. Jednou z verzí je nap°íklad zvý²ené mnoºství cirkulujících androgen· u m·º· £erné pleti, viz [6].
2.4.3 Rodinná anamnéza T°etím významným faktorem je rodinná anamnéza, tedy údaje o tom, zda bratr, otec, d¥de£ek nebo t°eba strýc trp¥li zhoubným nádorem prostaty. Muº·m s pozitivní rodinnou anamnézou je vhodné v¥novat v¥t²í pozornost a provád¥t pravidelný skreening jiº p°ed dosaºením 50 let, protoºe riziko výskytu karcinomu prostaty je u nich zvý²ené. Zárove¬ platí, ºe vy²²í riziko je p°i onemocn¥ní bratra neº otce, dal²í nár·st je v p°ípad¥, ºe k tomuto onemocn¥ní do²lo p°ed dosaºením 60. roku ºivota nebo pokud nemocí trp¥li jak bratr tak i otec, viz [7]. Ve stejné literatu°e se také uvádí moºná trochu p°ekvapiv¥j²í spojitost rakoviny prostaty s rakovinou prsu matky £i sestry.
2.4.4 Obezita Pokud stanovujeme nadváhu jen pomocí indexu BMI, tedy pom¥ru hmotnosti a kvadrátu vý²ky, jsou výsledky studií vztahu mezi obezitou a výskytem karcinomu prostaty rozporuplné. Z tohoto hlediska se jeví lep²í pouºívat indexu WHR, který je de nován jako pom¥r obvodu pasu a bok·. P°i pouºití tohoto indexu se ukazuje, ºe u muº·, kte°í mají hodnotu WHR v¥t²í neº 0,9 je trojnásobn¥ vy²²í riziko diagnostiky rakoviny prostaty. Dal²í spojitost lze pozorovat mezi hodnotou BMI a rizikem úmrtí na tuto chorobu, kdy u osob s BMI vy²²ím neº 30 kg/m2 je toto riziko vy²²í. [8]
12
2.4.5 Sexuální aktivita asto jsem se setkávala s otázkou, zda výskyt rakoviny prostaty n¥jak souvisí se sexuální aktivitou. áste£ná odpov¥¤ je uvedena v [6]. Výsledky studií jsou ale bohuºel nejednozna£né. První studie se zabývala rozdílem ve výskytu karcinomu u katolických kn¥z· a svobodných muº·. Tento rozdíl se ukázal jako výrazný a záv¥rem bylo, ºe vy²²í riziko vzniku zhoubného nádoru prostaty je u kn¥z·. Jiná studie ov²em naopak ukázala, ºe vy²²í riziko onemocn¥ní se vyskytuje u muº·, kte°í byli za sv·j ºivot vícekrát ºenatí, vyst°ídali více partner· nebo prod¥lali pohlavní chorobu. Na základ¥ t¥chto dvou studií tedy sice nelze vliv sexuální aktivity na vznik nádoru vylou£it, ale nejde také ani ur£it, jakou podobu tento vliv vlastn¥ má.
2.5
Diagnostika
T°emi základními diagnostickými prost°edky jsou sledování hladiny prostatického speci ckého antigenu (PSA) a z n¥j odvozených charakteristik, digitální rektální vy²et°ení (DRV) a transrektální sonogra e (TRUS).
2.5.1 Prostatický speci cký antigen PSA PSA a je antigen speci cký pro prostatickou tká¬, který se za£al sledovat zhruba v 80. letech a znamenal pom¥rn¥ velký pr·lom v diagnostice rakoviny prostaty. Antigen je ale produkován jak zdravou, tak i zhoubnou tkání a jeho zvý²ené hodnoty tedy nemusí být nutn¥ zap°í£in¥ny zhoubným novotvarem, ale m·ºou indikovat nap°íklad i benigní hyperplasii prostaty. I p°es tuto nevýhodu je PSA nejvýznamn¥j²ím diagnostickým prost°edkem, protoºe jeho hodnoty jsou pom¥rn¥ snadno, rychle a levn¥ získatelné z krve. P°edm¥tem dal²ího bádání tedy z·stává zjistit, jak vyuºít m¥°ení hladiny PSA nejen pro odhalení novotvaru na prostat¥, ale také k rozli²ení, jestli je tento nádor zhoubný nebo ne. Malé mnoºství PSA je obsaºeno i v t¥le zdravých muº· a za hrani£ní hodnotu se v¥t²inou povaºují 4 ng/ml. Tato hodnota se v²ak nyní zdá být pom¥rn¥ vysoká, protoºe se vyskytují mnohé p°ípady, ve kterých byl karcinom odhalen i u muº·, kte°í m¥li hladinu PSA niº²í neº 4 ng/ml, viz [9]. I nadále se tedy hledá optimální hrani£ní hodnota, která by pomohla identi kovat maximum zhoubných novotvar· a zárove¬ nevedla k velkému po£tu zbyte£ných biopsií. P°irozená hladina se zvy²uje s v¥kem a pro lep²í senzitivitu (pravd¥podobnost odhalení choroby u nemocného £lov¥ka) je dobré pouºívat r·zné referen£ní hodnoty pro r·zné v¥kové skupiny. Normální hodnota PSA se u muº· ve v¥ku mezi 40 a 49 lety pohybuje od 0 do 2,5 ng/ml, u muº· o deset let star²ích od 0 do 3,5 ng/ml, ve v¥ku od 60 do 69 let je horní hranice jiº 4,5 ng/ml a u muº· star²ích dokonce 13
6,5 ng/ml, viz [10]. Dal²í nevýhodou tohoto ukazatele je, ºe hladina PSA kolísá a také nap°íklad to, ºe u ku°ák· je systematicky o 0,1 ng/ml niº²í, viz [11]. Od PSA je odvozeno n¥kolik dal²ích veli£in, které by m¥ly napomáhat k p°esn¥j²í diagnostice a k rozli²ení benigních a maligních nádor·. Jedná se zejména o densitu (PSAD), velocitu (PSAV) a PSA index. PSA densita je dána jako pom¥r hladiny PSA a objemu prostaty. A£koliv n¥které studie nazna£ují, ºe pacienti s karcinomem mají mají hodnotu PSAD signi kantn¥ vy²²í neº pacienti s benigním nádorem, a stanovují hrani£ní hodnotu na 0,15, není vhodné se na PSAD p°íli² spoléhat, protoºe její senzitivita je velmi nízká a velká £ást karcinom· by tak nebyla odhalena. Pozornost se nyní ubírá spí²e k pom¥ru hladiny PSA ku objemu tzv. tranzitorní zóny prostaty (PSAT), který má p°i hrani£ní hodnot¥ 0,35 senzitivitu vy²²í, viz [12]. PSA velocita vyjad°uje absolutní rychlost zm¥ny hladiny PSA za rok. Pro její zji²t¥ní je doporu£eno provést t°i m¥°ení v pr·b¥hu jednoho a p·l nebo dvou let. Jako hrani£ní se standartn¥ uvádí hodnota 0,75 ng/ml/rok, tuto hodnotu je v²ak t°eba upravit v závislosti na v¥ku pacienta a jeho sou£asné hladin¥ PSA, p°esn¥ji u mlad²ích muº· ve v¥ku od 50 do 59 let sníºit kritické hodnoty PSA velocity a PSA na 0,4 ng/ml/rok a 2ng/ml a obecn¥ p°i PSA > 4 ng/ml je výnamnou hodnotou jiº 0,4 ng/ml/rok. PSA velocita se ukazuje být dobrým indikátorem zejména agresivních, rychle rostoucích nádor· nebo v p°ípad¥, kdy je PSA < 4,0 ng/ml, viz [11],[12] a [13]. Dal²ím odvozeným indikátorem je PSA index. Ten vyuºívá toho, ºe je jedna £ást PSA vázaná na alfa-1-antichymotrypsin, £ást zbylou pak nazýváme volné PSA (free PSA, fPSA). PSA index je dán jako pom¥r hladiny volného ku celkovému PSA vyjád°ený v procentech. A jelikoº pravd¥podobnost vzniku zhoubného novotvaru roste s podílem vázaného PSA, jsou rizikové zejména nízké hodnoty PSA indexu. Za významnou hodnotu se u pacient· s PSA v rozmezí od 4 do 10 ng/ml a normální nebo mírn¥ zv¥t²enou prostatou povaºuje 20 %-25 %. Dále se ukazuje, ºe je pom¥r volné frakce PSA dobrým pomocníkem p°i rozli²ení karcinomu a benigní hyperplazie prostaty. P°i zhoubném novotvaru je jeho hodnota totiº významn¥ niº²í. Pokud je ov²em objem prostaty v¥t²í neº 75 cm3 , nelze tento pom¥r p°i diagnostice pouºít, viz [12].
2.5.2 Digitální rektální vy²et°ení DRE Jedná se o rutinní a pom¥rn¥ nenákladné vy²et°ení, které je v²ak t°eba brát spí²e jako podp·rnou diagnostickou metodu, protoºe jeho vyhodnocení je pom¥rn¥ subjektivní a ne v²echny typy nádor· jsou touto metodou zachytitelné. Obecn¥ lze ale °íct, ºe za podez°elé lze povaºovat jakékoliv takto nalezené tvrdé loºisko, viz [14]. 14
2.5.3 Transrektální sonogra e TRUS I tuto metodu lze povaºovat za podp·rnou. Stejn¥ jako v p°ípad¥ vy²et°ení per rectum je vyhodnocení tohoto vy²et°ení velmi subjektivní a nelze jím spolehliv¥ odhalit v²echny typy novotvar·. P°istupuje se k ní tedy zejména v p°ípad¥ podez°ení vzniklého p°i vyhodnocování PSA nebo DRE, kdy lze s její pomocí ur£it velikost prostaty £i samotného novotvaru a také lokalizovat místo, ve kterém se provede biopsie, viz [14].
2.5.4 Biopsie Pokud n¥kterou z vý²e uvedených metod vznikne podez°ení na výskyt zhoubného novotvaru, je pacient odeslán na biopsii, která by m¥la vést k potvrzení nebo vyvrácení tohoto podez°ení. B¥hem tohoto procesu se z r·zných míst prostaty odebírá v¥t²inou ²est vzork· tkán¥ a ty jsou následn¥ histologicky vy²et°eny, aby se ur£ilo, jakého typu tato tkᬠje. Tato tzv. sextantová biopsie, má ov²em své nevýhody. Jelikoº se p°i ní odebere jen ²est vzork· z celé prostaty, nemusíme se do postiºené tkán¥ jednodu²e tre t a výsledky biopsie tak mohou být negativní, a£koliv se nádor v t¥le skute£n¥ vyskytuje. Toto riziko lze sníºit nap°. odb¥rem více neº ²esti vzork· £i p°esn¥j²í lokalizací odebíraných vzork·. Více viz [14]. Výstupem biopsie je informace o tom, zda se na prostat¥ pacienta vyskytuje zhoubný novotvar, který je pot°eba okamºit¥ lé£it, nebo t°eba n¥který z benigních nádor·. Pokud se ov²em podez°ení na výskyt t¥chto novotvar· nepotvrdí, a£koliv p°edchozí vy²et°ení sv¥d£ila o opaku, je pacient dále sledován a v mnoha p°ípadech po n¥jakém £ase znovu podroben odb¥ru vzork· tkán¥ tzv. rebiopsii.
2.6
Lifetime risk - odhad
Jak jiº bylo uvedeno d°íve, pojem lifetime risk, který bychom mohli nahradit £eským ekvivalentem riziko, v sob¥ skrývá pravd¥podobnost nastoupení n¥jakého jevu v pr·b¥hu celého ºivota jedince, která je pro v¥t²í p°ehlednost vyjád°ena v procentech. V této kapitole se budu zabývat zp·sobem odhadu tohoto ukazatele, který potom aplikuji na data získaná z Ústavu zdravotnických informací a statistiky R, ÚZIS (viz [15]). M¥jme
• jevy Ai , i = 0, 1, 2, . . . , n, které zna£í ºe i-letý £lov¥k onemocní • a jevy Bi , i = 0, 1, 2, . . . , n, zna£ící, ºe se £lov¥k doºil v¥ku i let.
15
S Potom jev A = ni=0 Ai reprezentuje skute£nost, ºe £lov¥k onemocní n¥kdy v pr·b¥hu celého jeho ºivota. Protoºe jsou jevy Ai disjunktní, pro pravd¥podobnost jevu A platí ! n n [ X P (A) = P Ai = P (Ai ). i=0
i=0
Nyní vyuºijeme toho, ºe Ai = Ai ∩ Bi a následn¥ vztahu pro podmín¥nou prav, který je uveden nap°íklad v [17]. M·ºeme tedy d¥podobnost P (A|B) = P P(A∩B) (B) psát n n X X P (A) = P (Ai ∩ Bi ) = P (Ai |Bi )P (Bi ). i=0
i=0
Známe-li údaje o v¥kovém sloºení obyvatelstva, po£tu nových p°ípad· choroby a po£tu úmrtí v d·sledku této choroby v daném roce, lze podmín¥nou pravd¥podobnost P (Ai |Bi ) odhadnout takto
P (Ai |Bi ) =
di , mi
kde di je po£et osob ve v¥ku i, kterým byla v daném roce diagnostikována rakovina a mi je celkový po£et osob ve v¥ku i ºijících v daném roce. Odhad pravd¥podobnosti P (Bi ) získáme z úmrtnostních tabulek a ozna£íme jej pi . Pravd¥podobnost, ºe £lov¥k n¥kdy v pr·b¥hu onemocní, tedy odhadneme ze vztahu n X di P (A) = pi . mi i=0
A hledaný odhad lifetime risku je pak
LR =
n X di pi · 100%. m i i=0
Obdobn¥ odhadneme i lifetime risk vztahující se k úmrtí. A£koliv ÚZIS uvádí ve²kerá data pro p¥tileté v¥kové intervaly, lze odhad provést podle vý²e uvedeného vztahu a to tak, ºe pro kaºdé i z daného intervalu, pouºijeme pr·m¥r uvád¥né hodnoty. Tedy nap°íklad pokud pro i = 0, 1, 2, 3, 4 je po£et muº· roven m, pak m0 = m1 = · · · = m4 = m/5. Z údaj· o po£tu v roce 2008 nov¥ diagnostikovaných p°ípad· rakoviny prostaty a za zjednodu²ujícího p°edpokladu, ºe lze ztotoºnit po£et nov¥ diagnostikovaných p°ípadu s po£tem nov¥ vzniklých zhoubných novotvar·, jsem odhadla, ºe riziko onemocn¥ní rakovinou prostaty je 9,23 % a riziko úmrtí v d·sledku této choroby je 2,57 %. Pro karcinomy plic a tlustého st°eva jsou rizika výskytu 7,66 % a 4,69 % a rizika úmrtí 16
6,52 % a 2,34 %. Dále jsem se zabývala odhadem lifetime risku, ºe muº b¥hem ºivota onemocní rakovinou libovolného typu a také rizikem, ºe bude rakovina p°í£inou jeho smrti. Z t¥chto úvah jsem vynechala zhoubné novotvary k·ºe, které se svojí povahou (incidence, mortalita) dosti odli²ují od zbylých typ· zhoubných novotvar·. Výsledné hodnoty jsou 49,37 % pro výskyt a 26,35 % pro úmrtí. Tento zp·sob odhadu ov²em obsahuje n¥kolik zjednodu²ení, jako je nap°íklad pouºívání pr·m¥rných hodnot v rámci p¥tiletých interval·, ale zejména to, ºe ·daj o po£tu muº· v daném v¥ku v sob¥ zahrnuje i muºe nemocné, jimº byla choroba diagnostikována v minulých letech a se kterými bychom jiº po£ítat nem¥li. Pokud bychom byli schopni tyto jedince vy£lenit, získali bychom hodnoty je²t¥ o n¥co vy²²í.
Poznámka 1. Lifetime risk je ukazatel vhodný nap°íklad ke hrubému srovnávání
jednotlivných diagnóz z hlediska jejich rizikovosti. V ºádném p°ípad¥ ale nelze konkrétní osob¥ tvrdit, ºe pravd¥podobnost toho, ºe onemocní n¥jakou chorobou je rovna uvád¥nému lifetime risku. Skute£ná pravd¥podobnost, dá-li se v·bec ur£it, je totiº pro kaºdého individuální, závislá na jeho genetických p°edpokladech, ºivotním stylu a spoust¥ dal²ích faktor·, které p°i výpo£tu lifetime risku v·bec nejsou zohledn¥ny. Snad jen pro úplnost jeden p°íklad. Odhad pravd¥podobnosti, ºe muº onemocní rakovinou plic je 7,66 %, ale opravdu by m¥la být tato pravd¥podobnost stejná pro neku°áka ºijícího na samot¥ v horách a muºe, který kaºdý ve£er prosedí v hospod¥ a nehne se bez své oblíbené krabi£ky cigaret?
Poznámka 2. Druhým úskalím tohoto ukazatele je, ºe a£koliv se jeho hodnota
vztahuje k dlouhému budoucímu období a m¥la by vyjad°ovat riziko, ºe se v pr·b¥hu celého ºivota £lov¥ka n¥co stane, je odhadován z dat minulých. Hodnota odhadu lifetime risku se tak m¥ní v závislosti na tom, jaké hodnoty pro jeho odhad pouºijeme. Tento problém lze demonstrovat nap°íklad srovnáním odhadu rizika, ºe muº onemocní rakovinou prostaty, získaný z údaj· z roku 1998, který je 6,44 % a jiº zmín¥ného odhadu z roku 2008 9,23 %. Pokud bychom tedy muºi v roce 1998 °ekli, ºe riziko, ºe onemocní rakovinou prostaty je 6,44 % a ten se toho tak vyd¥sil, ºe okamºit¥ eliminoval v²echny rizikové faktory, jak bychom mu potom po deseti letech vysv¥tlili, ºe i p°es tato opat°ení se riziko nesníºilo, jak by z°ejm¥ o£ekával, ale naopak zvý²ilo? Op¥t tedy vidíme, ºe vztahovat hodnotu lifetime risku na konkrétní osobu není p°íli² vhodné. Na druhou stranu m·ºe být p°ínosné sledovat jeho vývoj v £ase. V²echny odhady pro rakovinu prostaty, plic a tlustého st°eva a v²echny karcinomy celkov¥ jsou uloºeny na p°iloºeném CD pod názvem lifetime_risk.xls.
17
3 3.1
Základní poznatky o kategoriálních datech Kategoriální prom¥nná
Pokud m¥°ená veli£ina nenabývá libovolné £íselné hodnoty, ale jen n¥které z hodnot z kone£né ²kály kategorií, jedná se o kategoriální prom¥nnou. S tímto typem prom¥nné se lze setkat ve studiích z mnoha r·zných obor·, ve kterých mohou popisovat nap°. názory voli£·, p°ítomnost sledovaného znaku, úrove¬ vzd¥lání nebo t°eba po£et £len· domácnosti. Kategoriální prom¥nné lze rozd¥lit podle n¥kolika kritérií a podle tohoto rozd¥lení s nimi pak dále pracovat. První moºností je rozd¥lení na prom¥nné vysv¥tlující a vysv¥tlované. Jak jejich ozna£ení napovídá, vysv¥tlovanou prom¥nnou je ta, jejíº chování nás primárn¥ zajímá. Toto chování by m¥la pomoci odhalit práv¥ prom¥nná vysv¥tlující.
P°íklad 1. Sledujeme-li závislost výsledku biopsie na hladin¥ PSA v krvi a po£tu
p°edchozích biopsií, pracujeme se dv¥ma kategoriálními prom¥nnými. Zatímco po£et p°edchozích biopsií nabývající hodnot {0, 1, 2, 3, 4, 5, 6} je prom¥nnou vysv¥tlující, výsledek biopsie, který je roven jedné v p°ípad¥ výskytu choroby a jinak je nulový, je prom¥nnou vysv¥tlovanou. Zanedbáme-li chyby zp·sobené nep°esností m¥°ení a zaokrouhlováním, lze hladinu PSA m¥°it spojit¥, o kategoriální prom¥nnou se tedy v tomto p°ípad¥ nejedná. Podle uspo°ádání kategorií, rozli²ujeme prom¥nné nominální, ordinální a intervalové. Pokud neexistuje p°irozené uspo°ádání kategorií, jedná se o prom¥nnou nominální a pro práci s ní lze pouºít jen ty metody, které povaºují v²echny kategorie za rovnocenné. V p°ípad¥ ordinálních prom¥nných jiº jednotlivé kategorie uspo°ádat lze a tak je p°i jejich zpracování moºné pouºít takové metody, které toto uspo°ádání vyuºijí. Pokud navíc známe vzdálenosti mezi kategoriemi, hovo°íme o prom¥nné intervalové.
P°íklad 2. Vysv¥tlovaná prom¥nná z p°edchozího p°íkladu je prom¥nnou nomi-
nální, zatímco po£et p°edchozích biopsií prom¥nnou intervalovou. Pokud bychom navíc mezi vysv¥tlující prom¥nné za°adili mnoºství za den vykou°ených cigaret, m¥°ené na ²kále {ºádná, málo, hodn¥}, jednalo by se o prom¥nnou ordinální. Kategoriální prom¥nné lze rozd¥lit také podle toho, zda jsou to prom¥nné
kvantitativní nebo kvalitativní. Zatímco v²echny nominální prom¥nné jsou zárove¬ kvalitativní a podobn¥ intervalové jsou kvantitativní, u prom¥nných ordinálních toto rozd¥lení tak jednozna£né není a lze na n¥ pohlíºet jak jako na kvalitativní, tak jako na kvantitativní. Pokud jsou totiº jednotlivé kategorie popsány slovn¥ a mají tedy kvalitativní podobu, m·ºeme jim v souladu s uspo°ádáním p°i°adit skóre a získat tak prom¥nnou kvantitativní, se kterou se v modelech £asto snáze pracuje. 18
P°íklad 3. Pokud op¥t naváºeme na p°edcházející p°íklad, vyskytují se v n¥m
prom¥nné obou typ·. Výskyt choroby je veli£inou kvalitativní a hodnoty {0, 1}, pomocí nichº je popsána mají pouze zástupný charakter. Druhou kvalitativní veli£inou je slovn¥ vyjád°ené mnoºství cigaret. Pokud bychom ov²em t¥mto slovn¥ rozli²eným kategoriím p°i°adili hodnoty {0, 1, 2}, stalo by se veli£inou kvantitativní. Nadále v²ak toto mnoºství z·stává pouze ordinální veli£inou, nikoli intervalovou, protoºe p°esné po£ty denn¥ vykou°ených cigaret odpovídající jednotlivým kategoriím nejsou ur£eny. Poslední, tentokrát jiº intervalová kategoriální prom¥nná udávající po£et dosavadních biopsií, je op¥t prom¥nnou kvantitativní. Vzhledem k cíli práce, kterým je popsat vlivy nejr·zn¥j²ích veli£in, které mohou, ale nemusejí mít kategoriální podobu, na výsledek biopsie, kterým je bu¤ nula, pokud nebyla odhalena rakovina prostaty nebo jedni£ka, pokud ano, se budu dále v¥novat zejména problematice model· s alternativn¥ rozd¥lenou nominální vysv¥tlovanou prom¥nnou. Pro ucelenost textu ov²em uvedu i základní informace o binomickém a multinomickém rozd¥lení pravd¥podobnosti.
3.2
Binomické rozd¥lení
Náhodná veli£ina X má binomické rozd¥lení s parametry n ∈ N a π ∈ (0, 1) (X ∼ Bi(n, π)), pokud je pravd¥podobnost, ºe bude rovna k = 0, 1, . . . , n, dána vztahem n k P (X = k) = π (1 − π)n−k , k kde k m·ºe vyjad°ovat nap°íklad po£et úsp¥ch· v n nezávislých pokusech. Pro st°ední hodnotu rozd¥lení platí
E(X) =
n P k=0
=
n P k=1
= nπ
n P
kP (X = k) =
k=1 n! π k (1 (n−k)!(k−1)!
n−1 P j=0
n−1 j
k
n k
π k (1 − π)n−k =
n P k=1
− π)n−k = nπ
n P k=1
n−1 k−1
n! π k (1 − π)n−k = k (n−k)!k!
π k−1 (1 − π)n−k =
π j (1 − π)n−1−j = nπ(π + 1 − π)n−1 = nπ.
Obdobn¥ lze odvodit, ºe E[X(X − 1)] = n(n − 1)π 2 , coº vyuºijeme p°i výpo£tu rozptylu, pro který platí var(X) = E(X 2 ) − [E(X)]2 = E[X(X − 1)] + E(X) − [E(X)]2 =
= n(n − 1)π 2 + nπ − n2 π 2 = n2 π 2 − nπ 2 + nπ − n2 π 2 = nπ(1 − π). 19
Speciálním p°ípadem binomického rozd¥lení je pro n = 1 rozd¥lení alternativní. Náhodná veli£ina X ∼ Alt(π) v tomto p°ípad¥ m·ºe nabývat pouze dvou hodnot a to jedni£ky v p°ípad¥ usp¥chu (ten v na²em p°ípad¥ znamená pozitivním výsledek biopsie) a nuly v p°ípad¥ neúsp¥chu (negativní výsledek biopsie). Parametr π zna£í pravd¥podobnost úsp¥chu a lze tedy psát
P (X = 1) = π
P (X = 0) = 1 − π.
a
íslené charakteristiky alternativního rozd¥lení jsou
E(X) = 0 · (1 − π) + 1 · π = π , var(X) = E(X 2 ) − [E(X)]2 = 0 · (1 − π) + 1 · π − π 2 = π(1 − π). Pro alternativní rozd¥lení dále platí, ºe jsou-li X1 , . . . , Xn nezávislé náhodné P veli£iny, v²echny s rozd¥lením Alt(π), potom pro náhodnou veli£inu Y = ni=1 Xi platí Y ∼ Bi(n, π). Této vlastnosti lze vyuºít p°i výpo£tu charakteristik binomického rozd¥lení, které byly uvedeny v úvodu kapitoly
E(Y ) = E (
Pn
i=1 Xi ) =
Pn
i=1 E(Xi ) =
Pn
i=1
π = nπ ,
P P P var(Y ) = var ( ni=1 Xi ) = ni=1 var(Xi ) = ni=1 π(1 − π) = nπ(1 − π).
3.2.1 Odhad parametru π Pokud neznáme skute£nou hodnotu parametru binomického rozd¥lení π , je t°eba jej odhadnout z dat, která nesou informaci o po£tu úsp¥ch· y v n pokusech. P°edpis pro tento odhad získáme pomocí metody maximáln¥ v¥rohodného odhadu. Jak jiº název napovídá, touto metodou hledáme maximum tzv. v¥rohodnostní funkce, kterou v p°ípad¥ binomického parametru získáme ze vztahu n y P (Y = y) = π (1 − π)n−y . y Jelikoº £len ny nezávisí na hledaném parametru π , zredukuje se v¥rohodnostní funkce na tvar l(π) = π y (1 − π)n−y . Výsledek se nezm¥ní pokud místo v¥rohodnostní funkce budeme uvaºovat její logaritmus, se kterým se snáze pracuje:
L(π) = ln[l(π)] = ln[π y (1 − π)n−y ] = y ln π + (n − y) ln(1 − π) 20
Dále hledáme maximum b¥ºným zp·sobem, provedeme tedy první derivaci L(π) a poloºíme ji rovnu nule.
y n−y y − nπ ! δL = − = = 0. δπ π 1−π π(1 − π) Odtud získáme maximáln¥ v¥rohodný odhad parametru π ve tvaru
π ˆ=
y , n
který je zárove¬ odhadem nestranným, nebo´
E(ˆ π) = E
y n
=
1 1 E(y) = nπ = π. n n
Kovarian£ní matici maximáln¥ v¥rohodného odhadu vektorového parametru β obecn¥ získáme jako inverzi z informa£ní matice, jejíº prvek na pozici (j, k) je dán vztahem 2 δ L(β) . −E δβj δβk Jelikoº ale odhadujeme jen jednorozm¥rný parametr π , zjednodu²í se tento vztah do tvaru 2 n−y n−y = −E − πy2 − (1−π) −E δδ2Lπ = E πy2 + (1−π) = 2 2
=
1 E(y) π2
+
1 π
1 1−π
= n
+
n (1−π)2
−
1 E(y) (1−π)2
1−π+π = = n π(1−π)
=
nπ π2
+
n−nπ (1−π)2
=
n . π(1−π)
Odtud jiº získáme vztah pro rozptyl odhadu
n var(ˆ π) = π(1 − π)
−1 =
π(1 − π) . n
Takto jsme odvodili asymptotickou hodnotu rozptylu odhadu, tu lze v p°ípad¥ binomického rozd¥lení ur£it i p°esn¥ var(ˆ π ) = var
y n
=
1 1 π(1 − π) var(y) = 2 nπ(1 − π) = . 2 n n n
3.2.2 Testy odhadu πˆ K testování hypotézy H0 : π = π0 lze p°istupovat t°emi zp·soby, které v²echny vycházejí z toho, ºe maximáln¥ v¥rohodné odhady testovaných parametr· mají 21
pro velký rozsah souboru normální rozd¥lení, viz V¥ty 7.91 a 7.92 v [17], jedná se tedy o testy asymptotické. První moºností je Waldova statistika, pro kterou za platnosti H0 platí
zw =
π ˆ − π0 as. π ˆ − π0 =q ∼ N (0, 1). SE π ˆ (1−ˆ π) n
Druhou metodou je skórový test, zaloºený na sklonu a o£ekávaném zak°ivení L(π) v bod¥ π0 . Ozna£me δL(π) y − nπ0 u(π0 ) = = δπ π=π0 π0 (1 − π0 ) a
ι(π0 ) = −E
δ2L δ2π
= π=π0
n . π0 (1 − π0 )
Statistika skórového testu má pak tvar
u(π0 ) y − nπ0 π ˆ − π0 zs = p =p =q π0 (1−π0 ) ι(π0 ) nπ0 (1 − π0 ) n
a op¥t platí, ºe za platnosti H0 má tato statistika asymptoticky normované normální rozd¥lení. Rozdíl mezi ob¥ma metodami je pouze v bod¥ vy£íslení sm¥rodatné odchylky, kdy se jako lep²í jeví skórový test, který vy£ísluje v bod¥ π0 a u n¥hoº je rozd¥lení výb¥rové hodnoty testové statistiky zs blíºe N(0, 1), neº v p°ípad¥ výb¥rové hodnoty zw . Poslední metodou je test zaloºený na v¥rohodnostním pom¥ru. Ozna£me l0 maximum v¥rohodnostní funkce za platnosti H0 a l1 maximum bez ohledu na hypotézu, pak l0 Λ= ≤1 l1 a testová statistika má obecn¥ tvar
−2 ln Λ = −2 ln
l0 = −2(L0 − L1 ). l1
V p°ípad¥ binomického parametru π mají L0 a L1 tvar
L0 = y ln π0 + (n − y) ln(1 − π0 ),
L1 = y ln π ˆ + (n − y) ln(1 − π ˆ ).
22
Testovou statistiku pak lze zapsat takto 1−ˆ π −2(L0 − L1 ) = 2 y ln ππˆ0 + (n − y) ln 1−π 0
n−y . = 2 y ln nπy 0 + (n − y) ln n−nπ 0 Porovnáváme tedy vºdy po£ty pozorovaných úsp¥ch· resp. neúsp¥ch· s po£ty o£ekávanými za platnosti nulové hypotézy. Tato statistika má za platnosti H0 asymtoticky χ2 rozd¥lení o jednom stupni volnosti, který zna£í rozdíl mezi po£tem hledaných parametr· bez p°edpokladu na platnost H0 a za její platnosti.
3.2.3 Intervaly spolehlivosti Vyuºijeme-li asymptotických vlastností maximáln¥ v¥rohodných odhad·, m·ºeme ur£it i jejich kon den£ní intervaly. První metodou je úprava statistiky Waldova testu do podoby r π ˆ (1 − π ˆ) , π ˆ ± u1− α2 n kde u1− α2 je 1 − α2 kvantil rozd¥lení N(0, 1). Takto získaný odhad ov²em ²patn¥ pokrývá skute£nou hodnotu π v p°ípad¥, kdy je blízko krajních hodnot 0 nebo 1 a mohlo by se stát, ºe interval spolehlivosti bude obsahovat i hodnoty mimo interval h0, 1i. Navíc, abychom mohli pouºít aproximaci normálním rozd¥lením, vyºaduje pom¥rn¥ velký rozsah výb¥ru n. Ukazuje se, ºe lep²í chování má interval získaný pomocí skórového testu, jehoº krajní body najdeme vy°e²ením rovnic
π ˆ − π0 q
π0 (1−π0 ) n
a
π ˆ − π0 q
π0 (1−π0 ) n
= u1− α2
= −u1− α2 ,
kde u1− α2 je op¥t 1 − α2 kvantil rozd¥lení N(0, 1), který budeme pro jednoduchost zna£it jen u. Tyto rovnice lze zjednodu²it do tvaru q 2 u π ˆ n(1 − π ˆ ) + 14 u2 2ˆ πn + u π0 = ± 2(n + u2 ) n + u2 a meze intervalu tak p°ímo spo£ítat dosazením. Stejné meze pro volbu α = 0, 05 získáme i p°íkazem prop.test(x,n,conf.level=0.95,correct=F) v R. 23
Interval spolehlivosti lze získat také pomocí testu zaloºeného na v¥rohodnostním pom¥ru. Tentokrát je dán jako mnoºina takových π0 , pro která je P-value testu vet²í neº dané α. V R lze krajní hodnoty vypo£ítat p°íkazem uniroot(f, interval, ...,), viz následující p°íklad.
P°íklad 4. Chceme odhadnout jaká je pravd¥podobnost, ºe výsledek biopsie bude
pozitivní. Biopsii podstoupilo 987 muº· ve v¥ku od 45 do 80 let, kterým byla nam¥°ena hladina PSA v krvi men²í neº 20 ng/ml a jejichº objem prostaty byl men²í neº 150 cm3 . Výsledk· potvrzujících výskyt zhoubného nádoru prostaty bylo 230. Odhad pravd¥podobnosti: π ˆ=
230 = 0, 233 987
Odhad rozptylu: var c (ˆ π) =
0, 233(1 − 0, 233) = 1, 81 · 10−4 987
95% interval spolehlivosti: 1. Wald·v test
p 0, 233 ± 1, 96 1, 81 · 10−4 ,
kde 1,96 je 0, 975-kvantil N(0, 1). Výsledný kon den£ní interval je h0, 207; 0, 259i.
2. Skórový test p 2 · 0, 233 · 987 + 1, 962 1, 96 0, 233 · 987(1 − 0, 233) + 0, 25 · 1, 962 ± 2(987 + 1, 962 ) 987 + 1, 962
takto získáme interval h0, 208; 0, 260i.
3. Test zaloºený na v¥rohodnostním pom¥ru Hledáme °e²ení nerovnice −2[230 ln π + 757 ln(1 − π) − 230 ln 0, 23 − 757 ln(1 − 0, 23)] ≤ 3, 84,
kde 3, 84 je 0, 95 kvantil χ2 rozd¥lení o jednom stupni volnosti. Takto získáme interval h0, 207; 0, 260i.
e²ení tohoto p°íkladu je uloºeno na p°iloºeném CD pod názvem priklad_4.R. 24
3.3
Multinomické rozd¥lení
Jedná se o roz²í°ení binomického rozd¥lení na p°ípad, kdy má sledovaná prom¥nná více neº dv¥ kategorie. M¥jme náhodnou veli£inu Yij , která je rovna 1, pokud byl výsledek i-tého pokusu z j-té kategorie a 0 jinak. Pravd¥podobnost, ºe Yij = 1 ozna£me πj . Celý i-tý pokus m·ºeme zapsat do vektoru
Y i = (Yi1 , . . . , Yic ),
c X
Yij = 1,
j=1
kde c je Pcelkový po£et kategorií sledované prom¥nné. Vzhledem k tomu, ºe je sou£et j Yij pevn¥ daný, sta£í znát jen c−1 £len· vektoru a zbylý z nich dopo£ítat. Zopakujeme-li pokus n-krát, je celkový po£et pozorování v j-té kategorii dán vztahem n X Nj = Yij i=1
a vektor N = (N1 , . . . , Nc ) má multinomické rozd¥lení s parametry n, π1 , . . . , πc a jeho pravd¥podobnostní funkce je dána vztahem
P (N = n) = P (N1 = n1 , . . . , Nc = nc ) = pro
nj = 1, 2, . . . , n
∀j,
n! π n1 · · · πcnc n1 ! · · · nc ! 1
c X
nj = n.
j=1
Protoºe celkové po£ty pozorování v jednotlivých kategoriích mají binomické rozd¥lení s parametry n a πj , platí
E(Nj ) = nπj
a
var(Nj ) = nπj (1 − πj ).
íselné charakteristiky náhodného vektoru N odvodíme následovn¥:
E(N ) = (E(N1 ), . . . , E(Nc )) = (nπ1 , . . . , nπc ) Pro výpo£et varian£ní matice nejprve pouºijeme veli£inu Yij , pro kterou platí
• E(Yij ) = 1 · πj + 0 · (1 − πj ) = πj • E(Yij2 ) = 12 · πj + 02 · (1 − πj ) = πj • E(Yij Yik ) = 0 pro ∀j 6= k , protoºe v kaºdém pokusu je jen jedna z veli£in Yij rovna jedné a zbylé jsou nulové • var(Yij ) = E(Yij2 ) − [E(Yij )]2 = πj − πj2 25
• cov(Yij , Yik ) = E(Yij Yik ) − E(Yij )E(Yik ) = 0 − πj πk = −πj πk pro ∀j 6= k . A protoºe Nj je sou£et hodnot Yij , i = 1, . . . , n ze v²ech n nezávislých pokus·, je jeho rozptyl ! n n X X var(Nj ) = var Yij = var(Yij ) = n(πj − πj2 ) ∀j = 1, . . . , c i=1
i=1
a pro kovarianci platí n X
cov(Nj , Nk ) = cov
i=1
Yij ,
n X
! Yik
i=1
=
n X
cov(Yij , Yik ) = n(−πj πk ) = −nπj πk .
i=1
Varian£ní matice vektoru N = (N1 , . . . , Nc ), který pochází z multinomického rozd¥lení je potom nπ1 (1 − π1 ) −nπ1 π2 ··· −nπ1 πc .. .. . . −nπ2 π1 var(N ) = .. .. . −nπc−1 πc . −nπc π1 ··· −nπc πc−1 nπc (1 − πc ) coº lze vektorov¥ zapsat jako var(N ) = n[diag(π) − ππ 0 ], kde π = (π1 , . . . , πc ).
3.3.1 Odhad parametru πj P°edpis pro odhad multinomického parametru op¥t získáme metodou maximální v¥rohodnosti. V¥rohodnostní funkce má v tomto p°ípad¥ tvar
l(π) =
c Y
n πj j ,
πj ≥ 0,
j = 1, 2, . . . , c
j=1
c X
πj = 1
j=1
kde π tentokrát zna£í jen c−1 sloºkový vektor (π1 , π2 , . . . , πc−1 ) a πc = 1− Op¥t se bude snáze hledat maximum logaritmu v¥rohodnostní funkce
L(π) =
c X
Pc−1 i=1
πi .
nj ln πj .
j=1
Tuto funkci proto zderivujeme podle π a poloºíme rovnu nule. Nejprve je ale pot°eba ur£it
δπc δ(1 − π1 − · · · − πc−1 ) = = −1 δπj δπj 26
∀j = 1, . . . , c − 1
a
1 δπc 1 δ ln πc = =− δπj πc δπj πc
∀j = 1, . . . , c − 1.
Nyní jiº m·ºeme ur£it derivaci logaritmu v¥rohodnostní funkce
δL(π) δ(n1 ln π1 + · · · + nc ln πc ) nj nc ! = = − =0 δπj δπj πj πc
∀j = 1, . . . , c − 1.
Odtud
π ˆj nj = π ˆc nc Pc toto dosadíme do j=1 πj = 1 c X π ˆc nj j=1
nc
∀j = 1, . . . , c − 1,
c π ˆc π ˆc X nj = n = 1, = nc j−1 nc
odkud jiº získáme maximáln¥ v¥rohodný odhad multinomického parametru ve tvaru nc ∀j = 1, . . . , c. π ˆj = n Vektorov¥ lze odhad zapsat takto
π ˆ=
N . n
Nejprve ov¥°íme nestrannost odhadu:
E(ˆ πj ) = E
n j
n
1 = E n
n X
! yij
i=1
n
1X 1 = E(yij ) = nπj = πj . n i=1 n
A pro varian£ní matici platí N 1 1 1 var(ˆ π ) = var = 2 var(N ) = 2 n[diag(π) − ππ 0 ] = [diag(π) − ππ 0 ]. n n n n
3.3.2 Testy multinomického parametru První moºností jak testovat hypotézu
H0 :
πj = πj0
j = 1, 2, . . . , c
c X j=1
je Pearsonova statistika, která má tvar c X (nj − nπj0 )2 X = . nπ j0 j=1 2
27
πj0 = 1
Pokud navíc ozna£íme o£ekávané frekvence jako
µj = nπj0 , zjednodu²í se statistika do podoby c X (nj − µj )2 X = . µ j j=1 2
Za platnosti hypotézy H0 má tato statistika pro velké soubory p°ibliºn¥ χ2 rozd¥lení o c − 1 stupních volnosti. Hypotézu zamítáme v p°ípad¥, kdy platí X 2 ≥ χ2c−1 (1 − α), tedy 1 − α kvantil χ2 rozd¥lení o c − 1 stupních volnosti. Z vý²e uvedeného vztahu lze také snadno ur£it, která kategorie k tomuto zamítnutí nejvíce p°isp¥la. Pro snaz²í výpo£et lze ale tento vztah je²t¥ upravit
X2 =
c P j=1
=
c P j=1
(nj −µj )2 µj
n2j µj
=
c P j=1
n2j −2nj µj +µ2j µj
c P
− 2n + n =
j=1
n2j µj
c P
=
j=1
n2j µj
−2
c P j=1
nj +
c P
µj =
j=1
− n.
Test pomocí Pearsonovy statistiky lze v R provést p°íkazem chisq.test(N,π0 ). Druhou moºností je test zaloºený na v¥rohodnostním pom¥ru. Jádro v¥rohodnostní funkce je Y n l(π) = πj j . j
Jeho maximum je za platnosti H0
l0 =
Y
n
πj0j
j
a obecn¥
l1 =
Y nj nj n
j
.
Pom¥r t¥chto dvou v¥rohodností je
Q nj l0 j πj0 Λ = = Q nj nj l1 j n a hledaná statistika pak
Q 2
G = −2 ln Λ = −2 ln Q
j
n
j
X X πj0j nj nj = 2 n ln = 2 nj ln . j n nj j nπj0 µj n j j 28
I tato statistika má pro velká n χ2 rozd¥lení o c − 1 stupních volnosti a hypotézu H0 zamítáme v p°ípad¥ kdy X 2 ≥ χ2c−1 (1 − α). Statistiky X 2 a G2 tedy pocházejí ze stejného rozd¥lení a pro velké soubory jsou jejich hodnoty za platnosti H0 tém¥° totoºné. Pearsonova statistika ov²em k teoretickému rozd¥lení konverguje rychleji. Pro správné fungování statistiky G2 je navíc pot°eba zajistit, aby nc > 5.
P°íklad 5. Op¥t se budeme zabývat výsledky biopsií prostaty u muº· ve v¥ku
od 45 do 80 let, jejichº hladina PSA byla v rozmezí 0 aº 20 ng/ml a zm¥°ený objem prostaty byl men²í neº 150 ml. Tentokrát ov²em sledujeme £ty°i choroby a to karcinom, hyperplazii, dysplazii a chronický zán¥t prostaty. Pacient· s jednou z t¥chto £ty° chorob bylo celkem 870, po£ty odpovídající jednotlivým onemocn¥ním, jsou uvedeny v tabulce. Onemocn¥ní Karcinom Hyperplazie Dysplazie Chronický zán¥t Celkem Ozna£ení n1 n2 n3 n4 n Po£et 230 379 91 170 870 Pravd¥podobnost diagnostikování jednotlivých výsledk· vy²et°ení odhadneme pro kaºdou chorobu ze vztahu n π ˆj =
j
n
∀j,
takto zíkáme odhady π ˆ1 =
379 91 170 230 = 0, 26 π ˆ2 = = 0, 44 π ˆ3 = = 0, 1 π ˆ4 = = 0, 2, 870 870 870 870
tedy π ˆ = (0, 26; 0, 44; 0, 1; 0, 2).
Dále je t°eba odhadnout varian£ní matici tohoto odhadu. var(ˆ c π) =
1 (diag(ˆ π) n
−π ˆπ ˆ0) =
0, 26 0 0 0 0, 26 0 0, 44 0 0 0, 44 1 = 870 − 0 0 0, 1 0 0, 1 0 0 0 0, 2 0, 2
(0, 26; 0, 44; 0, 1; 0, 2)
2.24 · 10−4 −1.32 · 10−4 −3.18 · 10−5 −5.94 · 10−5 −1.32 · 10−4 2.83 · 10−4 −5.24 · 10−5 −9.78 · 10−5 = −3.18 · 10−5 −5.24 · 10−5 1.08 · 10−4 −2.35 · 10−5 −5.94 · 10−5 −9.78 · 10−5 −2.35 · 10−5 1.81 · 10−4
Nyní p°edpokládejme, ºe v²echny choroby mohou být diagnostikovány se stejnou pravd¥podobností, tedy H0 :
πj0 = 29
1 4
∀j.
Nejprve ur£íme o£ekávané frekvence a hypotézu následn¥ otestujeme pomocí Pearsonovy statistiky: µj = nπj0 = X2 =
4 P j=1
(nj −µj )2 µj
=
(230−217,5)2 217,5
870 = 217, 5 ∀j 4
+
(379−217,5)2 217,5
+
(91−217,5)2 217,5
+
(170−217,5)2 217,5
=
= 0, 72 + 119, 92 + 73, 57 + 10, 37 = 204, 58
Hypotézu H0 zámítáme, protoºe 204, 58 > χ23 (0, 95) = 7, 815. Z výpo£tu je navíc vid¥t, ºe nejvíc k tomuto zamítnutí p°isp¥l velký po£et diagnostikovaných hyperplazií a také malý po£et dysplazií. Záv¥r m·ºeme ov¥°it testem pom¥ru v¥rohodností, hodnota jeho statistiky je 2
G
= 2
4 P j=1
nj ln
nj µj
= 2 230 ln
230 217,5
+ 379 ln
379 217,5
+ 91 ln
91 217,5
+ 170 ln
170 217,5
=
= 2(25, 71 + 420, 95 − 158, 58 − 83, 78) = 204, 29
Protoºe jsou jsme m¥li k dispozici velký po£et pozorování a sledovali jsme jen £ty°i kategorie, jsou hodnoty obou statistik prakticky shodné a podmínka nc > 5 je zjevn¥ spln¥na. Op¥t tedy zamítáme hypotézu, ºe jsou v²echny výsledky biopsie stejn¥ pravd¥podobné. Ilustrovali jsme tedy, ºe oba testy vedly ke stejnému záv¥ru, v praxi pak v¥t²inou sta£í pouºít pouze jeden z nich. e²ení toho p°íkladu je uloºeno pod názvem priklad_5.R.
30
4
Kontingen£ní tabulky
Jednou z moºností jak zkoumat vztah mezi dv¥ma kategoriálními prom¥nnými je pouºití kontingen£ní tabulky. Tato tabulka vznikne, pozorujeme-li v n nezávislých pokusech hodnoty prom¥nné X , která má I kategorií a prom¥nné Y o J kategoriích. Na pozici (i, j) v kontingen£ní tabulce potom zapí²eme hodnotu nij , která je rovna po£tu pozorování, pro neº platilo X = i a Y = j . Výsledná tabulka bude vypadat takto X\Y 1 2 .. . I P
1 n11 n21 .. .
2 n12 n22 .. .
··· ··· ···
J n1J n2J .. .
P
nI1 n.1
nI2 n.2
··· ···
nIJ n.J
nI. n
n1. n2. .. .
V tabulce se také vyskytují hodnoty ni. resp. n.j , které ozna£ují °ádkové resp. sloupcové sou£ty. Jsou tedy dány vztahem X X ni. = nij n.j = nij . j
i
Zárove¬ musejí v²echny hodnoty v tabulce spl¬ovat X X XX ni. = n.j = nij = n. i
j
i
j
Vektor (n11 , . . . , n1J , . . . , nI1 , . . . , nIJ ) má multinomické rozd¥lení s parametry n a π = (π11 , . . . , π1J , . . . , πI1 , . . . , πIJ ). Op¥t m·ºeme tyto pravd¥podobnosti zapsat do tabulky X\Y 1 2 .. . I P
1 π11 π21 .. .
2 π12 π22 .. .
··· ··· ···
J π1J π2J .. .
P
πI1 π.1
πI2 π.2
··· ···
πIJ π.J
πI. 1
π1. π2. .. .
kde
πij = P (X = i, Y = j)
πi. = P (X = i)
i = 1, 2, . . . , I
j = 1, 2, . . . , J. 31
π.j = P (Y = j)
Pokud ov²em skute£né hodnoty paramatr· πij neznáme, je pot°eba je odhadnout ze vztahu ni. n.j nij π ˆi. = π ˆ.j = π ˆij = n n n pro i = 1, 2, . . . , I j = 1, 2, . . . , J. Nej£ast¥j²í otázkou, na kterou se za pomoci kontingen£ních tabulek snaºíme odpov¥d¥t, je nezávislost prom¥nných X a Y . Pokud X a Y skute£n¥ jsou nezávislé, pak platí πij = P (X = i, Y = j) = P (X = i)P (Y = j) = πi. π.j
∀i = 1, 2, . . . , I
∀j = 1, 2, . . . , J.
V tomto p°ípad¥ tedy sta£í odhadnout jen π ˆi. , i = 1, 2, . . . , I − 1 a π ˆ.j , j = 1, 2, . . . , J − 1 a zbylé parametry z nich dopo£ítat podle vztah·
π ˆI. = 1 −
I−1 X
π ˆi.
π ˆ.J = 1 −
i=1
J−1 X
π ˆ.j
π ˆij = π ˆi. π ˆ.j .
j=1
Pro testování hypotézy H0 , ºe X a Y jsou nezávislé, lze pouºít Pearsonovu χ2 statistiku I X J n n 2 X nij − i.n .j 2 X = , ni. n.j i=1 j=1
n
která má za platnosti H0 asymptoticky χ2 rozd¥lení o (I − 1)(J − 1) stupních volnosti. Pokud navíc o£ekávanou frekvenci ozna£íme ni. n.j , µ ˆij = n zjednodu²í se vztah do podoby I X J X (nij − µ ˆij )2 X = . µ ˆij i=1 j=1 2
Hypotézu o nezávislosti tedy zamítneme, pokud X 2 ≥ χ2(I−1)(J−1) (1 − α). Pro správné fungování tohoto testu je pot°eba zajistit dostate£ný po£et pozorování a také spln¥ní podmínky µ ˆij ≥ 5, ∀i, j . Tento test lze v R provést pomocí p°íkazu chisq.test() viz P°íklad 6. Dal²í moºností, jak lze v kontingen£ních tabulkách testovat nezávislost, je statistika G2 , která je výsledkem testu pom¥ru v¥rohodností pro multinomický parametr. Jádro v¥rohodnostní funkce je obecn¥ XX YY n l(π) = πijij ∀πij ≥ 0 πij = 1. i
j
i
32
j
Za platnosti H0
ni. n.j n2
π ˆij = a tedy
l0 =
Y Y (ni. n.j )nij i
n2nij
j
.
A bez p°edpokladu platnosti H0
nij n
π ˆij = a
Y Y nnijij
l1 =
i
nnij
j
.
Pom¥r v¥rohodnostních funkcí je
l0 Λ= = l1 Ve výpo£tu jsme vyuºili rovnosti
(ni. n.j)nij n2nij Q Q nnijij i j nnij
Q Q i
P P i
j
j
nij = n. Statistika G2 má tvar
n (ni. n.j) ij 2nij i j n nij Q Q nij i j nnij
Q Q
G2 = −2 ln Λ = −2 ln
.
nij µ ˆ ij i j nnij nij Q Q nij i j nnij
Q Q
= −2 ln
=
h Q Q Q Q nij i nij ˆij − ln i j nij = = −2 ln i j µ = −2
hP P
= −2
P P
= 2
i
i
P P i
j
j
nij ln µ ˆij −
P P i
i
j
nij ln nij =
nij (ln µ ˆij − ln nij ) = −2
P P i
µ ˆ
j
nij ln nijij =
n
j
nij ln µˆijij .
I tato statistika má asymptoticky χ2 rozd¥lení o (I-1)(J-1) stupních volnosti a hypotézu zamítáme v p°ípad¥, kdy je G2 ≥ χ2(I−1)(J−1) (1 − α). Pro její dobré fungování je op¥t t°eba zajistit dostate£ný po£et pozorování a také dodrºet nerovnost n ≥ 5. IJ Statistiky X 2 a G2 srovnáváme se stejnými hodnotami. Zárove¬ pro n¥ platí, ºe rozdíl (X 2 − G2 ) konverguje podle pravd¥podobnosti k nule. Metody tedy s rostoucím po£tem pozorování mají velmi podobné výsledky, o n¥co rychlej²í konvergenci k teoretickému rozd¥lení v²ak vykazuje X 2 . 33
P°íklad 6. Cílem je ov¥°it tvrzení uvedené v kapitole 2.5.1, ºe existuje závislost
mezi vý²í PSA indexu a výskytem zhoubného nádoru prostaty. K dispozici jsou údaje o PSA indexu a výsledku biopsie od 511 pacient·, kterým byla nam¥°ena hodnota PSA z intervalu h4; 10i ng/ml a jejichº objem prostaty byl men²í neº 75 cm3 . Diag.\Index v % <15 h15; 30i >30 Jiná 154 201 30 385 Karcinom 76 44 6 126 230 245 36 511 Odhad pravd¥podobností bez p°edpokladu nezávislosti Diag.\Index v % <15 h15; 30i Jiná 0,30 0,39 Karcinom 0,15 0,09 0,45 0,48
>30 0,06 0,75 0,01 0,25 0,07 1
Odhad pravd¥podobností za p°edpokladu nezávislosti Diag.\Index v % <15 h15; 30i Jiná 0,34 0,36 Karcinom 0,11 0,12 0,45 0,48
>30 0,05 0,75 0,02 0,25 0,07 1
O£ekáváné po£ty za p°edpokladu nezávislosti Diag.\Index v % <15 h15; 30i >30 Jiná 173,29 184,59 27,12 385 Karcinom 56,71 60,41 8,88 126 230 245 36 511 Hodnota statistiky X 2 X2 =
(154 − 173, 29)2 (6 − 8, 88)2 + ··· + = 15, 86 173, 29 8, 88
Hodnota statistiky G2 154 6 G = 2 154 ln + · · · + 6 ln = 15, 84 173, 29 8, 88 2
Protoºe je rozsah výb¥ru dosti velký jsou hodnoty statistik X 2 a G2 tém¥° totoºné a ob¥ vedou k zamítnutí hypotézy o nezávislosti PSA indexu a výsledku biopsie, protoºe jsou vy²²í neº χ22 (0, 95) = 5, 991. Z dat je také vid¥t, ºe po£et diagnostikovaných karcinom· je nep°ímo úm¥rný vý²i PSA indexu, coº potvrzuje 34
tvrzení kapitoly 2.5.1, ºe hodnoty PSA indexu niº²í neº stanovená hranice (20 %25 %) indikují zvý²enou pravd¥podobnost výskytu zhoubného nádoru prostaty. e²ení tohoto p°íkladu je na CD uloºeno pod názvem priklad_6. Dále je v kapitole 2.5.1 zmín¥no, ºe pokud má pacient objem prostaty vy²²í neº 75 ml, není hodnota PSA indexu pro predikování výsledku biopsie sm¥rodatná. I toto tvrzení lze ov¥°it obdobn¥ jako bylo uvedeno vý²e, tentokrát ov²em uvaºujeme jen 77 pacient· se zv¥t²enou prostatou, tedy objemem v¥t²ím neº 75 ml. Kontingen£ní tabulka má tentokrát tuto podobu Diag.\Index v % <15 h15; 30i >30 Jiná 11 45 14 70 Karcinom 2 3 2 7 13 48 16 77 a o£ekávané po£ty jsou Diag.\Index v % <15 h15; 30i >30 Jiná 11,82 43,64 14,55 70 Karcinom 1,18 4,36 1,45 7 13 48 16 77 Z této tabulky je vid¥t, ºe není spln¥na podmínka pro správné pouºití statistiky X 2 aby o£ekávaný po£et pozorování v kaºdé bu¬ce byl aspo¬ 5, a ke spln¥ní této podmínky nepom·ºe ani sníºení kategorií pro PSA index na dv¥. Statistiku G2 ale pouºít lze. 2 11 + · · · + 2 ln = 1, 25 G = 2 11 ln 11, 82 1, 45 2
A protoºe 1, 25 ≤ χ22 (0, 95) = 5, 991 nem·ºeme zamítnout hypotézu, ºe výsledek biopsie nezávisí na hodnot¥ PSA indexu a v tomto p°ípad¥ tedy nem·ºeme najít takovou hodnotu PSA indexu, jejíº p°ekro£ení by s dostate£nou jistotou vylou£ilo p°ítomnost karcinomu prostaty v t¥le pacienta. Tato druhá £ást p°íkladu je uloºena pod názvem priklad_6b.R. 4.1
Vztahy mezi prom¥nnými v tabulce 2x2
Nyní se budeme zabývat tabulkami, zachycujícimi sledování dvou prom¥nných, které mají vºdy jen dv¥ kategorie. ádky budou p°estavovat kategorie vysv¥tlující prom¥nné X . Sloupce pak zna£í kategorie vysv¥tlované prom¥nné Y , kterými 35
budou nej£ast¥ji úsp¥ch (zna£íme 1) a neúsp¥ch (zna£íme 0). Pravd¥podobnosti v tabulce 2x2 zna£íme P X\Y 1 0 1 π11 π12 π1. 2 π21 π22 π2. P π.1 π.2 1 Nebo m·ºeme zavést podmín¥né pravd¥podobnosti a získat tak tabulku P X\Y 1 0 1 π1|1 π2|1 1 2 π1|2 π2|2 1 kde π1|i je pravd¥podobnost úsp¥chu v i-té kategorii a pro pravd¥podobnost neúsp¥chu v téºe kategorii platí π2|i = 1 − π1|i , m·ºeme tedy zna£it π1 = π1|1 a π2 = π1|2 a tabulku podmín¥ných pravd¥podobností zjednodu²it do tvaru X\Y 1 2
1 π1 π2
0 1-π1 1-π2
P 1 1
Pro podmín¥né pravd¥podobnosti platí
• πj|i = P (Y = j|X = i) =
πij πi.
∀i, j ,
• pokud jsou X a Y nezávislé, pak πj|i = • π ˆj|i =
πij πi.
=
πi. π.j πi.
= π.j
∀i, j ,
nij . ni.
4.1.1 Senzitivita a speci cita testu Tabulka z které získáváme informaci o senzitivit¥ a speci cit¥ má v °ádcích prom¥nnou s kategoriemi úsp¥ch a neúsp¥ch, jejíº °ádkové sou£ty jsou pevn¥ dané, takovou prom¥nnou m·ºe být nap°íklad výskyt choroby u pacienta, ve sloupci jsou pak uvedeny výsledky testu, který m¥l tuto chorobu odhalit. Pracujeme tedy s tabulkou P výskyt choroby\výsledek testu pozitivní negativní ano π1|1 π2|1 1 ne π1|2 π2|2 1 Senzitivita je potom dána jako podmín¥ná pravd¥podobnost, ºe bude výsledek testu pozitivní v p°ípad¥, kdy pacient chorobou skute£n¥ trpí. V tabulce je senzitivita rovna π1|1 . Oproti tomu speci cita testu je podmín¥ná pravd¥podobnost negativního výsledku testu v p°ípad¥ zdravého pacienta. Speci cita je tedy π2|2 . 36
Test funguje ideáln¥ v p°ípad¥, kdy jsou ob¥ tyto podmín¥né pravd¥podobnosti rovny jedné. Senzitivitu a speci citu testu odhadujeme ze vztah·
π ˆ1|1 =
n11 n1.
a
π ˆ2|2 = 1 −
n21 n22 = . n2. n2.
ni1 p°edstavuje po£et pozitivních výsledk· testu v i-té skupin¥, který pochází z binomického rozd¥lení s parametry π1|i a ni ., pro rozptyl odhad· senzitivity a speci city tak platí π1|1 (1 − π1|1 ) n11 1 var(π1|1 ) = var = 2 var(n11 ) = n1. n1. n1. a
π1|2 (1 − π1|2 ) . n2. P°i velkém po£tu pozorování v kaºdé skupin¥ m·ºeme stavit Waldovy intervaly spolehlivosti q var(π2|2 ) = var(1 − π1|2 ) = var(π1|2 ) =
var(ˆ π1|1 )
π ˆ1|1 ± u1− α2
pro senzitivitu a
π ˆ2|2 ± u
1− α 2
q var(ˆ π2|2 )
pro speci citu. u1− α2 je op¥t 1 − α/2 kvantil rozd¥lení N(0, 1).
P°íklad 7. Chceme zjistit, jaká je senzitivita a speci cita diagnostiky zaloºené jen na hladin¥ PSA v krvi. Pokusíme se také nalézt hrani£ní hodnotu, která bude z t¥chto dvou hledisek optimální. K dispozici jsou údaje o vý²i PSA a výsledku biopsie, které pocházejí od 987 pacient· ve v¥ku od 45 do 80 let s objemem prostaty men²ím neº 150 ml a hladinou PSA do 20 ng/ml.
Nejprve uvaºujme nej£ast¥j²í mez 4 ng/ml, kontingen£ní tabulka má v tomto p°ípad¥ podobu P Biopsie\P SA ≥4 <4 pozitivní 177 53 230 negativní 530 227 757 177 Odhad senzitivity je 230 = 0, 77 a to znamená, ºe pokud bychom se rozhodovali jen na základ¥ hladiny PSA, odhalili bychom chorobu u 77 % nemocných. Speci cita takového rozhodování je 227 = 0, 3, 30 % v²ech zdravých pacient· by tedy bylo 757 dal²ího vy²et°ení u²et°eno, na druhou stranu 70 % zdravých by bylo biopsii podrobeno zbyte£n¥, coº p°edstavuje zbyte£n¥ vynaloºené náklady. Abychom v¥d¥li, jak p°esné jsme získali odhady, je t°eba ur£it jejich rozptyly a 95% intervaly spolehlivosti.
var c (π1|1 ) =
0, 77(1 − 0, 77) = 0, 00077 230
p 0, 77±1, 96 0, 00077 = h0, 716; 0, 824i 37
var c (π2|2 ) =
0, 7(1 − 0, 7) = 0, 00028 757
p 0, 3 ± 1, 96 0, 00028 = h0, 267; 0, 333i
Zejména vy²²í senzitivitu by m¥lo zajistit sníºení hrani£ní hodnoty na 2,5 ng/ml, £ímº se zm¥ní kontingen£ní tabulka Biopsie\P SA ≥2,5 <2,5 pozitivní 221 9 230 negativní 725 32 757 P
Senzitivita se opravdu zvý²ila, chorobu tentokrát odhalíme u 221 · 100 % = 96 % 230 nemocných. Stinnou stránkou sníºení hrani£ní hodnoty je ov²em sníºení speci32 = 0, 04 a u 96 % zdravých pacient· by tak bylo rozhodnuto city, ta je nyní jen 757 chybn¥. Op¥t spo£ítáme i rozptyly a intervalové odhady var c (π1|1 ) =
0, 96(1 − 0, 96) = 0, 00017 230
p 0, 96±1, 96 0, 00017 → h0, 934; 0, 986i
var c (π2|2 ) =
0, 96(1 − 0, 96) = 0, 00005 757
p 0, 04±1, 96 0, 0005 → h0, 026; 0, 054i
Z kon den£ních interval· je vid¥t, ºe pokud bychom sníºili hranici na 2,5 ng/ml, byli by k dal²ímu vy²et°ení posláni tém¥° v²ichni pacienti a samotné sledování hladiny PSA by tak ztratilo smysl. Nabízí se otázka, jak volit hrani£ní hodnotu tak, abychom v£as odhalili co nejv¥t²í po£et karcinom· a zárove¬ provád¥li co moºná nejmén¥ zbyte£ných biopsií. První moºností jak tuto hodnotu najít je maximalizace výrazu [λ · senzitivita + (1 − λ) · speci cita],
ve kterém λ ∈ h0, 1i p°edstavuje váhu, kterou p°i°azujeme senzitivit¥ testu. Pokud volíme λ = 0, 6 nebo v¥t²í, jeví se podle tohoto pravidla jako nejvýhodn¥j²í stanovení hranice PSA=2 ng/ml, která má podle na²ich dat 100% senzitivitu a speci citu 1,6 %. Pokud bychom oba ukazatele povaºovali za stejn¥ d·leºité a λ by tedy bylo rovno 0, 5, byla by optimální hrani£ní hodnota PSA=4,5 ng/ml, p°i které je senzitivita 71,7 % a speci cita 41,4 %. Druhou moºností jak nalézt hrani£ní hodnotu je metoda minimaxu, p°i které hledáme maximum p°es v²echny p°ípustné hodnoty PSA z minim dvojice senzitivitaspeci cita ur£ené pro kaºdou hodnotu PSA. Chceme tedy nalézt takovou hladinu PSA, p°i níº budeme dosahovat co nejvy²²ích hodnot senzitivity i speci ty. Touto metodou byla ur£ena mez 5,5 ng/ml, která má sice pom¥rn¥ vysokou speci citu (54,7 %), ale senzitivita je pouze 55 %. P°íklad je uloºen pod názvem priklad_7.R. 38
4.1.2 Rozdíl proporcí, relative risk Pokud pracujeme s tabulkou X\Y 1 2
1 π1 π2
0 1-π1 1-π2
P 1 1
m·ºeme vliv prom¥nné X posoudit pomocí rozdílu pravd¥podobností úsp¥chu v obou kategoriích X π1 − π2 , ten nabývá hodnot mezi -1 a 1 a nezávislosti Y na hodnot¥ X odpovídá situace, kdy je tento rozdíl roven 0. Rozdíl proporcí odhadujeme za pomoci vztahu
n11 n21 − , n1. n2.
π\ ˆ1 − π ˆ2 = 1 − π2 = π
kde ni1 p°edstavuje po£et úsp¥ch· v i-té skupin¥ a ni. celkový po£et pozorování v této skupin¥. A protoºe výb¥ry n11 jednotek z 1. skupiny a n21 jednotek z 2. skupiny povaºujeme za nezávislé výb¥ry z binomického rozd¥lení s parametry ni. a πi , který reprezentuje pravd¥podobnost úsp¥chu v i-té skupin¥, je rozptyl tohoto odhadu var(ˆ π1 − π ˆ2 ) = var(ˆ π1 ) + var(ˆ π2 ) − 2cov(ˆ π1 , π ˆ2 ) =
=
π1 (1−π1 ) n1.
+
π1 (1−π1 ) n1.
+
π2 (1−π2 ) n2.
−0=
π2 (1−π2 ) . n2.
P°i dostate£ném po£tu pozorování lze ur£it Wald·v interval spolehlivosti s π ˆ1 (1 − π ˆ1 ) π ˆ2 (1 − π ˆ2 ) (ˆ π1 − π ˆ2 ) ± u1− α2 + , n1. n2. ve kterém u1− α2 je 1 −
α 2
kvantil rozd¥lení N(0, 1).
Pokud jsou ov²em ob¥ podmín¥né pravd¥podobnosti π1 a π2 blízké hrani£ním hodnotám 0 nebo 1, není toto jejich srovnání p°íli² vhodné a je lep²í vzájemnou závislost X a Y posuzovat pomocí tzv. relative risku
RR =
π1 , π2
který s t¥mito krajními hodnotami pracuje lépe. Pokud jsou X a Y nezávislé, je relative risk roven 1. Relative risk odhadujeme ze vztahu
ˆ1 d=π RR = π ˆ2 39
n11 n1. n21 n2.
.
N¥kdy je ale výhodn¥j²í pracovat s jeho logaritmem, to znamená, ºe odhadujeme
\ ln RR = ln n11 + ln n2. − ln n21 − ln n1. . Abychom odvodili jeho rozptyl, rozd¥líme logaritmus na dva
ln
π ˆ1 ˆ 2 = ln π ˆ1 − lnπ π ˆ2
a odvodíme rozptyl výrazu ln π ˆ1 . Jak jiº bylo uvedeno, po£et úsp¥ch· v první skupin¥ n11 je náhodná veli£ina s binomickým rozd¥lením s parametry π1 a n1. a podle kapitoly 3.2.2 pro odhad π ˆ1 platí
π ˆ − π1 q1
as.
∼ N (0, 1),
π1 (1−π1 ) n1.
neboli
√
as.
n1. (ˆ π1 − π1 ) ∼ N (0, π1 (1 − π1 )).
Tohoto vztahu vyuºijeme v tzv. Delta metod¥, která je popsána nap°. v [15], kapitola 14.1.2. Principem této metody je, ºe pokud platí √ as. n(ˆ π − π) ∼ N (0, σ 2 ), pak pro diferencovatelnou funkci g(π) platí
√ as. n(g(ˆ π ) − g(π)) ∼ N
0, σ 2
δg δπ
2 ! .
V na²em p°ípad¥ je g(π) = ln(π). A protoºe δg(π) = π1 , získáme vý²e uvedený δπ vztah ve tvaru √ 1 as. n1. (ln(ˆ π1 ) − ln(π1 )) ∼ N 0, π1 (1 − π1 ) 2 . π1 Odtud je jiº vid¥t, ºe asymptotický rozptyl ln(ˆ π1 ) lze neformáln¥ zapsat jako
var(ln(ˆ π1 )) =
1 − π1 . n1. π1
Analogicky bychom odvodili i neformální zápis asymptotického rozptylu ln(ˆ π2 ) , který je 1 − π2 var(ln(ˆ π2 )) = . n2. π2 Nyní op¥t vyuºijeme toho, ºe výb¥ry n11 jednotek z 1. skupiny a n21 jednotek z druhé skupiny jsou nezávislé a kovariance mezi ln(ˆ π1 ) a ln(ˆ π2 ) je tak nulová. Hledaný asmytotický rozptyl lze op¥t neformáln¥ zapsat jako
1 − π1 1 − π2 \ var(ln RR) = + . n1. π1 n2. π2 40
Pokud máme k dispozici dostate£ný po£et pozorování, jsme op¥t schopni zkonstruovat interval spolehlivosti r 1−π ˆ1 1 − π ˆ2 \ + ln RR ± u1− α2 π ˆ1 n1. π ˆ2 n2. a z n¥j potom dopo£ítat i interval spolehlivosti pro samotný relative risk, dosazením krajních hodnot intervalu za x do výrazu ex .
P°íklad 8. Op¥t se budeme zabývat problematikou vztahu mezi výsledkem biopsie a hladinou PSA pacienta, kterou d¥líme podle hrani£ní hodnoty PSA=4 ng/ml. Op¥t máme k dispozici údaje o 987 pacientech ve v¥ku od 45 do 80 let, s hladinou PSA men²í neº 20 ng/ml a objemem prostaty men²ím neº 150 ml. Vycházíme z této tabulky: P PSA\Biopsie pozitivní negativní ≥4 177 530 707 <4 53 227 280 P 230 757 987 Podmín¥né pravd¥podobnosti se tentokrát vztahují k výsledku biopsie v závislosti na tom, do které kategorie PSA pacient spadá. Jedná se tedy o opa£nou situaci neº v p°edchozím p°íklad¥. Nyní jsou pravd¥podobnosti odhadnuty takto: PSA\Biopsie pozitivní negativní ≥4 0,25 0,75 <4 0,19 0,81 platí
P
1 1
Pro pravd¥podobnosti pozitivní biopsie v jednotlivých kategoriích π1 a π2 tedy π ˆ1 = P (bio. = poz.|P SA ≥ 4) = π ˆ1|1 = 0, 25 π ˆ2 = P (bio. = poz.|P SA < 4) = π ˆ1|2 = 0, 19
a odhad rozdílu proporcí resp. relative risku je π ˆ1 − π ˆ2 = 0, 25 − 0, 19 = 0, 06
d = 0, 25 = 1, 32. RR 0, 19
Pravd¥podobnost, ºe bude výsledek biopsie pozitivní, je v p°ípad¥ kdy je PSA v¥t²í neº 4 ng/ml o 0,06 v¥t²í neº pokud je PSA men²í neº 4 ng/ml, coº p°edstavuje nár·st o 32 %. Nyní spo£ítáme odhad rozptylu pro rozdíl proporcí var c (ˆ π1 − π ˆ2 ) =
0, 25(1 − 0, 25) 0, 19(1 − 0, 19) + = 0, 00081 707 280
a jeho interval spolehlivosti je p 0, 06 ± 1, 96 0, 00081 = h0, 00529; 0, 11685i. 41
Nejprve spo£ítáme rozptyl logaritmu relative risku \ var c (ln RR) =
1 − 0, 19 1 − 0, 25 + = 0, 0195, 0, 25 · 707 0, 19 · 280
kon den£ní interval logaritmu je potom p ln 1, 32 ± 1, 96 0, 0195 = h0, 0057; 0, 5535i
odtud jiº dopo£ítáme interval spolehlivosti pro relative risk he0,0057 ; e0,5535 i = h1, 0057; 1, 7394i.
Jelikoº interval spolehlivosti pro rozdíl proporcí neobsahuje nulu a kon den£ní interval relative risku naopak jedni£ku, m·ºeme na hladin¥ 0,95 zamítnout hypotézu, ºe výsledek biopsie nezávisí na hladin¥ PSA v krvi. P°íklad je uloºen pod názvem priklad_8.R.
4.1.3 Pom¥r ²ancí Poslední charekteristikou kontingen£ních tabulek 2x2, kterou se budeme zabývat, je pom¥r ²ancí neboli odds ratio (OR). Nejprve de nujme ²anci
Ωi =
πi , 1 − πi
kde πi je podmín¥ná pravd¥podobnost úsp¥chu v i-té kategorii. ance tedy vyjad°uje, kolikrát je v dané kategorii pravd¥podobn¥j²í úsp¥ch oproti neúsp¥chu. Odds ratio potom vyjad°uje, kolikrát je v¥t²í ²ance na úsp¥ch v první kategorii oproti ²anci v kategorii druhé a je dáno jako pom¥r ²ancí v t¥chto dvou kategoriích, tedy π1 Ω1 1−π1 θ= = π2 . Ω2 1−π2 Tento vztah lze ale zjednodu²it tak, abychom byli schopni vypo£ítat OR p°ímo ze sdruºených pravd¥podobností πij .
θ=
π1 1−π1 π2 1−π2
=
π11 π1. π21 π2.
42
π1. π12 π2. π22
=
π11 π22 π12 π21
Mezi základní vlastnosti OR pat°í:
• θ ≥ 0, • Pokud θ = 1, jsou veli£iny X a Y nezávislé a platí, ºe £ím více je θ vzdáleno od 1, tím siln¥j²í je mezi X a Y závislost. • Pokud 1 < θ < ∞, je v první kategorii v¥t²í ²ance na úsp¥ch neº v druhé, pokud ale 0 < θ < 1, pak je to naopak. • N¥kdy m·ºe být výhodn¥j²í pracovat s logaritmem OR, v tom p°ípad¥ ln θ ∈ (−∞, ∞) a nezávislosti odpovídá situace kdy ln θ = 0. • Výhodou OR je, ºe se jeho hodnota nezm¥ní, zm¥níme-li orientaci tabulky, není tedy pot°eba ur£ovat, která prom¥nná je vysv¥tlující a která vysv¥tlovaná. • Pro vztah mezi OR a relative riskem platí OR =
π1 1−π1 π2 1−π2
=
π 1 1 − π2 1 − π2 = RR · , π 2 1 − π1 1 − π1
ob¥ hodnoty si tedy pro velmi malé π1 i π2 jsou blízko. Odhad pom¥ru ²ancí získáme ze vztahu
π ˆ11 π ˆ22 θˆ = = π ˆ12 π ˆ21
n11 n n12 n
n22 n n21 n
=
n11 n22 n12 n21
a jeho hodnota se nezm¥ní pokud oba °ádky resp. sloupce tabulky vynásobíme nenulovou konstantou. Problém ale m·ºe nastat, pokud je n¥který z po£t· nij roven 0, odhad odds ratia je potom roven nule nebo nekone£nu, pokud jsou navíc nulové ob¥ hodonoty ve stejném °ádku nebo sloupci, není tento odhad v·bec de nován. Abychom se tomuto problému vyhnuli, m·ºeme ke kaºdému z po£t· nij p°i£íst 0,5. Druhou nevýhodou neupraveného odds ratia je jeho nesymetri£nost kolem 1, pokud totiº p°edpokládáme, ºe je jeho skute£ná hodnota rovna jedné, odhadem se od této skute£né hodnoty m·ºeme jen velmi málo odchýlit sm¥rem k nule, ale mnohem víc sm¥rem k nekone£nu. Z tohoto hlediska je výhodn¥j²í práce s logaritmem pom¥ru ²ancí, jehoº odhad je
ln θˆ = ln
n11 n22 = ln n11 − ln n12 − ln n21 + ln n22 . n12 n21
Nyní za pomoci následující v¥ty odvodíme rozptyl tohoto odhadu.
43
V¥ta 1. Nech´ α1 , . . . , αk z R spl¬ují, ºe α1 +· · ·+αk = 0 a |αk |+· · · |αk |>0. Nech´ dále X = (X1 , . . . , Xk )0 má multinomické rozd¥lení s parametry n, p1 , . . . , pk . Poloºme δ=
k X
αi ln pi ,
d=
i=1
k X
αi ln Xi ,
i=1
v u k uX α2 i , σ=t p i i=1
v u k uX α2 i s=t . X /n i i=1
Potom pro n → ∞ platí √ d−δ d n → N (0, 1), σ
√ d−δ d n → N (0, 1) s
D·kaz viz [17], v¥ta 12.6. Vektor X v tomto p°ípad¥ p°edstavují po£ty pozorování v jednotlivých skupinách (n11 , n12 , n21 , n22 ), koe cienty αi jsou v p°ípad¥ odhadu podílu ²ancí rovny postupn¥ 1, -1, -1 a 1 a ve v¥t¥ uvaºovaná prom¥nná d je práv¥ ln θˆ. Dosadíme do vztahu pro s s r √ 1 1 1 1 1 1 1 1 + + + . s = n11 + n12 + n21 + n22 = n n11 n12 n21 n22 n n n n Úpravou vztahu
√
d
n d−δ → N (0, 1) získáme s s2 d − δ → N (0, ) n d
a protoºe δ = ln p11 − ln p12 − ln p21 + ln p22 je nenáhodné, pro velké rozsahy n platí s2 var(d) = . n Po dosazení za d a s jiº získáme výsledný vztah pro rozptyl ln θˆ 1 1 1 1 n n11 + n12 + n21 + n22 1 1 1 1 ˆ = var(ln θ) = + + + . n n11 n12 n21 n22 P°i velkém rozsahu souboru m·ºeme také ur£it Wald·v interval spolehlivosti pro logaritmus OR, který bude mít tvar r 1 1 1 1 ln θˆ ± u1− α2 + + + , n11 n12 n21 n22 kde u1− α2 je kvantil normovaného normálního rozd¥lení. Interval spolehlivosti pro θ získáme dosazením spo£ítaných krajních bod· do exponentu výrazu ex . 44
P°íklad 9. Op¥t se budeme v¥novat problematice p°edchozích p°íklad·, tedy vztahu mezi vý²í PSA a výsledkem biopsie. Pro p°ipomenutí uvádím p·vodní tabulku a tabulku s odhadnutými podmín¥nými pravd¥podobnostmi. PSA\Biopsie pozitivní negativní ≥4 177 530 707 <4 53 227 280 P 230 757 987 P
PSA\Biopsie pozitivní negativní ≥4 0,25 0,75 <4 0,19 0,81
P
1 1
Nejprve odhadneme ²ance, ºe výsledek biopsie bude v jednotlivých skupinách pozitivní Ω1 =
0, 25 = 0, 33 0, 75
Odhad odds ratia je
Ω2 =
0, 19 = 0, 23. 0, 81
177 · 227 θˆ = = 1, 43, 53 · 530
u pacient· jejichº hladina PSA je v¥t²í neº 4 ng/ml je 1,43krát v¥t²í ²ance, ºe se potvrdí p°ítomnost karcinomu prostaty. Abychom v¥d¥li jak p°esný tento ukazatel je, musíme ur£it jeho kon den£ní interval. Hodnota logaritmu odds ratia je ln θˆ = ln 177 − ln 53 − ln 530 + ln 227 = 0, 36,
odhad jeho rozptylu je var c ln θˆ =
1 1 1 1 + + + = 0, 03 177 53 530 227
a jeho 95% kon den£ní interval je potom p 0, 36 ± 1, 96 0, 03 = h0, 02; 0, 7i.
Nyní jiº m·ºeme dopo£ítat interval spolehlivosti pro samotný pom¥r ²ancí he0,02 ; e0,7 i = h1, 02; 2, 1i.
Protoºe interval spolehlivosti pro odds ratio neobsahuje jedni£ku, m·ºeme i tentokrát zamítnout hypotézu, ºe výsledek biopsie nezávisí na hladin¥ PSA, a to na hladin¥ 1-α=0,95. e²ení tohoto p°íkladu je uloºeno pod názvem priklad_9.R.
45
5
Logistická regrese
Nyní se budeme zabývat situací, kdy je alternativn¥ rozd¥lená prom¥nná vysv¥tlována pomocí p veli£in, které mohou být jak kategoriální, tak i spojité. Pro nalezení vztahu mezi t¥mito prom¥nnými máme k dispozici n nezávislých m¥°ení vysv¥tlované prom¥nné Yi a k ní p°íslu²ejícího (p+1)-rozm¥rného vektoru vysv¥tlujících prom¥nných X i = (1, Xi1 , . . . , Xip )0 . Jak bylo odvozeno v kapitole 3.2, st°ední hodnota Yi je rovna pravd¥podobnosti úsp¥chu p°i dané hodnot¥ vysv¥tlujících prom¥nných π(Xi ), platí tedy rovnost
E(Yi ) = µi = π(Xi ). Dále také víme, ºe rozptyl vysv¥tlované prom¥nné je
var(Yi ) = π(Xi )(1 − π(Xi )) = µi (1 − µi ), závisí tedy na její st°ední hodnot¥. První moºností jak najít vztah mezi Yi a Xi je pouºití lineárního regresního modelu E(Yi ) = µi = π(Xi ) = X0i β, kde β = (β0 , β1 , ..., βp )0 je (p + 1)-rozm¥rný vektor neznámých parametr·. Takto ˆ , které by pro n¥které vektory x0 generobychom ale získali odhady parametr· β valy hodnoty E(Y0 ) mimo interval h0, 1i a tedy hodnoty nesmyslné. Model lineární regrese tedy pouºít lze, ale jen pro omezenou mnoºinu vysv¥tlujících vektor· x. Abychom se tomuto omezení vyhnuli, je pot°eba najít tzv. linkovou funkci g(µi ), která bude zobrazením intervalu h0, 1i na mnoºinu v²ech reálných £ísel. Potom jiº budeme moci pouºít model g(µi ) = X0i β . Zatímco v p°ípad¥ lineární regrese se jednalo o identický link (g(µi ) = µi ), logistická regrese pouºívá tzv. kanonický link, který je obecn¥ dán vztahem
g(µi ) = Q(θi ), ve kterém Q(θi ) p°edstavuje tzv. p°irozený parametr, jehoº podobu pro alternativní rozd¥lení nyní odvodíme. Alternativní rozd¥lení stejn¥ jako binomické nebo nap°íklad Poissonovo pat°í do t°ídy rozd¥lení exponenciálního typu, jejichº hustota se dá zapsat ve tvaru f (y, θ) = a(θ)b(y)eyQ(θ) ,
θ je obecn¥ vektor parametr· rozd¥lení a hledaný p°irozený parametr Q(θ) se nachází v exponentu funkce. Abychom na²li jeho hodnotu pro p°ípad alternativního rozd¥lení, musíme jeho hustotu p°evést do stejného tvaru y π π y 1−y = (1 − π)ey ln 1−π . f (y, π) = π (1 − π) = (1 − π) 1−π 46
Odtud jiº vidíme, ºe θ = π , a(θ) = (1 − π), b(y) = 1 a
Q(θ) = ln
π . 1−π
Logistická regrese tedy pracuje s modelem
ln
π = X0 β. 1−π
Výraz na levé stran¥ je vlastn¥ logaritmus ²ance na úsp¥ch v jednotlivých pokusech, jedná se o tzv. logit transformaci parametru π . Nyní je²t¥ ov¥°me, ºe takto opravdu získáme takové odhady π ˆ , které padnou 0 do intervalu h0, 1i. Ozna£me η(β) = X β , pak platí
ln
π = η(β), 1−π
odtud
π = eη(β) 1−π a snadnou úpravou dosp¥jeme ke vztahu π=
1 eη(β) = , 1 + eη(β) 1 + e−η(β)
ze kterého uº je vid¥t, ºe pro libovolnou hodnotu x opravdu získáme E(Y ) = π jen z intervalu h0, 1i.
5.1
Interpretace parametr·
Vyjd¥me z modelu
ln
π(X) = β0 + β1 X1 + · · · + βp Xp 1 − π(X)
a v¥nujme se nejprve absolutnímu £lenu β0 . Jeho vliv je z°ejmý pokud budou v²echny hodnoty vysv¥tlujících prom¥nných nulové, pak bude platit:
ln
π(X1 = 0, ..., Xp = 0) = β0 1 − π(X1 = 0, ..., Xp = 0)
parametr β0 je tedy roven logitu ²ance na úsp¥ch oproti neúsp¥chu v p°ípad¥, kdy jsou v²echny regresory nulové. Pro samotnou ²anci potom platí
π(X1 = 0, ..., Xp = 0) = e β0 . 1 − π(X1 = 0, ..., Xp = 0) 47
Ne vºdy má ale praktický smysl uvaºovat situaci, kdy jsou v²echna x1 , . . . , xp rovna nule, v tom p°ípad¥ je m·ºeme p°e²kálovat tak, aby nule odpovídala je˜ , se sloºkami jich pr·m¥rná hodnota. Pracovali bychom tedy s novým vektorem x x˜i = xi − x¯i pro i = 1, 2, . . . , p a absolutní £len by odpovídal logaritmu ²ance na úsp¥ch oproti neúsp¥chu v p°ípad¥, kdy by v²echny regresory byly na úrovni svých pr·m¥rných hodnot. Ne vºdy je v²ak nutné a také vhodné m¥nit v²echny prom¥nné, ale sta£í zm¥nit jen n¥které z nich a podle toho pak upravit i interpretaci parametru β0 . Pon¥kud odli²ná je interpretace parametr· βi , i = 1, 2, . . . , p, jak ukáºeme nap°íklad na parametru β1 . Za xujeme-li regresory Xi , i = 2, . . . , p na n¥jaké pevné hodnot¥ xi , i = 2, . . . , p, získáme rovnici
ln
π(X1 = x1 ) = β0 + β1 x1 + β2 x2 + · · · + βp xp , 1 − π(X1 = x1 )
pokud nyní x1 zv¥t²íme o jedni£ku a zbylé regresory ponecháme stejné, máme novou rovnici
ln
π(X1 = x1 + 1) = β0 + β1 x1 + β1 + β2 x2 + · · · + βp xp . 1 − π(X1 = x1 + 1)
Nyní od sebe tyto rovnice ode£teme
ln
π(X1 = x1 + 1) π(X1 = x1 ) − ln = β1 , 1 − π(X1 = x1 + 1) 1 − π(X1 = x1 ) ln
π(X1 =x1 +1) 1−π(X1 =x1 +1) π(X1 =x1 ) 1−π(X1 =x1 )
= β1 .
Parametr β1 tedy vyjad°uje logaritmus pom¥ru ²ancí v situaci, kdy hodnotu regresoru X1 zv¥t²íme o jednotku a hodnoty v²ech ostatních necháme nezm¥n¥né. Toto tvrzení lze roz²í°it i na zbylé parametry. Platí tedy, ºe parametr βi , i = 1, 2, . . . , p je dán jako logaritmus pom¥ru ²ancí p°i zm¥n¥ xi o jedni£ku a pevných hodnotách xj , j 6= i, neboli βi
e =
π(Xi =xi +1) 1−π(Xi =xi +1) π(Xi =xi ) 1−π(Xi =xi )
= OR(xi +1,xi ) .
Druhou moºnost jak interpretovat význam parametru βi , lze ukázat, pokud uvaºujeme jen jednu spojitou vyv¥tlující prom¥nnou X a pracujeme tedy s modelem π(x) ln = β0 + β1 x 1 − π(x) a pro π(x) platí
π(x) =
eβ0 +β1 x . 1 + eβ0 +β1 x 48
Pomocí parametru β1 jsme schopni ur£it, jak se zm¥ní pravd¥podobnost π(x) p°i malé zm¥n¥ regresoru x z hodnoty x0 . Tato zm¥na je dána jako tangens úhlu, který svírá te£na k pravd¥podobnostní funkci π(x) v bod¥ x0 s osou x, neboli derivace funkce π(x) podle x vy£íslená v bod¥ x0 . β +β x 1 e 0 δ β1 eβ0 +β1 x0 (1+eβ0 +β1 x0 )−β1 e2(β0 +β1 x0 ) δπ(x) 1+eβ0 +β1 x = = = δx δx (1+eβ0 +β1 x0 )2 x=x0 x=x0
=
β1 eβ0 +β1 x0 (1+eβ0 +β1 x0 )2
= β1 π(x0 ) 1+eβ01+β1 x0 = β1 π(x0 ) 1 −
eβ0 +β1 x0 1+eβ0 +β1 x0
=
= β1 π(x0 )(1 − π(x0 )). Tento výraz je nejvy²²í pokud je π(x0 ) = 21 a v okolí bodu x0 , pro který je tato rovnost spln¥na se tedy pravd¥podobnost π nejvíce m¥ní. Pro bod x0 , kde π(x0 ) = 21 , platí: 1
ln 1−2 1
= β 0 + β 1 x0
ln 1 0 x0
= β 0 + β 1 x0 = β 0 + β 1 x0 = − ββ10 .
2
5.2
Odhad parametru β
Regresní parametr β odhadujeme metodou maximální v¥rohodnosti, jak je uvedeno zejména v [15], pop°ípad¥ v [19]. Máme-li k dispozici ni nezávislých pozorování p°i stejné hodnot¥ xi = (1, xi1 , . . . , xip ), pro i = 1, . . . , N ve kterých je yi úsp¥ch·, jsou P náhodné veli£iny Y1 , . . . , YN nezávislé a mají binomické rozd¥lení s paramtery ni , N i=1 ni = n a π(xi ) a jejich sdruºená pravd¥podobnostní funkce je p°ímo úm¥rná QN QN π(xi ) yi yi ni −yi [1 − π(xi )]ni = π(x ) [1 − π(x )] = i i i=1 i=1 1−π(xi )
=
QN
i=1 e
PN
= e
i=1
ln
π(xi ) 1−π(xi )
yi ln
y
π(xi ) 1−π(xi )
i
QN
i=1 [1
QN
− π(xi )]ni =
i=1 [1
− π(xi )]ni .
Pp
PN
Nyní vyuºijeme znalosti toho, ºe π(xi ) • ln 1−π(x = i)
Pp
j=0 βj xij −
• 1−π(xi ) = π(xi )e
Pp
j=0
a βj xij
e =
PN
i=1
yi
Pp
j=0
Pp β x j=0 j ij Pp β x 1+e j=0 j ij
e
49
βj xij
e−
=e
Pp
j=0
j=0 (
βj xij
i=1
yi xij )βj
−1 Pp = 1 + e j=0 βj xij
a vztah zjednodu²íme a zaneseme do n¥j hledaný parametr β . Takto získáme v¥rohodostní funkci Pp
l(β) = e
PN
j=0 (
i=1
yi xij )βj
N Y
Pp
1+e
j=0
βj xij
−ni
.
i=1
Op¥t se snáze maximalizuje její logaritmus ! p N N X X X Pp yi xij βj − ni ln 1 + e j=0 βj xij , L(β) = j=0
i=1
i=1
který derivujeme podle β a jednotlivé parciální derivace poloºíme rovny nule. N
Pp
N
X δL(β) X e k=0 βk xik ! Pp = 0 j = 0, ..., p. = yi xij − ni xij β x δβj 1 + e k=0 k ik i=1 i=1 A protoºe rovnic
Pp ˆ x β k=0 k ik Pp ˆ x β 1+e k=0 k ik
e
N X
=π ˆi , odhadujeme vektorový parametr β z p+1 normálních
yi xij −
N X
i=1
ni π ˆi xij = 0
j = 0, 1, 2, . . . , p
i=1
nebo s vyuºitím rovnosti µ ˆ i = ni π ˆi N X i=1
yi xij −
N X
µ ˆi xij = 0
j = 0, 1, 2, . . . , p,
i=1
tyto rovnice pak lze zapsat i vektorov¥
X0 y = X0 µ ˆ, X je matice pozorování xij neboli matice plánu. Protoºe jsou normální rovnice v parametrech beta nelineární, musíme je °e²it itera£n¥ nap°íklad pomocí Newton-Raphsonovy metody, která bude p°edstavena v následující podkapitole.
5.2.1 Newton-Raphsonova metoda Metoda v kaºdé iteraci odhaduje parametry β, π, µ, které zna£íme β (t) , π (t) , µ(t) , index t = 0, 1, . . . ozna£uje po°adí iterace. Nejprve ur£eme gradient a hessián logaritmu v¥rohodnostní funkce. Pro gradient v t-té iteraci t = 0, 1, . . . platí (t) gj
N X δL(β) (t) = = yi − ni πi xij δβj β(t) i=1 50
j = 0, 1, . . . , p
coº je maticov¥
g (t) = X0 (y − µ(t) )
a pro prvky hessiánu potom (t) hab
N X δ 2 L(β) (t) (t) 1 − π = = − x x n π ia ib i i i δβa δβb β(t) i=1
Pouºijeme-li diagonální matici diag prvek
(t) ni πi (1
−
(t) πi ),
h
(t) ni πi (1
−
i
(t) πi )
a, b = 0, 1, . . . , p. , která má na pozici [i, i]
i = 1, . . . , N , lze hessián maticov¥ zapsat takto: h i (t) (t) −X0 diag ni πi (1 − πi ) X. (t)
ˆi s vyuºitím parametru β (t) a V obou vztazích se vyskytuje πi , coº je odhad π vztahu P (t) p j=0 βj xij e (t) P . πi = (t) p j=0 βj xij 1+e Nyní aproximujeme L(β) v bod¥ β (t) pomocí Taylorova rozvoje, ve kterém uvaºujeme jen £leny nejvý²e druhého °ádu
1 L(β) ≈ L(β (t) ) + g (t)0 (β − β (t) ) + (β − β (t) )0 H (t) (β − β (t) ). 2 Její derivace je potom
δL(β) = g (t) + H (t) (β − β (t) ) δβ a odtud jiº získáme odhad β (t+1) pro dal²í iteraci
β (t+1) = β (t) − H (t)−1 g (t) . Dosadíme-li za gradient a hessián, získáme kone£ný vztah pro výpo£et nového odhadu parametru
β
(t+1)
=β
(t)
h i −1 (t) (t) 0 + X diag ni πi (1 − πi ) X X 0 (y − µ(t) ). (t+1)
Na základ¥ β (t+1) potom ur£íme nový odhad pravd¥podobnosti πi . Je-li dob°e zvolen po£áte£ní odhad β (0) , konverguje metoda pom¥rn¥ rychle k maximáln¥ v¥rohodným odhad·m β a π ˆ . Itera£ní postup zastavíme v okamºiku, kdy se hod(t) (t) nota odhadu β resp. π m¥ní jiº jen velmi málo. 51
Druhou itera£ní metodou pro odhad parametr· zobecn¥ného lineárního modelu je Fisherova skórová metoda, která namísto hessiánu pracuje s jeho o£ekávanou hodnotou, tyto dv¥ matice jsou ale v p°ípad¥ logistické regrese shodné, protoºe hessián nezávisí na hodnotách sledované veli£iny Y . Fisherova skórová metoda tedy p°i výpo£tu parametr· logistické regrese splývá s metodou NewtonRaphsonovou. Newton-Raphsonova metoda spo£ívá ve zp°es¬ování po£áte£ního odhadu parametr· β (0) . Abychom tento odhad získali, nahlíºíme na kaºdý z parametr· (0) eβj xj βj , j = 0, 1, . . . , p odd¥len¥ a p°edpokládáme, ºe derivace funkce π(xj ) = 1+e β j xj , j = 1, . . . , p je nejv¥t²í v okolí st°edu intervalu, na kterém veli£inu xj m¥°íme. Tento bod je dán vztahem
min(xj ) + max(xj ) . 2 Dále p°edpokládejme, ºe úsp¥ch je v tomto bod¥ stejn¥ pravd¥podobný jako neúsp¥ch, platí tedy π(¯ xj ) = 12 a pro úhel α, který svírá te£na k funkci v tomto bod¥ s osou x, platí 1 . tg(α) = max(xj ) − min(xj ) Derivace funkce π(xj ) je rovna x¯j =
δπ(xj ) = βj π(xj )(1 − π(xj )), δxj jak bylo odvozeno v p°edchozí kapitole. V bod¥ x ¯j je tato derivace rovna tg(α) a platí 1 βj = . 4 max(xj ) − min(xj ) Odtud jiº ur£íme po£áte£ní odhad (0)
βj =
4 , max(xj ) − min(xj )
j = 1, 2, . . . , p.
Dále musíme zohlednit, zda pravd¥podobnost úsp¥chu s r·stem j-tého regresoru spí²e roste nebo spí²e klesá. Pokud roste, p°i°adíme parametru kladné znaménko a pokud pravd¥podobnost klesá, bude parametr záporný. (0)
Vý²e uvedeným zp·sobem ale nelze odhadnout parametr β0 . Jeho hodnotu získáme ze vztahu π (0) ln = β0 , 1−π kdy za π dosadíme podíl celkového po£tu úsp¥ch· a celkového po£tu pozorování. Takto jiº získáme celý po£áte£ní odhad β (0) , který pouºijeme p°i iniciaci NewtonRaphsonovy metody. Pokud metoda p°i startu z tohoto bodu nekonverguje, máme alespo¬ p°ibliºnou informaci o tom, kde se hledané odhady nacházejí a v souladu s touto informací po£áte£ní odhady zm¥níme. 52
5.2.2 Varian£ní matice odhadu parametr· Podobn¥ jako tomu bylo v p°ípad¥ binomického rozd¥lení, získáme varian£ní matici odhadu jako inverzi z informa£ní matice I , jejíº prvky jsou dány vztahem 2 δ L(β) pro a, b = 0, 1, . . . , p. −E δβa δβb V úvodu kapitoly 5.2 jiº bylo odvozeno, ºe N
Pp
N
X δL(β) X e k=1 βk xik Pp = yi xia − ni xia . k=1 βk xik δβa 1 + e i=1 i=1 Nyní je²t¥ provedeme druhou derivaci podle βb δL2 (β) δβa βb
= −
PN
= −
PN
Pp Pp Pp Pp β x β x β x β x k=1 k ik (1+e k=1 k ik )−e k=1 k ik e k=1 k ik 2 Pp βk xik 1+e k=1
e i=1 xia xib ni
i=1 ni xia xib
op¥t vyuºijeme toho, ºe
Pp β x k=1 k ik Pp β x 1+e k=1 k ik
e
Pp β x k=1 k ik Pp β x 1+e k=1 k ik
e
= πi a
=
1
Pp β x k=1 k ik
1+e
1
1+e
Pp β x k=1 k ik
= 1 − πi ,
N
X δL2 (β) =− xia xib ni πi (1 − πi ). δβa βb i=1 Protoºe druhá derivace nezávisí na hodnotách yi , je st°ední hodnota informa£ní matice stejná jako hodnota pozorovaná
I = X0 W X
W = diag[ni πi (1 − πi )]
a odhad varian£ní matice je
ˆ = (X0 W ˆ X)−1 . var( c β) Protoºe je informa£ní matice rovna záporn¥ vzatému hessiánu, lze varian£ní matici získat jiº v pr·b¥hu itera£ního hledání odhad· β . P°esnost jednotlivých odhad· βˆj potom ur£íme jako odmocninu z diagonálních prvk· varian£ní matice. Výhodné také je, ºe matice W je pozitivn¥ de nitní a v p°ípad¥ plné sloupcové hodnosti matice plánu X je pozitivn¥ de nitní i informa£ní matice I a tedy i hessián H . Ten pak lze snadno invertovat a pouºít v itera£ním procesu.
53
5.2.3 Intervaly spolehlivosti a testování parametru βj Máme-li k dispozici odhad varian£ní matice, m·ºeme ur£it interval spolehlivosti pro jednotlivé odhady βj . Jejich krajní body získáme ze vztahu q ˆ jj , ˆ c β) βj ± u1− α2 var( ve kterém u1− α2 p°edstavuje (1 − α/2)-kvantil rozd¥lení N(0, 1). Takto lze také testovat, zda je parametr βj nulový, coº by znamenalo, ºe hodnota vysv¥tlující prom¥nné nezávisí na regresoru xj . Waldova testová statistika má tvar
βˆj z=q ˆ jj var( c β) a za platnosti H0 : βj = 0 má z 2 asymptoticky χ2 rozd¥lení o jednom stupni volnosti. Pro testování lze pouºít také test zaloºený na v¥rohodnostním pom¥ru a skórový test, které byly p°edstaveny v kapitole 3.2.2 a jsou blíºe popsány v [15]. Dále m·ºeme ur£it kon den£ní interval pro logit(π ˆ ) a následn¥ i pro samotný parametr π ˆ . Pro rozptyl logitu platí π ˆ (x) ˆ = xvar( ˆ 0 = var(xβ) c β)x var ln 1−π ˆ (x) a pro krajní body jeho intervalu spolehlivosti potom q π ˆ (x) ˆ 0. α ln ± u1− 2 xvar( c β)x 1−π ˆ (x)
x v obou vztazích zna£í °ádkový vektor regresor·. Kon den£ní interval hπd , πu i pro π(x) získáme p°epo£tem krajních bod· p°edchozího intervalu, které ozna£íme logitd a logitu podle vztahu
πd =
5.3
elogitd 1 + elogitd
πu =
elogitu . 1 + elogitu
Hodnocení kvality modelu
5.3.1 Deviance Jednou z moºností jak ov¥°it, zda uvaºovaný model dob°e popisuje skute£ná data, je pouºití deviance. M¥jme nejprve model, který má práv¥ tolik parametr·, kolik je v datovém souboru r·zných hodnot vysv¥tlujících prom¥nných xi , pro 54
i = 1, 2, . . . , N . Odhad st°ední hodnoty vysv¥tlované prom¥nné p°i hodnot¥ vysv¥tlujích prom¥nných xi , který zna£íme µ ˆi , je v tomto modelu roven p°ímo hodnot¥ yi pro i = 1, 2, . . . , N , neboli µ ˆ = y. Tento model nazveme saturovaný a a£koliv je hodnota jeho v¥rohodnostní funkce L(y, y) nejvy²²í ze v²ech model·, není pro praktické pouºití v·bec vhodný, nebo´ neumoº¬uje ºádné zobecn¥ní. Saturovaný model lze ale pouºít k hodnocení jiného modelu a to práv¥ prost°ednictvím deviance, která je obecn¥ dána vztahem −2[L(ˆ µ, y) − L(y, y)], ve kterém L(ˆ µ, y) je maximum logaritmu v¥rohodnostní funkce sledovaného modelu. Deviance je tedy statistikou testu pom¥ru v¥rohodností, který byl uveden v kapitole 3.3.2, v n¥mº nulovou hypotézu p°edstavuje platnost modelu a alternativu platnost saturovaného modelu. Nyní odvo¤me tvar deviance pro p°ípad binomicky rozd¥lené závislé prom¥nné. Ozna£me θ˜ odhad parametru θ v saturovaném modelu, ve kterém také platí µ˜i = yi pro kaºdé i, odhady ve sledovaném modelu pak ozna£me θˆ a odhadnuté st°ední hodnoty pak µ ˆi . Protoºe logaritmus v¥rohodnostní funkce je pro rozd¥lení exponenciálního typu dán vztahem
L(β) =
X yi θi − b(θi )
+
a(φ)
i
X
c(yi , φ),
i
lze devianci zapsat ve tvaru
2
X yi θ˜i − b(θ˜i ) a(φ)
i
nebo
2
−2
X yi θˆi − b(θˆi ) i
a(φ)
X yi (θ˜i − θˆi ) − b(θ˜i ) + b(θˆi ) a(φ)
i
.
V p°ípad¥ binomického modelu platí yi π ˆi • θˆi = ln 1−ˆ a θ˜i = ln 1−y , πi i ˆ • b(θˆi ) = ln 1 + eθi = ln 1 +
π ˆi 1−ˆ πi
• b(θ˜i ) = − ln(1 − yi ) a • a(φ) =
1 , ni
55
= − ln(1 − π ˆi ) a obdobn¥
coº je odvozeno nap°íklad v [15]. Devianci pro binomický model získáme dosazením t¥chto hodnot do p°edchozího vztahu i P h yi π ˆi − ln + ln(1 − y ) − ln(1 − π ˆ ) = D = 2 i ni yi ln 1−y i i 1−ˆ πi i
= 2
P
yi ni yi ln 1−y −2 i
= 2
P
i yi ni yi ln nin−n −2 i yi
= 2
P
= 2
P
ni yi ln nniiπyˆii + 2
= 2
P
ni yi ln
= 2 = 2
i
i
i ni yi ln
i
P
i
π ˆi ni yi ln 1−ˆ +2 πi
P
i
ni yi ni −ni π ˆi ni −ni yi ni π ˆi
ni yi ni π ˆi
P
i
P
1−yi ni ln 1−ˆ = πi
i
ˆi iπ ni yi ln nin−n +2 ˆi iπ
+2
P
i
1−ˆ πi 1−yi
yi
P
i
1−yi 1−ˆ πi
+2
ni ln
P
ni yi ln nniiπyˆii + 2
P
πi ni (yi − 1) ln 1−ˆ = 1−yi
P
ni yi ln nniiπyˆii + 2
P
i
i
i
i
i (ni
1−yi ni ln 1−ˆ = πi
1−yi ni ln 1−ˆ = πi
P
i
i
1−yi ni ln 1−ˆ = πi
πi ni yi ln 1−ˆ +2 1−yi
P
=
−ni yi − ni yi ) ln nnii−n . ˆi iπ
Deviance tedy srovnává po£et pozorovaných úsp¥ch·, resp. neúsp¥ch· s po£ty odhadnutými. Vlastnosti deviance se li²í v závislosti na tom, zda jsou data seskupená nap°íklad tak, ºe jsou spojité regresory rozd¥leny na intervaly, nebo neseskupená. V prvním p°ípad¥ pracujeme s binomickými veli£inami Yi ∼ Bi(πi , ni ) P i = 1, 2, . . . , N , i ni = n, ve druhém s alternativn¥ rozd¥lenou veli£inou Yi ∼ Alt(πi ) i = 1, 2, . . . , n. Pouze v p°ípad¥ seskupených dat, má deviance asymptoticky χ2 rozd¥lení, jehoº stupn¥ volnosti jsou dány jako rozdíl v po£tu parametr· v obou modelech, tedy N − p. Devianci lze pouºít pro srovnání dvou model·, z nichº jeden M0 je podmodelem druhého M1 . Protoºe v podmodelu M0 odhadujeme mén¥ parametr·, není maximum logaritmu v¥rohodnostní funkce v¥t²í neº maximum pro sloºit¥j²í model M1 a pro deviance tak platí opa£ná nerovnost
D(y, µ ˆ 1 ) ≤ D(y, µ ˆ 0 ). Pokud bychom pomocí testu pom¥ru v¥rohodností testovali hypotézu, ºe je platný podmodel M0 , oproti alternativ¥, ºe platí M1 , zjednodu²ila by se testová statistika do tvaru D(y, µ ˆ 0 ) − D(y, µ ˆ1) a má op¥t asymptoticky χ2 rozd¥lení o stupních volnosti daných rozdílem v po£tu parametr· v obou modelech. Proti vhodnosti podmodelu sv¥d£í velké hodnoty 56
rozdílu deviancí a protoºe jiº nepot°ebujeme znát po£et parametr· odhadovaných v saturovaném modelu, lze takto testovat i modely s neseskupenými daty. Pomocí rozdílu deviancí lze testovat také významnost parametru βj , který se vyskytuje v bohat²ím modelu M1 , hypotéza H0 : βj = 0 je totiº ekvivalentní s hypotézou, ºe platí podmodel M0 , který se od modelu M1 li²í jen tím, ºe neobsahuje j-tý regresor.
5.3.2 Rezidua Pokud je deviance modelu tak velká, ºe zamítneme hypotézu o jeho vhodnosti, m·ºeme se podrobn¥ji zabývat tím, které sloºky deviance k tomu nejvíce p°isp¥ly. Devianci lze zapsat také ve tvaru X D(y, µ) ˆ = di , i
kde
yi ni − yi yi (θ˜i − θˆi ) − b(θ˜i ) + b(θˆi ) = 2 yi ln + (ni − yi ) ln . di = 2 a(φ) ni π ˆi ni − ni π ˆi M·ºeme tedy pracovat s rezidui ve tvaru p p di · sign(yi − µ ˆi ) = di · sign(yi − ni π ˆi ). Tato rezidua jsou v anglicky psané literatu°e ozna£eny jako deviance reziduals, nezahrnují ov²em variabilitu vzniklou odhadem pparametr· βj . Tu zohledníme ˆ ii , ve kterém h ˆ ii je odhad jejich standardizací, neboli vyd¥lením výrazem 1 − h diagonálního prvku tzv. hat matice
H = W 1/2 X(X0 W X)−1 X0 W 1/2 , kde W = diag[ni πi (1 − πi )] a X je matice plánu. Pro hodnocení kvality odhadu lze pouºít také Pearsonova rezidua, která jsou dána vztahem yi − ni π ˆi . ei = p ni π ˆi (1 − π ˆi ) Pokud se£teme mocniny Pearsonových reziduí, získáme Pearsonovu statisPN druhé 2 2 tiku X = i=1 ei uvedenou v kapitole 3.3.2. I tentokrát lze pracovat se standardizovanými Pearsonovými rezidui, která vzniknou stejn¥, jako tomu bylo v p°ípad¥ deviance residuals, tedy yi − ni π ˆi ei ri = p =q . ˆ ii ˆ 1−h ni π ˆi (1 − π ˆi )(1 − hii ) Pokud je model zvolen správn¥, má reziduum ri normované normální rozd¥lení, situace kdy |ri | > 2 zna£ní, ºe zvolený model odhaduje hodnoty πi (xi ) ²patn¥. 57
5.3.3 Modi kace koe cientu determinace Podobn¥ jako lze v p°ípad¥ lineární regrese posuzovat model pomocí koe cientu determinace, lze takto po ur£itých úpravách hodnotit kvalitu i u logistického modelu. První moºností je McFadden·v koe cient determinace, který je p°ímou analogií b¥ºného koe cientu, která jen místo sou£tu £tverc· pracuje s deviancemi. Tento koe cient je dán vztahem
RL2 = 1 −
D(y, µ) ˆ , D0
kde D0 je deviance tzv. nulového modelu, Pn yi ve kterém jsou v²echny st°ední hodnoty µi stejné a jejich odhad je roven i=1 n . P ˆi + (1 − yi ) ln(1 − µ ˆi )] = D0 = −2 i [yi ln µ
= −2
P
= −2
P
P
yi ln
i
j
yj
n
+ (1 − yi ) ln 1 −
P
j
yj
n
=
i
yi ln
P
j
yj − yi ln n + (1 − yi ) ln(n −
P
j
yj ) − (1 − yi ) ln n =
P P P P P = −2 i yi ln j yj + i (1 − yi ) ln(n − j yj ) − i ln n = P P P P = −2 j yj ) − n ln n . i yi ) ln(n − j yj + (n − i yi ln Podobn¥ jako v lineární regresi model i tentokrát tím lépe popisuje reálná data, £ím blíºe je tento koe cient jedné. Druhým koe cientem je Nagelkerk·v koe cient determinace D(y,µ)−D ˆ 0
2 RN
n R2 1−e = 2 = D0 Rmax 1 − e− n
V tomto koe cientu je R2 jedna z moºností zápisu b¥ºného koe cientu deter2 minace (viz [19]) a Rmax je hodnota tohoto koe cientu v p°ípad¥ saturovaného modelu, jehoº deviance je nulová.
5.3.4 Akaikeho informa£ní kritérium Jako modi kaci deviance lze chápat Akaikeho informa£ní kritérium (AIC), které ov²em navíc zohled¬uje po£et parametr· p v modelu a up°ednost¬uje tak jedno58
du²²í modely p°ed sloºit¥j²ími. Kritérium je dáno vztahem
AIC = −2[L(µ, ˆ y) − p]. Srovnáváme-li více model·, lze podle tohoto kritéria povaºovat za nejlep²í ten, kterému náleºí nejniº²í hodnota AIC.
59
6
P°íklad pouºití logistické regrese - diagnostika rakoviny prostaty
Nyní pouºijeme model logistické regrese k tomu, abychom popsali vliv jednotlivých ukazatel· na pravd¥podobnost toho, ºe histologické vy²et°ení vzork· prostaty odhalí výskyt rakoviny. K dispozici máme údaje o vy²et°eních, která prob¥hla ve Fakultní nemocnici Olomouc od za£átku £ervna roku 2006 do konce prosince roku 2010. Tento datový soubor obsahoval záznamy z celkem 1723 vy²et°ení, kdy kaºdý záznam sestával z t¥chto poloºek: identi kátor pacienta a jeho v¥k, datum vy²et°ení, informace o tom, zda se jedná o rebiopsii, pop°ípad¥ kolikátou, hladina PSA, hladina volného PSA (fPSA) nebo PSA index, celkový objem prostaty, po£et p°i biopsii odebraných vzork·, vy²et°ující léka°, výsledná diagnóza a informace o tom, zda biopsie odhalila výskyt karcinomu nebo ne. N¥které záznamy navíc obsahovaly informaci o rodinné anamnéze, PSA velocit¥, objemu tranzitorní zóny a v p°ípadech, kdy byla odhalena rakovina prostaty také hodnotu tzv. Gleasonova skóre. Nejprve bylo pot°eba tuto vstupní tabulku upravit tak, aby byla pouºitelná k samotné analýze. Tato úprava zahrnovala zejména vy£len¥ní vºdy jen prvního záznamu u kaºdého pacienta, pokud bychom totiº v souboru ponechali v²echna vy²et°ení pacient·, tedy i záznamy o rebiopsiích, byl by poru²en p°edpoklad nezávislosti pozorování. Dále bylo pot°eba odstranit odlehlá a nesmyslná pozorování a omezit n¥které vstupní prom¥nné jen na zvolené intervaly, které byly stanoveny jednak na základ¥ rozloºení hodnot v základním souboru, ale také na základ¥ zkoumané problematiky. Takto jsme získali soubor 987 pozorování, která zahrnovala první vy²et°ení muº· ve v¥ku od 45 do 80 let, kterým byla nam¥°ena hladina PSA men²í neº 20 ng/ml a jejichº objem prostaty nep°esáhl hranici 150 ml. Ze souboru byla zcela vynechána poloºka obsahují údaje o rodinné anamnéze, která byla vypln¥na jen u velmi malého po£tu pozorování a také PSA velocita, která se vztahuje k opakovaným vy²et°ením. asto také chyb¥l údaj o objemu tranzitorní zóny, tato volná místa sice byla nahrazena nulami, v dal²í analýze v²ak tranzitorní zóna pouºita nebyla. Naproti tomu chyb¥jící údaje o hladin¥ volného PSA resp. PSA indexu ²ly snadno dopo£ítat ze vztahu SA a v analýze tedy pouºity být mohly. f P SA = index · P SA resp. index = fPPSA
6.1
Hledání modelu
Vhodný model lze hledat jak vzestupn¥, postupným p°idáváním nových regresor· a jejich interakcí a testováním jejich významnosti, tak i sestupn¥ postupným zjednodu²ováním modelu a testováním podmodel· za pomoci deviancí. P°i odhadu parametr· model· jsem pouºila software R, který na základ¥ p°íkazu 60
model <- glm(y∼x1 +x2 +..., family='binomial') tyto odhady jak spo£ítá a ur£í jejich p°esnosti, tak i otestuje jejich významnost a vy£íslí dal²í s modelem spojené statistiky. Poloºka family='binomial' avizuje práci s binomicky rozd¥lenou závislou prom¥nnou a tedy pouºití logitové transformace. Výchozím modelem p°i vzestupném p°ístupu byl model bez interakcí
ln
π(x) = β0 + β1 · vek + β2 · objem + β3 · PSA, 1 − π(x)
ve kterém π(x) zna£í pravd¥podobnost, ºe bude výsledek biopsie p°i dané hodnot¥ vysv¥tlujících prom¥nných x pozitivní a potvrdí se tedy podez°ení na výskyt karcinomu prostaty. V tomto modelu se v²echny parametry významn¥ li²ily od nuly. K tomuto modelu jsem následn¥ postupn¥ p°idávala jak dvojnou interakci zárove¬ mezi v¥kem, objemem i hladinou PSA, tak i v²echny interakce jednoduché, ºádná z nich se v²ak neukázala jako statisticky významná. Pokra£ovala jsem tedy s modelem bez interakcí, který jsem dopl¬ovala o dal²í prom¥nné a to hodnotu volného PSA, PSA index, densitu, která je dána jako podíl PSA a celkového objemu a také po£et odebraných vzork·. Z t¥chto veli£in se jako významné jevily jen hladina volného PSA a PSA index, které ov²em svoji významnost ztrácí, jsou-li v modelu sou£asn¥. Vzhledem ke snaz²í interpretovatelnosti jsem do modelu p°idala volné PSA. ádná z interakcí s p·vodními t°emi veli£inami významná nebyla, rozhodla jsem se tedy dále pracovat jen s modelem
ln
π(x) = β0 + β1 · vek + β2 · objem + β3 · PSA + β4 · fPSA. 1 − π(x)
V²echny testované modely jsou uvedeny v souboru hledani_modelu.R. Sestupný p°ístup k hledání modelu lze demonstrovat p°i °e²ení problému, zda do modelu, který obsahuje prom¥nné v¥k, celkový objem prostaty, celková hladina PSA v krvi a mnoºství jeho volné frakce, zahrnout i interakce t¥chto veli£in. Tentokrát tedy vyjdeme z nejsloºit¥j²ího modelu obsahujícího trojnou interakci v²ech veli£in, v²echny interakce dvojné i jednoduché a samoz°ejm¥ i v²echny veli£iny samostatn¥. Tento model následn¥ zjednodu²íme odebráním trojné interakce a platnost podmodelu testujeme podle kapitoly 5.3.1. Rozdíl deviancí obou model·, které lze v R zobrazit p°íkazem model$deviance nebo vypo£ítat za pomoci skriptu deviance.R, srovnáváme s 0,95% kvantilem χ2 rozd¥lení o jednom stupni volnosti, protoºe odebrání trojné interakce znamená úbytek jednoho parametru. Hypotézu o platnosti podmodelu zamítáme, pokud je rozdíl deviancí v¥t²í neº uvaºovaný kvantil. V tomto p°ípad¥ byla deviance bohat²ího modelu rovna 973,60 a deviance podmodelu 977,87, pro jejich rozdíl tedy platilo:
977, 87 − 973, 60 = 4, 27 > χ21 (0, 95) = 3, 84. Z toho hlediska je tedy lep²ím modelem model bohat²í. Takový model je ale velmi t¥ºko interpretovatelný a mnohé z parametr·, které obsahuje nejsou významné. 61
Je tedy pot°eba pokra£ovat v jeho zjednodu²ování. Dal²í moºností je rozloºení trojné interakce jen na t°i dvojné, v tomto p°ípad¥ jiº odebíráme dva parametry, jeden, který odpovídal trojné interakce a jeden, který byl na míst¥ vypu²t¥né dvojné interakce, nyní tedy srovnáváme s kvantilem χ2 rozd¥lení o dvou stupních volnosti, který je roven 5,99. V²echny takto vzniklé podmodely, jsou spolu s jejich testovými statistikami a hodnotami p-value, kterou lze získat za pomoci p°íkazu 1-pchisq(x,df) uvedeny v Tabulce 6. Model 1 2 3a 3b 3c 3d 4a 4b 4c 5a 5b 6 7a 7b 7c 7d 7e 7f 8a 8b 8c 8d 8e 9a 9b 9c 9d 10a 10b 10c 11a 11b 12 13
v¥k*objem*PSA*fPSA v*o*p+v*o*f+v*p*f+o*p*f v*o*p+v*o*f+v*p*f v*o*p+v*o*f+o*p*f v*o*p+v*p*f+o*p*f v*o*f+v*p*f+o*p*f v*o*f+v*p*f+o*p v*p*f+o*p*f+v*o v*o*f+o*p*f+v*p v*p*f+o*p+o*f+v*o o*p*f+v*p+v*f+v*o o*p+o*f+p*f+v*p+v*f+v*o o*f+p*f+v*p+v*f+p*o o*p+p*f+o*f+v*f+v*o o*p+o*f+v*p+p*f+v*o v*p+o*f+p*f+v*f+v*o o*p+v*f+p*f+v*p+v*o o*p+o*f+v*o+v*p+v*f o*f+o*p+v*f+p*f o*p+p*f+o*f+v*o v*p+o*f+p*f+v*o o*p+p*f+v*p+v*o o*p+o*f+v*p+v*o v*p+o*p+p*f v*o+o*p+p*f v*o+v*p+p*f v*o+v*p+o*p+fPSA o*p+p*f+v¥k v*o+p*f v*o+o*p+fPSA v*o+PSA+fPSA p*f+v¥k+objem v¥k+objem+PSA+fPSA v¥k+objem+PSA
Deviance
df
973,6 977,87 981,2 978,32 978,05 977,87 981,23 978,14 978,4 981,44 978,43 981,56 983,01 981,81 981,57 981,9 981,57 982,38 983,4 981,87 981,93 981,58 982,42 983,63 981,92 982 982,46 983,68 982,46 982,62 982,79 983,94 984,4 988,53
971 972 973 973 973 973 974 974 974 975 975 976 977 977 977 977 977 977 978 978 978 978 978 979 979 979 979 980 980 980 981 981 982 983
Srov. modely 2-1 3a-1 3b-1 3c-1 3d-1 4a-3d 4b-3d 4c-3d 5a-4b 5b-4b 6-5b 7a-6 7b-6 7c-6 7d-6 7e-6 7f-6 8a-7c 8b-7c 8c-7c 8d-7c 8e-7c 9a-8d 9b-8d 9c-8d 9d-8d 10a-9b 10b-9b 10c-9b 11a-10b 11b-10b 12-11a 13-12
Rozdíl dev. 4,27 7,6 4,72 4,45 4,27 3,36 0,27 0,53 3,3 0,29 3,13 1,45 0,25 0,01 0,34 0,01 0,82 1,83 0,3 0,36 0,01 0,85 2,05 0,34 0,42 0,88 1,76 0,54 0,7 0,33 1,48 1,61 4,13
Rozdíl df 1 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
p-value
AIC
0,039 0,022 0,094 0,108 0,118 0,067 0,603 0,467 0,069 0,59 0,077 0,229 0,617 0,92 0,56 0,92 0,365 0,176 0,584 0,549 0,92 0,357 0,152 0,56 0,517 0,348 0,185 0,462 0,403 0,566 0,224 0,204 0,042
1005,6 1007,9 1009,2 1006,3 1006,1 1005,9 1007,2 1004,1 1004,4 1005,4 1002,4 1003,6 1003 1001,8 1001,6 1001,9 1001,6 1002,4 1001,4 999,87 999,93 999,58 1000,4 999,63 997,92 998 998,46 997,68 996,46 996,62 994,79 995,94 994,4 996,53
Tabulka 6: Zjednodu²ení model· a jejich srovnání prost°ednictvím deviance Jak je z tabulky vid¥t, ve zjednodu²ování pokra£ujeme a do okamºiku, kdy pracujeme op¥t jen s modelem 12
ln
π(x) = β0 + β1 · vek + β2 · objem + β3 · PSA + β4 · fPSA, 1 − π(x)
ve kterém sice jsou v²echny parametry významné, hodnota statistiky ov¥°ující významnost parametru p°íslu²ejícího volnému PSA je v²ak hrani£ní (P-val=0,053), 62
jak bude ukázáno pozd¥ji. Nabízí se tedy volné PSA z modelu vypustit a pracovat jen s modelem
ln
π(x) = β0 + β1 · vek + β2 · objem + β3 · PSA, 1 − π(x)
který je v tabulce uveden jako model 13. Otestujme tedy, zda toto zjednodu²ení vylep²í vlastnosti modelu. Deviance bohat²ího modelu 12 je 984,4, po£et stup¬· volnosti, který tomuto modelu p°íslu²í je dán jako rozdíl celkového po£tu pozorování a po£tu v modelu odhadovaných parametr·, tedy df=987-5=982. Deviance podmodelu 13 se zvý²ila na 988,53 p°i 983 stupních volnosti. Testová statistika má tentokrát podobu
988, 53 − 984, 4 = 4, 13 > χ21 (0, 95) = 3, 84. A£koliv tedy volné PSA není na hladin¥ významnosti 5 % významné, ze srovnání deviancí model· 12 a 13 vyplývá, ºe je lep²í jej v modelu ponechat. Nadále se tedy budeme zabývat jen modelem 12. V²echny modely jsou uvedeny v souboru hledani_modelu_sestup.R
6.2
Odhad parametr·
6.2.1 Po£áte£ní odhady Abychom mohli pomocí Newton-Raphsonovy metody ur£it odhady regresních parametr· β , musíme mít alespo¬ p°ibliºnou p°edstavu o tom, ve kterém bod¥ s iteracemi za£ít. Po£áte£ní odhady ur£íme pomocí kapitoly 5.2.1. Datový soubor se skládá celkem z n = 987 pozorování a rakovina byla potvrzena ve y = 230 p°ípadech, po£áte£ní odhad absolutního £lenu je tedy (0)
β0 = ln
1
230 987 230 − 987
= ln
230 = ln 0, 2987 = −1, 191. 987 − 230
Nyní ur£eme po£áte£ní odhady zbylých parametr·. Protoºe se v¥k pacient· v souboru pohybuje mezi 45 a 80 lety, bude po£áte£ní odhad parametru β1 (0)
β1 =
4 4 = = 0, 11. 80 − 45 35
Objem prostaty jsme jiº na za£átku práce s daty omezili shora hodnotou 150 ml. Vykleslíme-li si navíc histogram s hodnotami objemu v p°ípadech, kdy byl výsledek biopsie pozitivní, který je zobrazen na Obrázku 5, vidíme, ºe s rostoucím objemem po£et pozitivních výsledk· ubýval a to na rozdíl od v¥ku, u kterého 63
takový trend patrný není. Lze tedy p°edpokládat, ºe pravd¥podobnost pozitivní biopsie bude s rostoucí velikosti prostaty klesat. Pro po£áte£ní odhad parametru β2 tedy platí 4 4 (0) =− = −0, 027. β2 = − 150 − 0 150
Histogram of objem
Frequency
30
40
0
0
10
10
20
20
30
Frequency
40
50
50
60
60
70
Histogram of vek
50
55
60
65
70
75
80
20
40
60
80
100
vek
objem
Histogram of PSA
Histogram of fPSA
120
140
40
Frequency
30
0
0
10
20
20
Frequency
40
60
50
80
60
45
5
10
15
20
0
PSA
1
2
3
4
5
fPSA
Obrázek 5: Rozloºení vysv¥tlujících prom¥nných v p°ípadech pozitivních biopsií Dal²ím regresorem je hladina PSA, i tu jsme hned v úvodu omezili hodnotou 20 ng/ml. P°i pohledu na histogram zachycující rozloºení PSA se m·ºe zdát, ºe jsou jeho hodnoty op¥t více koncetrovány v niº²ích polohách. V tomto p°ípad¥ je ale p°í£inou spí²e ta skute£nost, ºe hodnoty PSA nad 10 ng/ml jsou hodnotami 64
jiº do jisté míry extrémními a tyto p°ípady jsou tedy v datovém souboru málo zastoupeny. Po£áte£ní odhad parametru β3 je (0)
β3 =
4 4 = = 0, 2. 20 − 0 20
Posledním regresorem je volné PSA, jehoº hladina se v datovém souboru pohybuje v intervalu od 0 do 7 ng/ml. Z histogramu je vid¥t, ºe s r·stem hladiny volné frakce PSA v krvi prudce ubývá p°ípad· pozitivní biopsie. Podobn¥ jako v p°ípad¥ celkového objemu tedy o£ekáváme nep°ímou závislost mezi pravd¥podobností pozitivní biopsie a volným PSA a po£áte£ní odhad posledního parametru tak bude (0)
β4 = −
4 4 = − = −0, 571. 7−0 7
Newton-Raphsonova metoda bude tedy za£ínat v bod¥
β (0) = (−1, 191; 0, 11; −0, 027; 0, 2; −0, 571)0 .
6.2.2 Iterace Newton-Raphsonovy metody Pro výpo£et odhad· parametru β v jednotlivých krocích i odhadu výsledného pouºijeme software R a skript, který je uloºen pod názvem newton_raphson.R. Jelikoº p°i startu ve spo£ítaném bod¥ β (0) = (−1, 191; 0, 11; −0, 027; 0, 2; −0, 571)0 algoritmus nekonverguje, je p°i jeho iniciaci pot°eba pouºít bod jiný, blízký bodu odhadnutému. Za£n¥me tedy v bod¥, který nep°edpokládá závislost výsledku biopsie na ºádném z regresor· v¥k, objem prostaty, PSA a fPSA
β (0) = (−1, 191; 0; 0; 0; 0)0 . Výsledky jednotlivých iterací jsou uvedeny v Tabulce 7: t (t) β0 (t) β1 (t) β2 (t) β3 (t) β4
0 -1,191 0 0 0 0
1 -3,875 0,050 -0,023 0,155 -0,315
2 -3,911 0,053 -0,027 0,144 -0,304
3 -3,917 0,053 -0,027 0,144 -0,305
4 -3,917 0,053 -0,027 0,144 -0,305
Tabulka 7: Odhady parametr· v jednotlivých iteracích Z tabulky je vid¥t, ºe mezi odhady z t°etí a £tvrté iterace jiº není prakticky ºádný rozdíl a algoritmus tedy m·ºeme ukon£it po £tvrté iteraci. Druhou 65
k g(t)
0 -0,05 292,07 -1660,6 187,49 -10,81
1 -18 -1207,64 -1081,9 -153,51 -24,72
2 -0,85 -53,94 -64,24 -6,51 -1,31
3 -0,01 -0,42 -0,52 -0,05 -0,01
4
-4,69·10 -3·10−5 -3,72·10−5 -3,49·10−6 -7,26·10−7
−7
Tabulka 8: Hodnoty gradientu v jednotlivých iteracích moºností jak algoritmus ukon£it, je sledovat vývoj gradientu; ten by m¥l být v bod¥ maximáln¥ v¥rohodného odhadu nulový. Tabulka 8 tento vývoj zachycuje. Ve £tvrté iteraci je gradient jiº tém¥° nulový, maximáln¥ v¥rohodný odhad parametr· je tedy
ˆ = (−3, 91651; 0, 05299; −0, 02737; 0, 14402; −0, 30509)0 β a výsledný model pak
ln
π(x) = −3, 917 + 0, 053 · vek − 0, 027 · objem + 0, 144 · PSA − 0, 305 · fPSA. 1 − π(x)
Ze £tvrté iterace také podle vztahu uvedeného v kapitole 5.2.2 získáme odhad varian£ní matice odhadu 0, 6395 −1, 0038 · 10−2 −2, 7766 · 10−4 −1, 5490 · 10−3 0, 0226 −4 −6 −5 −0, 0100 1, 7079 · 10 −7, 1210 · 10 −2, 3133 · 10 −0, 0003 −5 −6 ˆ = −0, 0003 −7, 1211 · 10−6 2, 2988 · 10 6, 8520 · 10 −0, 0003 var cβ −6 −4 −0, 0015 −2, 3133 · 10−5 6, 8520 · 10 7, 9389 · 10 −0, 0025 −4 −4 −3 0, 0226 −3, 0171 · 10 −2, 8440 · 10 −2, 4691 · 10 0, 0248 V Tabulce 9 je²t¥ srovnejme jak moc se po£áte£ní odhad li²il od odhadu kone£ného. Je vid¥t, ºe po£áte£ní odhady nebyly p°íli² vzdáleny od t¥ch maxiParametr β0 β1 β2 β3 β4
(0)
Po£áte£ní odhad βj -1,191 0,11 -0,027 0,2 -0,571
Max. v¥rohodný odhad βˆj -3,917 0,053 -0,027 0,144 -0,305
Tabulka 9: Srovnání po£áte£ního a maximáln¥ v¥rohodného odhadu parametru β 66
máln¥ v¥rohodných a podávaly tedy dobrou informaci o tom, kde se maximáln¥ v¥rohodné odhady vyskytují a kde za£ít itera£ní mechanismus. V²imn¥me si zejména po£áte£ního odhadu parametru β2 , který odpovídá celkovému objemu prostaty a který je prakticky roven maximáln¥ v¥rohodnému odhadu βˆ2 .
6.2.3 Intervaly spolehlivosti a testování významnosti parametr· Nejprve ur£eme Waldovy intervaly spolehlivosti parametr· βj , j = 0, 1, . . . , 4 na hladin¥ α = 0, 05. Obecn¥ jsou tyto intervaly dány vztahem q ˆ c βˆj . βj ± u(0,975) var Kon den£ní intervaly jednotlivých parametr· tedy jsou: √ β0 : h −3, 9165 ± 1, 96 0, 6395 i = h −5, 4839 ; −2, 3492 i
β1 : h
0, 053
p ± 1, 96 1, 7079 · 10−4 i = h
0, 0274
;
0, 0786
i
p β2 : h −0, 0274 ± 1, 96 2, 2988 · 10−5 i = h −0, 0368 ; −0, 0180 i β3 : h
0, 144
p ± 1, 96 7, 9389 · 10−4 i = h
β4 : h −0, 3051 ±
√ 1, 96 0, 0248
;
0, 1992
i
i = h −0, 6138 ;
0, 0037
i
0, 0888
V²imn¥me si, ºe interval spolehlivosti pro parametr β4 obsahuje nulu, to zna£í, ºe význam volného PSA v modelu není p°íli² velký. Ov¥°m¥ tedy významnost parametr· pomocí Waldovy statistiky. Testujeme postupn¥ hypotézy H0j : βj = 0, j = 0, 1, . . . , 4 a druhou mocninu testové statistiky srovnáváme s 95% kvantilem χ2 rozd¥lení o jednom stupni volnosti, který je roven 3,84. Hodnoty statistik, které vypo£ítáme pomocí vztahu uvedeného v kapitole 5.2.3 a p°íslu²né p-value, jsou uvedeny v Tabulce 10. Z tabulky je op¥t vid¥t, ºe volné PSA není p°íli² významné, v modelu jej ov²em ponecháme, protoºe p-value je pom¥rn¥ blízko uvaºované hladin¥ spolehlivosti 5 % a p°i testování model· v úvodu kapitoly se jako lep²í jevil ten, který volnou frakci PSA obsahoval. V²em zbývajícím parametr·m p°íslu²ejí velmi malé p-value a hypotézu o jejich nulovosti tedy zamítáme. Nyní ur£eme odhad pravd¥podobnosti, ºe biopsie potvrdí p°ítomnost rakoviny prostaty, p°i dané hodnot¥ regresor· x = (1, x1 , x2 , x3 , x4 ). Ten je pro libovolný vektor x roven
π ˆ (x) =
e−3,917+0,053x1 −0,027x2 +0,144x3 −0,305x4 . 1 + e−3,917+0,053x1 −0,027x2 +0,144x3 −0,305x4 67
Parametr β0 β1
β2 β3 β4
Testová statistika −3,9165 √ = −4, 898 0,6395 √ 0,053 = 4, 055 −4
1,7079·10 −0,0274 √
2,2988·10−5 √ 0,144
P-value 9, 699 · 10−7 5, 014 · 10−5
= −5, 717
1, 136 · 10−8
= 5, 111
3, 199 · 10−7
7,9389·10−4 −0,3051 √ = −1, 937 0,0248
5, 278 · 10−2
Tabulka 10: Testování parametr· pomocí Waldovy statistiky
Pro ilustraci nyní uvaºujme dva p°ípady, kdy máme dva muºe ve v¥ku 60ti let, kterým byl nam¥°en objem prostaty 50 ml a jejichº hladina volné frakce PSA byla 2 ng/ml. Oba p°ípady se li²í jen hladinou celkového PSA, ta je v prvním p°ípad¥ 5 ng/ml a v druhém 10 ng/ml. V prvním p°ípad¥ je vektor regresor·
x1 = (1, 60, 50, 5, 2) a v druhém
x1 = (1, 60, 50, 10, 2). Odhady obou pravd¥podobností jsou
π ˆ (x1 ) =
e−3,917+0,053·60−0,027·50+0,144·5−0,305·2 = 0, 1196618 1 + e−3,917+0,053·60−0,027·50+0,144·5−0,305·2
a
e−3,917+0,053·60−0,027·50+0,144·10−0,305·2 = 0, 2183084. 1 + e−3,917+0,053·60−0,027·50+0,144·10−0,305·2 Vidíme tedy, ºe zvý²ením hladiny PSA v krvi se pravd¥podobnost pozitivního výsledku tém¥° zdvojnásobila. π ˆ (x2 ) =
Abychom mohli stanovit 95% intervaly spolehlivosti pro tyto pravd¥podobnosti, musíme nejprve ur£it kon den£ní intervaly jejich logitových transformací, které jsou obecn¥ rovny q π ˆ (x) ˆ 0. ln ± 1, 96 xvar c βx 1−π ˆ (x) Tento interval je v prvním p°ípad¥ roven q 0, 1196618 ˆ 60, 50, 5, 2)0 = ln ± 1, 96 (1, 60, 50, 5, 2)var c β(1, 1 − 0, 1196618
= h−2, 395165; −1, 596110i
68
a v druhém
q 0, 2183084 ˆ 60, 50, 10, 2)0 = ± 1, 96 (1, 60, 50, 10, 2)var ln c β(1, 1 − 0, 2183084 = h−1, 5901734; −0, 9609301i. Krajní body interval· spolehlivosti pro samotné pravd¥podobnosti nyní získáme ex dosazením krajních bod· spo£ítaných interval· do výrazu 1+e x . Interval spolehlivosti pro π(x1 ) je tedy −2,395165 e e−1,596110 = h0, 0835422; 0, 1685261i ; 1 + e−2,395165 1 + e−1,596110 a pro π(x2 ) je
e−0,9609301 e−1,5901734 ; 1 + e−1,5901734 1 + e−0,9609301
= h0, 1693595; 0, 2766920i.
I intervaly spolehlivosti lze v R spo£ítat pomocí skriptu newton_raphson.R, který lze pouºít také pro testování významnosti parametr·.
6.3
Ov¥°ení kvality modelu
První moºností jak ohodnotit kvalitu modelu je pouºití deviance, která se p°i seskupených datech vypo£ítá podle vztahu uvedeného v kapitole 5.3.1. Protoºe ale pracujeme s daty neseskupenými, získáme devianci ze vztahu
D=2
987 X i=1
987
X 1 − yi yi yi ln + 2 (1 − yi ) ln . π ˆi 1−π ˆi i=1
Po dosazení v²ech pozorování získáme za pomoci vhodného softwaru hodnotu
D = 984, 4, která je uvedena i v Tabulce 6. Kvalitu modelu lze posuzovat i pomocí Akaikeho informa£ního kritéria, pro které je v kapitole 5.3.4 uveden vztah
AIC = −2[L(µ, ˆ y) − p], ten se prakticky rovná devianci modelu, zv¥t²ené o dvojnásobek po£tu odhadovaných parametr·. Hodnota tohoto kritéria je tedy pro uvaºovaný model
AIC = 984, 4 + 2 · 5 = 994, 4. 69
Jak deviance, tak i hodnota Akaikeho informa£ního kritéria se mohou zdát pom¥rn¥ velké, to je zp·sobeno zejména velkou redukcí v parametrech ve srovnání se saturovaným modelem, která je ale vzhledem k pouºitelnosti a snadné interpretovatelnosti modelu nutná. Z Tabulky 6 je navíc vid¥t, ºe uvaºovaný model bez interakcí má nejniº²í hodnotu AIC ze v²ech model· obsahujících v¥k, celkový objem, a hladinu volného a celkového PSA v krvi a z tohoto hlediska se tedy op¥t jeví jako optimální. V R lze devianci a AIC získat p°íkazem model$deviance resp. model$aic nebo pomocí skriptu deviance.R. Abychom mohli ur£it McFadden·v a Nagelkerk·v koe cient determinace, musíme nejprve spo£ítat devianci nulového modelu. D0 = −2 230 ln 230 + (987 − 230) ln(987 − 230) − 987 ln 987 = 1071, 7 McFaddn·v koe cient determinace je
RL2 = 1 − a Nagelkerk·v 2 RN
984,4−1071,7 987
= 0, 128. −1071,7 1 − e 987 I koe cienty determinace lze spo£ítat pomocí skriptu deviance.R, devianci nulového modelu navíc získáme i p°íkazem model$null.deviance.
6.4
=
1−e
984, 4 = 0, 081 1071, 7
Interpretace výsledk·
V p°edchozích kapitolách byl pro odhad pravd¥podobnosti pozitivního výsledku biopsie odvozen model
ln
π(x) = −3, 917 + 0, 053 · vek − 0, 027 · objem + 0, 144 · PSA − 0, 305 · fPSA, 1 − π(x)
pomocí deviance a Akaikeho informa£ního kritéria byla posouzena jeho kvalita a byla také otestována významnost jednotlivých parametr·. Nyní tedy vysv¥tleme, co hodnoty jednotlivých parametr· znamenají. Za£n¥me absolutním £lenem β0 . Jeho odhad je βˆ0 = −3, 917 a jak jiº bylo °e£eno v kapitole 5.1, je roven logaritmu ²ance na pozitivní biopsii oproti negativní v p°ípad¥, kdy jsou v²echny regresory nulové, samotná ²ance by pak byla rovna 0,0199 a negativní biopsie by tak byla 50krát pravd¥podobn¥j²í neº pozitivní. V tomto p°ípad¥ je v²ak p°edpoklad nulovosti regresor· nesmyslný. Nejen ºe bychom uvaºovali muºe nulového v¥ku s nulovou hladinou PSA a tedy i 70
fPSA, ale také s nulovým objemem prostaty a tedy bez moºnosti vzniku zhoubného novotvaru na tomto orgánu. Abychom tomuto paradoxu zabránili, m·ºeme pouºít zmín¥né p°e²kálování regresor· tak, ºe nula bude tentokrát odpovídat jejich pr·m¥rným hodnotám. P°i odhadu parametr· z takto upravených dat se zm¥ní jen hodnota parametru β0 , jehoº odhad je tentokrát roven β˜0 = −2, 205. Za pouºití tohoto odhadu m·ºeme ur£it ²anci na pozitivní biopsii oproti negativní v p°ípad¥ pacienta ve v¥ku 62,5 let, s objemem prostaty 76 ml, hladinou celkového PSA v krvi 10,775 ng/ml a volného PSA 3,514 ng/ml. Tato ²ance je
π ˆ (x) = e−2,205 = 0, 11, 1−π ˆ (x) negativní výsledek biopsie je tedy v tomto p°ípad¥ asi 9krát pravd¥podobn¥j²í n¥º pozitivní. Stejnou hodnotu bychom získali i kdybychom vektor regresor· x = (62, 3; 76; 10, 775; 3, 514)0 dosadili do p·vodního modelu, ve kterém by pro ²anci platilo
π ˆ (x) = e−3,917+0,053·62,5−0,027·76+0,144·10,775−0,305·3,514 = e−2,17667 = 0, 11. 1−π ˆ (x) Parametr β1 se interpetuje jako logaritmus pom¥ru ²ancí mezi pacienty se stejnými hodnotami objemu, PSA i fPSA, ale mezi nimiº je v¥kový rozdíl jeden rok. Protoºe je odhad parametru β1 kladný, s r·stem v¥ku poroste i ²ance na pozitivní biopsii a tedy pravd¥podobnost výskytu rakoviny prostaty, coº odpovídá tvrzení kapitoly 2.4.1, ºe v¥k je jedním z rizikových faktor· pro vznik tohoto druhu novotvaru. Odhad pom¥ru ²ancí p°i rozdílu jeden rok je
d = e0,053 = 1, 0544, OR ²ance, ºe bude výsledek biopsie pozitivní tedy s kaºdým rokem rozdílu vzroste o 5,44 %. Zajímá-li nás, p°i jakém v¥kovém rozdílu je tato ²ance dvojnásobná, °e²íme tuto rovnici en·0,053 = 2. Odtud vypo£ítáme, ºe se ²ance zdvojnásobí p°i rozdílu
n=
ln 2 . = 13 let. 0, 053
Obdobn¥ jako v p°ípad¥ odhadu parametru β1 , p°i interpretaci βˆ2 uvaºujeme dva pacienty se stejnými hodnotami regresor·, kdy jediný rozdíl je tentokrát v objemu prostaty, který se li²í o 1 ml. Jiº na základ¥ toho, ºe je odhad parametru β2 záporný, m·ºeme usuzovat, ºe s r·stem objemu prostaty bude klesat ²ance na pozitivní biopsii a s ní i pravd¥podobnost odhalení karcinomu prostaty. Odhad pom¥ru ²ancí pacienta s v¥t²ím objemem a pacienta s objemem men²ím je
d = e−0,027 = 0, 9734, OR 71
²ance tedy p°i nár·stu objemu o 1 ml klesne asi o 2,66 %. Tento záv¥r m·ºe mimo jiné souviset se zp·sobem provád¥ní biopsíí, popsaným v kapitole 2.5.4 a tím, ºe ve v¥t²ích prostatách je napadená tkᬠh·°e odhalitelná. Namísto zdvojnásobení se tentokrát m·ºeme zajímat o to, p°i jakém rozdílu mezi objemy obou pacient· bude ²ance na pozitivní biopsii pacienta s v¥t²í prostatou polovi£ní oproti ²anci druhého pacienta. Tentokrát tedy °e²íme tuto rovnici
1 en(−0,027) = , 2 jejímº °e²ením je
ln 12 n= = 25, 67 ml. −0, 027
ance tedy bude polovi£ní p°i rozdílu 25,67 ml. Odhad parametru β3 je op¥t kladný, coº zna£í, ºe vy²²í hladina PSA v krvi povede k vy²²í ²anci na pozitivní výsledek biopsie a tedy i pravd¥podobnosti diagnostikování rakoviny prostaty. Tento záv¥r odpovídá kapitole 2.5.1 a tomu, ºe zvý²ené hodnoty tohoto ukazatele jsou v praxi hojn¥ vyuºívaným indikátorem rizika výskytu karcinomu prostaty. Odhad pom¥ru ²ancí je
d = e0,144 = 1, 1549 OR a pacient, který se od jiného li²í jen hladinou PSA, kterou má o 1 ng/ml vy²²í, má zárove¬ i o 15,49 % vy²²í ²anci na pozitivní biopsii. Tato ²ance bude dvojnásobná p°i rozdílu ln 2 = 4, 81 ng/ml. n= 0, 144 Posledním parametrem je parametr β4 , který odpovídá volnému PSA. Jeho odhad je záporný a p°i vy²²ích hodnotách volné frakce PSA v krvi tak lze o£ekávat niº²í ²anci na pozitivní biopsii a tedy i pravd¥podobnost výskytu rakoviny prostaty. Jak bylo uvedeno v kapitole 2.5.1, zvý²ené riziko rakoviny souvisí s vysokými hodnotami PSA vázaného na alfa-1-antichymotrypsin a tedy nízkými hodnotami volného PSA. Jaký je vztah mezi ²ancemi na pozitivní biopsii dvou pacient·, z nichº jednomu byla nam¥°ena o 1 ng/ml vy²²í hladina volného PSA a jejichº hodnoty zbylých regresor· jsou shodné, udává odhad pom¥ru ²ancí
d = e−0,305 = 0, 7371, OR z kterého vidíme, ºe rozdíl ve volném PSA o velikosti 1 ng/ml, znamená sníºení ²ance o 26,29 %. ance bude polovi£ní, p°i rozdílu v hladinách volného PSA o velikosti ln 12 = 2, 27 ng/ml. n= −0, 305 72
P°i pohledu na absolutní hodnoty odhad· parametr· a následné vy£íslení jednotlivých pom¥r· ²ancí se m·ºe zdát, ºe pravd¥podobnost pozitivní biopsie nejvíce ovliv¬uje hladina volného PSA, protoºe jednotkový rozdíl mezi hladinami fPSA dvou pacient· s jinak stejnými parametry znamená, ºe ²ance na pozitivní biopsii pacienta s vy²²í hladinou bude zhruba o 26 % men²í neº ²ance pacienta s hladinou vy²²í. Je ale pot°eba si uv¥domit, ºe jednotlivé regresory nejsou m¥°eny na stejn¥ velkých intervalech a rozdíl o jednotku tak není pro v²echny rovnocenný. Pro lep²í porovnatelnost významu jednotlivých regresor· ur£eme, jaký bude pom¥r ²ancí p°i rozdílu, který odpovídá jedné setin¥ intervalu, na kterém jej m¥°íme. Výsledné hodnoty jsou uvedeny v Tabulce 11. Parametr β1 β2 β3 β4
Interval h45, 80i h0, 150i h0, 20i h0, 7i
Rozdíl xj 80−45 = 0, 35 100 150−0 = 1, 5 100 20−0 = 0, 2 100 7−0 = 0, 07 100
Pom¥r ²ancí % e0,35·0,053 = 1, 01872 1,87 % e1,5(−0,027) = 0, 96031 -3,97 % e0,2·0,144 = 1, 02922 2,92 % 0,07(−0,305) e = 0, 97887 -2,11 %
Tabulka 11: Srovnání pom¥r· ²ancí p°i rozdílu odpovídajícím 1 % intervalu Z tabulky je vid¥t, ºe nejv¥t²í zm¥nu pom¥ru ²ancí vyvolá 1% rozdíl v objemu, který odpovídá 1,5 ml. Naopak jako nejmén¥ významný se z tohoto hlediska jeví procentní rozdíl ve v¥ku.
6.5
Srovnání s modelem pro hyperplazii prostaty
Nej£ast¥j²í diagnózou ve zpracovávaných datech byla hyperplazie prostaty, odhadn¥me pro ni parametry ve stejném modelu, jako byl pouºit pro rakovinu prostaty, a oba odhady srovnejme. P°íkazem model<-glm(hyper∼vek+objem+PSA+fPSA,family='binomial') byl v R odhadnut model
ln
π(x) = 1, 928 − 0, 039 · vek + 0, 007 · objem − 0, 079 · PSA + 0, 201 · fPSA, 1 − π(x)
ve kterém byly v²echny parametry s vyjímkou β4 významn¥ nenulové. hyper je op¥t alternativn¥ rozd¥lená prom¥nná, která je rovna jedné v p°ípad¥, kdy biopsie potvrdila hyperplazii a nulová ve zbylých p°ípadech. Tento model lze získat také pomocí skriptu hyperplasie.R. Parametry modelu pro rakovinu prostaty a pro hyperplazii srovnává Tabulka 12.
73
Parametr β0 β1 β2 β3 β4
Rakovina -3,917 0,053 -0,027 0,144 -0,305
Hyperplazie 1,928 -0,039 0,007 -0,079 0,201
Tabulka 12: Srovnání odhadu parametr· v modelu s rakovinou prostaty a modelu s hyperplazií Jak je z tabulky vid¥t, odhady parametr· v jednotlivých modelech mají vºdy opa£né znaménko a zatímco pravd¥podobnost výskytu rakoviny prostaty je tím vy²²í, £ím vy²²í je v¥k pacienta, £i hladina PSA v krvi a £ím niº²í je objem prostaty, v p°ípad¥ pravd¥podobnosti výskytu hyperplazie je tomu naopak. V tomto sm¥ru lze tedy hyperplazii prostaty vnímat jako opak karcinomu.
74
7
Záv¥r
Úvod práce byl v¥nován problematice rakoviny prostaty, jejímu výskytu, úmrtnosti a také diagnostickým prost°edk·m. Záv¥r druhé kapitoly tvo°il odhad rizika, ºe muº onemocní rakovinou prostaty nebo ºe v d·sledku této choroby zem°e. Následovaly kapitoly 3 a 4, ve kterých byly odvozeny a na p°íkladech objasn¥ny základní pojmy, které byly následn¥ vyuºity v kapitole páté, v¥nované logistické regresi. está kapitola se jiº zabývala zpracováním záznam· z vy²et°ení na Urologické klinice Fakutní nemocnice Olomouc. Po úvodní úprav¥ datového souboru, která zahrnovala zejména odstran¥ní odlehlých £i nesmyslných pozorování, byl jak vzestupn¥ tak i sestupn¥ nalezen model, který zahrnoval jen vliv v¥ku pacienta, objemu jeho prostaty, hladiny PSA v krvi a mnoºství jeho volné frakce a v n¥mº nebyla ºádná z interakcí mezi t¥mito veli£inami významná. Práv¥ hladina volného PSA byla pon¥kud problematická z hlediska její významnosti a následného za°azení do modelu. Hodnota statistiky testující nenulovost odpovídajícího parametru totiº byla hrani£ní a sv¥d£ila spí²e pro jeho nevýznamnost. Na druhou stranu test podmodelu bez volného PSA pomocí rozdílu deviancí sv¥d£il spí²e pro jeho ponechání v modelu. Volné PSA tedy v modelu z·stalo. Jakmile jsme m¥li nalezen vhodný model, následovalo odhadování jeho parametr·. Pro iniciaci Newton-Raphsonovy metody byly nalezeny po£áte£ní (0) odhady βj , j = 0, 1, . . . , 4 a a£koliv se pozd¥ji ukázalo, ºe tyto odhady nebyly p°íli² vzdáleny od odhad· maximáln¥ v¥rohodných, metoda p°i startu z tohoto bodu nekonvergovala a bylo proto pot°eba za£ít v bod¥ jiném, který nep°edpokládal závislost pravd¥podobnosti pozitivní biopsie na ºádném z regresor· a v²echny jim p°íslu²né parametry tedy byly nulové. Nyní jiº metoda maximáln¥ v¥rohodné odhady parametr· nalezla. D·leºitou sou£ástí ²esté kapitoly je interpretace parametr·. Ta ukázala, ºe srovnáváme-li dva pacienty li²ící se pouze v¥kem resp. hladinou PSA v krvi, je pravd¥podobnost toho, ºe biopsie prokáºe výskyt karcinomu prostaty, vy²²í v p°ípad¥ pacienta s vy²²í hodnotou jednoho ze zmín¥ných regresor·. V p°ípad¥ objemu prostaty resp. volného PSA je tomu naopak. Ze srovnání vlivu, který má procentní rozdíl vºdy jednoho z regresor· na pom¥r ²ancí obou pacient· vyplynulo, ºe je tento vliv nejv¥t²í v p°ípad¥ objemu prostaty a naopak nejmen²í v p°ípad¥ v¥ku. Odhadnutý model byl také srovnán s modelem popisujícím vztah mezi benigní hyperplazií prostaty a jiº zmín¥nými veli£inami v¥k pacienta, objem prostaty a hladina volného a celkového PSA. V tomto srovnání se ukázalo, ºe rozdíl v hodnotách regresor· mezi dv¥ma pacienty má na pravd¥podobnost diagnostikování hyperplazie opa£ný efekt, neº je tomu v p°ípad¥ rakoviny. Tato pravd¥75
podobnost je tedy niº²í pro skupinu s vy²²ím v¥kem a hladinou PSA a naopak vy²²í ve skupin¥ s vy²²ím objemem prostaty a mnoºstvím volného PSA, kdy op¥t uvaºujeme, ºe se pacienti resp. skupiny li²í v hodnot¥ vºdy jen jednoho z regresor·. Pro mne osobn¥ byla práce velmi p°ínosná, jiº její první £ást vyºadující prostudování velkého mnoºství odborných £lánk·, mi pomohla lépe se orientovat nejen v problematice rakoviny prostaty, která pro mne do té doby byla oblastí zcela neznámou. Jelikoº byly prakticky v²echny výpo£ty provedeny v softwaru R, který je sice v praxi hojn¥ vyuºívaný, ale se kterým jsem dosud nepracovala, bylo pro mne cenným i kdyº místy náro£ným bonusem nau£it se v n¥m nejen provád¥t výpo£ty, ale také tvo°it jednoduché výpo£etní programy. V neposlední °ad¥ bylo velmi zajímavé seznámení se s metodou, se kterou jsem se navzdory tomu, ºe má ²iroké uplatn¥ní nejen v oblasti léka°ství, v pr·b¥hu studia setkala jen okrajov¥.
76
Literatura [1] Internetové stránky softwaru pro vizualizaci http://www.svod.cz [citováno 25. 11. 2011]
onkologických
dat
[2] Zdravotnická statistika Novotvary 2008 R, Ústav zdravotnických informací a statistiky R, Národní onkologický registr R, 2011, dostupné na http://www.uzis.cz/publikace/novotvary-2008 [citováno 1. 11. 2011] [3]
R., Ward, E., Brawley, O., Jemal, A.,, The Impact of Eliminating Socioeconomic and Racial Disparities on Premature Cancer Deaths, CA: A Cancer Journal for Clinicals, Vol. 61, Issue 4, 2011, 212-236 Siegel,
dostupné na http://onlinelibrary.wiley.com/doi/10.3322/caac.20121/pdf [citováno 13. 11. 2011] [4]
Altekruse, SF. a kol.
SEER Cancer Statistics Review, 1975-2007, Na-
tional Cancer Institute. Bethesda, MD, 2010 dostupné na http://seer.cancer.gov/csr/1975_2007/ [citováno 15. 11. 2011]
[5] Prostate Cancer Fact Sheet [online] dostupné na http://www.histoscanning.com/media/ les/Prostate_Cancer_ Fact_Sheet.pdf [citováno 15. 11. 2011] [6]
Bujdák, P., Cuninková, M., Karcinóm prostaty-trendy výskytu a rizikové faktory, Urologie pro praxi 4/2004, 169-171
dostupné na http://www.urologiepropraxi.cz/pdfs/uro/2004/04/07.pdf [citováno 5. 11. 2011]
[7]
Král, M., Vyhnánková, V., tudent, V., Bouchal, J.,
riziko karcinomu prostaty eská urologie 3/2010, 139-147
Genetické
dostupné na http://www.czechurol.cz/cislo_akt.php?cis=48 [citováno 5. 11. 2011] [8]
Houska, D., Haluzík, M., Vernerová, Z., Housová, J., Herá£ek, J.,
Obezita, adipocytokiny a karcinom prostaty, Urologie pro praxi 1/2007,
10-14 dostupné na http://www.urologiepropraxi.cz/pdfs/uro/2007/01/02.pdf [citováno 5. 11. 2011] [9]
Grepl, M., Král, M., Bure²ová, E., tudent, V.,
Agresivni karcinom
prostaty u pacient· s nízkým PSA, eská urologie 1/2010, 48-54
dostupné na http://www.czechurol.cz/cislo_akt.php?cis=46 [citováno 5. 11. 2011]
77
[10]
Oesterling JE, Jacobsen SJ, Chute CG, Guess HA, Girman CJ,
[11]
tudent, V., Grepl, M., Král, M., Hartmann, I.,
Serum prostate speci c antigen in a communitybased population of healthy men. Establishment of age speci c reference ranges, JAMA 1993;270(7):8604
Panser LA, Lieber MM,
Má vy²et°ení
PSA stále význam p°i vyhledávání karcinomu prostaty?, Urologie pro praxi 5/2006, 214-218 dostupné na http://www.urologiepropraxi.cz/pdfs/uro/2006/05/04.pdf [citováno 5. 11. 2011]
[12]
Pe²l, M., Záme£ník, L., Soukup, V., Dvo°á£ek, J., Prostatický speci cký antigen a odvozené parametry, Urologie pro praxi 2/2004, 59-63
dostupné na http://www.urologiepropraxi.cz/pdfs/uro/2004/02/05.pdf [citováno 5. 11. 2011] [13]
PSA-velocita a její význam pro v£asnou diagnostiku Ca prostaty, Urologie pro praxi 6/2008, Hrabec, M., tudent, V., Grepl, M., Král, M.,
309-312 dostupné na http://www.urologiepropraxi.cz/pdfs/uro/2008/06/06.pdf [citováno 5. 11. 2011] [14]
Diagnostika karcinomu prostaty-sou£asné moºnosti a limitace transrektální ultrazvukem vedené biopsie prostaty,
ermák,
A.,
Pacík,
D.,
Urologie pro praxi 4/2002, 142-149 dostupné na http://www.urologiepropraxi.cz/pdfs/uro/2002/04/03.pdf [citováno 5. 11. 2011]
Categorical Data Analysis, second edition, John Wiley &
[15]
Agresti, Alan,
[16]
Laura
[17]
Ji°í And¥l,
[18]
Pavla Kunderová,
[19]
Karel Zvára,
Sons, Inc., 2002
pany
S-Plus (and R) Categorical Data Analysis,
Thompson,
A.
Agresti's
Manual
to
Accom-
2005 dostupné http://www.karlin.m.cuni.cz/∼kulich/vyuka/Rdoc/SplusAgresti.pdf [citováno 5. 2.2012]
na
Základy matematické statistiky, 3. vydání, Matfyzpress, vyda-
vatelství matematicko-fyzikální fakulty univerzity Karlovy, 2011
Základy pravd¥podobnosti a matematické statistiky,
1. vydání, Univerzita Palackého v Olomouci, 2004
Regrese, 1. vydání, Matfyzpress, vydavatelství matematickofyzikální fakulty univerzity Karlovy, 2008
78
