Univerzita Karlova v Praze Matematicko-fyzikální fakulta
BAKALÁŘSKÁ PRÁCE
Jakub Petrásek
Aplikace zobecněného modelu lineární regrese v bankovnictví Katedra pravděpodobnosti a matematické statistiky
Vedoucí bakalářské práce: Mgr. Karel Vaníček Studijní program: Matematika Studijní plán: Ekonometrie
2006
Poděkování Rád bych poděkoval Mgr. Karlu Vaníčkovi za zajímavé téma, za cenné rady a podněty, věcné připomínky a za poskytnutá data.
Prohlášení: Prohlašuji, že jsem svou bakalářskou práci napsal samostatně a výhradně s použitím citovaných pramenů. Souhlasím se zapůjčováním práce a jejím zveřejňováním. V Praze dne 23. května 2006
Jakub Petrásek 2
Obsah 1 Úvod 1.1 Motivační příklad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Zobecněný lineární model 2.1 Základní popis . . . . . . . . . . . 2.2 Příklad . . . . . . . . . . . . . . . . 2.3 Momenty rozdělení exponenciálního 2.4 Odhad parametrů . . . . . . . . . . 2.5 Testování submodelu . . . . . . . .
5 5
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
8 8 9 9 10 11
3 Logistická regrese 3.1 Základní vlastnosti . . . . . . . . . . . . . . . . 3.2 Vztah logistické regrese k normálnímu rozdělení 3.3 Odhady parametrů pro model logistické regrese 3.4 Testování správnosti odhadu parametru β . . . 3.5 Diagnostika . . . . . . . . . . . . . . . . . . . . 3.6 Deviance a její použití . . . . . . . . . . . . . . 3.7 Predikční síla modelu . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
12 12 14 15 16 18 19 21
. . . . . . typu . . . . . .
. . . . .
. . . . .
. . . . .
4 Aplikace modelu logistické regrese 22 4.1 Základní statistické charakteristiky . . . . . . . . . . . . . . . . . . . 23 4.2 Hledání modelu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.3 Diagnostika . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5 Závěr
35
3
Název práce: Aplikace zobecněného modelu lineární regrese v bankovnictví Autor: Jakub Petrásek Katedra (ústav): Katedra pravděpodobnosti a matematické statistiky Vedoucí bakalářské práce: Mgr. Karel Vaníček E-mail vedoucího:
[email protected] Abstrakt: Tato práce se zabývá popisem, diagnostikou a aplikací zobecněných lineárních modelů. Soustředí se na model logistické regrese, který se hojně užívá v praxi. Interpretuje jeho parametry a uvádí vztah modelu logistické regrese k normálnímu rozdělení. Vedle odhadu těchto parametrů a užívaných testových statistik klade důraz na diagnostiku. Následně pojednává o omezeném užití testů míry dobré shody a navrhuje alternativní řešení. Popisuje postupy, které hodnotí predikční sílu. V závěru diskutované poznatky aplikuje na konkrétních datech. Klíčová slova: logistická regrese, linková funkce, deviance, rezidua
Title: Banking Application of Generalized Linear Regression Model Author: Jakub Petrásek Department: Department of Probability and Mathematical Statistics Supervisor: Mgr. Karel Vaníček Supervisor’s e-mail address:
[email protected] Abstract: This thesis deals with description, diagnostics and application of generalized linear models. It focuses on logistic regression model, which is much used in practice. After interpreting its parameters, the thesis shows the relation between the logistic regression model and normal distribution. Besides estimates of these parameters and types of inference, emphasis is laid on diagnostics. Subsequently, the reader is informed about limited use of goodness of fit tests and variant solutions are suggested. Procedures which assess predictive power are described. Finally, discussed principles are applied on actual data. Keywords: logistic regression, link function, deviance, reziduals 4
Kapitola 1 Úvod Zobecněné lineární modely popisují závislost středních hodnot náhodných veličin na veličinách nenáhodných. Termín „zobecněnéÿ značí, že náhodné veličiny se řídí rozdělením z rodiny exponenciálního typu. Do rodiny exponenciálního typu patří například alternativní, Poissonovo a normální rozdělení. Práce se soustředí na studium náhodných veličin s alternativním rozdělením. Model logistické regrese, speciální případ zobecněných lineárních modelů, je stále více využíván zejména v bankovnictví. Velmi názorným příkladem jeho aplikace je hodnocení rizikovosti poskytnutí úvěru. Na základě získaných údajů o klientovi se snažíme předpovědět pravděpodobnost, s jakou klient nesplatí, neboli nastane default. Odhad pravděpodobnosti defaultu je užitečný, neboť o většině bankou přijatého rizika se rozhoduje právě v čase poskytnutí úvěru. Neodhadne-li banka klienta správně, další její aktivity mohou vést jen ke snížení již vzniklé ztráty. Jelikož odhadujeme pravděpodobnost, hledáme funkci, jejíž obor hodnot náleží do intervalu [0, 1] při minimálním omezení na její definiční obor. Tomuto omezení vyhovuje funkce deexp(α+βx) finovaná vztahem 1+exp(α+βx) , kterou je určen logistický regresní model. Logistická regrese má řadu užitečných vlastností. Může obsahovat polynomy a jiné transformace a není omezena ani typem dat, lze ji aplikovat i na kategorické proměnné. Další výhodou je možnost bezprostřední interpretace odhadnutých parametrů.
1.1
Motivační příklad
Následující příklad ilustruje, kdy je vhodnější použít logistickou regresi namísto regrese lineární. Pracujeme s daty uvedenými v tabulce 1.1. Na data aplikujeme dva modely, lineární a logistický. Podívejme se na analýzu rozptylu1 a odhady parametrů.
1
V tabulkách zachováme značení, které uvádí statistický program R. V případě logistického modelu pro nás Deviance bude test poměrem věrohodností a Resid. Dev bude značit devianci
5
Věk Default ≤ 25 10 (25, 30] 33 (30, 35] 33 (35, 40] 34 (40, 45] 24 (45, 50] 19 (50, 55] 15 > 55 11
Spláceno 77 304 495 496 513 442 321 169
Výběrový průměr 0.11494253 0.09792285 0.06250000 0.06415094 0.04469274 0.04121475 0.04464286 0.06111111
Výběrový rozptyl 0.10291366 0.08859686 0.05870493 0.06014909 0.04277495 0.03960200 0.04277719 0.05769708
Tabulka 1.1: Seskupené údaje o četnosti výskytu defaultu v závislosti na věku klienta NULL vek
Df
Deviance
1
13.23
Resid. Df 2995 2994
Resid. Dev 1355.80 1342.58
P(> |Chi|) 0.000276
Tabulka 1.2: Analýza rozptylu pro případ logistické regrese
vek Residuals
Df Sum Sq Mean Sq 1 0.728 0.728 2994 167.577 0.056
F value 13.012
Pr(>F) 0.0003145
Tabulka 1.3: Analýza rozptylu pro případ lineární regrese Dostali jsme tedy dva modely pro odhad pravděpodobnosti defaultu, logistický a lineární, které jsou určeny vztahy (1.1) a (1.2). 1 1 + exp(1.502277 + 0.031524x) p(x) = 0.129823 − 0.001713x. p(x) =
(1.1) (1.2)
Při práci s modelem logistické regrese je třeba rozlišovat mezi pojmy šance a pravděp(x) . Z modelu (1.1) lze vyčíst, že s každým podobnost. Šanci definujeme jako poměr 1−p(x) −0.031524 rokem se šance defaultu snižují e =0.97-krát, ˙ tedy asi o 3 %. Model lineární regrese nám okamžitě odpoví, jak se mění pravděpodobnost defaultu s rostoucím věkem. S každým rokem se pravděpodobnost defaultu podle modelu (1.2) snižuje o 0.001713. Pokud srovnáme čtyřicetiletého klienta s klientem třicetiletým, vidíme z modelu (1.1), že šance defaultu pro staršího z nich klesla e10·(−0.031524) =0.73-krát, ˙ to znamená o 27 %. Na druhé straně pro třicetiletého klienta vychází odhad pravděpodobnosti defaultu 0.082 a 0.061 pro klienta čtyřicetiletého. Tedy pravděpodobnost defaultu klesla o necelých 16 %. Z grafu 1.1 je patrné, že křivka modelu logistické regrese je blíže skutečným hodnotám než křivka regrese lineární. Použití logistické regrese je vhodnější také vzhledem k faktu, že pokud bychom z dat chtěli učinit závěr pro širší skupinu obyvatel, křivka lineární regrese by se mohla dostat do záporných hodnot. Dosazením 6
Estimate Std. Error z value (Intercept) -1.502277 0.349166 -4.302 vek -0.031524 0.008805 -3.580
Pr(> |z|)) 1.69e-05 0.000343
Tabulka 1.4: Odhady parametrů pro případ logistické regrese Estimate Std. Error t value (Intercept) 0.129823 0.019902 6.523 vek -0.001713 0.000475 -3.607
Pr(> |t|) 8.05e-11 0.000315
Tabulka 1.5: Odhady parametrů pro případ lineární regrese
0.08
Logisticka regrese
0.06
Linearni regrese
0.04
Sance defaultu
0.10
0.12
do vztahu (1.2) zjístíme, že již pro klienta 76 let starého nastává default se zápornou pravděpodobností.
20
30
40
50
60
Vek
Obrázek 1.1: Grafické srovnání logistické a lineární regrese
7
Kapitola 2 Zobecněný lineární model 2.1
Základní popis
Stejně jako v lineárních modelech, také i v zobecněných lineárních modelech vystupují nezávislé náhodné veličiny Y1 , ..., Yn a daná matice X = (xij ) typu n × k. Rozdělení těchto náhodných veličin závisí na daných xij , které budeme nazývat vysvětlující veličiny (prediktory). Náhodné veličiny Yi budeme nazývat vysvětlované. První komponentu modelu tvoří nezávislé náhodné veličiny Y1 , ..., Yn , o kterých budeme předpokládat, že patří do rodiny exponenciálního typu s disperzním parametrem φ. To znamená, že vyhovují hustotě f (yi ; θi , φ) = exp {[yi θi − b(θi )] /a(φ) + c(yi , φ)} ,
(2.1)
kde θi se nazývá přirozený parametr. Po funkci b(θi ) chceme, aby byla ryze monotónní, dvakrát spojitě diferencovatelná a druhá derivace kladná. Druhou složku tvoří vektor {η1 , ..., ηn }, který je lineární funkcí vysvětlujících proměnných {x1 , ..., xn }. Závislost je tvaru ηi =
k X
i = 1, ..., n.
βj xij ,
j=1
Tyto dvě komponenty jsou propojeny takzvanou linkovou (spojovací) funkcí, která je diferencovatelná a ryze monotónní a platí pro ni vztah g(EYi ) = ηi . Zde je vhodné poznamenat, že v zobecněných lineárních modelech se nevyskytuje chybový člen jako u lineární regrese. Důvodem je, že levá strana předchozí rovnice je funkcí střední hodnoty veličiny Yi , nikoliv samotné veličiny Yi . Pokud platí, že g(EYi ) = θi , hovoříme o kanonické P linkové funkci. Pro přirozený parametr pak dostáváme přímou závislost θi = kj=1 βj xij . Z výše uvedeného vyplývá, že nalezení vhodné spojovací funkce je při aplikaci modelu podstatné. Nalezením její kanonické verze se navíc mnoho výpočtů výrazně zjednoduší, jak bude ukázáno později.
8
2.2
Příklad
Definované pojmy aplikujme na následujícím jednoduchém příkladu. Pro alternativní rozdělení označme p = P (Y = 1), pak P (Y = 0) = 1 − p, kde p ∈ (0, 1). Hustota je dána vztahem · ¸ p y 1−y f (y; p) = p (1 − p) = exp y log + log(1 − p) , y ∈ {0, 1} . 1−p
Alternativní rozdělení je definováno jen jedním parametrem p. Dále ze zápisu vyplývá, že toto rozdělení je exponenciálního typu a můžeme dle hustoty (2.1) označit p θ = log 1−p a b(θ) = − log(1 − p). Ukažme si další vztahy mezi jednotlivými parametry, později se nám to bude eθ hodit. Předně, úpravou okamžitě dostáváme, že p = 1+e θ . Dosazením získaného vyjádření parametru p do funkce b(θ) máme b(θ) = log(1 + eθ ). Podívejme se ještě na derivace funkce b(θ). b′ (θ) =
eθ = p, 1 + eθ
b′′ (θ) =
eθ = p(1 − p). (1 + eθ )2
Víme, že střední hodnota alternativního rozdělení je p a tedy příslušná kanonická p linková funkce je tvaru log 1−p . Dále varY = p(1 − p), tudíž rozptyl je funkcí střední hodnoty. Z prvních dvou derivací funkce b(θ) vidíme vztah mezi touto funkcí a prvními dvěma momenty 0-1 rozdělení. Zobecněné lineární modely se neomezují jen na binární vysvětlované proměnné. Dalšími příklady rozdělení spadajících do rodiny exponenciálního typu jsou Rozdělení vel. Y Alternativní Alternativní Multinomické Normální Poissonovo
2.3
Typ dat v matici X Spojitá i kategorická Spojitá i kategorická Spojitá i kategorická Spojitá Spojitá i kategorická
Linková funkce Probit Loglog Zobecněný logit Identická Log
Model Probitový model Log-log model Regrese Regrese Loglineární model
Momenty rozdělení exponenciálního typu
V této kapitole vycházíme z obecného popisu zobecněných lineárních modelů daných hustotou (2.1). Výpočtem prvních dvou momentů ukážeme některé vztahy mezi parametry rozdělení exponenciálního typu. K řešení našeho problému poslouží momentová vytvořující funkce. Označme K = {y ∈ R; f (y; θ, φ) > 0}. ½ ¾ Z yθ − b(θ) tY M (t) = Ee = exp {ty} exp + c(y, φ) dy a(φ) K ½ ¾ Z y [a(φ)t + θ] − b(θ) = exp + c(y, φ) dy a(φ) K 9
= exp
½
b [a(φ)t + θ] − b(θ) a(φ)
¾Z
exp
K
½
¾ y [a(φ)t + θ] − b [a(φ)t + θ] + c(y, φ) dy. a(φ)
Výraz je upraven na tvar, kde v integrálu vystupuje hustota exponenciálního typu, a proto je integrál roven 1. Teď již stačí jen zderivovat momentovou vytvořující funkci ¾ ½ b [a(φ)t + θ] − b(θ) ′ ′ b (a(φ)t + θ), M (t) = exp a(φ) ½ ¾ o b [a(φ)t + θ] − b(θ) n ′ 2 ′′ ′′ M (t) = exp {b [a(φ)t + θ]} + b [a(φ)t + θ] a(φ) . a(φ)
Dostáváme
EY varY
= M ′ (0) = b′ (θ), 2 = M ′′ (0) − [M ′ (0)] = b′′ (θ)a(φ).
(2.2) (2.3)
Vztah (2.2) a předpoklady kladené na funkci b(θ) umožňují, že kanonickou linkovou funkci snadno nalezneme a že je tvaru θ = (b′ )−1 (µ). Dále ze vztahu (2.3) dostáváme, že rozptyl je funkcí střední hodnoty. Dokonce pouhá znalost vztahu mezi střední hodnotou a rozptylem již jednoznačně určuje, o jaké rozdělení z rodiny exponenciálního typu se jedná ([1], str. 136). Obecně platí, že ani znalost všech momentů neurčuje jednoznačně rozdělení.
2.4
Odhad parametrů
Věrohodnostní rovnice Nechť Y1 , ..., YN je náhodný výběr z rozdělení s hustotou (2.1). Nazvěme ℓ(θi , φ) = Q N i=1 f (Yi ; θi , φ) věrohodnostní funkce. Maximálně věrohodným odhadem parametru β je hodnota, při které je maximalizována věrohodností funkce a tedy i její logaritmus L(β) = log(ℓ(β)) =
N X i=1
{[Yi θi − b(θi )] /a(φ) + c(Yi , φ)} .
(2.4)
Funkce je maximalizována pro hodnotu, kdy ∂L(θi , φ)/∂βj = 0 pro j = 0, ..., k. Použitím řetízkového pravidla dostaneme ∂L ∂θi ∂µi ∂ηi ∂L = ∂βj ∂θi ∂µi ∂ηi ∂βj
j = 1, ..., k.
(2.5)
Užitím vztahů (2.2) a (2.3) máme ∂µi /∂θi = ∂b′ (θi )/∂θi = b′′ (θi ) = var(Yi )/a(φ). Podle věty o derivaci inverzní funkce víme, že ∂µi /∂ηi = 1/g ′ (µi ). Dosazením do výrazu (2.5) získáváme soustavu rovnic ¸ N · ∂L(β) X Yi − b′ (θi ) 1 = x =0 ′ (µ ) ij ∂βj varY g i i i=1 10
j = 1, ..., k.
(2.6)
Výsledná soustava rovnic (2.6) se dá interpretovat také tak, že odhad parametru β metodou maximální věrohodnosti závisí na středních hodnotách a rozptylech náhodných veličin Y1 , ..., YN , nikoliv na jejich rozdělení. A navíc z předchozí kapitoly víme, že v případě rozdělení daného hustotou (2.1) je rozptyl funkcí střední hodnoty. Za použití kanonické linkové funkce se soustava (2.6) zjednoduší na tvar N X Yi − b′ (θi ) i=1
a(φ)
xij = 0
j = 1, ..., p.
Pokud je a(φ) konstantní, pak hledaný odhad vyhovuje maticovému zápisu X ′Y = X ′µ, ˆ kde X je matice N × p s hodnotami xij a µ ˆ značí odhad střední hodnoty. Tedy pro zobecněný lineární model s kanonickou linkovou funkcí vyhovují odhady soustavě normálních rovnic, neboli maximalizací věrohodnostní funkce docházíme ke stejným výsledkům, k jakým bychom se dostali aplikací metody nejmenších čtverců.
2.5
Testování submodelu
Pro ověření, zda je model vhodně použit, se nejčastěji užívá test poměrem věrohodností. Jeho princip si vysvětleme s použitím deviance. K tomuto účelu se zavádí nasycený (saturovaný) model. Takovýto model obsahuje odlišný parametr pro každé pozorování. Sám o sobě není použitelný, avšak může sloužit jako dobrá výchozí pozice pro testování správnosti ostatních modelů. Deviance Uvažujme saturovaný model, tedy model, který obsahuje odlišný parametr pro každé pozorování. Označme θ˜i odhad parametru θi pro saturovaný model. Z definice nasyceného modelu dostáváme, že µ ˜i = Yi pro každé i = 1, ..., N . θˆi a µ ˆi jsou příslušné odhady pro testovaný model. Dále LS značí maximální hodnotu logaritmické věrohodnostní funkce pro saturovaný model a L1 pro testovaný. Dle vztahu (2.4) a možným zjednodušením a(φ) = φ/ωi dostaneme škálovanou devianci D(Y, µ)/φ ˆ = −2(L1 − LS ) = 2
N X i=1
h i ωi Yi (θ˜i − θˆi ) − b(θ˜i ) + b(θˆi ) /φ.
Nazvěme dále D(Y, µ) ˆ deviance. Z definice je zřejmé, že čím jednodušší model testujeme, tím větší má devianci. Deviance určuje přiléhavost testovaného modelu na naše data. Aplikaci deviance lze zobecnit na srovnávání dvou nesaturovaných do sebe vnořených modelů, čímž získáme test poměrem věrohodností −2(L1 − L2 ) = −2(L1 − LS ) − [−2(L2 − LS )] = D(Y ; µ ˆ 1 ) − D(Y ; µ ˆ 2 ). Tato statistika má za platnosti testovaného podmodelu asymptoticky rozdělení χ2n−k , kde n − k je rozdíl nezávislých parametrů v testovaných modelech. To znamená, že test poměrem pravděpodobností pro dva do sebe vnořené modely je jednoduše rozdílem jejich škálovaných deviancí. 11
Kapitola 3 Logistická regrese Logistická regrese je speciálním případem zobecněných lineárních modelů a představuje užitečný model při vyhodnocování závislosti binárních veličin na spojitých a kategorických prediktorech. Proč nepoužít jednodušší lineární regresi? Jedním z důvodů je to, že v praxi se mnohem častěji vyskytuje právě nelineární vztah mezi podmíněnou pravděpodobností náhodné vysvětlované veličiny P (Y = 1|x) = p(x) a nenáhodnými prediktory. Změna x má menší vliv na funkci p(x), která je blízko 0 nebo 1. Tato důležitá vlastnost se projevuje typickým tvarem křivky do S, viz obr. 3.1.
3.1
Základní vlastnosti
Mějme náhodnou veličinu Y nabývající hodnot 0 a 1. Jelikož E [Y |x] = 1 · P(Y = 1|x) + 0 · P(Y = 0|x) = P(Y = 1|x), sledujeme podmíněnou pravděpobnost P (Y = 1|x) = p(x). V logistické regresi platí vztah exp(α + βx) p(x) = . (3.1) 1 + exp(α + βx) Ze vztahu (3.1) vyplývá, že u logistické regrese se oproti obyčejné lineární regresi nesetkáme s nutným omezením vysvětlujících proměnných, neboť logistická funkce zobrazuje celé R na interval (0, 1). Zkusme najít linkovou funkci. Z příkladu 2.2 již víme, že p(x) logit[p(x)] = log = α + βx. 1 − p(x)
Takto je definována pro logistickou regresi nejdůležitějsí funkce tzv. LOGIT, neboli logistická funkce. Vidíme, že jde dokonce o kanonickou spojovací funkci. Z omezení p(x) ∈ (0, 1) dostáváme, že logistická funkce je dobře definována. Poslední rovnice se také dá interpretovat jako základní předpoklad modelu logistické regrese, tj. že logaritmus šance je lineární funkcí vysvětlujících proměnných.
12
Alternativní linkové funkce pro binární vysvětlovanou proměnnou Pro přehlednost se omezme na model logistické regrese s jednou vysvětlující proměnnou. Při používání logistické regrese transformujeme data funkcí p(x) =
exp(α + βx) . 1 + exp(α + βx)
Naproti tomu probitový model využívá distribuční funkci normálního rozdělení, platí tedy p(x) = Φ(α + βx). Log-log model předpokládá vztah p(x) = exp(− exp(α + βx)).
NORMALNI ROZDELENI
LOG−LOG FUNKCE
0.8 0.6 0.2 0.0
0.0
0.2
0.2
0.4
0.4
0.4
0.6
0.6
0.8
0.8
1.0
LOGISTICKA FUNKCE
−3
−2
−1
0
1
2
3
−3
−2
−1
0
1
2
3
−3
−2
−1
0
1
2
3
Obrázek 3.1: Srovnání dvou linkových funkcí podobných logistické Interpretace parametru β 1. Analytický pohled: Kladné resp. záporné znaménko nám říká, zda je pravděpodobnost p(x) rostoucí resp. klesající. Sklon p(x) je dán velikostí |β|. Pokud zderivujeme p(x) dle x β {exp(α + βx) [1 + exp(α + βx)] − [exp(α + βx)]2 } ∂p (x) = ∂x [exp(α + βx)]2 = βp(x)[1 − p(x)],
snadno nahlédneme symetrii, křivka p(x) se blíží k jedné pod stejným úhlem jako k nule a dosahuje maximálního sklonu pro p(x) = 1/2. Pro tuto hodnotu je argument v exponenciále roven nule, které se nabývá pro x = −α/β. 2. Statistický pohled: Hodnota eβ vyjadřuje poměr šancí v modelu pro prediktory x + 1 a x. Platí p(x+1) eα+β(x+1) 1−p(x+1) = = eβ . p(x) eα+βx 1−p(x)
13
3.2
Vztah logistické regrese k normálnímu rozdělení
Předpokládejme, že náhodná veličina X se řídí normálním rozdělením N(µj , σ 2 ) za podmínky Y = j, j ∈ {0, 1}. Zajímá nás podmíněná pravděpodobnost P (Y = 1|X = x). Použijme Bayesovu větu P (X = x|Y = 1)P (Y = 1) = P (X = x|Y = 1)P (Y = 1) + P (X = x|Y = 0)P (Y = 0) i h −(x−µ1 )2 √1 exp P (Y = 1) 2σ 2 σ 2π i h i o n h = −(x−µ1 )2 −(x−µ0 )2 √1 P (Y = 0) + exp P (Y = 1) exp 2σ 2 2σ 2 σ 2π
P (Y = 1|X = x) =
Zlomek upravíme a rozšíříme výrazem exp dostáváme P (Y = 1|X = x) =
h
1 −(x−µ0 )2 2σ 2
exp
h
i
P (Y = 0)
(µ1 −µ0 ) x σ2
1 + exp
h
+
(µ1 −µ0 ) x σ2
µ20 −µ21 2σ 2
+
,
=1) + log PP (Y (Y =0)
µ20 −µ21 2σ 2
i
=1) + log PP (Y (Y =0)
i.
0) Výsledkem je logistická funkce s parametrem β = (µ1σ−µ (srovnejme s rovnicí 2 (3.1)). Z příkladu je patrné, že model logistické regrese odpovídá právě té transformaci, kterou hledáme při studiu dat s podmíněným normálním rozdělením. V případě, že by náhodná veličina X v závislosti na Y měnila také svůj rozptyl, pak stejným výpočtem získáváme, že tato podmíněná pravděpodobnost také vyhovuje modelu logistické regrese, avšak s kvadratickým členem x2 . Platí h ³ ´ ³ ´ i µ1 µ0 1 1 1 2 exp 2 σ2 − σ2 x + σ2 − σ2 x + K 1 0 h ³0 ´ ³1 ´ i, P (Y = 1|X = x) = 1 + exp 12 σ12 − σ12 x2 + σµ12 − σµ02 x + K 0 1 1 0 ³ 2 ´ 2 µ µ (Y =1) kde K = 2σ02 − 2σ12 + log σσ10 + log PP (Y . =0) 0
1
Podle ([3], str. 501) platí, že pokud by náhodná veličina X v závislosti na Y měla gamma rozdělení Ga(aj , pj ), pak je vhodné zavést kromě lineárního členu ještě člen log(x). Zkusme toto tvrzení dokázat. Stejnými úpravami jako v předchozím dostáváme p
P (Y = 1|X = x) = h
a11 −a1 x p1 −1 e x P (Y Γ(p1 ) p a0 0 Γ(p0 )
e−a0 x xp0 −1 P (Y = 0) + a
p1
p a1 1 Γ(p1 )
= 1)
e−a1 x xp1 −1 P (Y = 1) i p a0 0 P (Y =1) − log Γ(p0 ) + log P (Y =0)
=
exp (a0 − a1 ) x + pp01 log(x) + log Γ(p1 1 ) h i. = p p a 1 a 0 =1) 1 + exp (a0 − a1 ) x + pp10 log(x) + log Γ(p1 1 ) − log Γ(p0 0 ) + log PP (Y (Y =0) 14
3.3
Odhady parametrů pro model logistické regrese
V sekci 2.4 jsme určili odhad neznámých parametrů pro zobecněný lineární model. Vzhledem k jeho důležitosti si odhad odvoďme konkrétně pro model logistické regrese. Věrohodnostní rovnice Náhodné veličiny Yi , i = 1, ..., n, které vystupují v logistické regresi jako vysvětlované proměnné, mají alternativní rozdělení a předpokládáme u nich nezávislost. Pro zdůraznění závislosti střední hodnoty na vysvětlujících proměnných budeme místo pi psát p(xi ). To znamená, že platí P (Yi = 1) = p(xi ), P (Yi = 0) = 1 − p(xi ). Nechť yi značí počet úspěchů v bodě xi . Pak náhodná veličina odpovídající pozorování yi má binomické rozdělení, jelikož součet náhodných veličin s alternativním rozdělením je roven náhodné veličině s rozdělením binomickým. Označme ji Zi , Zi ∼ Bi(ni , ni p(xi )), kde ni je počet pozorování v bodě xi a náhodná veličina Zi =
ki X
Yij ,
j=1
odpovídá pozorování v bodě xi pro i = 1, ..., N . Celkem máme N nezávislých náhodných veličin s binomickým rozdělením. Platí ℓ(β) =
N µ ¶ Y ni i=1
yi
(p(xi ))yi (1 − p(xi ))ni −yi .
(3.2)
Opět pracujme s logaritmem funkce ℓ(β). Označme b hodnotu, při které je maximalizována L(β) = log ℓ(β). Pro větší přehlednost označme absolutní člen α jako β0 a přidejme do matice vysvětlujících proměnných umělý sloupec jedniček. Před dalším odvozováním je vhodné si celý výraz upravit a použít vztahů p
X p(xi ) log βk xik , = 1 − p(xi ) k=0 1 − p(xi ) =
1 Pp . 1 + exp( k=0 βk xik )
Užitím vztahů (3.3) a (3.4) se logaritmus výrazu (3.2) zjednoduší na tvar ( µ ¶ ¸) · p N X X 1 ni Pp . βk xik + ni log + yi log L(β) = 1 + exp( y i k=0 βk xik ) i=1 k=0 15
(3.3) (3.4)
Nutná podmínka pro hledání extrémů funkce L(β) dává rovnici P ¸ N · exp( pk=0 βk xik ) ∂L(β) X Pp 0= = , yi xij − ni xij ∂βj 1 + exp( k=0 βk xik ) i=1
j = 0, ..., p.
Dostáváme nelineární rovnici pro odhad parametru β (srovnejme se vztahem (2.4)) N X i=1
[yi xij − ni xij pˆ(xi )] = 0,
j = 0, ..., p.
(3.5)
Výsledkem je, že maximálně věrohodný odhad b splňuje nelineární rovnici X ′y = X ′
3.4
exp(b′ x) . 1 + exp(b′ x)
Testování správnosti odhadu parametru β
Metodou maximální věrohodnosti jsme získali b odhad parametru β. K rozhodnutí o správnosti našeho odhadu nám poslouží vhodný test. K jeho sestavení a zejména k intervalovým odhadům potřebujeme druhou derivaci logaritmické věrohodnostní funkce. Platí # " P N X exp( pj=0 βj xij ) ∂L(β) ∂ P Jab (β) = − = yi xij − ni xij =− ∂βa ∂βb ∂βb 1 + exp( pj=0 βj xij ) i=1 P N X xia xib ni exp( pj=0 βj xij ) = h i2 = Pp i=1 1 + exp( j=0 βj xij ) =
N X i=1
xia xib ni p(xi )(1 − p(xi )).
Dospěli jsme ke vztahu, který již neobsahuje informace o náhodných veličinách Z1 , ..., ZN . Tuto vlastnost mají všechny zobecněné lineární modely, ve kterých je použita kanonická linková funkce ([1], str. 148-149). Nyní použijme druhou derivaci k odhadu rozdělení parametru b. Za podmínky splnění předpokladů uvedených v ([2], str. 149-151) platí · ¸ √ 1 d n (b − β) → N 0, . J(β) Úpravou dostaneme
b ∼ N (β, E(X ′ diag [ni p(xi )(1 − p(xi ))] X)−1 ).
(3.6)
Pro názornost uvažujme model s jednou nezávisle proměnnou danou vztahem (3.1) a test hypotézy H0 : β1 = β10 = 0, H1 : β1 6= β10 = 0, 16
(3.7) (3.8)
tedy nezávislost na prediktoru. Waldův test: Známe odhad parametru β1 i jeho rozptyl v11 , kde užitím (3.6) po dosazení odhadů, jak je uvedeno v ([2], str. 177), je )−1 ( n X . ni x2i [ˆ p(xi )(1 − pˆ(xi ))] v11 = i=1
Odhadnutý parametr b1 má za platnosti hypotézy H0 asymptoticky normální rozdělení N (0, v11 ). Testujeme tedy pomocí statistiky b1 S=√ , v11
která má za platnosti hypotézy H0 asymptoticky rozdělení N (0, 1). Poměr věrohodností: Je založen na poměru ℓ1 /ℓ0 , kde ℓ0 značí maximální hodnotu věrohodnostní funkce za platnosti hypotézy H0 a ℓ1 je maximální hodnotou věrohodnostní funkce přes všechny možné hodnoty odhadovaného parametru β. V modelu se užívá asymptoticky χ2 chování výrazu −2[log ℓℓ01 ]. Předpokládáme opět test hypotézy (3.7) proti alternativě (3.8). Dostáváme statistiku −2[log
ℓ0 ] > χ21 (α), ℓ1
pak zamítáme hypotézu H0 ve prospěch alternativy H1 na hladině α. A tudíž −2[log
ℓ0 ] ≤ χ21 (α) ℓ1
je intervalový odhad parametru β o spolehlivosti 1 − α. Test poměrem věrohodností využívá informaci o chování věrohodnostní funkce v odhadnutém bodě b1 a také v bodě β10 . Waldův test v sobě zahrnuje informaci o chování věrohodnostní funkce pouze v odhadnutém bodě b1 , a to vzdálenost od β10 a křivost věrohodnostní funkce (z rozptylu). V sekci 2.5 jsme si test poměrem věrohodností odvodili s pomocí deviance, využijme toho tedy v našem výpočtu. Vyjádřeme příslušnou devianci pro náš model logistické regrese. Používáme značení ze sekce 3.3. Dle definice D(Y ; µ) ˆ = −2 [L1 − LS ] , po dosazení získáme výraz ¸ N · X yi yi − 2 yi log pˆi + (ni − yi ) log(1 − pˆi ) − yi log − (ni − yi ) log(1 − ) ni ni i=1 ¸ N · X yi ni − yi = 2 yi log . + (ni − yi ) log ni pˆi ni − ni pˆi i=1
Došli jsme tedy k interpretaci 2
N X i=1
empirické četnosti × log 17
empirické četnosti . teoretické četnosti
(3.9)
3.5
Diagnostika
Zjistíme-li pomocí zvoleného statistického testu, že náš model není vhodný pro pozorovaná data, nastupuje další fáze rozhodovacího postupu - hledání, kde k tomuto selhání dochází, neboli diagnostika. Diagnostikou se rozumí zejména analýza reziduí. Rezidua Jedná se o srovnání odhadů získaných s pomocí modelu a nasledovaných dat. Při zachování označení z části 3.3 definujme Pearsonovo reziduum yi − ni pˆi rip = p , ni pˆi (1 − pˆi )
i = 1, ..., N.
(3.10)
Dostali jsme další statistiku, totiž Pearsonovu χ2 , neboť podle ([2], str. 270) 2
χ =
N · X (yi − ni pˆi )2 i=1
(ni − yi − ni (1 − pˆi ))2 + ni (1 − pˆi )
ni pˆi
¸
a úpravou získáváme 2
χ =
N · X (yi − ni pˆi )2 i=1
ni pˆi
¸ X N (yi − ni pˆi )2 = (rip )2 . + ni (1 − pˆi ) i=1
Názorné vyjádření Pearsonovu χ2 statistiky pro srovnání s výrazem (3.9) lze v našem případě zapsat následnovně N X (empirické četnosti − teoretické četnosti)2 . teoretické četnosti i=1
Pokud by ni pˆi byla skutečně střední hodnota náhodné veličiny Zi , pak bychom podle centrální limitní věty pro stejně rozdělené náhodné veličiny ([4], str. 101) věděli, že rip mají asymptoticky normované normální rozdělení pro ni → ∞ a tedy χ2 ∼ χ2N −1−k . Avšak ni pˆi je odhad vypočtený pomocí informací o yi , proto má náhodná veličina (3.10) menší rozptyl než normované normální rozdělení. Častěji se proto používá standardizované Pearsonovo reziduum ris = p
rip ˆ ii 1−h
,
i = 1, ..., N,
ˆ ii je i-tý diagonální prvek matice H. kde h c 1/2 X(X ′W c X)−1X ′W c 1/2 , H=W
c je diagonální matice N × N s prvky wii = ni pˆi (1 − pˆi ). Standardizované PearsoW novo reziduum má také asymptoticky normované normální rozdělení. V tomto kontextu pojmem „asymptotickyÿ chápeme, že pro každou hodnotu vysvětlujících proměnných roste počet pozorování do nekonečna. S jeho použitím lze sestavit vhodný test hypotézy a detekovat tak pozorování, která nevyhovují předpokladům. 18
Další možný přístup představuje určení deviančních reziduí definovaných jako p riD = sgn(yi − ni pˆi ) di , i = 1, ..., N,
kde
·
¸ yi ni − yi di = 2 yi log , + (ni − yi ) log ni pˆi ni − ni pˆi
i = 1, ..., N,
je výraz, se kterým jsme se již setkali na konci předchozí kapitoly, proto tedy devianční rezidua. Stejně jako u Pearsonových reziduí, je i zde vhodné definovat standardizované devianční reziduum. riD∗ = p
riD ˆ ii 1−h
.
Pro asymptotické chování deviančních reziduí platí stejné závěry jako pro chování reziduí Pearsonových. Velmi názornou aplikací jsou grafy reziduí v závislosti na vysvětlujících proměnných. Pokud by ovšem vysvětlující proměnné byly spojité, pak by zřejmě téměř každé pozorování bylo v jiném bodě. V takovém případě je třeba postupovat obezřetně a devianci i Pearsonovu χ2 statistiku (označovány jako testy míry dobré shody) používat jen omezeně. Uveďme demonstrační příklad, který mimo jiné osvětlí, proč jsme si v sekci 3.3 mohli dovolit nasčítat úspěchy ve shodných bodech.
3.6
Deviance a její použití
Mějme pro každé i = 1, ..., N dáno ni ∈ N a náhodný výběr {Yi1 , ..., Yini } z alterPN nativního rozdělení v bodě xi . Označme n = i=1 ni . Uvažujme logistický model logit(ˆ pi ) = α + βxi , kde pˆi = P (Yij = 1). Odvoďme hodnotu věrohodnostní funkce: • pracujeme-li s daty jako s výběrem z binomického rozdělení, máme # " µ ¶ X ni ni N X X ni Yij ) log(1 − pˆi ) , Yij log pˆi + (ni − + Lb = log Pni Y ij j=1 j=1 j=1 i=1 • z alternativního pak La =
ni N X X i=1 j=1
=
[Yij log pˆi + (1 − Yij ) log(1 − pˆi )]
"n N i X X i=1
j=1
Yij log pˆi + (ni −
ni X j=1
#
Yij ) log(1 − pˆi ) .
Tedy odhad parametru β je shodný pro obě věrohodnostní funkce. A stejné je to s testem poměrem pravděpodobností, avšak deviance již stejná být nemusí. Uvažujeme-li saturovaný model, pro první alternativu dostáváme věrohodnostní funkci s celkem N parametry (každému pozorování odpovídá jeden parametr). Provedeme-li totéž, 19
ale pracujeme-li s daty jako s výběrem z 0-1 rozdělení, máme pro nasycený model věrohodnostní funkci s n parametry. Nakonec se podívejme na model, kdy ni = 1 pro i = 1, ..., n. Ukažme, jak potom deviance závisí na náhodných veličinách Yi . Vzhledem k faktu, že náhodná veličina Yi může nabývat jen hodnot 0 nebo 1, je D(p; Y ) = −2 = −2 = −2
N X i=1
N X
[Yi log pˆi − (1 − Yi ) log(1 − pˆi ) − Yi log Yi − (1 − Yi ) log(1 − Yi )] [Yi log pˆi + (1 − Yi ) log(1 − pˆi )]
i=1 N · X i=1
¸ pˆi Yi log + log(1 − pˆi ) . 1 − pˆi
P pˆi = pj=1 bj xij a že pro odhad pˆ v našem případě podle (3.5) Již víme, že log 1−ˆ pi platí N X (Yi xij − pˆi xij ) = 0, j = 1, ..., p. i=1
Proto
N µ X i=1
pˆi pˆi − pˆi log Yi log 1 − pˆi 1 − pˆi p
=
X j=1
bj
N X i=1
¶
=
N X i=1
(Yi − pˆi )
p X
bj xij
j=1
(Yi xij − pˆi xij ) = 0.
Záměna pořadí sčítání je v tomto případě bez potíží, neboť oba indexy, do kterých sčítáme, jsou konečné. To znamená, že můžeme devianci vyjádřit vztahem N · X D(p; Y ) = −2 pˆi log i=1
¸ pˆi + log(1 − pˆi ) . 1 − pˆi
V tomto případě tedy deviance závisí jen na odhadnuté pravděpodobnosti, nikoliv přímo na našich pozorováních. Zřejmě bude shodná pro jakékoliv různé soubory dat, které dávají stejnou hodnotu odhadnuté pravděpodobnosti. Avšak na srovnávání míry shody s daty jednotlivých modelů aplikovaných na stejná data použita být může. Řešením tohoto problému je seskupení dat. Pak budeme moci hovořit o asymptotickém chování jednotlivých reziduí, deviance i Pearsonovy χ2 statistiky. Nicméně, může se nám stát, že model bude obsahovat příliš mnoho prediktorů a z různých důvodů nebude možné data seskupit. V tomto případě Hosmer a Lemeshow, jak je uvedeno v ([1] str.176-177), radí postupovat následovně: 1. Data rozdělíme do deseti „shodně velkýchÿ skupin podle předpovězené pravděpodobnosti pˆ(xi ). 20
2. Nechť Yij resp. pˆij značí výsledek pokusu resp. odhadnutou pravděpodobnost pro j-té pozorování v i-té skupině, kde i = 1, ..., 10, j = 1, ..., ni . Pak statistika ³P ´2 Pni ni 10 Y − p ˆ X ij ij j=1 j=1 ³P ´h ³P ´ i ni ni p ˆ 1 − p ˆ /n i=1 i j=1 ij j=1 ij má asymptoticky χ28 .
3.7
Predikční síla modelu
Jedním z měřítek vypovídací kvality modelu je zobecněný koeficient determinace. Jde o další aplikaci věrohodnostní funkce. My použijeme Nagelkerkeho koeficient determinace, odvozen v ([6], str. 125-126), definován jako 2 RN =
1 − [ℓ0 /ℓ1 ]2/n
1 − [ℓ0 /ℓS ]2/n
=
1 − exp(2(L0 − L1 )/n) . 1 − exp(2(L0 − LS )/n)
Vrátíme-li se k věrohodnostní funkci dané vztahem (3.2), vidíme, že ℓ1 je součinem složek, které udávají kvalitu předpovědi v každém bodě. Proto se tento koeficient přiblíží nule, kdykoliv je alespoň jedno pozorování chybně předpovězeno. Jeho nevýhodou je tedy poměrně složitá interpretace a také dost často malá hodnota koeficientu i při aplikaci na celkově kvalitní modely, neboli vysoká citlivost na chybná pozorování. Nakonec uvedeme velmi názorný test kvality modelu, totiž ROC (receiver operating characteristic) křivku. Uvažujme dj , j = 1, ..., k dělení intervalu [0, 1]. Definujme Yˆi = 1, pokud pˆi > dj , jinak Yˆi = 0 pro nějaké j. Pojmenujme senzitivita = P (Yˆ = 1|Y = 1) a specificita = P (Yˆ = 0|Y = 0) P P
(v našem případě) =
Yˆi Yi i=1 Yi
n
=
Pi=1 n
n (1−Yˆi )(1−Yi ) i=1 Pn i=1 (1−Yi )
.
ROC křivkou je pak graf senzitivity v závislosti na (1-specificitě). ROC křivky jsou vhodné ke srovnávání predikční síly modelů, přičemž za relativní ukazatel predikční síly modelů lze považovat veličinu c = plocha pod křivkou. Případ, kdy c = 0.5, signalizuje jen nahodilé odhady narozdíl od modelů s c ≥ 0.8, které považujeme za velmi kvalitní.
21
Kapitola 4 Aplikace modelu logistické regrese Pracujeme s anonymními demografickými daty o 2996 subjektech (klientech), o nichž je znám výsledek splácení v jednoletém horizontu. O každém subjektu máme následujích deset údajů, s jejichž pomocí se pokusíme předpovědět selhání (default) klienta. sex edu depend housing vek term marstat loan income def
pohlaví dosažené vzdělání počet vyživovaných osob typ bydlení klienta věk požadovaná doba splácení úvěru stav (ženatý, svobodný, . . .) výše úvěru (v Kč) příjem, jaký je uveden v daňovém přiznání klienta (v Kč) vysvětlovaná proměnná, 0 - default nenastal, 1 - default (tzn. subjekt 30 a více dní po splatnosti nesplatil dluh) Tabulka 4.1: Přehled dat
default . Symbolem X|Y označme podmíněné Definujme defaultní míru poměrem spláceno rozdělení veličiny X za podmínky Y. Učiníme následující kroky: • Určíme základní charakteristiky datového souboru. • Bude-li to vhodné, provedeme transformaci. Pro samotné hledání modelu je možno použít jeden ze dvou přístupů, vzestupný a sestupný výběr, popsány v ([6] str.111). K tomuto účelu použijeme proceduru step. Procedura step využívá Akaikeho informační kritérium (AIC), kde je užita věrohodnostní funkce a současně zohledněn počet parametrů v modelu. AIC = −2 (maximum log. věrohodnostní funkce − počet parametrů v modelu) 22
• Zanalyzujeme rezidua a aplikujeme teoreticky probrané testy míry dobré shody. Pro srovnání použijeme také probitovou linkovou funkci. Nakonec zhodnotíme modely také z hlediska predikční síly pomocí zobecněného koeficientu determinace a ROC křivky.
4.1
Základní statistické charakteristiky
K dispozici máme celkem devět prediktorů. Rozdělme je do tří základních skupin dle následujících charakteristik: • Intervalové - lze je přirozeně uspořádat a zároveň určit numerickou vzdálenost mezi každými dvěma pozorováními. V našem případě: depend, vek, term, loan a income. • Kategorické – Nominální - nedají se přirozeně uspořádat. Z našich dat do této skupiny spadá sex, housing a marstat. – Ordinální - lze je přirozeně uspořádat, nelze stanovit vzdálenost dvou různých kategorií. Jako ordinální můžeme uvažovat edu. Kategorické proměnné Neuvedeno Odborné Střední všeobecné Vysokoškolské Základní
Default 1 70 93 12 3
Spláceno Výběrový průměr Výběrový rozptyl 6 0.14285714 0.14285714 913 0.07121058 0.06620699 1387 0.06283784 0.05892906 462 0.02531646 0.02472770 49 0.05769231 0.05542986
Tabulka 4.2: Základní vlastnosti vysvětlující proměnné edu
Družstevní byt Jiné Pronajatý dům U rodičů Vlastní dům
Default 12 6 35 20 106
Spláceno 226 59 423 213 1896
Výběrový průměr Výběrový rozptyl 0.05042017 0.04807999 0.09230769 0.08509615 0.07641921 0.07073376 0.08583691 0.07880716 0.05294705 0.05016872
Tabulka 4.3: Základní vlastnosti vysvětlující proměnné housing
23
Default Druh (Družka) 6 Nevyplněno 2 Rozvedený (á) 21 Svobodný (á) 35 Vdovec (Vdova) 6 Ženatý (á) 109
Spláceno Výběrový průměr Výběrový rozptyl 43 0.12244898 0.10969388 6 0.25000000 0.21428571 408 0.04895105 0.04666362 380 0.08433735 0.07741109 49 0.10909091 0.09898990 1931 0.05343137 0.05060127
Tabulka 4.4: Základní vlastnosti vysvětlující proměnné marstat
Female Male Celkem
Default 34 145 179
Spláceno 682 2135 2817
Výběrový průměr Výběrový rozptyl 0.04748603 0.04529437 0.06359649 0.05957811 0.05974633 0.05619546
Tabulka 4.5: Základní vlastnosti vysvětlující proměnné sex Intervalové proměnné Základní data týkající se vysvětlující proměnné vek jsou již uvedena v úvodním motivačním příkladu. Dle histogramu 4.1 nelze předem zamítnout normalitu rozdělení vek|def. Pokud by se skutečně podmíněné rozdělení řídilo normálním rozdělením, pak dle bychom sekce 3.2 do modelu zanesli prediktor vek ve tvaru polynomu druhého stupně. O základních statistických charakteristikách prediktorů loan a income si uděláme představu vyhotovením boxplotů a histogramů. Z boxplotů 4.2 je vidět, že jsou data poměrně rozptýlena. U proměnné income je rozptyl způsoben přítomností právnických osob v datech. Není proto účelné zkoumat přímou závislost poměru defaultu na těchto prediktorech. Vzhledem k malému počtu vzorků, kde jsou uvedeny vysoké příjmy, bychom mohli dojít k chybným závěrům. Naším úkolem je nyní data rozdělit do skupin tak, abychom získali grafy, které budou mít maximální vypovídací schopnost. Pokud bychom se omezili jen na interval hodnot prediktoru income, pro která máme dostatek dat, ztratili bychom podstatnou informaci o zbytku klientů. Jestliže bychom dokázali, že rozdělení income|def resp. loan|def se řídi gamma rozdělením, zanesli bychom do modelu dle sekce 3.2 příslušně transformované prediktory. Z histogramů 4.3 však nelze rozhodnout, zda skutečně jde o gamma rozdělení. V takovém případě je vhodné data transformovat za účelem normalizace podmíněných rozdělení. Použijme transformaci logaritmem. Obrázek 4.4 ukazuje, jak cennou informaci tímto postupem získáme. Histogram transformovaného prediktoru income by naznačoval normální rozdělení, pokud by v našich datech nebylo velké množství klientů s nulovým příjmem. Eliminujeme-li tyto klienty z našich dat, připravíme se o příliš mnoho pozorování. Transformovaná data již lze seskupit a vyčíst přímou závislost defaultní míry na prediktorech income resp. loan.
24
10
20
30
40
50
60
70
Vek
Obrázek 4.1: Histogram proměnné vek a porovnání rozdělení věku klientů, u kterých nastal default (tečkovaně) a u kterých default nenastal
Vyse pujcky
0 e+00
2 e+07
4 e+07
6 e+07
2 e+05 4 e+05 6 e+05 8 e+05 1 e+06
Prijem
Obrázek 4.2: Boxploty pro vysvětlující proměnné income a loan
25
0 e+00
2 e+07
4 e+07
6 e+07
0 e+00
2 e+05
Prijem
4 e+05
6 e+05
8 e+05
1 e+06
Vyse pujcky
Obrázek 4.3: Histogramy proměnných income a loan a porovnání rozdělení příjmu resp. výše půjčky klientů, u kterých nastal default (tečkovaně) a u kterých default nenastal
0
5
10
15
10
log(Prijem)
11
12
13
14
log(Vyse pujcky)
Obrázek 4.4: Histogramy transformovaných proměnných income a loan a porovnání rozdělení příjmu resp. výše půjčky klientů, u kterých nastal default (tečkovaně) a u kterých default nenastal
26
4.2
Hledání modelu
Příprava dat Před vlastním hledáním modelu provedeme vhodnou transformaci dat, kterou lze aplikovat jen na spojité proměnné, v našem případě se jedná o prediktory vek, term, depend, income a loan.
−2.4 −2.6 −2.8 −3.0
Log(def/nedef)
−2.2
−2.0
Logaritmus sanci defaultu
25
30
35
40
45
50
55
Vek
Obrázek 4.5: Závislost logaritmu defaultní míry na proměnné vek Graf 4.5 naznačuje kvadratickou závislost mezi prediktorem vek a logaritmem defaultní míry. Náš odhad dle histogramů byl tedy správný. Nejprve se podíváme na „chudýÿ model pouze s jedním prediktorem vek. Estimate Std. Error z value (Intercept) -1.502277 0.349166 -4.302 vek -0.031524 0.008805 -3.580
Pr(> |z|) 1.69e-05 0.000343
Tabulka 4.6: Odhady parametrů pro model s jediným prediktorem vek Z tabulky 4.6 pro odhad parametrů je patrné, že prediktor vek hraje významnou roli. Test poměrem věrohodností dává dokonce hladinu testu 0.000276, což je ještě silnější důkaz, než jaký poskytuje Waldův test. Analyzujme, zda je vhodné do modelu zanést prediktor vek také v kvadratickém tvaru, jak naznačuje předchozí graf. Tabulka 4.7 dokazuje, že přidání prediktoru vek v kvadratickém tvaru, vede ke značnému zlepšení kvality modelu. Na hladině pět procent zamítáme hypotézu, že vek2 nemá vliv na vyšší kvalitu modelu. S pomocí testu poměrem věrohodnostní se ukázalo, že takováto transformace je nejvhodnější. Následující grafy s intervalovými odhady také poměrně dobře zdůvodňují použití prediktoru vek ve tvaru polynomu 27
vek poly(vek,2)
Resid. Df 2994 2993
Resid. Dev 1342.58 1336.63
Df
Deviance
P(> |Chi|)
1
5.95
0.01
Tabulka 4.7: Srovnání dvou modelů testem poměrem věrohodností druhého stupně. Body v grafech jsou skutečné hodnoty poměru defaultu, které jsou uvedeny v tabulce v motivačním příkladu.
0.12 0.10 0.08
Pravdepodobnost defaultu
0.04
0.06
0.12 0.10 0.08 0.06 0.04
Pravdepodobnost defaultu
0.14
Model:def~poly(vek, 2)
0.14
Model:def~vek
20
30
40
50
60
20
Vek
30
40
50
60
Vek
Obrázek 4.6: Grafické srovnání modelů s prediktorem vek Prediktor term nabývá 14-ti hodnot a téměř pro polovinu z nich máme méně než tři pozorování. Není proto vhodné zkoumat přímou závislost defaultní míry na tomto prediktoru. Lepším přístupem je nejprve data seskupit do čtyř skupin vždy po jednom roce a pak konstruovat příslušný graf. Na základě grafu 4.7, vzhledem k malému počtu bodů, nelze zamítnout linearitu. Do modelu zaneseme prediktor term v lineárním tvaru. V případě prediktoru depend již nejsme omezeni množstvím různých hodnot, jak tomu bylo u předchozích intervalových proměnných. Transformace při práci s takovouto proměnnou, vzhledem k malému počtu různých hodnot, nemá příliš smysl a tak si vysvětleme jiný přístup. Při práci se spojitou proměnnou můžeme data rozdělit do skupin a uvažovat tuto proměnnou jako kategorickou. Každý interval pak lze považovat za jednu kategorii. Je také možné sloučit kategorie se shodnou defaultní mírou. Na základě grafu 4.7 se nabízí možnost zakódovat proměnnou depend jako kategorickou. Ukázalo se, že tato volba vede ke ztrátě signifikantní interakce, tedy do modelu zaneseme depend v lineárním tvaru. Z grafů 4.8 plyne, že mezi prediktory income resp. loan a def není výrazná závislost. Bylo by však chybou je do modelu nezařadit, mohli bychom se připravit o interakci zvyšující kvalitu modelu. 28
Logaritmus sanci defaultu
−1.7 −1.9
−1.8
Log(def/nedef)
−0.96 −0.98
−2.0
−1.02
−1.00
Log(def/nedef)
−0.94
−0.92
Logaritmus sanci defaultu
25
30
35
40
45
50
55
60
0
1
Doba splatnosti
2
3
4
5
Pocet vyzivovanych osob
Obrázek 4.7: Závislost logaritmu defaultní míry na proměnných term a depend
Logaritmus sanci defaultu
−1.9 −2.1
−2.0
Log(def/nedef)
−1.4 −1.6
−2.2
−1.8
Log(def/nedef)
−1.2
−1.8
−1.0
−1.7
Logaritmus sanci defaultu
0
5
10
15
11.0
log(Prijem)
11.5
12.0
12.5
13.0
13.5
log(Vyse pujcky)
Obrázek 4.8: Závislost logaritmu defaultní míry na proměnných income a loan
29
Zkusme si data ještě trochu zjednodušit. Prediktor edu lze přirozeně kódovat jako ordinální.
0.20
Model: Default~poly(vek, 2) + dosazene vzdelani
0.15 0.10
Stredni vseobecne Odborne Vysokoskolske Zakladni
0.05
Pravdepodobnost defaultu
Neuvedeno
20
30
40
50
60
Vek
Obrázek 4.9: Závislost pravděpodobnosti defaultu na proměnné edu Nejvhodnější asi bude • klienty s vysokoškolským vzděláním kódovat číslem 1, • klienty, kteří kolonku nevyplnili, číslem 3, • zbylé číslem 2. poly(vek,2)+edu.ord poly(vek,2)+edu
Resid. Df 2992 2989
Resid. Dev 1326.85 1326.15
Df
Deviance
P(> |Chi|)
3
0.71
0.87
Tabulka 4.8: Srovnání modelu s ordinální a kategorickou proměnnou edu testem poměrem věrohodností Test poměrem věrohodností dokazuje, že model nepokazíme, pokud budeme vzdělání považovat za ordinální. Mohli bychom se ještě podívat na grafy interakcí, my se ale spoléháme, že procedura step významné interakce odhalí.
30
Celkové modely 1. Nejprve užijme vzestupný výběr na náš základní model obsahující jen transformovaný prediktor vek. Step
Df Deviance NA NA + edu.ord -1 9.772400 + term -1 6.635964 + poly(vek, 2):edu.ord -2 6.255951
Resid. Df 2993 2992 2991 2989
Resid. Dev 1336.625 1326.853 1320.217 1313.961
AIC 1342.625 1334.853 1330.217 1327.961
Tabulka 4.9: Vzestupný výběr pomocí procedury step 2. Druhým možným přístupem, který statistici preferují, je sestupný výběr. Vzhledem k velkému počtu prediktorů můžeme sestupný výběr aplikovat maximálně na „bohatýÿ model obsahující všechny prediktory i s interakcemi do prvního řádu. Je dobré si uvědomit, že pokud by sestupný výběr objevil nějaké signifikantní interakce, následující postup je nalezne také. 3. Nakonec na data aplikujeme výběr oběma směry. Výchozím modelem bude model obsahující všechny prediktory. Step
Df NA + housing:income -4 + depend:marstat -5 + poly(vek, 2):edu.ord -2 - loan 1 - sex 1
Deviance Resid. Df NA 2978 16.27439631 2974 16.85658130 2969 6.05634541 2967 0.01845206 2968 0.97934056 2969
Resid. Dev 1306.233 1289.959 1273.102 1267.046 1267.064 1268.043
AIC 1342.233 1333.959 1327.102 1325.046 1323.064 1322.043
Tabulka 4.10: Výběr oběma směry pomocí procedury step Tento model již vypadá velmi kvalitně. Ukázalo se, že ač se zdálo vhodnější do modelu nezanést prediktor income, v konečném modelu se objevil. Důvodem je silná interakce housing:income. Dále je na první pohled zřejmé, že pokud bychom použili jen vzestupný výběr, připravili bychom se o několik cenných interakcí. Výběr se totiž ke zkoumání interakcí ani nedostal, protože zamítl přidání samotných prediktorů. Podíváme se ještě na další diagnostiky, které naši teorii o horší kvalitě prvního modelu jen potvrdí.
31
NULL poly(vek, 2) depend housing term marstat income edu.ord housing:income depend:marstat poly(vek, 2):edu.ord
Df
Deviance
2 1 4 1 5 1 1 4 5 2
19.18 0.71 3.42 7.66 8.60 0.03 8.88 16.30 17.14 5.85
Resid. Df 2995 2993 2992 2988 2987 2982 2981 2980 2976 2971 2969
Resid. Dev 1355.80 1336.63 1335.91 1332.50 1324.84 1316.24 1316.21 1307.33 1291.03 1273.89 1268.04
P(> |Chi|) 6.846e-05 0.40 0.49 0.01 0.13 0.87 2.878e-03 2.642e-03 4.248e-03 0.05
Tabulka 4.11: Analýza rozptylu pro model získaný výběrem oběma směry užitím procedury step
4.3
Diagnostika
Pracujeme se dvěma modely. Označme si je podle výběru, jakým jsme je nalezli M.for a M.both. Zkusme ještě na model M.both použít alternativní linkovou funkci probit. Následně na uvedené modely aplikujme analýzu reziduí a testy probrané v části 3.5. Pro srovnání jejich kvality nám poslouží Hosmer-Lemeshow statistika, zobecněný koeficient determinace, plocha pod ROC křivkou a pro úplnost uveďme ještě Pearsonovu χ2 statistiku. χ2 Resid. Dev M.both 3003.833 1268.043 M.for 2955.892 1313.961 M.probit 3025.321 1269.348
H-L df 7.705413 8 3.172338 8 4.748138 8
P-value 0.462762 0.923081 0.784125
Nagel.R2 c 0.079309 0.68 0.038103 0.632 0.078147 0.681
Tabulka 4.12: Základní testy kvality Analýza reziduí nepotvrdila porušení předpokladů. Rezidua, až na několik málo výjimek, jsou soustředěna kolem dvou přímek. Na základě Hosmer-Lemeshow statistiky nelze zamítnout dobrou shodu s daty ani u jednoho modelu. Podle Nagelkerkeho kritéria vychází jako lepší model M.both, jak jsme již tušili vzhledem k rozdílné devianci v tabulkách 4.9 a 4.11. Je zajímavé, že pokud bychom modely posuzovali podle Pearsonovy χ2 statistiky, dostali bychom přesně opačný výsledek, než jaký dává deviance. My jsme model hledali s pomocí deviance, a tak i teď použijme toto kritérium, tedy model M.both preferujeme před modelem M.probit. Model s probitovou linkovou funkcí se téměř ve všech směrech ukázal jako méně kvalitní. Nepřisuzujme to ale kvalitě metody, jde spíše o to, že data jsme transformovali na model logistické regrese nikoliv probitové. 32
Analyza standardizovanych rezidui pro model M.both
xx
xx
x
x
2 −1
2
2141 x
1
6
x xx x
4
x
1645 x
0
10 8
x
x
800 x
3
2141 x
1645 x
Standardizovana deviancni rezidua
12
800 x
0
Standardizovana pearsonova rezidua
Analyza standardizovanych rezidui pro model M.both
0
500
1000
1500
2000
2500
3000
0
500
1000
Index
1500
2000
2500
3000
Index
Obrázek 4.10: Analýza reziduí pro model M.both
x
Analyza standardizovanych rezidui pro model M.for
x
6
x x
x
−1
2
4
xx
2
x
1
x
0
Standardizovana deviancni rezidua
1645 x
0
Standardizovana pearsonova rezidua
3
10
800 x
8
Analyza standardizovanych rezidui pro model M.for
0
500
1000
1500
2000
2500
3000
0
Index
500
1000
1500 Index
Obrázek 4.11: Analýza reziduí pro model M.for
33
2000
2500
3000
Č. poz. χ2 Resid. Dev - 800 2940.365 1257.932 - 1645 2879.881 1248.087 - 2141 2824.129 1235.678
H-L df P-value 4.712761 8 0.787789 5.506431 8 0.702327 8.31802 8 0.403038
Nagel.R2 c 0.083568 0.684549 0.087603 0.688125 0.093948 0.696184
Tabulka 4.13: Vliv odlehlých pozorování v modelu M.both Tabulka 4.13 ilustruje, jak se mění kvalita modelu při postupném odebírání odlehlých pozorování. Výrazná změna nastala u Nagelkerkeho koeficientu determinace. Jeho hodnota se oproti původní zvýšila téměř o 20% a to po eliminaci pouhých 3 pozorování z 2996-ti, což potvrzuje naši teorii o vysoké citlivosti tohoto koeficientu. Markantní změna veličiny χ2 nám poskytuje možné vysvětlení její vysoké hodnoty v případě modelu M.both uvedené v tabulce 4.12. Na její velikost mají vliv především odlehlá pozorování, neboť v Pearsonově χ2 statistice vystupují rozdíly empirických četností a jejich teoretických četností umocněné na druhou. Pro srovnání, při eliminaci pozorování 800 a 1645, avšak pro model M.for, dostaneme hodnotu Pearsonovy χ2 statistiky 2900.587. Tato hodnota je výrazně vyšší než v modelu M.both po stejné úpravě dat (viz tabulka 4.13). Poznamenejme ještě, že není pravdou, že jsme eliminovali jen ta pozorování, která vykazují vysoké absolutní hodnoty Pearsonových reziduí pouze v modelu M.both. Pozorování číslo 800 a 1645 mají nejvyšší absolutní hodnoty těchto reziduí i pro model M.for, jak je ukázáno v grafu 4.11. Graf 4.12 dokazuje vyšší diskrimační sílu modelu M.both.
0.6 0.4
M.both 0.68 M.for 0.648
0.0
0.2
Senzitivita
0.8
1.0
Srovnani modelu M.both a M.for pomoci ROC krivky
0.0
0.2
0.4
0.6
0.8
1.0
1−Specificita
Obrázek 4.12: Srovnání predikční síly modelů pomocí ROC křivky
34
Kapitola 5 Závěr Tato bakalářská práce se zabývala zobecněnými lineárními modely. Popsala tyto modely a uvedla postup odhadu parametrů. Definovala pojem deviance a popsala její možné užití. Dále se pak práce zaměřila na model logistické regrese. Detailně popsala vlastnosti logistické funkce. Odvodila odhad metodou maximální věrohodnosti, srovnala dva používané testy nulovosti parametrů a uvedla příslušné intervalové odhady. Ovšem největší důraz kladla na diagnostiku. Vedle popisu reziduí a příslušných testů míry dobré shody ukázala, za jakých okolností se tyto dají aplikovat. Následně se také zmínila o charakteristikách, které hodnotí predikční sílu modelu. Poslední část práce byla věnována aplikaci teoretických poznatků na konkrétní data. S pomocí statistického software R uvedla základní charakteristiky prediktorů, nalezla vhodné modely a provedla diagnostiku. Při volbě prediktorů ukázala výhodnost alternativního grafického přístupu. Aplikovala standardní testy míry dobré shody a vyhodnotila predikční sílu nalezených modelů. Na základě obdržených výsledků navrhla vhodný model. Nakonec naznačila, jak by se dal model vylepšit eliminací odlehlých pozorování. Námětem k dalšímu studiu by mohla být hlubší analýza možností transformace spojitých vysvětlujících proměnných. Bylo by také zajímavé zjistit, do jaké míry je vhodná eliminace odlehlých pozorování. Eliminace odlehlých pozorování vede ke zvýšení kvality modelu, avšak na úkor ztráty informace získané z těchto pozorování.
35
Literatura [1] Agresti A. (1990): Categorical Data Analysis, John Wiley & Sons. [2] Anděl J. (2005): Základy matematické statistiky, Univerzita Karlova v Praze, Matematicko-fyzikální fakulta. [3] Cook R.D., Weisberg S. (1999): Applied Regression Including Computing and Graphics, John Wiley & Sons. [4] Lachout P. (2004): Teorie pravděpodobnosti, učební text pro posluchače MFF UK. [5] Marhoun P. (2005): Skóringové a klasifikační metody v bankovnictví, diplomová práce, MFF UK Praha. [6] Zvára K. (2003): R a Regrese, učební text pro studenty, KPMS MFF UK Praha.
36
