Analýza přežití Cumulative Proportion Surviving (Kaplan-Meier) Complete Censored 0,9 1,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 CumulativeProportionSurviving Petr Škarka 348736 Group 0, Group 1, 0 5 10 15 20 25 30 35 40 Time -0,1 0,0 0,1 Analýza přežití (Survival analysis; SA) •Soubor statistických procedur analyzujících chování NV doba přežití •Doba přežití udává čas od doby počátku sledování do události •Jednoduše zkoumáme dobu do nastání nějaké události •Událostí je obecně jakékoliv selhání (smrt, default,…) •Čas měřen v hodinách, dnech, měsících, kvartálech,…•Čas měřen v hodinách, dnech, měsících, kvartálech,… •Aplikace: •Demografie •Odchod od rodičů,recidiva vězně po propuštění, rozchod páru či rozvod,… •(Bio)medicína •Doba do propuknuti nákazy virem HIV u narkomanů, remise u onkolog. pacientů,… •Technika •Porucha přístroje,…•Porucha přístroje,… •Marketing •Míra návratnosti po zavedení výrobku na trh, požadavek na reklamaci po jeho prodeji •Finance •Kdy klient přestane splácet úvěr •Jaká je doba do zaplacení dlužného pojistného •Modelování úmrtnosti pro účely finančního rozhodování •Atd. Proč analýza přežití a ne jiné metody? •Využití neúplné informace, kdy nenastala sledovaná událost – cenzorovaná pozorování •Sledování času do selhání objektu (po uzavření úvěru, podepsání pojistné smlouvy,…) •Doba trvání studie < doba do selhání… neznáme skutečnou dobu přežití •Objekt přestal vykazovat činnost… známe datum poslední kontaktu •Vyloučení objektu z nějakých důvodů nesouvisejících se studií… stále cenná informace•Vyloučení objektu z nějakých důvodů nesouvisejících se studií… stále cenná informace •Více dat pro modelování doby přežití a následné (finanční) rozhodování •Ideální by byl dataset bez cenzorovaných dat → nereálný předpoklad •Studie = pojistná doba, doba trvání úvěru,... •5 sledovaných subjektů: •A – selhání (zde smrt) •B – přežil dobu studie (a zanedlouho•B – přežil dobu studie (a zanedlouho zemřel (cenzorování)) •C – přežil dobu studie (cenzorování) •D – cenzorování – vyňat se studie (např. subjekt přejel vlak nebo přestal reagovat, což nesouvisí s naší studií) •E – cenzorování – nevíme, kdy přesně zemřel, ale víme, že zemřel Cenzorování a jiné důvody pro analýzu přežití •Cenzorování •Zprava – víme pouze to, že skutečná doba do události je větší než pozorovaná doba (konec studie) •Zleva – víme jen, že doba přežití není větší než pozorovaná doba (HIV a narkomani) •Intervalové – událost nastala v intervalu (viz subjekt E)•Intervalové – událost nastala v intervalu (viz subjekt E) •Neinformativní– důvody pro cenzorování nesouvisí s rizikem výskytu události;subjekty cenzorované v čase ti jsou stále v daném časem součástí „přežívajících“ subjektů •Informativní – souvislost s rizikem nastání události; vychyluje odhady estimátorů •„Usekávání“ (truncation) •Další důvody upřednostnění metod SA: •Vysvětlující proměnná je čas •Data přežívání nemívají symetrické rozdělení •Zešikmené rozložení v histogramu – posunuté těžiště •U rizikové funkce (definice později) často proložení přes Weibull, log-normal,… Označení a terminologie doba přežití T Spojitá NV udávající čas, jež uplyne od počátku sledování po selhání, příp. cenzorování. Realizaci T značíme t a nabývá nezáporných hodnot. Rozložení prstí popsáno hustotou pravděpodobností f(t) a distribuční funkcí F(t). indikace selhání δ Alternativní NV, kde 1 = selhání a 0 = cenzorování, tj. přežití doby studie neboindikace selhání δ Alternativní NV, kde 1 = selhání a 0 = cenzorování, tj. přežití doby studie nebo vyloučení ze sledovaného souboru (vlak, ztracený kontakt,…). funkce přežití S(t) S(t) = P(T > t), tedy vyjadřuje pravděpodobnost přežití okamžiku t. Teoretická funkce přežití je spojitá a nerostoucí. Většinou uvažujeme S(0) = 1 a S(t) → 0 pro t → inf. Opačný jev distribuční funkce S(t) = P(T > t) = 1 – P(T ≤ t) = 1 – F(t). odhad Ŝ(t) Data nesbíráme nekonečně dlouho, pak v čase i může platit S(ti) ≥ 0.odhad Ŝ(t) Data nesbíráme nekonečně dlouho, pak v čase i může platit S(ti) ≥ 0. riziková funkce h(t) Platí, že Nevyjadřuje pravděpodobnost selhání v čase t, ale podmíněnou míru rizika selhání v čase t při dožití se času t („tendence selhat po přežití okamžiku t“). Nezáporná (P) a hodnoty závisí na jednotkách času, v nichž měříme. dt tTdttTtP th dt )/( lim)( 0 ≥+<≤ = → •Jak funkce přežití, tak riziková funkce popisují chování NV T. Interpretace se ale liší. S(t) jako pravděpodobnost přežití okamžiku t (P neselhání), h(t) naopak podmíněným selháním. •Jednu lze odvodit z druhé; po vyjádření: a •Znalost průběhů je vhodná k posouzení pravděpodobnosti přežití t (S(t)), popř. podmíněného t tS tS th ∂ ∂ ⋅−= )( )( 1 )(       −= ∫ t duuhtS 0 )(exp)( •Znalost průběhů je vhodná k posouzení pravděpodobnosti přežití t (S(t)), popř. podmíněného rizika selhání v t (h(t)). •Cíle SA: I. Odhad a interpretace funkce přežití a rizikové funkce II. Srovnání funkce přežití a rizikové funkce ve více skupinách III. Modelování vztahu mezi časem přežití a vysvětlujícími proměnnými Ukázka distribučních a rizikových funkcí pro různá rozdělení* Normální rozdělení Exponenciální rozdělení Weibullovo rozdělení *a různé parametry Kaplan-Meier estimátor •S(t) není obvykle známá → odhad z dat •Bez cenzorovaných údajů (klasický přístup): P(T > ti) = (počet objektů s T > ti / počet všech objektů souboru) •Smysl SA je však právě ve využití cenzorovaných dat•Smysl SA je však právě ve využití cenzorovaných dat •Nejčastěji užívaný odhad S(t) je Kaplan-Meierův •Uvažujme okamžiky, kdy došlo k selhání: t0 < t1 < t2 < … < tk. •Lze odvodit: pro j = 1,…,k •Odvození: )/()()( 1 jjjj tTtTPtStS ≥>⋅= − )/()()/()( )()()()( 111 1 jjjjjj jjjjjj tTtTPtStTtTPtTP tTtTPtTtTPtTPtS ≥>⋅=>>⋅>= =>∧>=>∧≥=>= −−− − •Vztah je rekurentní, pak: )/()()/()( 111 jjjjjj tTtTPtStTtTPtTP ≥>⋅=>>⋅>= −−− )/(1)/( )/()/()()( 1 22110 ii j i jj j tTtTPtTtTP tTtTPtTtTPtStS ≥>∏⋅=≥>⋅ ⋅≥>⋅≥>⋅= = K K Kaplan-Meier estimátor •K-M odhaduje podílem , kde: •ni = počet objektů, které jsou v čase ti v riziku (tj. platí pro ně T ≥ ti) •mi = počet objektů, které selhaly přesně v čase ti •Čím větší rozsah výběru n a zároveň čím více selhání, tím přesnější výsledky )/( ii tTtTP ≥> i ii n mn − •Čím větší rozsah výběru n a zároveň čím více selhání, tím přesnější výsledky •U SA je především důležitý počet selhání •Mnoho cenzorovaných dat → moc se toho nedozvíme •Pro j = 1,…,k lze pak odhad pomocí K-M napsat jako •Pro t ∈ [ ti ; ti+1 ) klademe Ŝ(t) = Ŝ(ti) •Cenzorovaná pozorování nezpůsobují skok v odhadu – navyšují rozsah výběru ve smyslu i ii j i j n mn tS − ∏= =1 )( ^ ∈ •Cenzorovaná pozorování nezpůsobují skok v odhadu – navyšují rozsah výběru ve smyslu počtu objektů v riziku (ni) •(při cenzorování zleva trochu komplikovanější) •K-M odhad fce přežití je neparametrický ML odhad Kaplan-Meier estimátor •Jsou-li data necenzorovaná, pak K-M vede ke stejným hodnotám jako klasický přístup •Necenzorovaná data: nj = nj-1 – mj-1 n mn n mn mn mn mn mn n mn n mn tS iijjjjii j i j − = − = − − ⋅⋅ − − ⋅ − = − ∏= −− = 2211 1 )( K ^ •Na obrázku vpravo výpočet Ŝ(t) •Sloupec qj obsahuje cenzorovaná pozorování •Lze si všimnout, že v tj nejsou součástí mj , tzn. stále je využíváme jako součást nj. nnmnmnnn jji i −− −− = 111111 1 •Ŝ(tj) je statistika → var(Ŝ(tj))? •Odhad pomocí Greenwoodovy formule •Pro fixní t (t ∈ [ ti ; ti+1 )) má statistika přibližně N rozložení Ŝ(t) ~ N(S(t), var (S(t)) •Pak lze odvodni 100(1-α)%-ní asymptotický IS pro skutečnou hodnotu S(t) •IS je symetrický kolem hodnoty Ŝ(t) ∑= −⋅ = j i iii i j mnn m tStS 1 ^ 2 ^ )( )())(var( Grafický výstup (STATISTICA) Logrank test •Obecně test pro posouzení rozdílu mezi dvěma skupinami za určitých předpokladů (zde OK) •Postaveno na χ² statistice •Rozdíl mezi pozorovanými (O) a očekávanými (E) četnostmi •Funkční hodnoty funkce přežití mohou být pro různé skupiny rozdílné*•Funkční hodnoty funkce přežití mohou být pro různé skupiny rozdílné* •Je tato rozdílnost statisticky významná? •H0: neexistuje rozdíl mezi objekty ze skupin × H1: doba přežití se mezi skupinami liší •Nezamítnutí H0 •Rozdíl být nemusí, •Ale může – problém křížení funkcí − Cumulative Proportion Surviving (Kaplan-Meier) Complete Censored 0,6 0,7 0,8 0,9 1,0 CumulativeProportionSurviving •Formálně: , kde LR ≈ χ² (1) (více viz stanford.edu/~kcobb/hrp262/lecture2.ppt) *různé odlišnosti mezi skupinami – různá síla logrank testu → vážené logrank testy )var( )( 22 2 22 EO EO LR − − = Group 1, Group 0, 0 5 10 15 20 25 30 35 40 Time -0,1 0,0 0,1 0,2 0,3 0,4 0,5 CumulativeProportionSurviving Coxův* model proporciálního rizika •Regresní model v SA; na závislou proměnnou doba přežití působí několik regresorů •Alternativa k parametrickému regresnímu modelu AFT (accelerated failure time) •Nejčastěji ve tvaru vycházejícím z rizikové funkce h(t), lze jej přepsat i do tvaru s S(t) •Coxův model: pro t ≥ 0 •X = (X1,…, Xp)‘ je vektor prediktorů •β = (β1,…, βp)‘ je vektor neznámých parametrů, které odhadujeme •h0(t) je základní riziko (baseline hazard function). Neznáme ji, plní roli interceptu. Pro (X1,…, Xp) = (0,…,0) je h(t,X) = h0(t). { }β β ´exp)()(),( 00 1 XthethXth p i ii X ⋅= ∑ ⋅= = p • výraz nezávislý na čase! *Sir David Roxbee Cox (1924 – ), britský statistik •Nonparametric estimation from incomplete observations, #2 WoS (#1 K-M) ∑ = p i ii X e 1 β Coxův model •Coxův model je semiparametrický •Základní riziko h0(t) není specifikovaná funkce – o jejím rozložení nic nepředpokládáme •Parametrická regresní složka •Proč je Coxův model populární?•Proč je Coxův model populární? •exp{X‘β} vždy kladné, pak 0 ≤ h(t,X) < inf •Dobře aproximuje parametrické modely •(β1,…, βp) a poměry rizik (viz dále) lze odhadnout i při nespecifikovaném riziku h0(t) !!! •Lze s ním získat odhady zákl. rizika h0(t), rizikové funkce h(t,X) a funkce přežití S(t,X) •Interpretace: •Př.) datový soubor určitého pojistného kmene (selhání je přestání placení pojistného), u jehož•Př.) datový soubor určitého pojistného kmene (selhání je přestání placení pojistného), u jehož objektů sledujeme 5 regresorů (pohlaví, věk, vzdělání, plat, rodinný stav). Zajímá nás, jak se liší riziko selhání v čase t liší u SŠ vzdělaných žen ve věku 50-55 let s platem z intervalu [15k;20k) lišící se v rodinném stavu (0 = svobodná, 1 = vdaná). Poměr rizik lze pak zapsat jako: • říká, kolikrát vzroste riziko selhání u vdané (h1), viz { } { } )( 55443322110 * 55443322110 2 * 1 5 * 55 exp)( exp)( ),( ),( xx e xxxxxth xxxxxth xth xth HR −⋅ = ⋅+⋅+⋅+⋅+⋅⋅ ⋅+⋅+⋅+⋅+⋅⋅ == β βββββ βββββ ),(),( 2 )( 1 5 * 55 xthexth xx ⋅= −⋅β)( 5 * 55 xx e −⋅β Coxův model •Číslo nazýváme relativní riziko vzhledem k prediktoru xi. •Nezmíněné prediktory jsou fixní •Obecně lze podíl rizikových funkcí, hazard ratio (HR), zapsat jako )( * iii xx e −⋅β •Výše uvedený vztah platí pro libovolné rozdíly v hodnotách prediktorů (o počtu 1 až p) •Odhady a získáváme pomocí parciální věrohodnostní funkce •Hustota doby přežití v Coxově modelu není specifikovaná parametrické funkce → úpravy       −== ∑= )(exp ),( ),( 1 * * i p i ii xx xth xth HR β i ^ β i e ^ β •Hustota doby přežití v Coxově modelu není specifikovaná parametrické funkce → úpravy •Intervaly spolehlivosti pro dané odhady statistik pak pomocí Waldovy statistiky •Lze kromě parametrických β1,…, βp odhadnout i neparametrickou složku h0(t)? •Ano, Gehan-Breslow test •Pak lze kreslit odhady funkce přežití pro specificky nastavené hodnoty regresorů •Např. jak vypadá funkce přežití (defaultu) pro objekty „muž – [20;25) let – ZŠ“. Coxův model - předpoklady •Žádné na specifikaci základního rizika h0(t), ALE za to předpoklad proporciálního rizika – relativní riziko HR se vzhledem k jednotlivým regresorům v čase t nemění. Pak pro dva objekty, které se liší hodnotou regresorů platí . ),( ),( * konst xth xth HR == •Na obrázku vpravo porušení proporcionality (překřížení rizik. funkcí) → Coxův model není vhodný •Jak ověříme proporcionalitu? •Graficky – na ose Y logaritmy kumulativních rizikových ),( xth •Graficky – na ose Y logaritmy kumulativních rizikových funkcí (lze získat přes Kaplan-Meiera) pro každý regresor a na ose X čas nebo log času → musí být rovnoběžné •Rezidua – vhodná jsou Schoenfeldova rezidua pro každý z regresorů Aplikace ve financích •Metody credit scoringu = jestli dojde k defaultu × SA = kdy k defaultu dojde •Využítí → kvantifikace profitability klienta •Výstupem je odhad funkce v čase •Navíc na rozdíl od credit scoringu SA zohledňuje cenzorovaná data (v CS zřejmě zmenšení tréningové matice vypuštěním těchto dat, čímž ale přícházíme o data, se kterými se v realitětréningové matice vypuštěním těchto dat, čímž ale přícházíme o data, se kterými se v realitě také setkáváme) Survival analysis in credit scoring (2009) •Dataset k bankovním úvěrům •Skupiny proměnných •Popis dlužníků, typ úvěru, definice defaultu (ihned/90 dní/180 dní po splatnosti)•Popis dlužníků, typ úvěru, definice defaultu (ihned/90 dní/180 dní po splatnosti) •Úpravy dat (outliers, kategorizace proměnných, korelující kovariáty,…) Survival analysis in credit scoring (2009) – average client •Průměrný klient datasetu + (M 58 %, Ž 42 %): •Ukázka odhadů parametru β tří statisticky nejvýznamnějších proměnných: •Funkce přežití s IS (α = 5 %) a riziková funkce průměrného klienta: Survival analysis in credit scoring (2009) •Efekt vzdělání (nalevo) a splátek úvěru z BÚ u úvěrové či externí instituce (napravo) na h(t): •Ploty rizikové funkce dle definic defaultu: Ihned 90 dnů po splatnosti 180 dnů po splatnosti Zdroje •Mgr. Maria Králová, PhD. – materiály k předmětu Statistika 3 (BPM_STA3), ESF MU. •Kristin Sainani, PhD. – materiály k předmětu Longitudinal Data Analysis, dostupné z: •Mgr. Michaela Stehlíková – diplomová práce Analýza přežití a populační dynamika, PřF MU.•Mgr. Michaela Stehlíková – diplomová práce Analýza přežití a populační dynamika, PřF MU. Dostupné z: •Mgr. Martin Hrba – prezentace Odhady aktuárských předpokladů pomocí analýzy přežití, dostupné z: •Pazdera, Rychnovský a Zahradník - projekt Survival analysis in credit scoring, MFF UK. Dostupné z: