Modelování heterogenity ročních příjmů českých domácností Modeling Heterogeinity in the Czech Household Incomes Marie Forbelská Abstract: The distribution of income in most populations is heterogeneous, with several modes and highly skewed, with a long right-hand-side tail and high density at the lower percentiles. Mixture-model-based clustering has become a popular approach of modeling heterogeinity for its statistical properties and the implementation simplicity of the EM algorithm. We focused on the partitions of household incomes into homogeneous subgroups via mixture models. Key words: Household income, Finite mixture models, EM algorithm, Generalized lambda distribution
1. Úvod Analýza rozdělení příjmů domácností je důležitým ukazatelem pro posuzování životní úrovně obyvatelstva. Příspěvek je zaměřen na stochastické modelování heterogenity ročních příjmů domácností pomocí konečných směsí. Podrobně jsou popsány postupy při vytváření stochastických modelů s využitím normálních směsí pro logaritmy ročních příjmů domácností. Pozornost je také věnována jinému typu směsí, a to GLD směsím, kde GLD značí zobecněné lambda rozdělení. Na závěr je zmíněn mnohem obecnější přístup založený na GLM modelech.
2. Modely konečných směsí Roční příjmy domácností mají rozdělení, které je silně zešikmené a je vícemodální, jak je názorně vidět z obrázků 1a, 1b, 1c, kde jsou vykresleny jádrové odhady hustot dané vzorcem 1 n x − xi fˆn ( x) = , ∑ K nh i=1 h
kde x1,...,xn jsou pozorování, K je tzv. jádro, h je vyhlazovací parametr (více lze najít například v monografii Silverman, 1978, [11], popř. v práci Horová, Zelinka, 2000, [4]).
Obrázek 1a: Jádrový odhad hustoty ročních příjmů českých domácností v roce 2005 (zdroj: EU SILC 2005)
Obrázek 1b: Jádrový odhad hustoty ročních příjmů českých domácností v roce 2006 (zdroj: EU SILC 2006)
Obrázek 1c: Jádrový odhad hustoty ročních příjmů českých domácností v roce 2007 (zdroj: EU SILC 2007)
Za těchto okolností se přímo nabízí modelovat rozdělení ročních příjmů domácností pomocí konečných směsí. Jestliže náhodná veličina X má hustotu tvaru f ( x) = π 1 f1 ( x) + + π k f k ( x) , kde
π j > 0,
j = 1,, k ,
π1 + + π k = 1
a f j (x) jsou hustoty, pak řekneme, že hustota náhodné veličiny X je konečnou směsí k hustot. Pravděpodobnosti π j ( j = 1,, k ) se nazývají váhy směsí (proporce směsí nebo apriorní pravděpodobnosti). Více lze najít v monografii McLachlan and Peel (2000, viz [7]). Abychom mohli rozčlenit směs f ( x; Ψ) = π 1 f1 ( x; θ1 ) + + π k f k ( x; θ k ) , kde Ψ = (π 1 ,, π k −1 , θ1 ,, θ k ) jsou neznámé parametry, je třeba nejprve určit typ rozdělení a následně odhadnout vektor neznámých parametrů Ψ . Standardní metodou odhadu neznámých parametrů na základě náhodného výběru x1,...,xn je tzv. EM algoritmus (Dempster et al., 1977, viz [1]), který ve dvou krocích, E a M, nabízí maximálně věrohodné odhady. Klasifikace, tj. určení příslušnosti jedince či objektu charakterizovaného veličinou x k jedné z konečného počtu tříd, se děje pomocí aposteriorních pravděpodobností
ω j ( xi ) =
π j f j( xi ; θ j ) f ( xi ; Ψ)
,
a to tak, že jedinec či objekt je přiřazen do té třídy, pro kterou je aposteriorní pravděpodobnost maximální.
S ohledem na tvar rozdělení ročních příjmů domácností lze postupovat dvojím způsobem a) buď modelovat hustotu příjmů jako směs logaritmicko normálních rozdělení b) nebo hustotu logaritmu příjmů modelovat jako směs normálních rozdělení, neboť náhodná veličina X má logaritmicko normální rozdělení LN ( µ , σ 2 ) s hustotou tvaru (log( x) − µ )2 2 f ( x) = exp− x > 0, µ ∈ R, σ > 0 , 2 2 2σ x 2πσ 1
pokud náhodná veličina Y = log(X ) má normální rozdělení N ( µ , σ 2 ) . Pro modelování ročních příjmů domácností je výhodnější zvolit postup b). V dalším tedy budeme pracovat už jen s logaritmy ročních příjmů domácností. Obecně, jestliže chceme použít konečné směsi hustot, musíme nejprve určit parametr k počtu komponent. Pro tento účel je opět výhodné využít jádrové odhady hustot.
Obrázek 2a: Jádrový odhad hustoty logaritmů ročních příjmů českých domácností v roce 2005 (zdroj: EU SILC 2005)
Obrázek 2b: Jádrový odhad hustoty logaritmů ročních příjmů českých domácností v roce 2006 (zdroj: EU SILC 2006)
Obrázek 2c: Jádrový odhad hustoty logaritmů ročních příjmů českých domácností v roce 2007 (zdroj: EU SILC 2007)
Ze tvarů odhadnutých hustot je zřejmé, že směs bude obsahovat nejméně tři komponenty. Výpočet provedeme pomocí EM algoritmu. Grafické znázornění výsledků je na obrázcích 3a, 3b, 3c a konkrétní maximálně věrohodné odhady vektoru neznámých parametrů
(
ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ2 ˆ2 Ψ MLE = (π 1 ,, π k −1 , θ1 ,, θ k ) = π 1 ,, π k −1 , µ1 ,, µ k , σ 1 ,, σ k
)
jsou pro k=3 uvedeny v tabulce č. 1.
Obrázek 3a: Výsledky EM algoritmu pro 3 - složkovou směs logaritmů ročních příjmů českých domácností v roce 2005(zdroj: EU SILC 2005). Zelená přerušovaná čára značí jádrový odhad hustoty.
Obrázek 3b: Výsledky EM algoritmu pro 3 - složkovou směs logaritmů ročních příjmů českých domácností v roce 2006(zdroj: EU SILC 2006). Zelená přerušovaná čára značí jádrový odhad hustoty.
Obrázek 3c: Výsledky EM algoritmu pro 3 - složkovou směs logaritmů ročních příjmů českých domácností v roce 2007(zdroj: EU SILC 2007). Zelená přerušovaná čára značí jádrový odhad hustoty.
Tabulka 1: Výsledné maximálně věrohodné odhady parametrů ve 3–složkové směsi logaritmů ročních příjmů českých domácností (zdroj: EU SILC 2005, 2006, 2007). Rok
2005
2006
2007
Komponenta
Proporce πj
Střední hodnota μj
Rozptyl σ2j
1
0.281
11.855
0.200
2
0.400
12.581
0.150
3
0.320
12.309
0.542
1
0.094
11.501
0.017
2
0.623
12.468
0.194
3
0.282
12.351
0.566
1
0.094
11.574
0.015
2
0.632
12.560
0.196
3
0.274
12.400
0.588
Není bez zajímavosti, že ve všech třech případech jedna komponenta, jejíž podíl je kolem jedné třetiny, má vysokou variabilitu a v podstatě pokrývá všechny hodnoty logaritmů ročních příjmů domácností. Tuto skupinu domácností lze také interpretovat jako nerozlišitelnou z hlediska výšky ročních příjmů. Podíl takto nerozlišitelných domácností je příliš vysoký. Navíc z obrázků 3a, 3b a 3c je zřejmé, že pouze 3 komponenty k modelování logaritmů ročních příjmů českých domácností v letech 2005, 2006 a 2007 nejsou postačující. Chceme-li najít nějaký optimální počet složek směsi, nabízí se velmi jednoduchý postup. Provedeme rozčlenění směsi postupně pro jednotlivá k a vybereme takové, které maximalizuje logaritmus sdružené hustoty. Ovšem tento postup nevede k cíli, neboť s přidáním každé další složky se hodnota logaritmu sdružené hustoty neustále zvyšuje. Proto se pro výběr optimálního počtu komponent používají různá penalizovaná kritéria, například BIC kritérium (Bayesian Information Criterion) BIC = 2log(maximized likelihood) – m log(n), kde m je počet neznámých parametrů (v našem případě m=3k-1) a n je počet pozorování (více např. Schwarz ,1978, viz [10]), popř. McLachlan and Peel, 2000, viz [7]).
Obrázek 4. Hodnoty BIC kritéria pro logaritmy ročních příjmů v letech 2005, 2006, 2007.
Na obrázku 4 jsou znázorněny hodnoty BIC kritéria pro různé volby parametru k (k=1,...,9). Funkce Mclust z balíku mclust programovacího jazyka R (viz [2] a [9]) navrhuje volbu k=6 pro všechny roky. Grafické výsledky jsou opět prezentovány na obrázcích 5a, 5b, ˆ 5c a maximálně věrohodné odhady Ψ MLE vektoru neznámých parametrů Ψ pro k=6 jsou uvedeny v tabulce č. 2.
Obrázek 5a: Výsledky EM algoritmu pro 6ti - složkovou směs logaritmů ročních příjmů českých domácností v roce 2005 (zdroj: EU SILC 2005). Zelená přerušovaná čára značí jádrový odhad hustoty.
Obrázek 5b: Výsledky EM algoritmu pro 6ti - složkovou směs logaritmů ročních příjmů českých domácností v roce 2006 (zdroj: EU SILC 2006). Zelená přerušovaná čára značí jádrový odhad hustoty.
Obrázek 5c: Výsledky EM algoritmu pro 6ti - složkovou směs logaritmů ročních příjmů českých domácností v roce 2007 (zdroj: EU SILC 2007). Zelená přerušovaná čára značí jádrový odhad hustoty.
Tabulka 2: Maximálně věrohodné odhady parametrů v 6ti–složkové směsi logaritmů ročních příjmů českých domácností (zdroj: EU SILC 2005, 2006, 2007). Rok
2005
2006
2007
Komponenta
Proporce πj
Střední hodnota μj
Rozptyl σ2j
1
0.130
11.828
0.375
2
0.116
11.506
0.020
3
0.191
12.114
0.020
4
0.159
12.524
0.035
5 6
0.232 0.171
12.843 12.400
0.107 0.610
1
0.153
11.542
0.030
2
0.056
12.154
0.007
3
0.164
12.128
0.034
4
0.185
12.596
0.042
5
0.250
12.839
0.127
6
0.192
12.336
0.666
1
0.165
11.627
0.032
2
0.107
12.151
0.014
3
0.118
12.319
0.021
4
0.190
12.694
0.043
5
0.267
12.929
0.144
6
0.154
12.334
0.736
Podíváme-li se podrobněji na obrázky 5a, 5b, 5c a tabulku č. 2, vidíme, že opět existuje složka, jejíž variabilita je výrazně větší. Tentokrát však tyto nerozlišitelné složky tvoří 15ti až 20ti procentní podíl. I když 6ti-složková směs věrně popisuje celkové rozdělení logaritmů příjmů domácnosti, z hlediska interpretace není příliš vhodná. Vraťme se proto znovu k obrázku č. 4 s hodnotami BIC kritéria pro jednotlivá k. Vidíme, že pro roky 2005 a 2006 stačí uvažovat pouze směs se 4 komponentami. Grafické výsledky EM algoritmu pro tyto dva roky jsou ˆ znázorněny na obrázcích 6a a 6b, maximálně věrohodné odhady Ψ MLE vektoru neznámých parametrů Ψ jsou uvedeny v tabulce č. 3. Pro dokreslení jsou uvedeny i výsledky roku 2007 (na obrázku 6c a odhady parametrů v tabulce č. 3), i když z hlediska BIC kritéria se směs se 4 komponentami jeví jako výrazně horší.
Obrázek 6a: Výsledky EM algoritmu pro 4 - složkovou směs logaritmů ročních příjmů českých domácností v roce 2005 (zdroj: EU SILC 2005). Zelená přerušovaná čára značí jádrový odhad hustoty.
Obrázek 6b: Výsledky EM algoritmu pro 4 - složkovou směs logaritmů ročních příjmů českých domácností v roce 2006 (zdroj: EU SILC 2006). Zelená přerušovaná čára značí jádrový odhad hustoty.
Obrázek 6c: Výsledky EM algoritmu pro 4 - složkovou směs logaritmů ročních příjmů českých domácností v roce 2007 (zdroj: EU SILC 2007). Zelená přerušovaná čára značí jádrový odhad hustoty.
Tabulka 3: Maximálně věrohodné odhady parametrů ve 4 - složkové směsi logaritmů ročních příjmů českých domácností (zdroj: EU SILC 2005, 2006, 2007). Rok
2005
2006
2007
Komponenta
Proporce πj
Střední hodnota μj
Rozptyl σ2j
1
0.120
11.495
0.020
2
0.104
12.103
0.013
3
0.505
12.564
0.163
4
0.271
12.204
0.631
1
0.138
11.543
0.029
2
0.112
12.129
0.013
3
0.458
12.618
0.132
4
0.292
12.374
0.567
1
0.310
11.790
0.169
2
0.103
12.242
0.012
3
0.316
12.718
0.077
4
0.271
12.872
0.290
Modelujeme-li logaritmy ročních příjmů pomocí 4 komponent, opět v roce 2005 a 2006 dostaneme složku, která má vysokou variabilitu a jeví se z pohledu logaritmů ročních příjmů jako nerozlišitelná a tvoří opět dosti vysoký podíl (27 a 29 procent). V roce 2007 tuto složku nenalezneme, což je ovšem na úkor modelování nejnižší příjmové skupiny.
Chceme-li modelovat heterogenitu ročních příjmů domácností, nejsme odkázáni pouze na logaritmicko normální či normální rozdělení, ale existuje celá řada dalších systémů, které dokáží velmi efektivně modelovat rozdělení podobného typu. Jako příklad můžeme uvést Pearsonův či Johnsonův systém křivek. Nevýhodou těchto systémů je ovšem nesnadná interpretace výsledných křivek. V posledních desetiletích se začal využívat další systém rozdělení, a to systém založený na tzv. zobecněném lambda rozdělení (GLD rozdělení), který se prosadil mimo jiné také díky snadné interpretaci parametrů a možnosti efektivně generovat pseudonáhodná čísla při Monte Carlo studiích. GLD rozdělení je zobecněním Tukeova symetrického lambda rozdělení definovaného pomocí kvantilové funkce takto u λ − (1 − u λ ) λ Q(u; λ ) = u) log( 1− u
λ≠0 pro
λ =0
0 ≤ u ≤ 1.
Kvantilová funkce je definovaná vztahem
Q(u ) = inf {x ∈ R : F ( x) ≥ u} , kde F (x) je distribuční funkce.
a) rovnoměrné rozdělení
b) přibližně Cauchyovo rozdělení
c) přibližně normální rozdělení
Obrázek 7: Ukázky Tukeova symetrického lambda rozdělení při různé volbě parametru λ .
Zobecněné lambda rozdělení (GLD rozdělení) je pak definováno pomocí čtyř parametrů, přičemž se používá dvojí parametrizace Ramberg and Schmeiser (1974, viz [9]) – RS GLD rozdělení Freimer, Mudholkar, Kollia, Lin (1988, viz [3]) – FMKL GLD rozdělení RS GLD rozdělení je definováno pomocí kvantilové funkce takto Q(u; λ1 , λ2 , λ3 , λ4 ) = λ1 +
u λ3 − (1 − u λ4 )
λ2
pro
0 ≤ u ≤ 1.
Pokud chceme explicitně vyjádřit hustotu RS GLD rozdělení, použijeme vzorec f (Q(u ) ) =
1 . Q' (u )
Pak hustota má pak tvar f RS −GLD ( x) =
λ3u
λ3 −1
λ2 , + λ4 (1 − u ) λ −1 4
Obrázek 8: Ukázky RS GLD rozdělení se stejnými momenty.
kde
λ1 je parametr polohy a platí λ1 ∈ R , λ2 je parametr měřítka, λ2 > 0 , λ3 , λ4 jsou parametry tvaru a jejich parametrický prostor je tvořen 6 oblastmi S1 = {(λ3 , λ4 ) : λ3 ≤ −1, λ4 ≥ 1} S 2 = {(λ3 , λ4 ) : λ3 ≥ 1, λ4 ≤ −1} S3 = {(λ3 , λ4 ) : λ3 ≥ 0, λ4 ≥ 0} S 4 = {(λ3 , λ4 ) : λ3 ≤ 0, λ4 ≤ 0} (1 − λ3 )1−λ λ (λ4 − 1) λ −1 < − 3 S5 = (λ3 , λ4 ) : − 1 < λ3 ≤ 0, λ4 > 1, λ −λ (λ4 − λ3 ) λ4 1−λ (1 − λ4 ) λ (λ3 − 1) λ −1 < − 4 S 6 = (λ3 , λ4 ) : λ3 > 1, − 1 < λ4 ≤ 0, λ −λ (λ3 − λ4 ) λ3 3
4
4
3
4
3
3
4
Aby se zjednodušil parametrický prostor pro λ3 a λ4 byla navržena jiná parametrizace a FMKL GLD rozdělení je pomocí kvantilové funkce definováno takto
Q(u; λ1 , λ2 , λ3 , λ4 ) = λ1 +
1 u λ3 − 1 (1 − u ) λ4 − 1 − λ2 λ3 λ4
pro
0 ≤ u ≤ 1, λ2 > 0.
Pokud λ3 = 0 , pak výraz u λ3 − 1
λ3 je nahrazen výrazem log(u ) , obdobně pokud λ4 = 0 , pak (1 − u ) λ4 − 1
λ4 .
Obrázek 9: Ukázky FMKL GLD rozdělení
je nahrazeno log(1 − u ) .
Na příkladu ročních příjmů domácností důchodců bez ekonomicky aktivních členů budeme (s využitím balíčku GLDEX v programovacích prostředí R, viz [12]) demonstrovat modelování dvousložkové FMKL GLD směsi pomocí EM algoritmu nejprve pro původní data a taká pro logaritmovaná data, viz obrázky 10a, 10b, 10c a 11a, 11b, 11c.
Obrázek 10a: Výsledky EM algoritmu pro 2 - složkovou FMKL GLD směs ročních příjmů domácností důchodců bez ekonomicky aktivních členů v roce 2005 (zdroj: EU SILC 2005)
Obrázek 10b: Výsledky EM algoritmu pro 2 - složkovou FMKL GLD směs ročních příjmů domácností důchodců bez ekonomicky aktivních členů v roce 2006 (zdroj: EU SILC 2006)
Obrázek 10c: Výsledky EM algoritmu pro 2 - složkovou FMKL GLD směs ročních příjmů domácností důchodců bez ekonomicky aktivních členů v roce 2007 (zdroj: EU SILC 2007)
Obrázek 11a: Výsledky EM algoritmu pro 2 - složkovou FMKL GLD směs logaritmů ročních příjmů domácností důchodců bez ekonomicky aktivních členů v roce 2005 (zdroj: EU SILC 2005)
Obrázek 11b: Výsledky EM algoritmu pro 2 - složkovou FMKL GLD směs logaritmů ročních příjmů domácností důchodců bez ekonomicky aktivních členů v roce 2006 (zdroj: EU SILC 2006)
Obrázek 11c: Výsledky EM algoritmu pro 2 - složkovou FMKL GLD směs logaritmů ročních příjmů domácností důchodců bez ekonomicky aktivních členů v roce 2007 (zdroj: EU SILC 2007)
Roční příjmy domácností důchodců bez ekonomicky aktivních členů jsou jasnou směsí dvou skupin, a to skupiny, kde hlavou domácností je žena a skupiny, kde hlavou domácností je muž (rozlišení hodnot je provedeno modrou a zelenou barvou). Z hlediska grafické interpretace výsledků je výhodnější modelovat logaritmy ročních příjmů domácností.
Na závěr se jen velmi krátce zmiňme o možnosti modelovat roční příjmy domácností pomocí konečných směsí regresních GLM modelů. Tento přístup však předpokládá hluboké znalosti vnitřní struktury modelu, neboť nepopisuje rozdělení ročních příjmů osamoceně, ale s ohledem na další, například sociální, demografické či geografické charakteristiky domácností. Směsi GLM modelů předpokládají, že podmíněné hustoty f1 ,, f k jsou exponenciálního typu, tj. lze je napsat ve tvaru f j ( y;θ j ) = exp{a( y )b(θ j ) + c(θ j ) + d ( y )}, kde a (⋅) , b(⋅) , c(⋅) a d (⋅) jsou známé funkce a parametr θ j = θ j ( x1 ,, xm ) závisí na nějakých regresorech x1 , , x m . neznámé pravděpodobnosti π 1 ,, π k opět závisejí na obecně jiných regresorech u1 ,, u p , tj. π j = π j (u1 ,, u p ) . V této GLM směsi se pomocí ryze monotonních linkovacích funkcí g1 a g 2 modelují podmíněné střední hodnoty
µ j ( x1 ,, xm ) = g1−1 (β1 x1 + + β m xm ) a pravděpodobnosti
π j (u1 ,, u p ) = g 2−1 (α1u1 + + α p u p ),
kde linkovací funkcí g 2 může být logit linkovací funkce: g 2 = log 1−ππ , probit linkovací funkce: g 2 = Φ −1 (π ) , kde Φ −1 značí kvantilovou funkci N (0,1) , log-log linkovací funkce: g 2 = − log(− log(−π )) , komplementární log-log linkovací funkce: g 2 = − log(− log(1 − π )) . Připomeňme, že mezi rozdělení exponenciálního typu patří ze spojitých například normální a gama rozdělení, z diskrétních například alternativní, binomické, Poissonovo, negativně binomické rozdělení. Více podrobností o GLM modelech lze najít např. v práci McCullagh a Nelder (1994, viz [6]), a regresní směsi GLM modelů jsou podrobně popsané v monografii McLachlan a Peel (2000, viz [7]).
3. Závěr Shrneme-li předchozí úvahy, je vidět, že pomocí EM algoritmu dokážeme díky vhodně zvolenému modelu konečných směsí provádět bayesovskou klasifikaci, a to na základě aposteriorních pravděpodobností. Následně je pak možné vypracovat podrobnou analýzu struktury jednotlivých komponent směsi, a to jak z hlediska sociálního složení domácností, tak i demografického či geografického. Konečné směsi jsou tedy vhodným stochastickým nástrojem pro provádění tzv. klasifikace bez učitele, někdy se také mluví o stochastické klastrové analýze (stochastic cluster analysis, popř. mixture-model-based clustering) .
4. Literatura [1] DEMPSTER, A. P., LAIRD, N. M. RUBIN, D. B.: Likelihood from Incomplete Data via the EM Algorithm. In Journal of the Royal Statistical Society. Series B (Methodological) 39 (1), pp. 1–38, 1977. [2] FRALEY, C., RAFTERY, A. E.: MCLUST: Normal Mixture Modeling and ModelBased Clustering. R package version 3.0-0; 2006. [3] FREIMER, M., MUDHOLKAR, G.S, KOLLIA, G, LIN, C.T. A study of the generalized Tukey lambda family. In Communications in Statistics – Theory and Methods, 17, pp. 3547-3567, 1988. [4] HOROVÁ, I., ZELINKA, J. Contribution to the bandwidth choice for kernel density estimates. In Computational Statistics, Springer, 22, 1, pp. 31-47, 2007. [5] JOHNSON, N. L., KOTZ, S, BALAKRISHAN, N.: Continuous univariate distributions, Vol. 1, 2nd edition.New York: Wiley & Sons, 1994. [6] MCCULLAGH, P., NELDER, J.A.: Generalized Linear Models. Chapman and Hall, London 1994. [7] MCLACHLAN, G. J. , PEEL, D.: Finite mixture models. New York: Wiley & Sons, 2000. [8] R Development Core Team: R: A language and environment for statistical computing. R. Foundation for Statistical Computing, Vienna, Austria. 2008. URL http://www.Rproject.org [9] RAMBERG, J.S, SCHMEISER, B.W. : An approximate method for generating asymmetric random variables, In Communications of the Associaion for Computing Machinery, 17, pp. 78-82, 1974. [10] SCHWARTZ, G.: Estimating the Dimension of a Model. In The Annals of Statistics, 6 (2), pp. 461-464, 1978. [11] SILVERMAN, B. W.: Density Estimation for Statistics and Data Analysis. Chapman and Hall, New York, 1986. [12] SU, S.: GLDEX: Fitting Single and Mixture of Generalized Lambda Distributions (RS and FMKL) Using Discretized and Maximum Likelihood Methods. R package version 1.0.3.4, 2007.
Kontakt: Marie Forbelská, RNDr., PhD. Ústav matematiky a statistiky, Přírodovědecká fakulta Masarykovy univerzity Kotlářská 2 611 37 Brno Česká republika
[email protected]
