Osnova přednášky „Binární logistická regrese“ 1. Motivace 2. Odvození modelu 3. Kódování proměnných 3.1. Případ dvou kategorií 3.2. Případ aspoň tří kategorií 4. Význam parametrů 5. Odhady parametrů 6. Intervaly spolehlivosti 6.1. Interval spolehlivosti pro regresní parametr 6.2. Interval spolehlivosti pro podíl šancí 7. Dílčí testy významnosti regresních parametrů 7.1. Waldův test 7.2. Test poměrem věrohodnosti 7.3. Skórový test 8. Test významnosti modelu jako celku 9. Výstavba modelu 9.1. Principy výstavby modelu 9.2. Výběr podmnožiny vysvětlujících proměnných 10. Hodnocení modelu z různých hledisek 10.1. Testy dobré shody 10.2. Koeficienty determinace 10.3. Informační kritéria 10.4. Klasifikační tabulka 10.5. ROC křivka 11. Příklad Binární logistická regrese 1. Motivace Tato metoda umožňuje odhad pravděpodobnosti nastoupení nějakého jevu (zapíšeme ho pomocí náhodné veličiny Y jako {Y = 1}) pomocí k známých regresorů X1, …, Xk, které mohou být jak spojitého, tak kategoriálního typu. Byla vytvořena v 60. letech 20. století. Použití v praxi: v medicíně, veličina Y popisuje přítomnost či nepřítomnost nějaké choroby; v bankovnictví, veličina Y popisuje splácení či nesplácení úvěru; v marketingových kampaních, veličina Y popisuje odezvu na reklamu nějakého výrobku; v pojišťovnictví, veličina Y popisuje uplatnění či neuplatnění pojistného nároku … 2. Odvození modelu Uvažme závisle proměnnou náhodnou veličinu Y, která nabývá hodnoty 1 s pravděpodobností ϑ a hodnoty 0 s pravděpodobností ϑ−1 , tj. ( )ϑA~Y a její pravděpodobnostní funkce má tvar: ( ) ( ) 0,1ypro jinak0 1 y 1y =    ϑ−ϑ =π ϑ− . Jev {Y = 1} často interpretujeme jako úspěch, jev {Y = 0} jako neúspěch. Předpokládejme, že máme k vysvětlujících proměnných X1, …, Xk. Označme ( )T k1 X,,X,1 K=X vektor vysvětlujících proměnných s absolutním členem a ( )T k10 ,,, βββ= Kβ vektor regresních koeficientů. Pro predikci hodnot veličiny Y nelze použít lineární regresní model tvaru βXT kk110 XXY =β++β+β= K , protože na levé straně jsou jen 0 a 1, zatímco pravá strana může nabývat jakoukoli reálnou hodnotu. Budeme tedy modelovat nikoliv hodnoty veličiny Y, ale pravděpodobnost úspěchu ϑ (tj. střední hodnotu veličiny Y) (za předpokladu, že známe hodnoty vektoru vysvětlujících proměnných). Pokud bychom využili model ( ) βXxX T 1YP === , mohlo by se stát, že některé predikované pravděpodobnosti úspěchu (při daném x) by ležely vně intervalu (0,1). Tento problém lze částečně řešit zavedením šance: ( ) ( ) ( ) ( ) ( )xX xX xX xX x ==− == = == == =ω 1YP1 1YP 0YP 1YP . Šance vyjadřuje, kolikrát je při daném x vyšší pravděpodobnost úspěchu než neúspěchu. Nabývá hodnot z intervalu (0,∞). Nyní je zapotřebí vzájemně jednoznačně transformovat interval (0,∞) na interval (-∞,∞). K tomuto účelu použijeme přirozený logaritmus šance: ( ) ( ) ( )xX xX x ==− == =ω 1YP1 1YP lnln (jde o tzv. logitovou transformaci pravděpodobnosti úspěchu za předpokladu X = x, zkráceně logit). Logaritmickou šanci již můžeme modelovat pomocí lineárního regresního modelu: ( ) ( ) βX xX xX T 1YP1 1YP ln = ==− == . Odtud můžeme vyjádřit šanci ( ) βXT x e=ω a dále podmíněnou pravděpodobnost úspěchu: ( ) βXβX βX TT T xX − + = + === e1 1 e1 e 1YP resp. podmíněnou pravděpodobnost neúspěchu: ( ) ( ) βXβX βX TT T xXxX e1 1 e1 e 11YP10YP + = + −===−=== , celkem ( ) y1y e1 1 1 e1 1 yYP − −−       + −      + === βXβX TT xX pro y = 0, 1. Tímto vztahem tedy modelujeme pravděpodobnost úspěchu či neúspěchu v závislosti na realizacích x1, …, xk. Upozornění: pravděpodobnost úspěchu, šance na úspěch a logit úspěchu jsou tři různé způsoby vyjádření téhož v tom smyslu, že jsou na sebe vzájemně převoditelné. Pro interpretaci jsou vhodnější pravděpodobnosti a šance než logity. 3. Kódování kategoriálních proměnných 3.1. Případ dvou kategorií V tomto případě kategorie vysvětlující proměnné X kódujeme nejčastěji pomocí 0 a 1. Např. proměnná X udává pohlaví pacienta. Zvolíme X = 0 pro ženu a X = 1 pro muže. 3.2. Případ aspoň tří kategorií Vysvětlující proměnná X má r ≥ 3 kategorií, např. X udává barvu vlasů a má tři kategorie: černá, blond, rezavá. 3.2.1. Kódování pomocí referenční kategorie: zavedeme r-1 nezávislých indikátorových proměnných (indikátorů) Z1, …, Zr-1. Každý z těchto indikátorů odpovídá jedné kategorii proměnné X, vynechaná kategorie je referenční, kódujeme ji pomocí 0 a můžeme si ji sami zvolit.    = jinak0 vlasyčernépro1 Z1 ,    = jinak0 vlasyblondpro1 Z2 . Vyjádřeno tabulkou: černé vlasy blond vlasy rezavé vlasy Z1 1 0 0 Z2 0 1 0 V tomto případě roli referenční kategorie hrají rezavé vlasy, jsou kódovány 0. 3.2.2. Kódování pomocí sigma omezené parametrizace: na rozdíl od předešlého případu se referenční kategorie kóduje pomocí -1.    = jinak0 vlasyčernépro1 Z1 ,    = jinak0 vlasyblondpro1 Z2 . Vyjádřeno tabulkou: černé vlasy blond vlasy rezavé vlasy Z1 1 0 -1 Z2 0 1 -1 Součet v každém řádku tabulky je 0. 3.2.3. Kódování přeparametrizovaného modelu: v tomto případě máme r závislých indikátorů Z1, …, Zr.    = jinak0 vlasyčernépro1 Z1 ,    = jinak0 vlasyblondpro1 Z2 ,    = jinak0 vlasyezavérpro1 Z3 . Vyjádřeno tabulkou: černé vlasy blond vlasy rezavé vlasy Z1 1 0 0 Z2 0 1 0 Z3 0 0 1 4. Význam parametrů Ze vztahu pro logit ( ) ( ) kk110 T XX 1YP1 1YP ln β++β+β== ==− == KβX xX xX plyne: parametr β0 udává velikost logitu pro nulové hodnoty všech vysvětlujících proměnných. Je-li β0 = 0, je logaritmus šance 0, tedy šance = 1 neboli pravděpodobnost úspěchu je 0,5. Pro β0 > 0 je šance na úspěch větší než 0,5 a pro β0 < 0 je šance na úspěch menší než 0,5. Pro interpretaci parametrů βj zavedeme podíl šancí: ( ) ( ) ( ) j e x,,x,,x x,,1x,,x xOR kj1 kj1 j β == ω +ω = K KK KK . Jednotková změna j-té vysvětlující proměnné znamená v průměru j e β násobnou změnu šance na úspěch, zůstanou-li všechny ostatní vysvětlující proměnné stejné. U kategoriálních proměnných závisí interpretace parametrů na způsobu, jakým kódujeme kategorie. Parametry u indikátorových proměnných vyjadřují po odlogaritmování příslušné násobky šance na úspěch v referenční kategorii. 5. Odhady parametrů Pro odhad parametrů v logistickém regresním modelu musíme mít k dispozici n > k nezávislých pozorování y1, …, yn závisle proměnné veličiny a příslušných regresorů xi1, …, xik, i = 1, …, n. Tato pozorování získáme na n objektech. V logistickém regresním modelu nelze kvůli charakteru závisle proměnné veličiny Y použít metodu nejmenších čtverců. Odhady parametrů hledáme metodou maximální věrohodnosti. Logaritmická věrohodnostní funkce pro model binární logistické regrese má tvar: ( ) ( ) ( ) ( ) ( )[ ]∑ ∑∏ = = − − − − = − − − +−−= + = + = n 1i n 1i T ii y1 n 1i y1 i i i i i i i i e1ln1y e1 e ln e1 e ln; βx βx βx βx βx T T T T T βxxβl . Řešením systému věrohodnostních rovnic ( ) 0 ; j i = β∂ ∂ xβl , j = 0, 1, …, k získáme odhady k10 ˆ,,ˆ,ˆ βββ K . Toto řešení nelze obecně nalézt v algebraickém tvaru, proto se hledá numericky, např. použitím Newtonovy – Raphsonovy iterační metody. Pro každý objekt pak můžeme vypočítat odhad pravděpodobnosti úspěchu či neúspěchu: ( ) ( ) βx βx T i T i xX ˆ y1ˆ iiii e1 e yYPˆ i − − − + === . Tomuto odhadu se říká skóre a značí se i ˆϑ . 6. Intervaly spolehlivosti Směrodatná chyba odhadu regresního parametru βj je ( ) jjj vˆse =β , kde vjj je j-tý diagonální prvek matice ( ) 1 lj 2 ;ˆ ˆ −         β∂β∂ ∂ −= xβ V l . Hodnoty ( )j ˆse β větší než 2 indikují numerické problémy, např. multikolonearitu mezi vysvětlujícími proměnnými nebo u kategoriálních proměnných nulové zastoupení objektů v některé kategorii. Toto upozornění se však nevztahuje na odhad směrodatné chyby odhadu β0. 6.1. Interval spolehlivosti pro regresní parametr 100(1-α)% asymptotický interval spolehlivosti pro regresní parametr βj má meze: ( ) ( )( )2/1jj2/1jj uˆseˆ,uˆseˆ α−α− β+ββ−β , j = 0, 1, …, k 6.2. Interval spolehlivosti pro podíl šancí Jestliže se j-tá vysvětlující proměnná zvětší o ∆ (a ostatní vysvětlující proměnné se nezmění), pak podíl šancí je j e β∆ a 100(1-α)% asymptotický interval spolehlivosti pro podíl šancí je ( ) ( ) ( )2/1jj2/1jj uˆseˆuˆseˆ e,e α−α− β+ββ−β , j = 0, 1, …, k . 7. Dílčí testy významnosti regresních parametrů Na hladině významnosti α testujeme hypotézu 0:H j0 =β proti 0:H j1 ≠β , j = 1, 2, …, k. Nulová hypotéza tvrdí, že j-tá vysvětlující proměnná je v modelu zbytečná. 7.1. Waldův test Testová statistika ( )j j 0 ˆse ˆ T β β = se za platnosti H0 asymptoticky řídí rozložením N(0,1) (tedy statistika ( ) 2 j j2 0 ˆse ˆ T         β β = se za platnosti H0 asymptoticky řídí rozložením χ2 (1)). Nulovou hypotézu zamítáme na asymptotické hladině významnosti α, když 2/10 uT α−≥ resp. ( )1T 1 22 0 α−χ≥ . 7.2. Test poměrem věrohodnosti Zavedeme několik nových pojmů. Saturovaný model S … má odlišný parametr pro každé pozorování (má tedy n parametrů a sám o sobě není použitelný). Každý model je jeho podmodelem. Jeho logaritmickou věrohodnostní funkci označme Sl . Námi zkoumaný model M … zahrnuje k regresorů X1, …, Xk. Jeho logaritmickou věrohodnostní funkci označme Ml . Přiléhavost modelu M k datům lze posoudit pomocí deviance ( )MSM 2D ll −= . Přiléhavější model než saturovaný model však neexistuje, proto Sl = 0 a tudíž MM 2D l−= neboli ( ) ( )[ ]∑= − +−−= n 1i ˆT iiM i e1lnˆ1yD βx T βx . Čím je model méně přiléhavý, tím je jeho deviance vyšší. Model Mj … zahrnuje k-1 regresorů X1, …, Xj-1, Xj+1, …, Xk. Jeho devianci označme jMD Pro test hypotézy, že j-tá vysvětlující proměnná je v modelu zbytečná, použijeme testovou statistiku MM0 DDT j −= , která se za platnosti H0 asymptoticky řídí rozložením χ2 (1). Nulovou hypotézu zamítáme na asymptotické hladině významnosti α, když ( )1T 1 2 0 α−χ≥ . Upozornění: Test poměrem věrohodnosti lze použít i v případě, kdy porovnáváme model M a jeho podmodel M ~ . Předpokládáme, že model M má k1 regresorů a model M ~ má k2 < k1 regresorů. Testujeme hypotézu, že dostačující je model M ~ , tedy k1 - k2 regresorů je nulových proti alternativě, že aspoň jeden je nenulový. Testová statistika MM ~0 DDT −= se za platnosti H0 asymptoticky řídí rozložením χ2 (k1 - k2). Nulovou hypotézu zamítáme na asymptotické hladině významnosti α, když ( )211 2 0 kkT −χ≥ α− . 7.3. Skórový test Testová statistika ( ) ( ) 0 22 2 0 0 L E L T =β =β       β∂ β∂ − β∂ β∂ = se za platnosti H0 asymptoticky řídí rozložením χ2 (1). Nulovou hypotézu zamítáme na asymptotické hladině významnosti α, když ( )1T 1 2 0 α−χ≥ . 8. Test významnosti modelu jako celku Na hladině významnosti α testujeme hypotézu 0:H k10 =β==β K proti H1: aspoň jeden parametr 0j ≠β . Nulová hypotéza tvrdí, že dostačující je tzv. nulový model obsahující pouze parametr β0, tj. model ( ) 0 e1 1 1YP β xX − + === . Označme D0 devianci nulového modelu a DM devianci našeho modelu s k regresory. Testová statistika T0 = D0 - DM se za platnosti H0 asymptoticky řídí rozložením χ2 (k). Nulovou hypotézu zamítáme na asymptotické hladině významnosti α, když ( )kT 1 2 0 α−χ≥ . 9. Výstavba modelu binární logistické regrese Při výstavbě modelu rozhodujeme, které z vysvětlujících proměnných X1, …, Xk jsou důležité pro vysvětlení proměnné Y. Nejprve je vhodné se zabývat významností modelu jako celku, tj. otestovat, zda je vůbec možné z daných k proměnných vytvořit model, který bude lepší než model konstanty. Použijeme test poměrem věrohodnosti, jehož testová statistika je rozdílem deviancí nulového modelu a modelu s k regresory. Není-li na dané hladině významnosti nulová hypotéza zamítnuta, nemá smysl se modelem zabývat. 9.1. Principy výstavby modelu Jde o vytvoření takového modelu, který bude obsahovat co nejmenší množství proměnných (resp. jejich kombinací) a přitom bude ještě dostatečně dobře vysvětlovat zkoumaná data. Na začátku celého procesu se doporučuje prozkoumat vztah mezi vysvětlovanou veličinou a každou vysvětlující veličinou zvlášť. U spojitých proměnných použijeme dvouvýběrový test. U kategoriálních proměnných provedeme test nezávislosti v kontingenční tabulce. Je-li některá četnost v kontingenční tabulce nulová, budou konečné výstupy modelu s takovou proměnnou obsahovat nesmyslné hodnoty. Poměr šancí totiž bude kvůli dosazení nuly do vzorce buď nula nebo nekonečno. Této situaci zabráníme, když logicky sloučíme varianty této proměnné nebo – je-li to možné – variantu s nulovou četností vyloučíme. 9.2. Výběr podmnožiny vysvětlujících proměnných a) Ruční postup 1. krok: Provedeme jednorozměrnou analýzu pro všechny vysvětlující proměnné a do modelu zahrneme ty, pro které test významnosti poskytne p-hodnotu menší než 0,25. (Ignoruje se původní k-rozměrná struktura dat.) 2. krok: Postupně vynecháváme proměnné, které jsou vysoce nevýznamné. 3. krok: Do modelu zahrneme ty interakce mezi proměnnými, které mají věcný smysl a jejichž p-hodnota je nejvýše 0,05. Může se stát, že při zahrnutí interakce do modelu se jedna ze vstupních proměnných stane nevýznamnou. Je lépe ji v modelu ponechat. b) Automatizovaný postup Používají se krokové (stepwise) metody – dopředná nebo zpětná. 10. Hodnocení vytvořeného modelu z různých hledisek 10.1. Testy dobré shody Nulová hypotéza tvrdí, že naměřené a predikované hodnoty se neliší. Jsou založeny na porovnání naměřených hodnot yi a odhadnutých skóre i ˆϑ , i = 1, 2, …, n. Pearsonův χ2 test Testová statistika má tvar ( ) ( )∑= ϑ−ϑ ϑ− =χ n 1i ii 2 ii2 ˆ1ˆ ˆy a za platnosti nulové hypotézy se asymptoticky řídí rozložením ( )1kn2 −−χ . Hypotézu o shodě dat s modelem tedy zamítáme na asymptotické hladině významnosti α, když ( )1kn1 22 −−χ≥χ α− . Hosmerův – Lemeshowův test Tento test je preferován pro spojité nezávisle proměnné a vyžaduje dostatečně rozsáhlý datový soubor. Datový soubor je uspořádán vzestupně podle skórů i ˆϑ a je rozdělen do G (zpravidla G = 10, musí být splněna podmínka G > k + 1) přibližně stejně velkých skupin. V každé z těchto skupin se zjišťuje skutečný a očekávaný počet objektů, pro něž Y = 1 resp. Y = 0. Označme ng rozsah g-té skupiny, ∑= = gn 1i ig1 yo resp. ( )∑= −= gn 1i ig0 y1o skutečný počet objektů v gté skupině, pro něž Y = 1 resp. Y = 0. Analogicky ∑= ϑ= gn 1i ig1 ˆe resp. ( )∑= ϑ−= gn 1i ig0 ˆ1e je očekávaný počet objektů v g-té skupině, pro něž Y = 1 resp. Y = 0. Testová statistika ( ) ∑∑= = − = 1 0k G 1g kg 2 kgkg 0 e eo T se za platnosti nulové hypotézy asymptoticky řídí rozložením ( )2G2 −χ . Nulovou hypotézu zamítáme na asymptotické hladině významnosti α, když ( )2GT 1 2 0 −χ≥ α− . Pro korektní použití H-L testu je nutné, aby všechny teoretické četnosti byly větší než 1 a většina z nich musí být větší než 5. 10.2. Koeficienty determinace Tyto koeficienty porovnávají nulový model s deviancí D0 a náš model s devianci DM. McFaddenův koeficient determinace: 0 M2 MF D D 1R −= . Tento koeficient se získá prostým dosazením deviance místo příslušných součtů čtverců do vztahu pro koeficient determinace u lineární regrese. Coxové – Snellův koeficient determinace: ( ) n/DD2 CS 0M e1R − −= Nevýhodou tohoto koeficientu je, že nemůže překročit hodnotu n/D0 e1 − − , je tedy vždy menší než 1, což ztěžuje interpretaci. Nagelkerkův koeficient determinace: ( ) n/D n/DD 2 N 0 0M e1 e1 R − − − − = Nagelkerkův koeficient vznikne z Coxové – Snellova koeficientu vydělením maximální možnou hodnotou n/D0 e1 − − . Čím je posuzovaný model M více vzdálen od nulového modelu, tím jsou koeficienty determinace vyšší. 10.3. Informační kritéria Informační kritéria slouží k porovnání modelů (vytvořených na týchž datech) s různým počtem regresorů. S rostoucím počtem regresorů roste i hodnota logaritmické věrohodnostní funkce (a tím i „důvěryhodnost“ modelu), na druhé straně však velký počet regresorů nemusí být vždy vhodný, např. z ekonomického hlediska. Informační kritéria jsou navržena tak, aby penalizovala velký počet regresorů. Za lepší je považován model, který poskytuje nižší hodnotu informačního kritéria. Označme l hodnotu logaritmické věrohodnostní funkce nějakého modelu, který má k regresorů a byl vytvořen na základě datového souboru rozsahu n. Akaikeovo informační kritérium: k22AIC +−= l Bayesovo informační kritérium: nlnk2BIC +−= l 10.4. Klasifikační tabulka Do klasifikační tabulky zaznamenáváme počty správně a nesprávně zařazených objektů: skutečnostpredikce Y = 1 Y = 0 celkem Y = 1 a b a+b Y = 0 c d c+d celkem a+c b+d n Na hlavní diagonále je tedy počet objektů, které model správně predikoval. Relativní četnost správně predikovaných objektů je n da + . Abychom mohli sestavit tuto klasifikační tabulku, musíme stanovit tzv. dělicí bod C pro odhadnutou pravděpodobnost úspěchu i ˆϑ , i = 1, …, n. Můžeme volit jakoukoli hodnotu z intervalu (0, 1), zpravidla však predikujeme Y = 1 pro 5,0ˆ i ≥ϑ a Y = 0 pro 5,0ˆ i <ϑ , tedy C = 0,5. 10.5. ROC křivka Pomocí klasifikační tabulky lze odhadnout senzitivitu a specificitu daného klasifikačního procesu (v našem případě logistického regresního modelu). Senzitivita … pravděpodobnost, že objekt, u něhož nastal úspěch, byl správně zařazen mezi úspěšné objekty. Odhad senzitivity: ca a TPF + = (true positive fraction – relativní četnost správně klasifikovaných úspěšných objektů). Specificita … pravděpodobnost, že objekt, u něhož nastal neúspěch, byl správně zařazen mezi neúspěšné objekty. Odhad specificity: db d TNF + = (true negative fraction – relativní četnost správně klasifikovaných neúspěšných objektů). Pro různé hodnoty dělicího bodu C dosáhneme různé hodnoty odhadů senzitivity a specificity. Graficky se jejich vztah reprezentuje právě pomocí ROC křivky. Na vodorovnou osu vynášíme hodnoty 1 – TNF (tj. 1 – odhad specificity) a na svislou osu hodnoty TPF (tj. odhady senzitivity). Pro ideální model má ROC křivka tvar lomené čáry procházející body [0;0], [0;1] a [1;1]. Pro náhodný model ROC křivka kopíruje úsečku spojující body [0;0] a [1;1]. ROC křivka pro reálný model by měla být pokud možno co nejblíže ke křivce pro ideální model. Ukázka ROC křivek pro ideální model (a), reálný model (b), náhodný model (c) (Obrázek je převzat z bakalářské práce Filipa Zlámala Logistická regrese v R) Velikost A plochy AUC pod ROC křivkou je nejběžnější kvantitativní index popisující ROC křivku. Predikční schopnost modelu hodnotíme pomocí tabulky: A hodnocení 0,9 - 1 výborně 0,8 – 0,9 velmi dobře 0,7 – 0,8 dobře 0,6 - 0,7 dostatečně 0,5 – 0,6 nedostatečně 11. Příklad V souboru 50 rodin byly zjišťovány tyto údaje: - zda v posledních dvou letech rodina navštívila jistou rekreační oblast (veličina ID, nabývá hodnoty 0 pro odpověď „ne“, hodnoty 1 pro odpověď „ano“) - roční příjem v tisících dolarů (veličina X1) - postoj k cestování (veličina X2, devítibodová škála, 1 = naprosto odmítavý, 9 = veskrze kladný) - význam přičítaný rodinné dovolené (veličina X3, devítibodová škála, 1 = nejnižší, 9 = nejvyšší) - počet členů rodiny (veličina X4) - věk nejstaršího člena rodiny (veličina X5). a) Proveďte jednorozměrnou analýzu dat. b) Sestavte model binární logistické regrese, který umožní odhadnout pravděpodobnost návštěvy dané rekreační oblasti v závislosti na vhodně vybraných veličinách X1, …, X5. Odhadněte regresní parametry a vypočtěte odhady podílů šancí Sestrojte 95% intervaly spolehlivosti pro regresní parametry a pro podíly šancí. Na hladině významnosti 0,05 proveďte dílčí testy významnosti. c) Na hladině významnosti 0,05 proveďte celkový test významnosti. d) Proveďte hodnocení kvality modelu Řešení: Ad a) Testy normality proměnných X1, …, X5 v daných dvou skupinách: Pro skupinu rodin, které danou rekreační oblast nenavštěvují: Statistiky – Základní statistiky/tabulky – Select cases – ID=0 – OK – Tabulky četností – Proměnné X1 až X5 – OK – Normalita – zaškrtneme S-W test – Testy normality Testy normality (dovolena.sta) Zhrnout podmínku: ID=0 Proměnná N W p X1: roční příjem v tisících dolarů X2: postoj k cestování (škála 9 bodů) X3: význam rodinné dovolené (škála 9 bodů) X4: počet členů rodiny X5: věk nejstaršího člena 29 0,940188 0,101411 29 0,964071 0,412187 29 0,964432 0,420319 29 0,917696 0,026668 29 0,944508 0,131598 Pro skupinu rodin, které danou rekreační oblast navštěvují: Statistiky – Základní statistiky/tabulky – Select cases – ID=1 – OK – Tabulky četností – Proměnné X1 až X5 – OK – Normalita – zaškrtneme S-W test – Testy normality Testy normality (dovolena.sta) Zhrnout podmínku: ID=1 Proměnná N W p X1: roční příjem v tisících dolarů X2: postoj k cestování (škála 9 bodů) X3: význam rodinné dovolené (škála 9 bodů) X4: počet členů rodiny X5: věk nejstaršího člena 21 0,935874 0,180430 21 0,930271 0,139382 21 0,934717 0,171087 21 0,928224 0,126815 21 0,967589 0,679311 Na hladině významnosti 0,05 zamítáme hypotézu o normalitě u veličiny X4 ve skupině rodin, které danou rekreační oblast nenavštěvují. Dvouvýběrové t-testy: Statistiky – Základní statistiky a tabulky – t-test, nezávislé, dle skupin – OK, Proměnné –Závislé proměnné X1 až X5, Grupovací proměnná ID – OK. Po kliknutí na tlačítko Výpočet dostaneme tabulku t-testy; grupováno: ID (dovolena.sta) Skup. 1: návštěva ne Skup. 2: návštěva ano Proměnná Průměr návštěva ne Průměr návštěva ano t sv p Poč.plat návštěva ne Poč.plat. návštěva ano Sm.odch. návštěva ne Sm.odch. návštěva ano F-poměr Rozptyly p Rozptyly X1 X2 X3 X4 X5 42,84483 59,76190 -7,40751 48 0,000000 29 21 7,013894 9,142783 1,699176 0,193069 4,24138 5,14286 -1,89805 48 0,063712 29 21 1,661651 1,651839 1,011916 0,995884 4,27586 5,76190 -3,15623 48 0,002760 29 21 1,623412 1,670472 1,058816 0,872933 3,72414 4,33333 -1,73042 48 0,089980 29 21 1,130630 1,354006 1,434168 0,372786 46,93103 53,61905 -3,01289 48 0,004122 29 21 7,568407 7,990471 1,114643 0,776989 Na hladině významnosti 0,05 zamítáme hypotézu o shodě středních hodnot pouze pro proměnné X1, X3 a X5. Do modelu logistické regrese však zahrneme všechny proměnné, protože i u proměnných X2 a X4 jsou odpovídající p-hodnoty menší než 0,25. Ad b) Provedení logistické regrese Statistiky – Pokročilé lineární/nelineární modely – Zobecněné lineární/nelineární modely – Logitový model – OK – Proměnné – Závislé ID, Spojité nezáv. prom. X1, …, X5 – OK – na záložce Detaily zaškrtneme Parametrizace Ref. – OK. V okně výsledků zvolíme tlačítko Odhady a ve stromové struktuře v levé části Pracovního sešitu vybereme tabulku ID – Odhady parametrů. ID - Odhady parametrů (dovolena.sta) Rozdělení : BINOMICKÉ, Linkující funkce: LOGIT Modelovaná pravděpodobnost, že ID = návštěva ano Efekt Úroveň Efekt Sloupec Odhad Standard chyba Wald. Stat. Dolní LS 95,0% Horní LS 95,0% p Abs.člen X1 X2 X3 X4 X5 Měřítko 1 -34,7469 12,50700 7,718363 -59,2602 -10,2336 0,005466 2 0,3203 0,11219 8,148696 0,1004 0,5401 0,004309 3 0,8072 0,44522 3,287311 -0,0654 1,6799 0,069817 4 0,7156 0,38627 3,432061 -0,0415 1,4727 0,063942 5 -0,0885 0,51445 0,029594 -1,0968 0,9198 0,863416 6 0,2404 0,10208 5,545822 0,0403 0,4405 0,018525 1,0000 0,00000 1,0000 1,0000 Ve sloupci Odhad vidíme odhady regresních parametrů, dále směrodatné chyby těchto odhadů, hodnoty Waldových statistik pro test nevýznamnosti regresních parametrů, meze 95% intervalů spolehlivosti pro regresní parametry a p-hodnoty pro test nevýznamnosti regresních parametrů. (Řádek Měřítko nehraje roli.) Pouze proměnná X4 je vysoce nevýznamná, ostatní proměnné v modelu ponecháme. Nová tabulka: ID - Odhady parametrů (dovolena.sta) Rozdělení : BINOMICKÉ, Linkující funkce: LOGIT Modelovaná pravděpodobnost, že ID = návštěva ano Efekt Úroveň Efekt Sloupec Odhad Standard chyba Wald. Stat. Dolní LS 95,0% Horní LS 95,0% p Abs.člen X1 X2 X3 X5 Měřítko 1 -35,0954 12,43292 7,968074 -59,4634 -10,7273 0,004761 2 0,3185 0,11222 8,057523 0,0986 0,5385 0,004532 3 0,8193 0,44140 3,445366 -0,0458 1,6844 0,063429 4 0,7208 0,38601 3,486900 -0,0358 1,4774 0,061856 5 0,2405 0,10140 5,624647 0,0417 0,4392 0,017710 1,0000 0,00000 1,0000 1,0000 Pravděpodobnost, že rodina patří do skupiny rodin, které navštěvují danou rekreační oblast, je vyjádřena rovnicí ( ) 5321 x2405,0x7208,0x8193,0x3185,00954,3555332211 e1 1 xXxXxXxX/1IDP −−⋅−⋅− + ==∧=∧=∧== Např. : Parametr 0,3185 znamená, že při vzrůstu proměnné X1 o jednotku (tj. průměrný roční příjem se zvýší o 1000 dolarů) vzroste šance na zařazení do skupiny rodin navštěvujících danou rekreační oblast v průměru na 375,1e 3185,0 = krát. Model zařadí rodinu do skupiny s ID = 1, pokud odhadnutá pravděpodobnost je aspoň 0,5. Dále se podíváme na tabulku ID – Poměry šancí. ID - Poměry šancí (dovolena.sta) Rozdělení : BINOMICKÉ, Linkující funkce: LOGIT Modelovaná pravděpodobnost, že ID = návštěva ano Efekt Úroveň Efekt Sloupec Šance Poměr Dolní LS 95,0% Horní LS 95,0% p Abs.člen X1 X2 X3 X5 Měřítko 1 2 1,375118 1,155175 1,595062 0,004532 3 2,268945 1,403815 3,134075 0,063429 4 2,056078 1,299518 2,812638 0,061856 5 1,271861 1,073123 1,470600 0,017710 1,000000 Zde jsou uvedeny odhady podílů šancí na návštěvu dané rekreační oblasti, 95% intervaly spolehlivosti pro teoretické podíly šancí a p-hodnoty pro testy hypotéz, že teoretické podíly šancí je 1, tj. pro test hypotézy, že jednotlivé regresory neovlivňuje návštěvu oblasti. Povšimněte si, že p-hodnoty pro testy těchto hypotéz jsou totožné s p-hodnotami pro testy nevýznamnosti regresních parametrů - viz tabulka ID – Odhady parametrů. Ad c) Celkový test významnosti provedeme takto: Návrat do GLZ-Výsledky – Kval. proložení – vybereme tabulku Testování globální nulové hypotézy. Testování glonální nulové hypotézy: BETA=0 (dovolena.sta) Rozdělení : BINOMICKÉ, Linkující funkce: LOGIT Modelovaná pravděpodobnost, žeID = návštěva ano (Vzorek pro analýzu) Chí-kvadrát SV p Poměr věrohodnos Skóre Wald. 47,114384 4 0,000000 30,871760 4 0,000003 8,738087 4 0,067990 Zajímá nás test poměrem věrohodnosti. Protože jeho p-hodnota je velmi blízká nule, hypotézu o nevýznamnosti modelu zamítáme na hladině významnosti 0,05. Ad d) Hodnocení kvality modelu provedeme několika způsoby. Ve výstupní tabulce GLZ – Výsledky vybereme Kvalita proložení a zvolíme tabulku ID – Statistiky kvality modelu. ID - Statistiky kvality modelu (dovolena.sta) Rozdělení : BINOMICKÉ, Linkující funkce: LOGIT Modelovaná pravděpodobnost, žeID = návštěva ano (Vzorek pro analýzu) SV Stat. Stat/sv Odchylka Deviance v měřít Pearsonovo Chi2 Scaled P. Chi2 AIC BIC Cox-Snell R2 Nagelkerke R2 Log-věrohodnost 45 20,914816 0,464774 45 20,914816 0,464774 45 19,602100 0,435602 45 19,602100 0,435602 30,914816 40,474931 0,610265 0,820812 -10,457408 Pearsonův 2 χ test dobré shody (nulová hypotéza tvrdí, že naměřené a predikované hodnoty se neliší): na řádku Pearsonovo Chi2 je uvedena hodnota testové statistiky, v našem případě 19,6. Kritický obor ( ) ) )∞=∞χ= ;6562,61,45W 95,0 2 , tedy nulovou hypotézu nezamítáme na hladině významnosti 0,05. Nagelkerkův koeficient nabývá hodnoty 0,8208, což je hodnota dosti blízká 1 a svědčí o dobré kvalitě našeho modelu. Výsledky Hosmerova-Lemeshowova testu jsou uvedeny v tabulce ID – Kvalita proložení. ID - Kvalita proložení: Hosmer-Lemeshow Test (dovolena.sta) Rozdělení : BINOMICKÉ, Linkující funkce: LOGIT Hosmer Lemeshow = 3,8441, p hodn. = 0,870904 Odezva Skupi1a Skupi2a Skupi3a Skupi4a Skupi5a Skupi6a Skupi7a Skupi8a Skupi9a Skupi10 Row Tot. 0: Pozorov. Očekáv. 1: Pozorov. Očekáv. Vš. skup. 5,00 5,00 5,00 5,00 4,00 2,00 3,00 0,00 0,00 0,00 29,0 4,99 4,97 4,92 4,80 4,35 3,15 1,57 0,24 0,01 0,00 0,00 0,00 0,00 0,00 1,00 3,00 2,00 5,00 5,00 5,00 21,0 0,01 0,03 0,08 0,20 0,65 1,85 3,43 4,76 4,99 5,00 5,00 5,00 5,00 5,00 5,00 5,00 5,00 5,00 5,00 5,00 50,0 Hosmerův _ Lemeshowův test poskytl p-hodnotu 0,8709, tedy na hladině významnosti 0,05 nelze zamítnout hypotézu, že pozorované a modelované hodnoty se neliší. Klasifikační tabulku najdeme v systému STATISTICA takto: Ve výstupní tabulce GLZ – Výsledky vybereme záložku Rezid. 1 – Klasif & odds ratio. Klasifikace případů (dovolena.sta) Odds ratio: 52,000000 Log odds ratio: 3,951244 Předpovězená: návštěva ano Předpovězená: návštěva ne Procento správných Pozorované: návštěva ano Pozorované: návštěva ne 18 3 85,7142857 3 26 89,6551724 Z 21 rodin, které sledovanou rekreační oblast navštívily, model správně klasifikoval 18 rodin, tj. 85,7 %. Z 29 rodin, které sledovanou rekreační oblast nenavštívily, model správně klasifikoval 26 rodin, tj. 89,7 %. Celková úspěšnost správné klasifikace je tedy %88 50 2618 = + . ROC křivku získáme takto: Ve výstupní tabulce GLZ – Výsledky vybereme záložku Rezid. 1 – ROC křivka ROC křivka Oblast: 0.97 -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 1-Specificita -0,2 0,0 0,2 0,4 0,6 0,8 1,0 1,2 Citlivost Vidíme, že ROC křivka se blíží ideálnímu tvaru a plocha pod ní je AUC = 0,97, tedy predikční schopnost modelu je výborná. Shrnutí: Vytvořený model logistické regrese je statisticky významný na hladině významnosti 0,05 a není v rozporu s danými daty. Nagelkerkova koeficientu determinace nabývá hodnoty 0,82. Úspěšnost správné klasifikace je 88 %. Plocha AUC pod ROC křivkou je 0,97. Lze soudit, že pravděpodobnost návštěvy dané rekreační oblasti lze uspokojivě vysvětlit působením čtyř sledovaných proměnných (roční příjem v tisících dolarů, postoj k cestování, význam přičítaný rodinné dovolené a věk nejstaršího člena rodiny).