π = P(y bekövetkezik)

Biomatematika (SZIE ÁOTK, 2011. tavasz)

1

A logit modell (=logisztikus regresszió) Ha a függő változó (y ) dichotom (=két lehetséges értéke van, pl. túlélés-halál, siker-kudarc stb.), akkor általában azt feltételezzük, hogy a magyarázó változók az eredmény bekövetkezési valószínűségét befolyásolják, ezért inkább a

π = P(y bekövetkezik) valószínűséget tekintjük függő változónak. Példa: Függ-e a tanulásra fordított időtől annak valószínűsége, hogy sikerül a vizsga statisztikából? Függ-e attól is, hogy vizsga előtt iszunk kávét? Gyengét vagy erőset? Pozitív vagy negatív az összefüggés? Fellép-e interakció a magyarázó változók között? A tanulásra fordított idő folytonos kovariáns, a kávéivás pedig faktor (3 szinttel: nem iszunk, gyengét iszunk, erőset iszunk). 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


2

A kínálkozó legegyszerűbb modell, a lineáris regresszió alkalmazhatatlan, mert a bekövetkezési valószínűség becsült értékei nem mindig fognak 0 és 1 közé esni. A logit modell alapgondolata, hogy a valószínűség helyett egy olyan – a valószínűséggel egyenértékű – mérőszámot használunk, – az ún. logit-et – amelynek értékei nem korlátozódnak a [0, 1] tartományra.

π

→

O = π / (1- π) →

valószínűség oddsz

l = ln( π / (1- π) )

logit

Visszafelé:

l

→

O = exp(l)

→

π = exp(l) / (1+exp(l))

13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


Ismétlés:

0.2

0 (szinte) lehetetlen

0 (szinte) lehetetlen

0.5 azonos eséllyel igen vagy nem

Oddsz

Logit

1 (szinte) biztos

3

1 azonos eséllyel igen vagy nem –1.386

–∞ (szinte) lehetetlen

0.75

Valószínűség

0.25

3

∞ (szinte) biztos

1.098

0 azonos eséllyel igen vagy nem

Házi feladat: ellenőrizzék a számításokat! 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt

∞

(szinte) biztos


4

A logit transzformáció egyértelmű megfeleltetést teremt a [0, 1] és a [–, +] tartományok között.    oddsz 5 logit

logit ( )  ln   1  

valószínűség 0 0.001 0.01 0.1 0.25 0.5 0.75   

logit – –6.907 –4.595 –2.197 –1.099 0 1.099   

4 3 2 1 probab.

0 -1 0 -2 -3 -4 -5


0,5

1


5

A regressziós egyenlet logit(π) = 0 + 1x1 + 2x2 + …, itt π = P(y bekövetkezik). Az egyenletben az x-ek folytonosak vagy dichotomok. Dichotom magyarázó változó esetén az értékeket 0/1-gyel kódoljuk. A 0 lesz a „referencia-csoport”, az 1 a „vizsgált csoport” kódja. Például ha az egyik magyarázó változó „A beteg manduláját korábban eltávolították-e? 0-nem, 1-igen” akkor a mandulaműtéten átesetteket hasonlítjuk a többiekhez, mint referenciacsoporthoz. A kapott regressziós együtthatóból kiolvashatjuk, hogy befolyásolja-e (és hogyan, mennyire) az y bekövetkezésének esélyét a mandula eltávolítása.

 > 0: növeli,

 < 0: csökkenti,

 = 0: nem befolyásolja



6

Ha olyan kategóriás magyarázó változónk van, amelynek k > 2 kategóriája van, akkor abból (k –1) darab x-et készítünk. Például ha a magyarázó változó Milyen gyakran fogyaszt alkoholt?

(A) soha (B) évente 1-2-szer (C) havonta 1-2-szer (D) hetente többször akkor ebből három x-et képezünk:

„dummy változók”

válasz x1

x2

x3

A

0

0

0

B

1

0

0

C

0

1

0

D

0

0

1

és ezzel a (B), (C) és (D) csoportokat hasonlítjuk az (A)-hoz mint referenciacsoporthoz.

A kapott három regressziós együttható az y bekövetkezésének esélyéről informál a (B), (C) és (D) csoportban az (A)-hoz képest! 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


7

Mit olvashatunk ki a regressziós együtthatókból? A regressziós együtthatóból esélyhányadost (=odds ratio, OR ) számolhatunk

OR  e



amelynek jelentése a) ha az x magyarázó változó folytonos (pl. testtömeg): az x egy egységgel való növekedése hányszorosára növeli az y bekövetkezésének oddszát átlagosan! (feltéve,

hogy mindig ugyanannyira) b) ha az x magyarázó változó dichotom a vizsgált csoportban az y bekövetkezésének oddsza hányszorosa a referenciacsoportbelinek c) ha a magyarázó változó több mint 2 kategóriás a szóban forgó csoportban az y bekövetkezésének oddsza hányszorosa a referenciacsoportbelinek 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


8

Példa: Mondjuk, a tanulással töltött időhöz tartozó regressziós együttható becslése és a hozzá tartozó SE, illetve p-érték a következők: b = 0.71, SE = 0.23, p = 0.002. Ekkor az esélyhányados

OR = eb = 2.7180.71 = 2.03, vagyis 1 nappal több tanulás az y bekövetkezésének (a sikeres vizsgának) az oddszát 2.03-szorosára növeli. 95%-os konfidenciaintervallum az oddsz növekedésére: eb  1.96SE = 2. 7180.71  1.960.23 = (1.30, 3.19) A p-érték pedig azt mondja, hogy a sikeres vizsga esélye függ a tanulással töltött időtől (szignifikáns, bizonyító erejű). 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


9

Példa: Mondjuk, a mandulaműtétre, mint magyarázó változóra az alábbi regressziós együtthatót kaptuk:

b = –0.73, SE = 0.31, p = 0.0112. Ennek alapján OR = eb = 2.718–0.73 = 0.48, vagyis azoknál, akik korábban átestek mandulaműtéten, az y bekövetkezésének (pl. torokgyulladásnak) az oddsza kisebb (0.48-szorosa a referenciacsoportbelinek). A p-érték alapján a különbség szignifikáns. Konfidenciaintervallumot is adhatnánk, ugyanúgy, mint az előbb.



Példa: Mondjuk, a vizsga előtti kávéivásra az alábbi együtthatókat kaptuk: bgyenge kávé = 0.31, SE = 0.15, p = 0.039, berős kávé = -0.43, SE = 0.33, p = 0.193,

10

regressziós

Ezek alapján

OR gyenge kávé = 2.7180.31 = 1.36, OR erős kávé = 2.718-0.43 = 0.65, vagyis egy gyenge kávé elfogyasztása a sikeres vizsga oddszát 1.36-szorosára növeli a kávét nem iváshoz képest. Az oddsz szignifikánsan nő, mert p < 0.05. Az erős kávé 0.65-szeresre csökkenti az oddsz-t a kávét nem iváshoz képest. A p-érték 0.193, a csökkentés nem szignifikáns. 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


11

Ha két vizsgált csoport közötti esélyhányados érdekel, akkor azt így kaphatjuk:

ORi  j

ORi ( bi  b j )  e OR j

Például az erős és gyenge kávét ivók csoportja között: 2.718(-0.43–0.31) = 2.718-0.74 = 0.48, vagyis az erős kávét ivók esélye 48%-a a gyenge kávét ivók esélyének.

Ezek fiktív megállapítások, illusztrációk, nem támasztja alá kísérlet! 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


12

Kitekintés: az általánosított lineáris modell Láttuk, hogy a logisztikus regresszió a lineáris regresszióhoz képest abban különbözik, hogy a baloldalon nem maga az y, hanem egy függvénye, logit( P(y bekövetkezik) ) áll. (Nem láttuk, de az is különbség, hogy az  véletlen komponensnek most nem normális eloszlásúnak kell lennie, hanem másmilyennek.) A jobboldal egyébként ugyanaz, az x -ek egy lineáris kifejezése. Az általánosított lineáris modellben megválaszthatjuk, hogy  az y helyett mely függvényét használjuk (pl. logit, probit, log),  az  milyen eloszlású (pl. binomiális, Poisson stb).

R-ben: family

„link függvény” 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


13

A logisztikus regressziós modellezés menete azonos a lineáris regressziós modellezésével. Főbb különbségek:  A modell illesztése nem a legkisebb négyzetek módszerével, hanem a maximum likelihood módszerrel történik. Ennek lényege, hogy olyan modellt illesztünk, melynek valószínűsége – az ún. likelihood – a megfigyelt adatok (a minta) mellett a lehető legnagyobb.  A magyarázó erő számszerűsítésére az R2 helyett más mérőszámokat használnak (az R2 interpretációja problémás).  A modell diagnosztikája bonyolult, rendszerint statisztikussal való konzultációt igényel.

Bánjunk takarékosan a magyarázó változókkal! Ha sok irreleváns változót teszünk a modellbe, nem sok jóra számíthatunk! 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


14

Példa: Egérembriók fejlődése során a blastocysta állapot elérésének valószínűségét vizsgálták in vitro körülmények között. Összesen 94 embrió szerepelt a kísérletben. A magyarázó változók az egycellástól a kétcellás állapotig tartó idő, a kétcellástól a háromcellásig tartó idő, valamint az, hogy az embrió időközben fragmentálódott-e (károsodott-e). Logisztikus regresszióval modellezzük a blastocysta állapot elérésének valószínűségét. Minden embrióra feljegyezzük:  elérte-e a blastocysta állapotot (blastocyst: 1, ha elérte, 0, ha nem – dichotom vagy más néven dummy változó),  a kétcellás állapotig tartó időt (time2C, folytonos),  a kétcellástól a háromcellásig tartó időt (D3C2C, folytonos),  fragmentálódott-e (fragmented: 1, ha igen, 0, ha nem, dichotom). 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


15

A célváltozó a blastocyst, folytonos kovariánsok a time2C és a D3C2C magyarázó változók, fragmented pedig kétszintű faktor (amit azonban rögtön dummy változó formájúra írtunk át). Az adatokat táblázatba rendeztük, melyben minden embriónak egy sor és minden változónak egy oszlop felel meg: embryo id time2C D3C2C fragmented blastocyst 1 1556 1254 0 1 2 1355 1164 0 1 3 1436 1605 0 0 4 1355 1295 0 1 5 1496 1705 0 0 6 1536 1374 1 0 7 1294 1221 1 1 8 1424 1302 0 1 … … … … … 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


16

A regressziós egyenlet logit( P(blastocyst=1) ) = 0 + 1  time2C + 2  D3C2C + 3  fragmented

Valójában két regressziós egyenletről van szó, ugyanúgy, mint az ANCOVA modellnél (vö. 11. előadás): logit( P(blastocyst=1) ) = 0 + 3  0 + 1  time2C + 2  D3C2C

a nem fragmentálódott embriókra, logit( P(blastocyst=1) ) = 0 + 3

+ 1  time2C + 2  D3C2C

a fragmentálódottakra. Az egyenlet skálája logaritmusos, ezért a blastocysta oddsza exp(3)-szeres a fragmentálódott embriókra a nem fragmentálódottakhoz képest. (Azt gondoljuk, hogy az oddsz csökken, ezért 3 –ra negatív értéket várunk.) 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


17

Vizsgáljuk meg a kapcsolatot először grafikusan! Graphs >> XY conditioning plot nem érte el elérte 1100

nem fragmentálódott

D3C2C

2000

1800

1600

1400

1200

1200

1300

1400

1500

1300

1400

1500

fragmentálódott

2200

1100

1200

1600

1600

Minél rövidebb a 2C-ig, illetve a 2C-től a 3C-ig terjedő időszak, annál nagyobb a blastocysta valószínűsége. A fragmentálódás csökkenti a valószínűséget, bár ezt nehéz a grafikonról megítélni, mert kevés a fragmentálódott embrió.

time2C



18

A modell illesztéséhez válasszuk a Statistics >> Fit models >> Generalized linear model funkciót. A párbeszéd-panelt az alábbi módon állítsuk be:



19

Az eredmény: Estimate Std. Error z value Pr(>|z|) (Intercept) 60.752853

16.139927

3.764 0.000167 ***

time2C

-0.022677

0.007294

-3.109 0.001878 **

D3C2C

-0.018798

0.004987

-3.769 0.000164 ***

fragmented

-3.413947

1.281424

-2.664 0.007718 **

Amit kaptunk, formailag pontosan ugyanolyan, mint az ANCOVA vagy a lineáris regresszió eredménytáblája. A különbség az értelmezésben van. Ha a 2C állapot eléréséig tartó idő (time2C) egy perccel több, akkor a blastocysta állapot elérésének oddsza exp(-0.022677)≈0.978szeresére változik. Reálisabb és könnyebben elképzelhető, ha a 10 perccel több idő hatását számítjuk ki: ekkor az oddsz exp(10(-0.022677))≈0.80-szoros. A csökkenés szignifikáns, p=0.001878. 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


20

Ha a 2C-től 3C-ig tartó idő (D3C2C) 10 perccel megnő – miközben a többi magyarázó változó ugyanaz marad –, akkor a blastocysta állapot elérésének oddsza exp(-0.18798)≈0.83-szorosára csökken. Ha az embrió a fejlődése közben fragmentálódik, akkor a blastocysta állapot elérésének oddsza drasztikusan, exp(-3.413947)≈0.03-szeresére csökken. Vegyük azért figyelembe, hogy a standard hiba nagy, a csökkenés 95%-os konfidenciaintervalluma exp(-3.41 ± 1.961.28)≈[0.003, 0.41]. Ezért, bár a csökkenés szignifikáns, elképzelhető, hogy populációs mértéke csak 0.41-szeres. Az általánosított lineáris modellek diagnosztizálása jóval nehezebb, mint a lineáris regresszió – ANOVA – ANCOVA modelleké. Az R lehetőséget ad ugyan a grafikus modellvizsgálatra, de ez inkább csak a kiugró értékek, torzító pontok kiszűrésére elegendő. 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt


21

A menüből válasszuk a Models >> Graphs >> Basic diagnostic plots funkciót: A reziduumok szórása nem kell, hogy homogén legyen, eloszlása sem szükségképpen normális.

15 59

-20 -15 -10

-5

0

5

Normal Q-Q 2

17

0

0

2

17

-2

Residuals

Residuals vs Fitted

-2

Std. deviance resid.

glm(blastocysta ~ time2C + D3C2C + fragmented)

59

10

-2

0

Predicted values

0

1

2

5

10

Residuals vs Leverage

0 2 4

1.0

1759 15

17 93

15

-4

Std. Pearson resid.

Scale-Location

-5

-1

Theoretical Quantiles

0.0

Std. deviance resid.

Predicted values

-20 -15 -10

15

Cook's distance 0.00

0.10

0.20

Leverage

1 0.5

0.5 1

0.30

Néhány kiugró érték előfordul a bal felső reziduumplot-on. A jobb alsó ábráról látszik, hogy ezek nem torzító pontok.



22

Van azonban egy egyszerű módszer, amivel könnyen vizsgálható a logisztikus regresszió illeszkedése. Csoportosítsuk a 2C-ig és a 2C-3C-ig tartó időszakokat mondjuk 300 percenként, adjuk össze a csoportokban megfigyelt blastocystás embriók számát (megfigyelt gyakoriság) és a blastocysta modell szerinti valószínűségeit (várt gyakoriság). Ha e két érték minden csoportban jó egyezést mutat, akkor a modell jól illeszkedik az adatokhoz. Illesztett valószínűségek: > embryo$expected

= fitted(glmMod)

Csoportosítás: > embryo$Bintime2C = with(embryo, round(time2C/300)*300) > embryo$BinD3C2C

= with(embryo, round(D3C2C/300)*300)



23

Összegzés, megfigyelt és várt gyakoriságok: > aggregate(embryo[,c("blastocyst","expected")], by=list(Bintime2C=embryo$Bintime2C, BinD3C2C=embryo$BinD3C2C, fragmented=embryo$fragmented), function(x) round(sum(x),2)) Bintime2C BinD3C2C fragmented blastocyst expected 1 1200 1200 0 16 15.99 2 1500 1200 0 24 24.52 3 1200 1500 0 8 7.61 4 1500 1500 0 14 14.40 5 1200 1800 0 0 0.34 6 1500 1800 0 1 0.15 7 1200 2100 0 0 0.00 8 1200 1200 1 1 1.00 9 1500 1200 1 1 1.01 10 1200 1500 1 4 3.98 11 1500 1500 1 0 0.01 12 1500 2100 1 0 0.00 13. előadás (május 17.) Reiczigel Jenő, Lang Zsolt

Az illeszkedés minden csoportban jó!


24

Érdemes a magyarázó változók szóródási diagramján ábrázolni a blastocysta valószínűségének szintvonalait – ezek egyenesek lesznek, ami kifejezi a modell lineáris jellegét – library(graphics); contour(...):



25

A statisztika megalapozója, tudománnyá szervezője, a kurzus során tanultak nagy részének felfedezője Sir Ronald Aymler Fisher angol statisztikus, matematikus és biológus volt (1890-1962).


π = P(y bekövetkezik)

Recommend Documents