Analýza faktorů ovlivňujících délku doby nezaměstnanosti využitím metod analýzy přežití Jan Popelka Doktorand oboru Statistika Abstrakt: Tento článek věnuje pozornost analýze přežití aplikované na problém nezaměstnanosti. Data získaná z úřadu práce v Příbrami se týkají registrovaných uchazečů o zaměstnání. Jsou nesymetricky rozdělena a cenzorována, což jsou dva z hlavních důvodů, proč byly použity právě postupy analýzy přežití. Součástí analýzy je volba vhodného semiparametrického modelu, odhad jeho parametrů a odpovídajících poměrů intenzit, jejich interpretace a diagnostika odhadnutého modelu i jednotlivých parametrů. Na základě získaných dat je odhadnut průběh základní funkce přežití a z ní jsou pak odvozeny konkrétní funkce přežití pro vybrané skupiny uchazečů o práci. Jejich průběh je graficky znázorněn pro přehlednější srovnání vybraných skupin. Klíčová slova: cenzorovaná data, Coxův proporcionální model, intenzitní funkce, poměr intenzit, věrohodnostní funkce, parciální věrohodnostnífunkce, funkce přežití.
1 Úvod Pojem analýza přežití je využíván k popisu takových dat, která se vztahují k určitému přesně vymezenému období, jehož konec je stanoven konkrétní událostí. Název pochází z oblasti lékařských výzkumů, kde je často sledována právě doba přežití pacientů s určitou diagnózou a událostí ukončující sledování pacienta bývá často jeho smrt. Analýza přežití tak představuje nástroj, který odpovídá na otázku, zda a jakým způsobem závisí doba přežití jednotlivce nebo skupiny jednotlivců stejných vlastností na jednom nebo více sledovaných faktorech. Takovými faktory jsou většinou různé druhy ordinovaných léků, operativní a léčebné postupy, biologické charakteristiky pacientů, jejich zdravotní stav atd. Jedním z cílů analýzy je odhalit, které z možných faktorů připadajících v úvahu skutečně na dobu přežití působí a ovlivňují tak pravděpodobnost, že určitá událost (nejčastěji právě zmiňované úmrtí) nastane v konkrétním čase, za podmínky, že sledovaná událost do této doby nenastala. Jsou dva důvody proč data o přežití není vhodné analyzovat standardními analytickými metodami. Zaprvé jsou tato data ve většině případů rozdělena nesymetricky, převládá kladné zešikmení. Není tedy vhodné využívat analytické nástroje založené na předpokladu normality rozdělení základního souboru. Zadruhé bývají taková data velmi často cenzorovaná. U mnoha sledovaných subjektů nenastane očekávaná událost před koncem experimentu, pacient přežívá a je v lepším případě vyléčen, nebo prostě není možné zjistit, zda a kdy sledovaná událost nastala. To 186
Vědecký seminář doktorandů FIS – březen 2004
proto, že pacient z výzkumu například odejde nebo se odstěhuje a není možné jej již dále sledovat. Kromě využití v medicíně se pro analýzu přežití nabízejí i jiné oblasti aplikace. Je to například analýza faktorů působících na životnost výrobků a nebo, jak se snaží ukázat tento článek, sledování faktorů ovlivňujících dobu nezaměstnanosti. Data použitá v tomto článku byla získána v rámci grantu IGA Vysoké školy ekonomické s názvem "Analýza faktorů ovlivňujících dobu do znovuzaměstnání v ČR".1 Data pocházejí z Úřadu práce v Příbrami. Soubor obsahuje informace o uchazečích o práci, kteří byli na úřadu vedeni v lednu roku 2002. Z celkového počtu 597 uchazečů bylo 422 evidence úřadu během sledovaného období vyřazeno, tzn. nalezli novou práci. 175 pozorování je zprava cenzorováno. Tito uchazeči nedokázali do konce studie v červnu 2003 získat zaměstnání a v evidenci úřadu práce zůstali. Sledovanými faktory, jejichž vliv byl analyzován, jsou věk, vzdělání a pohlaví uchazečů.
2 Semiparametrický regresní model Rozdělení doby přežití může být popsáno dvěma způsoby. Prostřednictvím konkrétní známé hustotní funkce (tzv. parametrický regresní model) nebo pomocí intenzitního poměru (semiparametrický regresní model) v případech, kdy tvar distribuční a hustotní funkce rozdělení doby přežití není znám. Intenzitní poměr se zároveň uplatňuje i ve studiích, kdy je úkolem porovnat šance na přežití mezi vybranými skupinami. Intenzitní funkce vyjadřuje pravděpodobnost, že očekávaná událost (smrt, znovuzaměstnání) nastane v čase t za podmínky, že do tohoto času nenastala. Neboli:
P (t ≤ T < t + δ t T ≥ t ) δt
h(t ) = lim δ t→0
(1)
Regresní model intenzitní funkce s vektorem vysvětlující proměnných x a vektorem neznámých parametrů β má následující tvar: h(t , x, β ) = h0 (t ) r (x, β) ,
(2)
kde funkce h0(t) vyjadřuje změny intenzitní funkce závisející na době přežití. Tato složka je nazývána základní intenzitní funkcí (baseline hazard function). Funkce r ( x, β ) pak zachycuje působení vysvětlujících proměnných. Intenzitní poměr ψ vyjadřuje, kolikrát vyšší je šance na znovuzaměstnání jedince s hodnotami vysvětlujících proměnných definovaných vektorem x1 oproti jedinci s hodnotami vysvětlujících proměnných definovaných vektorem x0. Počítán je následovně: h(t , x1 , β) h0 (t )r ( x1 , β) r ( x1 , β) ψ (t , x1 , x0 ) = = = . (3) h(t , x 0 , β ) h0 (t ) r ( x0 , β ) r ( x 0 , β) Jak plyne z výše uvedené úpravy, je intenzitní poměr závislý pouze na funkci r ( x, β ) a konkrétní tvar základní intenzitní funkce h0(t) nemusí být pro potřeby výpočtu vůbec 1
Grant číslo IG 410043
Vědecký seminář doktorandů FIS – březen 2004
187
znám. Parametry těchto dvou funkcí mohou být tedy odhadovány odděleně. V praxi tak pro odhad intenzitních poměrů postačí pouze hodnoty vysvětlujících proměnných a znalost rozdělení doby přežití tak není nutná. Výhoda tohoto přístupu spočívá ve skutečnosti, že konkrétní rozdělení doby přežití není často vůbec známo a tvar distribuční a hustotní funkce musí být dodatečně zjišťován dalšími analytickými nástroji. Další výhodou je i jednoduchá interpretace odhadnutých parametrů semiparametrického modelu, resp. od nich odvozených poměrů intenzit. Konkrétní tvar funkce r ( x, β ) navrhl Cox jako r ( x, β ) = exp(x β ) . Uvedený model je často nazýván "Coxův proporcionální rizikový model" nebo zjednodušeně "Proporcionální rizikový model". Konkrétní tvar rizikové funkce uvádí vzorec (4). T
h(t , x, β ) = h0 (t ) exp(x β ) A tvar intenzitního poměru vzorec (5) T
(4)
ψ (t , x1 , x0 ) = exp[(x1 − x 0 ) β] (5) K odhadu parametrů semiparametrického regresního modelu je využívána metoda maximalizace věrohodnostní funkce. V případech, kdy je model plně specifikován (je známo rozdělení doby přežití) má věrohodnostní funkce tvar T
n
l(β ) = ∏ i =1
{[ h(t , x , β)]
ci
i
i
× [ S (t i , x i , β ) ]
1− 2 ci
}
(6)
Indikátor cenzorování ci nabývá hodnoty nula pokud je i-tá doba přežití cenzorovaná zprava. V ostatních případech je jeho hodnota jedna. Funkce S(ti ,xi ,β) je funkcí přežití a vyjadřuje pravděpodobnost, že doba přežití i-tého jedince bude stejná nebo delší než doba ti , neboli S(ti,xi ,β) = P(T≥ti ). Z výpočetního hlediska je jednodušší nahradit věrohodnostní funkci jejím logaritmem. Odhad parametrů se pak provede maximalizací funkce: n
{
L ( β ) = ∑ ci ln [ h0 (ti )] + ci xi β + (1 − 2ci )e i =1
T
T
xi β
}
ln [ S0 (t i ) ] .
(7)
Jak již bylo v této kapitole zmíněno, v mnoha případech ovšem není znám konkrétní tvar rozdělení doby přežití. Není tak známa ani funkce přežití S(ti,xi,β). Cox proto navrhl věrohodnostní funkci závislou pouze na vysvětlujících proměnných tzv. parciální věrohodnostní funkci. V případě, kdy se v modelu nenacházejí opakovaná data, je její tvar následující:
xβ l (β) = ∏ e i =1 n
T i
∑ e .2 j ∈R ( t ) ci
T
xjβ
(8)
(i)
Cox předpokládal, že parametry odhadnuté pomocí parciální věrohodnostní funkce budou mít stejné rozdělení jako parametry získané maximalizací plné věrohodnostní funkce. Matematický důkaz tohoto předpokladu poskytli v roce 1993 Andersen, Brogan, Gill a Keiding [1].
2
Součet v čitateli je pro skupinu všech jedinců, kteří v daném čase t(i) práci stále hledali, označeno jako R(t (i)).
188
Vědecký seminář doktorandů FIS – březen 2004
Pokud je tedy odhadován vektor p neznámých parametrů β = ( β1 , β 2 ,..., β p ) , je pak T
řešeno p následujících rovnic, jedna pro každý parametr: ∂L (β )
n
= ∑ ci xik −
∑
x jk exp( x j β ) T
∑
exp(x j β ) . T
(9) ∂β k i =1 j ∈R ( t ) j ∈R ( t ) Pokud se v modelu data opakují, je nutné provést modifikace parciální věrohodnostní funkce. Parametry mohou být odhadnuty prostřednictvím přesného vyjádření navrženého Kalbfleischem a Prenticem [10], nebo pomocí vybrané aproximace (Breslow [2], Efron [4], Cox [3]), které jsou využívány především díky menší výpočetní náročnosti. i
i
3 Diagnostika odhadnutého regresního modelu Pro potřeby testování statistické významnosti odhadnutých parametrů a sestrojení jejich intervalů spolehlivosti je nutné získat odhady jejich směrodatných chyb sˆ( βˆ ) . Veškeré potřebné informace jsou obsaženy v informační matici o rozměrech p x p ∂ L (β ) 2
I (β) = −
∂β
2
.
(10)
ˆ (βˆ ) = I (β ) −1 . Diagonální prvky Kovariační matice je inverzí matice informační Var kovariační matice jsou odhady rozptylů odhadnutých parametrů, jejichž odmocněním se vypočtou hledané odhady směrodatných chyby sˆ( βˆ ) . Za předpokladu, že odhady parametrů jsou normálně rozděleny se střední hodnotou β a směrodatnou odchylkou sˆ( βˆ ) , lze 100(1-α)% interval spolehlivosti vypočítat podle jednoduchého tvaru βˆk ± u1−α / 2 sˆ( βˆk ) . Odhadnuté parametry a vypočtené intervaly spolehlivosti parametrů Coxova modelu jsou ovšem velmi špatně interpretovatelné. Je proto vhodnější dopočíst a interpretovat ) odhadnuté intenzitní poměry ψˆ = exp β , odhady jejich směrodatných chyb sˆ(ψˆ ) = ψˆ sˆ( βˆ ) a následně i jejich intervaly spolehlivosti ψˆ ± u sˆ(ψˆ ) . 1−α / 2
K posouzení statistické významnosti odhadnutých parametrů se nejčastěji používá Waldův test. Waldova statistika má tvar z = βˆ sˆ( βˆ ) . Při platnosti hypotézy o nevýznamnosti parametru má normální rozdělení s jedním stupněm volnosti. Některé statistické programy (např. SAS) počítají Waldovu statistiku v upravením tvaru 2 2 z = [ βˆ sˆ( βˆ )] , který má chí-kvadrát rozdělení s jedním stupněm volnosti. K testování statistické významnosti odhadnutého modelu se nejčastěji využívají test věrohodnostním poměrem, Waldův test a test skórů. Test věrohodnostním poměrem. Tento test je velmi jednoduchý. Porovnává hodnoty logaritmů parciálních věrohodnostních funkcí modelu bez vysvětlujících proměnných s modelem, do kterého Vědecký seminář doktorandů FIS – březen 2004
189
bylo zahrnuto p vysvětlujících proměnných. Za podmínky platnosti nulové hypotézy že veškeré odhadnuté regresní parametry jsou rovny nule, tedy jsou statisticky nevýznamné, má testové kritérium G = 2[ L (βˆ ) − L (0)] chí-kvadrát rozdělení s p stupni volnosti.3 Waldův test Vícerozměrné Waldovo testové kritérium již vyžaduje maticové početní operace. Má T ˆ ˆ . Za podmínky platnosti nulové hypotézy, že veškeré odhadnuté regresní tvar βˆ I(β)β parametry jsou statisticky nevýznamné, má opět chí-kvadrát rozdělení s p stupni volnosti. Test skórů Stejně jako u vícerozměrného Waldova testu je testové kritérium dáno součinem matic u ( 0 ) [I ( 0 )] u ( 0 ) , kde u ( 0 ) = u ( β ) β = 0 je tzv. vektor skórů a I ( 0 ) = I ( β ) β = 0 je informační matice odvozená od nulového vektoru parametrů. Za stejných podmínek jaké platí pro oba výše uvedené testy má chí-kvadrát rozdělení s p stupni volnosti. T
−1
4 Volba vhodného modelu Velmi důležitou součástí analýzy přežití je i volba vhodného modelu. V mnoha případech totiž není riziková funkce ovlivněna jen jedním faktorem.Výběr vhodných vysvětlujících proměnných a vyřazení těch, které jsou v modelu přebytečné je vhodné provádět na základě porovnání alternativních modelů. Jako vhodná statistika pro porovnávání je v literatuře doporučována hodnota věrohodnostní funkce modelu. Právě hodnota věrohodnostní funkce totiž v sobě obsahuje informaci o všech datech obsažených v modelu. Ve statistických packetech je tato statistika často počítána upraveně ve formě −2 log Lˆ . Její hodnota je vždy kladná a obecně platí, čím je nižší, tím je daný model vhodnější. Touto statistikou je však možné pouze porovnávat modely založené na stejných datových souborech. Její hodnota se totiž se změnou rozsahu datového souboru mění. Na podobném principu je založeno Akiakeho informační kritérium AIC = −2 log Lˆ + α q . V tomto tvaru je α předem definovaná konstanta, jejíž hodnota se pohybuje většinou v rozmezí 2 až 6 a q je počet parametrů modelu.
5 Nezaměstnanost na Příbramsku Analýza doby potřebné ke znovuzaměstnání uchazeče o práci vychází ze souboru o 597 nezaměstnaných evidovaných úřadem práce v Příbrami v období mezi lednem 2002 a červnem 2003. Z celkového počtu 597 uchazečů bylo 311 žen a 286 mužů. 422 uchazečů bylo z evidence úřadu během sledovaného období vyřazeno, tzn. nalezli novou práci, byli odvedeni nebo nastoupili náhradní vojenskou službu, odešli na mateřskou dovolenou, do důchodu, přestěhovali se nebo zemřeli. 175 pozorování je zprava cenzorováno. Tito uchazeči nedokázali do konce studie v červnu 2003 získat zaměstnání a v evidenci úřadu práce zůstali i nadále, nebo o nich neexistují žádné další 3
Na matematické detaily testů upozorňuje [1].
190
Vědecký seminář doktorandů FIS – březen 2004
informace. Sledovanými faktory, jejichž vliv na dobu potřebnou ke znovuzaměstnání byl analyzován, jsou věk, vzdělání a pohlaví uchazečů. Věk je jedinou spojitou proměnnou. Vzdělání a pohlaví jsou binární proměnné. Protože u proměnné vzdělání byly rozlišeny čtyři stupně, musely být zavedeny následující umělé proměnné: EDU_2 - nabývá hodnoty 1 pro středoškolské vzdělání bez maturity, jinak hodnoty 0 EDU_3 - nabývá hodnoty 1 pro středoškolské vzdělání s maturitou, jinak hodnoty 0 EDU_4 - nabývá hodnoty 1 pro vysokoškolské vzdělání, jinak hodnoty 0 Pokud ani jedna z výše uvedených proměnných nenabývá hodnoty 1, jedná se pak o uchazeče se základním vzděláním. Proměnná SEX_M nabývá hodnoty 1 pro muže a 0 pro ženu. Pro potřeby zvolení nejvhodnějšího modelu bylo odhadnuto pět alternativních modelů s různými kombinacemi v úvahu připadajících proměnných. Kritériem pro volbu nejvhodnějšího modelu je záporná hodnota dvojnásobku logaritmu věrohodnostní funkce −2 log Lˆ a hodnota Akiakeho testového kritéria tak, jak je nabízí program SAS.
Tabulka 1 Porovnání alternativních modelů Model
Proměnné v modelu
−2 log Lˆ
AIC
1 2 3 4 5 6
žádné AGE AGE+AGE^2 AGE+AGE^2+ SEX_M AGE+AGE^2+EDU_2+EDU_3+EDU_4 AGE+AGE^2+EDU_2+EDU_3+EDU_4+SEX_M
4284.118 4273.778 4258.692 4258.684 4241.168 4241.142
4284.118 4275.778 4262.692 4264.684 4251.168 4253.142
V porovnání s modelem bez vysvětlujících proměnných (v tabulce 1 označen jako „žádné“) je jasné, že alespoň jedna z uvažovaných proměnných má na dobu potřebnou ke znovuzaměstnání vliv. Potvrzují to obě sledované statistiky, které jsou u prvního navrženého modelu nejvyšší.
Tabulka 2 Srovnání alternativních modelů pomocí testu věrohodnostním poměrem Porovnávané modely 2 vs 1 3 vs 2 4 vs 3 5 vs 3 6 vs 5
G 10,34 15,086 0,008 17,524 0,026
Df 1 1 1 3 1
p-value 0,001 0,000 0,929 0,001 0,872
Přidáním umělé proměnné AGE^2, která je odvozena jako druhá mocnina věku (AGE), bylo zohledněno zjištění, že vliv věku není lineární. Důkaz o tomto tvrzení podává Jarošová v pracích [7] a [8]. Tato nově vložená proměnná zohledňuje skutečnost, že osoby velmi mladé nebo naopak v pokročilém věku, mají šance na získání nového zaměstnání nižší, než uchazeči ve věku středním.
Vědecký seminář doktorandů FIS – březen 2004
191
Vložení věku a jeho druhé mocniny je pouze jednou z možností jak do modelu nelinearitu věku zahrnout. Alternativní modely, včetně modelu využívajících splinů, jsou uvedeny v [8]. Model s proměnnou AGE^2 (označen jako model 3) se podle sledovaných kritérií ukazuje jako vhodnější, což potvrzuje i test věrohodnostním poměrem (viz. tabulka 2). Stejně tak se hodnota −2 log Lˆ statisticky významně snižuje se zahrnutím vlivu vzdělání, jak opět dokazuje test věrohodnostním poměrem v tabulce 2. Překvapivě však žádné zlepšení nepřináší zahrnutí proměnné pohlaví. Ani model 4 ani model 6 totiž v porovnání s alternativními modely nepřinášejí statisticky významnou změnu ukazatele −2 log Lˆ . Přesto existuje předpoklad, že ženy mají horší možnosti znovuzaměstnání než muži a proto bude i faktor pohlaví do výsledného modelu zahrnut. To zda je tento předpoklad správný či nikoliv může potvrdit nebo vyvrátit analýza založená na rozsáhlejším souboru z celé České republiky, která bude v rámci uvedeného grantu také zpracována. Jako nejvhodnější byl tedy zvolen model s pořadovým číslem 6. Hodnota obou sledovaných kritérií je v jeho případě ze všech uvažovaných modelů nejnižší. Pomocí programu SAS byly odhadnuty parametry zvoleného modelu. Protože se v souboru vyskytují opakovaná data, bylo nutné provést odhad pomocí modifikované parciální věrohodnostní funkce, tak jak ji navrhli Kalbfleisch a Prentice. Výstup programu SAS je v tabulkách 3 a 4.
Tabulka 3 Odhady parametrů modelu Variable D Parameter F Estimate
AGE AGE^2 SEX_M EDU_2 EDU_3 EDU_4
1 1 1 1 1 1
Standard Error
ChiSquare
0.09556 0.03236 8.7217 -0.00149 0.0004417 11.4382 0.01627 0.09951 0.0267 0.53016 0.15113 12.3053 0.61198 0.16030 14.5752 0.44136 0.32190 1.8800
Pr>ChiS Hazard q Ratio
0.0031 0.0007 0.8701 0.0005 0.0001 0.1703
1.100 0.999 1.016 1.699 1.844 1.555
95% Hazard Ratio Confidence Limits 1.033 1.172 0.998 0.999 0.836 1.235 1.264 2.285 1.347 2.525 0.827 2.922
Tabulka 4 Testování významnosti odhadnutého modelu Testing Global Null Hypothesis: BETA=0 Test Chi-Square DF Pr > ChiSq Likelihood Ratio 42.9759 6 <.0001 Score 38.8386 6 <.0001 Wald 37.8585 6 <.0001 Všechny tři testy uvedené v tabulce 4 potvrzují statistickou významnost navrženého modelu. Jinými slovy je alespoň jeden z odhadnutých parametrů statisticky významný. 192
Vědecký seminář doktorandů FIS – březen 2004
Odhady jednotlivých parametrů pak znázorňuje tabulka 3. Jak již bylo zmíněno kapitole 3 nabízejí statistické packety kromě odhadů parametrů βˆ a jejich směrodatných ) chyb sˆ( βˆ ) i hodnoty intenzitních poměrů vypočtených jako ψˆ = exp β a jejich intervalů spolehlivosti. Důvodem je právě jejich snadná interpretovatelnost. Podle Waldova testu (sloupce Chi-Square a Pr>ChiSq) se potvrzuje to, co již odhalilo srovnávání alternativních modelů. Vliv pohlaví se ukazuje jako statisticky nevýznamný. Jeho p-value je 0,8701. Na základě analyzovaných dat se tedy nepodařilo prokázat, že by existoval statisticky významný rozdíl mezi ženami a muži v šancích na znovuzískání zaměstnání. Stejně tak se jako nevýznamný ukazuje i vliv vysokoškolského vzdělání v porovnání se vzděláním základním. Také tento závěr je nutné ověřit na rozsáhlejším datovém souboru, protože neodpovídá předpokladu, že šance s rostoucím vzděláním rostou. Naopak parametr druhé mocniny věku statisticky významný je a znovu tak potvrzuje vhodnost zahrnutí této dodatečné proměnné. Interpretace odhadnutých intenzitních poměrů je tedy následující. S každým dalším dosaženým rokem věku nezaměstnaného roste příležitost k získání práce 1,1 krát. S pravděpodobností 95% se tento poměr bude pohybovat v rozmezí 1,03 až 1,17. Oproti uchazeči se základním vzděláním má středoškolák bez maturity 1,70 krát a středoškolák s maturitou 1,84 krát větší šanci, že získá nové zaměstnání. Analogicky je možné interpretovat i odhadnuté intervaly spolehlivosti. Jak již bylo zmíněno dříve, obsahuje analyzovaný soubor opakovaná data. K odhadu parametrů byl použit přesný tvar modifikované parciální věrohodnostní funkce. Statistický packet SAS však nabízí i odhad parametrů založený na Breslowově a Efronově aproximaci parciální věrohodnostní funkce. Srovnání odhadů s použitím alternativních metod je uvedeno v tabulce 5.
Tabulka 5 Porovnání poměrů intenzit odhadnutých aproximacemi parciální věrohodnostní funkce Metoda Přesné vyjádření Beslowova aproximace Efronova aproximace
AGE 1.100
AGE^2 0.999
Poměr intenzit SEX_M EDU_2 1.016 1.699
EDU_3 1.844
EDU_4 1.555
1.100 1.100
0.999 0.999
1.016 1.016
1.840 1.844
1.552 1.555
1.695 1.699
Z tabulky vyplývá, že rozdíly mezi odhady získanými exaktním přístupem a Efronovou aproximací jsou minimální. Výpočetně nejjednodušší aproximace Breslowova přináší odhady více odchýlené.
6 Odhad funkce přežití Výše uvedený odhad parametrů proporcionálního modelu vychází z předpokladu, že rozdělení doby přežití není známo. Přesto je na základě získaných dat možné odhadnout jak tvar funkce přežití S(t,xi ,β), tak i intenzitní funkce h(t,xi ,β) i-tého jedince nebo skupiny jedinců stejných vlastností. Pro připomenutí vyjadřuje funkce přežití Vědecký seminář doktorandů FIS – březen 2004
193
pravděpodobnost, že doba nezaměstnanosti i-tého jedince bude stejná nebo delší než doba t, neboli S(t,xi ,β) = P(T≥t). T Odhadnutá intenzitní funkce i-tého jedince hˆ(t , x , βˆ ) = hˆ (t ) exp(x βˆ ) je závislá jednak i
o
i
na odhadnutých parametrech proporcionálního modelu, jednak na odhadu základní intenzitní funkce hˆ (t ) , který založili Kalbfleisch a Prentice [9] na metodě maximální o
věrohodnosti. Odhad základní rizikové funkce v čase t(j) vychází z tvaru hˆ0 (t ( j ) ) = 1 − ξˆ j , kde ξˆ j je řešením následující rovnice:
∑
l ∈D ( t ( j )
T exp( xl βˆ )
ˆ exp( x βˆ ) ) 1− ξ j T l
=
∑
T exp( xl βˆ ) pro j = 1,2, ... ,r
(11)
l ∈R ( t ( j ) )
Výše uvedená rovnice vychází z předpokladu, že doby do znovuzaměstnání byly seřazeny vzestupně, takže t(1)
k
přežití jako Sˆ0 (t ) = ∏ ξˆj a analogicky s odhadem rizikové funkce i samotné funkce j =1
exp( x β ) přežití Sˆ (t , xi , βˆ ) = [ Sˆo (t )] . T l
Odhad základní funkce přežití je součástí nabídky statistického paketu SAS. Její číselný výstup překračuje možnosti tohoto článku. Přehlednější je její grafické zpracování, vyhotovené pomocí programu EXCEL. Obrázek 1 zobrazuje průběh tří odhadnutých funkcí přežití pro muže se základním vzděláním a věkem 24 let (hodnota dolního kvartilu sledovaného souboru), 33 let (hodnota mediánu) a 46 let (hodnota horního kvartilu). Graficky je tak zobrazen předpoklad o nelinearitě věku, který byl do modelu dodatečně zahrnut. Muž ve věku 33 let tak má ve srovnání s ostatními nejlepší vyhlídky na znovuzískání zaměstnání. Pravděpodobnost, že zůstane v evidenci úřadu práce déle než je doba přežití t, je v každém okamžiku sledovaného období nejnižší. Nejhorší vyhlídky má naopak muž 46-letý. Pravděpodobnost, že zůstane bez zaměstnání déle než je doba přežití t, je v každém okamžiku nejvyšší.
194
Vědecký seminář doktorandů FIS – březen 2004
1
24 let 33 let 46 let
Odhad funkce přežití
0,9 0,8 0,7 0,6 0,5 0,4 0,3 0
50
100
150
200
250
300
350
400
450
500
550
Doba přežití
Obr. 1 Odhad funkce přežití pro muže se základním vzděláním ve věku 24, 33 a 46 let. Obdobným způsobem je na obrázku 2 srovnán vliv různých stupňů dokončeného vzdělání. Jak již vyplynulo z odhadů intenzitních poměrů v kapitole 5, jsou vyhlídky na znovuzaměstnání nejhorší u uchazečů se základním vzděláním. Nejlepší pak u uchazečů s maturitou. 1
Středoškolské bez maturity Středoškolské s maturitou Základní
0,9 Odhad funkce přežití
0,8 0,7 0,6
Vysokoškolské
0,5 0,4 0,3 0,2 0,1 0 0
100
200
300
400
500
Doba přežití
Obr. 2 Odhad funkce přežití pro muže ve věku 33 se vzděláním základním, středoškolským bez maturity, středoškolským s maturitou a vysokoškolským. Vědecký seminář doktorandů FIS – březen 2004
195
7 Závěr Uvedený článek přináší další z možných oblastí aplikace nástrojů analýzy přežití. Problematika nezaměstnanosti je v České republice stále více aktuální a sledování faktorů působících na tento jev nabývá na významu. Analýza kvantifikuje vliv jen malého množství vybraných faktorů, které jsou v souvislosti s dobou nezaměstnanosti asi nejvíce zmiňovány. Je možné ovšem vyhodnotit i vliv faktorů méně významných, pokud budou k dispozici vhodná data. Pro další studium je třeba soustředit se na problém nelinearity věku v souvislosti s volbou nejvhodnějšího modelu. Dále pak prozkoumat působení pohlaví a také vysokoškolského vzdělání, jejichž vliv se na základě sledovaného datového souboru ukázal jako statisticky nevýznamný. Mezi další faktory, které by bylo vhodné do analýzy zařadit, určitě patří i vliv regionu, ve kterém se nezaměstnaný o práci uchází, protože rozdíly v tomto směru jsou v České republice dost významné.
Literatura [1] ANDERSEN, P.K., BORGAN, O., GILL, R.D., KEIDING, N.: Statistical Models Based on Counting Processes, Springer–Verlag, N.Y. 1993 [2] BRESLOW, N.: Covariance Analysis of Survival Data under the Proportional Hazards Model, International Statistical Review 1974, č.43 [3] COX, D.R.: Regression Models and Life Tables, Journal of the Royal Statistical Society, Series B 1972, č.34 [4] EFRON, B.: The Efficiency of Cox’s Likelihood Function for Censored Data, Journal of the American Statistical Association 1977, č.72 [5] ESSER, M., POPELKA, J.: Analysis of Factors Influencing Time of Unemployment Using Survival Time Analysis, Zborník 12. medzinárodného seminára Výpočtová štatistika, SŠDS, Bratislava 2003 [6] HOSMER, D.W., LEMESHOW, S.: Applied Survival Analysis, J.Wiley & Sons, N.Y. 1999 [7] JAROŠOVÁ, E.: Analysis of Interval Censored Data, Universita Mateja Bela, Banská Bystrica 2003 [8] JAROŠOVÁ, E.: Exploring the Functional Form of Covariates in Cox Model, Zborník 12. medzinárodného seminára Výpočtová štatistika, SŠDS, Bratislava 2003 [9] KALBFLEISCH, J.D., PRENTICE, R.L.: Marginal Likelihoods Based on Cox’s Regression and Life Table Model, Biometrika 1973, č.60 196
Vědecký seminář doktorandů FIS – březen 2004
[10] KALBFLEISCH, J.D., PRENTICE, R.L.: The Statistical Analysis of Failure Time Data, Wiley, N.Y. 1980 [11] THERENEAU, T.M., GRAUBSH, P.M.: Modeling Survival Data: Extending The Cox Model, Springer–Verlag, N.Y. 2000
Summary ANALYSIS OF FACTORS INFLUENCING TIME OF UNEMPLOYMENT USING SURVIVAL TIME ANALYSIS Survival time analysis approach is used to examine factors influencing the hazard ratios and the length of unemployment. Analysis is based on data aquired from the Labour office in Příbram. Cox proportional model for right censored data is fitted to obtain the hazard ratio estimates. More alternative models are compared to choose the apropriate one. Tests of model and parameters significance are evaluated. Survivorship function is estimated.
Vědecký seminář doktorandů FIS – březen 2004
197