--- title: "Ekonometrie v R" author: "Štěpán Mikula" documentclass: article output: html_document: theme: cerulean toc: yes toc_float: yes pdf_document: default fontsize: 10pt classoption: a4paper --- # Ekonometrie v R {#kap:ekonometrie} Skupina balíků sdružená v `tidyverse` představuje kompaktní sadu nástrojů pro práci s daty. Sdílejí datové struktury (tabulky v tidy formátu) i logiku ovládání. Bohužel podobná unifikace chybí u nástrojů, které výzkumník potřebuje pro ekonometrickou (statistickou) analýzu dat, která následuje po zpracování dat. Ekonometrické nástroje jsou v R rozptýleny v mnoha balících různé kvality, které nejsou žádným způsobem koordinovány nebo unifikovány. To velmi znepříjemňuje ekonometrickou analýzu dat v R. I přesto R poskytuje všechny základní a mnoho pokročilých nástrojů a v zásadě pokrývá všechny běžné potřeby výzkumníků. Účelem této kapitoly není vysvětlovat teorii ekonometrické analýzy nebo diskutovat vhodnost nasazení jednotlivých metod nebo designu regresní analýzy (identifikační strategie), ale seznámit čtenáře se základní logikou fungování ekonometrických nástrojů v R. Právě základní logika fungování je totiž něco, co drtivá většina ekonometrických nástrojů v R sdílí. Tuto základní logiku lze popsat pomocí obrázku XXX: Na počátku ekonometrické analýzy má výzkumník data a představu o tom, jak svět funguje (o data generujícím procesu). Pro naprostou většinu ekonometrických nástrojů (a všechny základní) je vhodné připravit si data do podoby tabulky v tidy formátu. Představa o fungování světa (tj. o vztazích mezi proměnnými v datové tabulce) se formalizuje do podoby modelu. Jeho podoba se předává R v podobě speciální datové třídy `formula`. Data a specifikovaný model jsou nezbytnými vstupy pro estimátor. Estimátor je typicky funkce, která provede odhad parametrů modelu a vytvoří výstupní objekt. Tento objekt obsahuje typicky pouze samotný odhad parametrů modelu, tabulku dat použitých pro odhad, rezidua, vyrovnané hodnoty a další doplňující údaje. Typicky neobsahuje žádné testy, robustní směrodatné chyby a podobně. Analýza a zobrazení odhadu (výsledků) se v R typicky provádí pomocí specializovaných funkcí, pro které je výstup estimační funkce vstupem. V této kapitole se naučíte: - Specifikovat model ve třídě `formula` - Použít základní estimační funkce - Provést základní diagnostické testy odhadu - Vypočítat robustní směrodatné chyby - Provádět hromadně odhady více modelových specifikací na stejném vzorku dat - Provádět hromadně odhady jedné modelové specifikace na různých vzorcích dat - Exportovat výsledky do přehledných tabulek ## R a ostatní R představuje spolu s Pythonem špičkový software ohledně práce s daty. Tomu odpovídá i to, že právě tyto dva jazyky jsou mezi datovými analytiky pravděpodobně nejvíce používané. Pozice R ohledně ekonometrické analýzy však zdaleka není tak jasná. V ekonomii se rozhodně nejedná o nejpoužívanější nástroj. Tím je stále bezesporu Stata, která se logikou svého fungování velmi podobá Gretlu. (Gretl je navržen podle Stata.) Stata pracuje nad jednou tabulkou dat, která je nahraná v paměti. V rámci každého odhadu specifikujete například i typ robustních chyb, které chcete použít. Své rozhodnutí však později nemůžete změnit. Výsledkem takového návrhu je velmi obtížná práce s daty ve Stata (např. slučování dvou tabulek). Nutnost specifikovat všechny parametry již v rámci odhadu modelu může být také velmi omezující v případě odhadů, které jsou časově náročné. Stata také dává oproti R výzkumníkovi často menší svobodu v detailním nastavení jednotlivých procedur. V neposlední řadě je Stata také velmi nákladná (levné verze jsou omezené tak, že nejsou v reálném boji použitelné). Omezené možnosti jsou ovšem lehce vyvažovány jednoduchostí ovládání. Hlavní výhodou Stata je však implementace obrovského množství nejrůznějších estimátorů, které mají velmi konzistentní rozhraní. ## Specifikace modelu pomocí `formula` Fungování ekonometrických nástrojů můžeme ilustrovat na příkladu datasetu s výškou, váhou, věkem, sportem, a pohlavím sportovců, kteří se zúčastnili olympijských her v Londýně v roce 2012 (dostupný v `VGAMdata::oly12`): ```r VGAMdata::oly12 %>% mutate( Sport = str_extract(Sport, "[:alnum:]*") ) %>% select( Height, Weight, Sex, Age, Sport ) %>% as_tibble %>% drop_na() %>% mutate_if(is.factor,as.character) -> oly12 print(oly12, n=5) ``` ``` ## # A tibble: 9,038 x 5 ## Height Weight Sex Age Sport ## ## 1 1.7 60 M 23 Judo ## 2 1.93 125 M 33 Athletics ## 3 1.87 76 M 30 Athletics ## 4 1.78 85 F 26 Athletics ## 5 1.82 80 M 27 Handball ## # … with 9,033 more rows ``` Nejprve nás bude zajímat vztah mezi váhou a výškou sportovců. Ten si můžeme vykreslit pomocí jednoduchého bodového grafu: Zdá se, že existuje lineární vztah mezi váhou a výškou sportovců, který můžeme popsat rovnicí: $$Weight = \alpha + \beta Height + \varepsilon$$ kde $\alpha$ a $\beta$ jsou neznámé parametry, které chceme odhadnout a $\varepsilon$ je náhodná složka. Tuto rovnici si můžeme představit jako přímku, která je proložena mrakem pozorování: ``` ## `geom_smooth()` using formula 'y ~ x' ``` Rovnici však musíme zapsat tak, aby byla srozumitelná i pro R. K tomu slouží objekty třídy `formula`. Jejich formulace následuje syntax: ```r LHS ~ RHS ``` Tedy výrazy na levé straně (LHS) $=$ (zastoupené znakem `~`) výrazy na pravé straně (RHS). Rovnici vysvětlující váhu bychom tedy mohli popsat jako: ```r model <- Weight ~ Height class(model) ``` ``` ## [1] "formula" ``` V ekonomii je závislá (vysvětlovaná) proměnná (tedy proměnná na levé straně rovnice) funkcí více než jedné nezávislé (vysvětlující) proměnné. Pokud bychom chtěli vysvětlovat váhu nejen výškou, ale i věkem, potom bychom rozšířili specifikaci následujícím způsobem: ```r model <- Weight ~ Height + Age ``` Dodatečná proměnné jsou přidávány pomocí `+`. Za povšimnutí stojí, že při použití takovéto specifikace by byl model odhadnut s úrovňovou konstantou $\alpha$ -- i když není explicitně v rovnici přítomna. Explicitně se naopak musí zadat její nepřítomnost připojením `-1` nebo `0`: ```r model <- Weight ~ Height - 1 model <- Weight ~ Height + 0 ``` Jména proměnných v rovnici musí odpovídat jménům sloupců v datové tabulce. Zásadní výhodou R je, že do rovnice je možné zadat i transformace dat: ```r model <- Weight ~ log(Height) model <- Weight ~ Height >= 1.8 ``` První varianta například odpovídá specifikaci: $$Weight = \alpha + \beta \log(Height) + \varepsilon$$ Ve druhé je váha vysvětlována dummy (umělou proměnnou). Výraz `Height >= 1.8` vytvoří logickou proměnnou, která je při volání estimační funkce transformována na vektor s hodnotami 1 a 0. **Není tedy potřeba vytvářet nové transformované sloupce v datové tabulce.** To je velký rozdíl mezi R a Statou, Gretlem,... (Jakkoliv je v některých případech vytvoření transformovaných proměnných do tabulky doporučeníhodné.) V rovnici může být přímo přítomna celá řada funkcí -- s výjimkou těch, jejichž interpretace by nebyla jasná. Pokud bychom například chtěli odhadnout modelovou specifikaci $$Weight = \alpha + \beta (Height + Age) + \varepsilon$$ potom je nutná sdělit R, že má nejdříve sečíst `Height` a `Age` a pro vysvětlení `Weight` použít až výsledný součet. Pro tyto účely se používá funkce `I()`. Ta sděluje, že se mají nejprve provést operace definované uvnitř funkce a pro odhad použít až výsledek. Výše uvedený příklad by se tak R tlumočil jako ```r model <- Weight ~ I(Height + Age) ``` Třída formula umožňuje jednoduše jednoduše definovat i interakční členy pomocí operátorů dvojtečka (`:`) a hvězdička (`*`): ```r model <- Weight ~ Height:Sex ``` Chápe R jako $$Weight = \alpha + \beta Height\times Sex_{male} + \gamma Height\times Sex_{female} + \varepsilon$$ Tato rovnice vyjadřuje možný odlišný mezní efekt výšky na váhu u mužů a žen. Výsledné proložení mrakem pozorování by se tak v bodovém grafu výše vlastně rozpadlo na dvě přímky s odlišným sklonem (parametry $\beta$ a $\gamma$) a shodnou úrovňovou konstantou $\alpha$. Další možností jak popsat interakci proměnných je `*`: ```r model <- Weight ~ Height*Sex ``` Tato `formula` odpovídá jiné rovnici: $$Weight = \alpha + \beta Height + \epsilon Sex_{male} + \gamma Height\times Sex_{male} + \varepsilon$$ Tato specifikace umožňuje odlišnost regresní křivky pro muže a ženy jak ve sklonu, tak v úrovňové konstantě. Zvláštní význam má i symbol stříšky (`^`). Rovnici $$Weight = \alpha + \beta_1 Height + \beta_2 Age + \beta_3 Sex_{male} + \gamma_1 Height\times Sex_{male} + \gamma_2 Height\times Age + \gamma_3 Sex_{male}\times Age + \varepsilon$$ můžeme zapsat jako: ```r model <- Weight ~ (Height + Age + Sex)^2 ``` Výraz se přeloží tak, že výsledná rovnice obsahuje jednotlivé proměnné a také jejich interakce do, v tomto případě, druhého řádu. Pro třídu formula jsou implementovány i některé zajímavé metody. Velmi užitečná je například funkce `update()` s velmi prostou syntaxí: ```r update(old,new,...) ``` Kde `old` je původní rovnice a `new` pravidla, která stanovují, jak se původní rovnice má upravit. ```r model <- Weight ~ Height + Age print(model) ``` ``` ## Weight ~ Height + Age ``` Takto můžeme zachovat všechny proměnné -- reprezentované pomocí `.` a jen některou přidat nebo odebrat: ```r model %>% update(. ~ . -Age) ``` ``` ## Weight ~ Height ``` Nebo je možné nahradit celou stranu rovnice: ```r model %>% update(. ~ Age) ``` ``` ## Weight ~ Age ``` ## Odhad modelu Při volání estimační funkce se dohromady spojí všechny tři komponenty: modelová specifikace, data a estimátor. Základní estimační funkcí je `lm()`, určená pro odhad lineárních modelů. Funkce přijímá celou řadu parametrů: ```r lm(formula, data, subset, weights, na.action, method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, ...) ``` Mezi základní patří: - `formula` je objekt třídy "formula", který obsahuje symbolický popis rovnice, jejíž parametry se mají odhadnout. - `data` je tabulka (`data.frame`) která obsahuje data, která se mají použít. Jména sloupců v tabulce musí odpovídat proměnným použitým v parametru `formula`. Datová tabulka nemusí být specifikována. Pokud je tento parametr prázdný, potom R interpretuje jména v modelu jako jména proměnných z prostředí, ze kterého je funkce `lm()` volána. - `subset` je nepovinný parametr, který určuje, které pozorování (řádky) se mají použít pro odhad - `weights` je vektor vah. Pokud není zadán je pro odhad použit obyčejný OLS estimátor, ve kterém mají všechna pozorování stejnou váhu. V případě zadání vah je vráceno WLS. Struktura a forma těchto základních vstupů je typicky shodná právě s `lm()`. ### Vytvoření matice plánu Pro odhad použijeme model v následující specifikaci: ```r model <- Weight ~ Height + Age + I(Age^2) + Sex ``` Proměnné `Weight`, `Height` a `Age` jsou numerické spojité proměnné. Proměnná `Sex` je zcela jiný případ. Jsou to kategoriální proměnná, se dvěma úrovněmi (`M` -- muž, `F` -- žena). Taková proměnná nemůže vstoupit do regrese přímo. Je nutné ji transformovat na matici nul a jedniček. Transformovaná matice bude mít vždy o jeden sloupec méně, než kolik má transformovaná proměnná úrovní. Estimační funkce si všechny vstupní faktory transformuje sama pomocí interního volání funkce `model.matrix()`. V rámci volání tento funkce se provedou i požadované transformace dat. V případě našeho modelu vytvoření druhé mocniny věku. Výsledkem volání `model.matrix()` je tedy matice plánu: ```r model.matrix(model, data = oly12) %>% head() ``` ``` ## (Intercept) Height Age I(Age^2) SexM ## 1 1 1.70 23 529 1 ## 2 1 1.93 33 1089 1 ## 3 1 1.87 30 900 1 ## 4 1 1.78 26 676 0 ## 5 1 1.82 27 729 1 ## 6 1 1.82 30 900 0 ``` Funkce `model.matrix()` provede požadované transformace numerických proměnných (přidán sloupec `I(Age^2)`) a pokud najde jinou než numerickou proměnnou (`factor`) provede jejich transformaci do matice binárních proměnných. Při této transformaci postupuje následovně. Prohlédne si faktor a určí referenční úroveň (*level*). Ta nebude ve výsledné matici obsažena, pokud by byla, model by nešel odhadnout. Pro každou zbývající úroveň vytvoří vektor, který nabývá hodnoty 1 pro pozorování s danou úrovní faktoru. V příkladu výše si `model.frame()` vybral jako referenční úroveň hodnotu `F` a vytvořil vektor `SexM`, který nabývá hodnoty 1 pro muže a 0 ve všech ostatních případech. Pokud by byla vstupní kategoriální proměnná třídy `character`, potom by ji `model.frame()`nejprve transformoval na `factor` a pokračoval tak, jak je popsáno výše. Toto chování může být motivem pro to, aby uživatel sám provedl konverzi na faktor. Třída faktor totiž umožňuje nastavit, která úroveň má být použita jako referenční. Toho můžeme docílit pomocí několika cest. První variantou je konverze do `factor` s pomocí `base::factor()`. Tato funkce umožňuje uživateli zadat vektor možných úrovní. První z nich se potom nastaví jako referenční: ```r factor(Sex$oly12, c("M","F")) ``` Tento postup lze aplikovat i na proměnné, které již jsou `factor`. Nevýhodou je nutnost zadávat všechny úrovně. Podobnou variantou je použití funkce `forcats::as_factor()`. Ta má vlastnost, že pořadí úrovní stanoví podle podle pořadí jejich výskytu v transformovaném vektoru. Problém této funkce je, že pracuje zvláštně s chybějícími hodnotami. (To se pravděpodobně bude během vývoje balíku ještě měnit.) Další principiálně odlišnou variantou je nastavení referenční úrovně pomocí funkce `stats::relevel()` (alternativně `forcats::fct_relevel()`). Ta umožňuje prosté nastavení referenční úrovně: ```r factor(Sex$oly12) %>% relevel("M") ``` Po změně referenční úrovně bude výsledek následující: ```r oly12 %>% mutate( Sex = factor(Sex) %>% relevel("M") ) %>% model.matrix(model,.) %>% head() ``` ``` ## (Intercept) Height Age I(Age^2) SexF ## 1 1 1.70 23 529 0 ## 2 1 1.93 33 1089 0 ## 3 1 1.87 30 900 0 ## 4 1 1.78 26 676 1 ## 5 1 1.82 27 729 0 ## 6 1 1.82 30 900 1 ``` Zejména v rámci průzkumu dat je zajímavá varianta provést konverzi na faktor v rámci modelu (`formula`) -- všimněte si, že může být aplikována i na numerickou proměnnou: ```r Weight ~ Height + factor(Age) + Sex ``` Konverze numerických hodnot do faktoru je poměrně častá zvláště u dotazníkových dat, kde čísla mohou pouze kódovat hodnoty kategoriální proměnné. Dalším případem je například vytvoření fixních efektů pro roky. Funkce `model.matrix()`, kterou si estimační funkce volá interně tedy vytvoří matici plánu, která je následně použité při samotném odhadu parametrů modelu. Odhadová procedura si liší podle použité estimační funkce. V případě `lm()` je to QR dekompozice (https://en.wikipedia.org/wiki/QR_decomposition), při použití `glm()` jsou to například různé numerické optimalizační algoritmy (defaultně IWLS:https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares). ### Výstupy estimační funkce ```r lm(model, data = oly12) -> est_model class(est_model) ``` ``` ## [1] "lm" ``` Funkce `lm()` vrací S3 objekt třídy `lm`. Jiné estimační funkce vrací jiné objekty, které jsou ovšem často potomky `lm`. Objekt třídy `lm` má následující strukturu (viz `?lm`): ```r str(est_model, max.level = 1) ``` ``` ## List of 13 ## $ coefficients : Named num [1:5] -1.06e+02 9.46e+01 4.10e-01 -4.09e-03 5.59 ## ..- attr(*, "names")= chr [1:5] "(Intercept)" "Height" "Age" "I(Age^2)" ... ## $ residuals : Named num [1:9038] -7.71 33.72 -9.15 14.68 0.11 ... ## ..- attr(*, "names")= chr [1:9038] "1" "2" "3" "4" ... ## $ effects : Named num [1:9038] -6926.6 -1161.4 -96.3 -15.7 -221.1 ... ## ..- attr(*, "names")= chr [1:9038] "(Intercept)" "Height" "Age" "I(Age^2)" ... ## $ rank : int 5 ## $ fitted.values: Named num [1:9038] 67.7 91.3 85.2 70.3 79.9 ... ## ..- attr(*, "names")= chr [1:9038] "1" "2" "3" "4" ... ## $ assign : int [1:5] 0 1 2 3 4 ## $ qr :List of 5 ## ..- attr(*, "class")= chr "qr" ## $ df.residual : int 9033 ## $ contrasts :List of 1 ## $ xlevels :List of 1 ## $ call : language lm(formula = model, data = oly12) ## $ terms :Classes 'terms', 'formula' language Weight ~ Height + Age + I(Age^2) + Sex ## .. ..- attr(*, "variables")= language list(Weight, Height, Age, I(Age^2), Sex) ## .. ..- attr(*, "factors")= int [1:5, 1:4] 0 1 0 0 0 0 0 1 0 0 ... ## .. .. ..- attr(*, "dimnames")=List of 2 ## .. ..- attr(*, "term.labels")= chr [1:4] "Height" "Age" "I(Age^2)" "Sex" ## .. ..- attr(*, "order")= int [1:4] 1 1 1 1 ## .. ..- attr(*, "intercept")= int 1 ## .. ..- attr(*, "response")= int 1 ## .. ..- attr(*, ".Environment")= ## .. ..- attr(*, "predvars")= language list(Weight, Height, Age, I(Age^2), Sex) ## .. ..- attr(*, "dataClasses")= Named chr [1:5] "numeric" "numeric" "numeric" "numeric" ... ## .. .. ..- attr(*, "names")= chr [1:5] "Weight" "Height" "Age" "I(Age^2)" ... ## $ model :'data.frame': 9038 obs. of 5 variables: ## ..- attr(*, "terms")=Classes 'terms', 'formula' language Weight ~ Height + Age + I(Age^2) + Sex ## .. .. ..- attr(*, "variables")= language list(Weight, Height, Age, I(Age^2), Sex) ## .. .. ..- attr(*, "factors")= int [1:5, 1:4] 0 1 0 0 0 0 0 1 0 0 ... ## .. .. .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..- attr(*, "term.labels")= chr [1:4] "Height" "Age" "I(Age^2)" "Sex" ## .. .. ..- attr(*, "order")= int [1:4] 1 1 1 1 ## .. .. ..- attr(*, "intercept")= int 1 ## .. .. ..- attr(*, "response")= int 1 ## .. .. ..- attr(*, ".Environment")= ## .. .. ..- attr(*, "predvars")= language list(Weight, Height, Age, I(Age^2), Sex) ## .. .. ..- attr(*, "dataClasses")= Named chr [1:5] "numeric" "numeric" "numeric" "numeric" ... ## .. .. .. ..- attr(*, "names")= chr [1:5] "Weight" "Height" "Age" "I(Age^2)" ... ## - attr(*, "class")= chr "lm" ``` Důležité jsou zejména následující položky: - `coefficients`...pojmenovaný vektor odhadnutých koeficientů ```r est_model$coefficients ``` ``` ## (Intercept) Height Age I(Age^2) SexM ## -1.060127e+02 9.463456e+01 4.096565e-01 -4.094392e-03 5.591553e+00 ``` - `residuals`...vektor reziduí -- tj. odhadů náhodné složky $\varepsilon$ ```r est_model$residuals %>% head ``` ``` ## 1 2 3 4 5 6 ## -7.7137568 33.7165890 -9.1502080 14.6799373 0.1103484 -1.8269273 ``` - `fitted.values`...vyrovnané hodnoty ```r est_model$fitted.values %>% head ``` ``` ## 1 2 3 4 5 6 ## 67.71376 91.28341 85.15021 70.32006 79.88965 74.82693 ``` - `model`...modifikovanou tabulku dat, která byla skutečně použita při odhadu (tj, například s vypuštěnými nepoužitými proměnnými nebo neúplnými řádky). ```r est_model$fitted.values %>% head ``` ``` ## 1 2 3 4 5 6 ## 67.71376 91.28341 85.15021 70.32006 79.88965 74.82693 ``` Za povšimnutí stojí, že objekt obsahuje skutečně pouze odhad -- data, odhadnuté koeficienty, rezidua, vyrovnané hodnoty a dodatečné informace. Neobsahuje žádné dodatečné výpočty jako například diagnostické testy. Pro objekty třídy `lm` existuje metoda generické funkce `print()`. Její výstupy však nejsou nijak zvlášť informativní: ```r print(est_model) ``` ``` ## ## Call: ## lm(formula = model, data = oly12) ## ## Coefficients: ## (Intercept) Height Age I(Age^2) SexM ## -1.060e+02 9.463e+01 4.097e-01 -4.094e-03 5.592e+00 ``` Obsahuje jen použité volání a odhadnuté hodnoty koeficientů. Více se dozvíme pomocí funkce `summary()`: ```r summary(est_model) ``` ``` ## ## Call: ## lm(formula = model, data = oly12) ## ## Residuals: ## Min 1Q Median 3Q Max ## -56.334 -5.717 -1.336 3.756 135.665 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -1.060e+02 2.385e+00 -44.444 < 2e-16 *** ## Height 9.463e+01 1.141e+00 82.918 < 2e-16 *** ## Age 4.097e-01 1.096e-01 3.738 0.000187 *** ## I(Age^2) -4.094e-03 1.819e-03 -2.250 0.024442 * ## SexM 5.592e+00 2.560e-01 21.840 < 2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 10.13 on 9033 degrees of freedom ## Multiple R-squared: 0.6031, Adjusted R-squared: 0.6029 ## F-statistic: 3431 on 4 and 9033 DF, p-value: < 2.2e-16 ``` V tomto případě vidíme kromě odhadů koeficientů i směrodatné chyby a výsledky testů statistické významnosti koeficientů (v tomto případě t-testů). Výstup je také doplněn několika dalšími údaji jako je výsledek F-testu a (adjustovaný) koeficient determinace. Funkce `summary()` poskytuje rychlý pohled na výsledky, nicméně ten může být zavádějící. Funkce totiž vrací pouze standardní výstupy. Nijak například nereflektuje možné zkreslení odhadu směrodatných odchylek kvůli heteroskedasticitě a podobně. Všimněte si, že `summary()` používá méně obvyklé -- a přísnější -- hvězdičkové značení! ### Práce s výsledkem odhadů V R najdete řadu specializovaných funkcí, které dokážou z výsledného modelu získávat, nebo na jeho základě dopočítávat užitečné informace. Tabulka obsahuje jejich základní přehled: |Funkce |Popis | |------------------|----------------------------------------------------------| |`coefficients()`, `coefs()` | Vrací odhadnuté koeficienty -- obsah `.$coefficients`| |`residuals()`, `resid()` |Vrací vektor reziduí -- pro objekty `lm` je výstup identický s `.$residuals`| |`fitted.values()`|Vrací očekávané hodnoty -- obsah `.$fitted.values`| |`predict()`|Umožňuje dopočítat predikce pro jiná data, než na kterých byl model odhadnut| |`plot()`|Vrací diagnostické grafy vytvořené pomocí základní grafiky (ne v `ggplot2`)| |`confint()`|Vrací konfidenční intervaly pro odhady koeficientů| |`deviance()`|Vrací MSE -- průměr čtverců reziduí| |`vcov()`|Vrací odhad kovarianční matice| |`logLik()`|Vrací log-likelihood| |`AIC()`|Vrací Akaikovo informační krítérium (AIC)| |`BIC()`|Vrací Schwarzovo Bayesian kritérium (BIC)| |`nobs()`|Vrací počet pozorování| Pro některé účely je vhodné seznámit se i se sofistikovanějšími funkcemi. Pro účely diagnostiky nebo okometrického posouzení kvality vyrovnání může být užitečné výsledky odhadu vizualizovat. Výsledky jsou však schované ve zvláštním objektu a ne v tabulce, kterou potřebujeme pro použití s `ggplot2`. Pro transformaci různých objektů do tabulky existuje v `ggplot2` funkce `fortify()`. Stejnou funkcionalitu poskytuje i balík **broom**. Pro objekty třídy `lm` je ekvivalentem funkce `broom::augment()`: ```r library(broom) augment(est_model) %>% as_tibble %>% print(n = 5) ``` ``` ## # A tibble: 9,038 x 12 ## Weight Height Age I.Age.2. Sex .fitted .se.fit .resid .hat .sigma ## > ## 1 60 1.7 23 529 M 67.7 0.216 -7.71 4.53e-4 10.1 ## 2 125 1.93 33 1089 M 91.3 0.221 33.7 4.74e-4 10.1 ## 3 76 1.87 30 900 M 85.2 0.175 -9.15 2.99e-4 10.1 ## 4 85 1.78 26 676 F 70.3 0.186 14.7 3.36e-4 10.1 ## 5 80 1.82 27 729 M 79.9 0.155 0.110 2.35e-4 10.1 ## # … with 9,033 more rows, and 2 more variables: .cooksd , .std.resid ``` Pozor, **broom** navrací `tibble`. Proto je užitečné výstup do `tibble` konvertovat. Výsledek je již jednoduše použitelný pro vizualizaci pomocí `ggplot2`. Následující obrázek ukazuje korelaci pozorovaných (*observed*) a vyrovnaných (*fitted*) hodnot. Do obrázku je přidána pomocí `geom_abline()` i osa kvadrantu na která se pozorované a vyrovnané hodnoty rovnají. ```r est_model %>% augment() %>% ggplot( aes(x = Weight, y = .fitted) ) + geom_point(alpha = 0.05) + geom_abline(slope = 1, intercept = 0, color = "red") + scale_x_continuous(name = "Observed values") + scale_y_continuous(name = "Fitted values") ``` Balík **broom** obsahuje i další užitečné funkce: `tidy()` a `glance()`. Funkce `tidy()` vrací tabulku s odhady parametrů, směrodatných odchylek, testových statistik a p-hodnot: ```r tidy(est_model) ``` ``` ## # A tibble: 5 x 5 ## term estimate std.error statistic p.value ## ## 1 (Intercept) -106. 2.39 -44.4 0. ## 2 Height 94.6 1.14 82.9 0. ## 3 Age 0.410 0.110 3.74 1.87e- 4 ## 4 I(Age^2) -0.00409 0.00182 -2.25 2.44e- 2 ## 5 SexM 5.59 0.256 21.8 4.38e-103 ``` Funkce `glance()` vrací opět tabulku. Jejím obsahem ovšem nejsou údaje týkající se jednotlivých proměnných, ale modelu jako celku: ```r glance(est_model) ``` ``` ## # A tibble: 1 x 11 ## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC ## ## 1 0.603 0.603 10.1 3431. 0 5 -33746. 67504. 67546. ## # … with 2 more variables: deviance , df.residual ``` Tyto funkce jsou mimořádně užitečné při strojovém zpracovávání velkého množství odhadnutých modelů. Funkce z balíku **broom** přitom neslouží jenom pro převádění objektů `lm` do tabulky. Jejich záběr je mnohem obecnější a obsahují metody pro mnohem více tříd objektů. Zvláštní pozornost se vyplatí věnovat funkci `predict()`. Ty v základní konfiguraci vrací prostě vyrovnané hodnoty: ```r all.equal(predict(est_model),fitted(est_model)) ``` ``` ## [1] TRUE ``` Při mnoha aplikacích nás spíše zajímá schopnost modelu predikovat out-of-sample. Tedy predikovat hodnoty vysvětlované proměnné na jiném vzorku dat, než který byl použit k odhadu parametrů. I pro tento účel se hodí funkce `predict()`. Nejprve si tabulku `oly12` rozdělíme na dvě dílčí tabulky `oly12_a` a `oly12_b`: ```r idx <- sample(nrow(oly12), round(nrow(oly12)/2)) oly12[idx,] -> oly12_a oly12[-idx,] -> oly12_b ``` Pro odhad parametrů použijeme tabulku `oly12_a`: ```r oly12_a %>% lm(model, data=.) -> est_model_a ``` Na základě odhadnutého modelu vypočítáme očekávané hodnoty `Weight` pro pozorování ze vzorku `oly12_b`: ```r predict(est_model_a, newdata = oly12_b) -> prediction_oly12_b ``` Výsledkem je vektor predikovaných hodnot. Podívejme se, jak si model vedl při predikci out-of-sample: ```r oly12_b %>% bind_cols(.fitted = prediction_oly12_b) %>% mutate( sample = "B" ) -> oly12_b est_model_a %>% augment() %>% mutate( sample = "A" ) %>% bind_rows(oly12_b) %>% ggplot( aes( x = Weight, y = .fitted ) ) + geom_point( alpha = 0.1 ) + geom_abline(slope = 1, intercept = 0, color = "red") + facet_wrap("sample", labeller = as_labeller( c("A"="In-sample predictions","B"="Out-of-sample predictions") )) ``` Obrázky zhruba ukazují, že odhadnutý model si vede dobře i v out-of-sample predikcích. tedy predikcích vypočítaných na jiných datech, než jaké byly použity pro odhad modelu. Vzhledem k tomu, že vzorek byl rozdělen náhodně to asi není žádné velké překvapení. ### Další estimační funkce a kde je najít Následující tabulka představuje pár estimačních funkcí dostupných v R. |Model/estimátor | Funkce | Poznámka | |------------------|---------------------|-------------------------| | Logit/probit | `glm()` | Modely diskrétní volby | | Poisson | `glm()` | Count data | | Poisson pseudo-ML (PPML) | `glm()` | Např. gravitační modely | | Lineární panelová data | `plm::plm()` | pooled, FE, RE,...| | Arellano-Bond | `plm::pgmmm()` | Dynamické panely | | Common Correlated Effects | `plm::pcce()` | Např. gravitační modely | | Poisson, logit, ordered probit,... (na panelu) | `pglm::pglm()` | | | Tobit | `AER::tobit()` | | | 2SLS | `AER::ivreg()` | | | Heckit (Tobit 2) |`sampleSelection::heckit()` | | ...a mnoho, mnoho dalších. ## Diagnostika ### Vlivná pozorování Populárním způsobem, jak najít pozorování, která významně ovlivňují výsledky modelu, je použití leave-one-out strategií. V nich se sleduje efekt vyloučení jednoho pozorování na odhady koeficientů, popřípadě na vyrovnaná pozorování. V R jsou implementovány v `influence.measures()`: ```r influence.measures(est_model) -> infl_obs infl_obs %>% class ``` ``` ## [1] "infl" ``` `influence.measures()` najednou volá celou řadu diagnostických testů -- mj. `dfbetas()` (vliv jednoho pozorování na odhadnuté koeficienty) a `dffits()` (vliv jednoho pozorování na vyrovnané hodnoty). Výstupní objekt třídy `infl` obsahuje dva prvky (matice). `.$infmat` obsahuje hodnoty jednotlivých statistik a `.$is.inf` rozhodnutí (o podezření), zda se jedná o vlivné pozorování. Viz obrázek: ```r augment(est_model) %>% bind_cols(.,as.data.frame(infl_obs$is.inf)) %>% ggplot( aes(x=Weight, y=.fitted) ) + geom_point( aes( color = dffit ), alpha = 0.2 ) + xlab("Observed values") + ylab("Fitted values") + scale_color_discrete(name="Influential\nobservation\n(DFFIT)") ``` Hodnoty, u kterých se predikce velmi odlišuje od pozorovaných hodnot jsou podle statistiky DFFIT jsou navrženy jako vlivné. Tyto výsledky naznačují, že se v datech může dít něco zvláštního, nebo že nemáme k dispozici dostatečný počet pozorování z určité části prostoru -- v tomto případě například u sportovců s vyšší hmotností. ### Normalita reziduí Prvním způsobem, jak ověřit normalitu reziduí je pohledem na jejich rozdělení. Pomoci může například Q-Q plot, který porovnává kvantily normálního rozdělení a standardizovaných reziduí. V případě, že je rozdělení reziduí normální, by pozorování měla ležet na ose kvadrantu: ```r augment(est_model) %>% ggplot( aes(sample = .std.resid) ) + geom_qq() + geom_abline( slope = 1, intercept = 0, color = "red" ) ``` Ani jeden okometrický test pro normalitu příliš nehovoří. Nicméně je potřeba nepodlehnout zdání. R poskytuje v balíku `normtest` celou baterii statistických testů. Například můžeme zkusit otestovat normalitu pomocí Jarque-Bera testu: ```r library(normtest) residuals(est_model) %>% jb.norm.test() -> jb class(jb) ``` ``` ## [1] "htest" ``` ```r print(jb) ``` ``` ## ## Jarque-Bera test for normality ## ## data: . ## JB = 82199, p-value < 2.2e-16 ``` Výstupem `jb.norm.test()` je objekt typu `htest`. Tato třída je oblíbená pro výsledky statistických testů. Pro tuto třídu existují metody v balíku `broom`, např.: ```r tidy(jb) ``` ``` ## # A tibble: 1 x 3 ## statistic p.value method ## ## 1 82199. 0 Jarque-Bera test for normality ``` Nulovou hypotézou Jarque-Bera testu je normalita, kterou tento test zamítá. ### Testy linearity Z testů linearity je dostupný například Rainbow test, který srovnává kvalitu vyrovnání v modelech odhadnutých "uprostřed" pozorování a na celém vzorku. V R je Rainbow test dostupný v balíku `lmtest` ve funkci `raintest`: ```r raintest(formula, fraction = 0.5, order.by = NULL, center = NULL, data=list()) ``` Parametr `fraction` udává velikost samplu pro odhad modelu "uprostřed". V parametru `order.by` se udává proměnná, podél které se mají pozorování řadit. Parametr `formula` může být naplněn i odhadnutým `lm` modelem. V tom případě není nutné zabývat se parametrem `data`. ```r library(lmtest) raintest(model, order.by = ~ Height, data = oly12) ``` ``` ## ## Rainbow test ## ## data: model ## Rain = 1.1287, df1 = 4519, df2 = 4514, p-value = 2.397e-05 ``` Rain test v tomto případě zamítá nulovou hypotézu. To znamená, že fit na neomezeném vzorku je významně horší, než na podmnožině "prostředních" pozorování. Tento výsledek indikuje možný problém s linearitou. ### Heteroskedasticita a odhad robustních chyb Pokud je porušena homoskedasticita -- tzn. konstantní rozptyl náhodné složky -- potom není OLS estimátor BLUE (Best linear unbiased estimator), nicméně odhady parametrů jsou stále nevychýlené (unbiased). To bohužel neplatí pro odhady jejich rozptylu -- ty naopak vychýlené jsou a v důsledku jsou nedůvěryhodné i testy statistické významnosti parametrů. ```r est_model %>% augment() %>% arrange(Height) %>% ggplot( aes( x = Height, y = .resid ) ) + geom_point( alpha = 0.1 ) + ylab("Residuals") ``` V reálném světě je podmínka homoskedasticity porušena velmi často. Jestli tomu tak je můžeme testovat pomocí mnoha testů. Jako příklad mohou sloužit příbuzné testy `bptest()` z balíku `lmtest` a `ncvTest` z balíku `car`: ```r lmtest::bptest(est_model) ``` ``` ## ## studentized Breusch-Pagan test ## ## data: est_model ## BP = 73.288, df = 4, p-value = 4.585e-15 ``` ```r car::ncvTest(est_model) ``` ``` ## Non-constant Variance Score Test ## Variance formula: ~ fitted.values ## Chisquare = 535.0632, Df = 1, p = < 2.22e-16 ``` Oba testy testují nulovou hypotézu o homoskedasticitě oba ji suverénně zamítají -- p-hodnoty obou testů jsou v podstatě nerozlišitelné od nuly. Odhady směrodatných chyb lze korigovat použitím tzv. HC *(heteroskedasticity consistent)* estimátorů kovarianční matice. Ty jsou implementovány například v balících `car` (`car::hccm`), ale zejména v balíku `sandwich`. Balík `sandwich` obsahuje celou řadu nástrojů pro odhad kovarianční matice pro průřezová i panelová data. Bližší pozornost budeme věnovat dvěma. Základní funkci `vcovHC()` a také relativně nově implementované funkci pro výpočet klastrovaných robustních chyb `vcovCL()`. Obě dvě funkce mají podobnou syntax: ```r vcovHC(x, type = c("HC3", "const", "HC", "HC0", "HC1", "HC2", "HC4", "HC4m", "HC5"), omega = NULL, sandwich = TRUE, ...) ``` ```r vcovCL(x, cluster = NULL, type = NULL, sandwich = TRUE, fix = FALSE, ...) ``` Vstupem funkce je odhadnutý model (`x`). V parametru `type` se potom specifikuje způsob korekce heteroskedasticity (viz `?vcovHC()`). V případě funkce `vcovCL()` je navíc nutné specifikovat proměnnou, nebo v případě více dimenzí proměnné, podle kterých se mají standardní chyby klastrovat. Klastrování se typicky používá v situaci, kdy určitá pozorování na sobě nejsou nezávislá. Klasický příklad je prospěch dětí z jedné třídy. Podstata dat sportovců naznačuje, že i v tomto případě je klastrování namístě. Jako specifické skupiny můžeme chápat sportovce věnující se jedné disciplině. Po `vcovCL()` tedy bude chtít spočítat kovarianční maticí klastrovanou podél proměnné `Sport`. Tu musíme dodat jako samostatný vektor. V objektu totiž není zahrnuta: ```r library(sandwich) vcovCL(est_model, cluster = oly12$Sport) ``` ``` ## (Intercept) Height Age I(Age^2) SexM ## (Intercept) 144.744530958 -87.73269188 0.544857828 -7.405519e-03 9.080987791 ## Height -87.732691881 57.87470724 -0.914600493 1.437195e-02 -6.275686441 ## Age 0.544857828 -0.91460049 0.079432487 -1.394717e-03 0.137949915 ## I(Age^2) -0.007405519 0.01437195 -0.001394717 2.566568e-05 -0.002251259 ## SexM 9.080987791 -6.27568644 0.137949915 -2.251259e-03 0.979011051 ``` Výsledkem je odhad kovarianční matice. Zpravidla nás ovšem více zajímá odhad významnosti parametrů -- provedení statistických testů s použitím robustní kovarianční matice jako vstupu. Pro tento účel použijeme funkci `coeftest()` z balíku `lmtest`: ```r coeftest(x, vcov. = NULL, df = NULL, ...) ``` Prvním vstupem funkce `coeftest()` je odhadnutý model. (Funkce má implementovány metody pro modely třídy `lm` a `glm`.) Druhým parametrem je kovarianční matice nebo funkce, která se má použít pro její výpočet. Pokud je tento parametr ponechán prázdný, potom `coeftest()` použije nekorigovanou kovarianční matici: ```r library(lmtest) coeftest(est_model) ``` ``` ## ## t test of coefficients: ## ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -1.0601e+02 2.3853e+00 -44.4438 < 2.2e-16 *** ## Height 9.4635e+01 1.1413e+00 82.9183 < 2.2e-16 *** ## Age 4.0966e-01 1.0960e-01 3.7378 0.0001868 *** ## I(Age^2) -4.0944e-03 1.8193e-03 -2.2505 0.0244416 * ## SexM 5.5916e+00 2.5603e-01 21.8397 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` Zadaní jména funkce jako parametru `vcov.` je vhodné tehdy, pokud nepotřebujete provádět žádné další nastavení parametrů výpočtu kovarianční matice: ```r coeftest(est_model, vcov. = vcovHC) ``` ``` ## ## t test of coefficients: ## ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -1.0601e+02 2.4526e+00 -43.2250 < 2.2e-16 *** ## Height 9.4635e+01 1.2120e+00 78.0841 < 2.2e-16 *** ## Age 4.0966e-01 9.2268e-02 4.4398 9.109e-06 *** ## I(Age^2) -4.0944e-03 1.5164e-03 -2.7002 0.006944 ** ## SexM 5.5916e+00 2.4404e-01 22.9123 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` V případě klastrovaných chyb však potřebuje takové nastavení provést -- minimálně potřebujeme zadat vektor podél kterého se má klastrovat. Do `vcov.` v tomto případě potřebujeme zadat celou kovarianční matici -- tedy výstup funkce `vcovCL`: ```r coeftest(est_model, vcov. = vcovCL(est_model, cluster = oly12$Sport)) ``` ``` ## ## t test of coefficients: ## ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -1.0601e+02 1.2031e+01 -8.8116 < 2.2e-16 *** ## Height 9.4635e+01 7.6075e+00 12.4396 < 2.2e-16 *** ## Age 4.0966e-01 2.8184e-01 1.4535 0.1461 ## I(Age^2) -4.0944e-03 5.0661e-03 -0.8082 0.4190 ## SexM 5.5916e+00 9.8945e-01 5.6512 1.642e-08 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` Klastrovat lze podél více proměnných. V tomto případě je potřeba do parametru `cluster` vložit list vektorů nebo `data.frame`. Nyní zkusíme cvičně klastrovat nejen podle sportovního odvětví, ale i podle pohlaví: ```r coeftest(est_model, vcov. = vcovCL(est_model, cluster = list(oly12$Sport, oly12$Sex))) ``` ``` ## ## t test of coefficients: ## ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -1.0601e+02 2.6150e+01 -4.0540 5.076e-05 *** ## Height 9.4635e+01 1.2569e+01 7.5295 5.584e-14 *** ## Age 4.0966e-01 3.3167e-01 1.2351 0.2168 ## I(Age^2) -4.0944e-03 4.2600e-03 -0.9611 0.3365 ## SexM 5.5916e+00 1.1367e+00 4.9189 8.857e-07 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` Výstupem funkce `coeftest()` je objekt speciální třídy `coeftest`. Pokud ho chceme převést na tabulku můžeme použít funkci `broom::tidy()`: ```r coeftest(est_model, vcov. = vcovCL(est_model, cluster = list(oly12$Sport, oly12$Sex))) %>% tidy() ``` ``` ## # A tibble: 5 x 5 ## term estimate std.error statistic p.value ## ## 1 (Intercept) -106. 26.1 -4.05 5.08e- 5 ## 2 Height 94.6 12.6 7.53 5.58e-14 ## 3 Age 0.410 0.332 1.24 2.17e- 1 ## 4 I(Age^2) -0.00409 0.00426 -0.961 3.37e- 1 ## 5 SexM 5.59 1.14 4.92 8.86e- 7 ``` ### Dignostika graficky V mnoha případech je užitečné se na diagnostiku podívat prostřednictvím vizualizací. Tuto možnost nabízí například funkce `plot()`, která má pro tento účel implementovanou metodu. Její fungování můžete vyzkoušet zadáním `plot(est_model)`. Tento postup logicky pracuje se základní grafikou. Existují však i řešení, které podobně uživatelsky přítulně vykreslí diagnostické grafy s pomocí `ggplot2`. Jedním z nich je funkce `ggplot2::autoplot()`, která využívá služeb balíku **ggfortify**: ```r library(ggfortify) autoplot(est_model) ``` ``` ## Warning: `arrange_()` is deprecated as of dplyr 0.7.0. ## Please use `arrange()` instead. ## See vignette('programming') for more help ## This warning is displayed once every 8 hours. ## Call `lifecycle::last_warnings()` to see where this warning was generated. ``` ## Odhad více modelů Často můžeme chtít odhadnout model v různých specifikacích, nebo na různých vzorcích dat. ### Odhad různých specifikací na stejných datech Častá úloha je odhad různých modelových specifikací na stejném vzorku pozorování. Pro tyto účely doporučuji následující postup: 1. Vytvořit si `list` se všemi rovnicemi (specifikacemi) 2. Odhadnout model pro všechny specifikace pomocí `lapply()` Nejprve vytvoříme různé modelové specifikace s pomocí funkce `update()`: ```r list( # Baseline model model, # Přidána interakce Sex*Height model %>% update(.~. + Sex*Height), # Přidán fixní efekt na zemi původu atleta model %>% update(.~. + factor(Country)) ) -> Models ``` Nyní pomocí `lapply()` zavoláme `lm()` na všechny prvky listu `Models` -- tedy na všechny modelové specifikace: ```r lapply(Models, lm, data = VGAMdata::oly12) -> Spec_models lapply(Models, function(x) lm(x, data = VGAMdata::oly12)) -> Spec_models ``` Oba výše uvedené způsoby použití `lapply` jsou v pořádku. Z nějakého důvodu však balík `stargazer`, který budeme používat pro tvorbu výstupních tabulek, neumí pracovat s listem odhadnutých modelů, který byl vytvořen prvním způsobem. Výsledné objekty jsou nyní uložené jako položky listu `Spec_models`. K nim můžeme přistupovat pomocí subsetování: ```r Spec_models[[2]] %>% summary ``` ``` ## ## Call: ## lm(formula = x, data = VGAMdata::oly12) ## ## Residuals: ## Min 1Q Median 3Q Max ## -59.122 -5.646 -1.235 3.622 135.458 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -78.890968 3.402356 -23.187 < 2e-16 *** ## Height 78.901728 1.814233 43.490 < 2e-16 *** ## Age 0.392204 0.108875 3.602 0.000317 *** ## I(Age^2) -0.003842 0.001807 -2.126 0.033527 * ## SexM -39.474467 4.065283 -9.710 < 2e-16 *** ## Height:SexM 25.695023 2.313338 11.107 < 2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 10.06 on 9032 degrees of freedom ## (1346 observations deleted due to missingness) ## Multiple R-squared: 0.6084, Adjusted R-squared: 0.6082 ## F-statistic: 2807 on 5 and 9032 DF, p-value: < 2.2e-16 ``` ### Odhad modelu na dílčích datech Základní model můžeme odhadnout na různých vzorcích dat -- v tomto případě pro různé sporty (proměnná `Sport`). Pro tuto operaci využijeme aparát funkcí z balíku **dplyr**. Tabulku `oly12` rozdělíme pomocí `group_by()` do skupin podle hodnoty proměnné `Sport`. Následně pomocí `do()` odhadneme model pro každou skupinu zvlášť. ```r oly12 %>% group_by(Sport) %>% do( est_model = lm(model, data = .) ) -> Sport_models ``` Výsledná proměnná je tabulka (tibble) se dvěma sloupci: grupující proměnnou `Sport` a `est_model`, který obsahuje odhadnuté modely: ```r Sport_models %>% print ``` ``` ## # A tibble: 26 x 2 ## # Rowwise: ## Sport est_model ## ## 1 Archery ## 2 Athletics ## 3 Badminton ## 4 Basketball ## 5 Beach ## 6 Canoe ## 7 Cycling ## 8 Diving ## 9 Equestrian ## 10 Fencing ## # … with 16 more rows ``` S modely uloženými ve sloupci `est_model` můžeme dále pracovat. Následující kód například opět za využití funkce `do()` a `glance()` vytvoří tabulku modelů seřazenou sestupně podle adjustovaného koeficientu determinace: ```r Sport_models %>% do(glance(.$est_model)) %>% bind_cols(Sport_models,.) %>% arrange(desc(adj.r.squared)) %>% print() ``` ``` ## # A tibble: 26 x 13 ## # Rowwise: ## Sport est_model r.squared adj.r.squared sigma statistic p.value df ## ## 1 Tria… 0.879 0.874 2.98 179. 4.32e- 44 5 ## 2 Mode… 0.878 0.870 3.81 110. 3.90e- 27 5 ## 3 Beach 0.874 0.869 4.56 153. 9.52e- 39 5 ## 4 Hand… 0.850 0.848 6.01 423. 9.22e-122 5 ## 5 Hock… 0.837 0.835 4.29 485. 2.43e-147 5 ## 6 Bask… 0.833 0.830 6.50 315. 6.52e- 97 5 ## 7 Tenn… 0.819 0.815 4.49 171. 4.96e- 55 5 ## 8 Rowi… 0.815 0.813 5.84 557. 5.13e-184 5 ## 9 Swim… 0.806 0.805 5.04 887. 2.58e-302 5 ## 10 Canoe 0.797 0.794 5.08 293. 3.55e-102 5 ## # … with 16 more rows, and 5 more variables: logLik , AIC , ## # BIC , deviance , df.residual ``` Je vidět, že nejvyšší $\bar{R}^2$ má model odhadnutý na vzorku triatlonistů. ## Tvorba pěkně formátovaných výsledků s balíkem `stargazer` Pro porovnání výsledků odhadů je užitečné srovnat si je vedle sebe do tabulky, jakou znáte z empirických článků. Tento úkol umí splnit balík Marka Hlaváče `stargazer`. ```r library(stargazer) ``` Stargazer dokáže vytvářet a do ASCII, HTML a LaTeXu exportovat: - regresní tabulky - tabulky popisných statistik - exportovat vstupní `data.frame` bez další transformace (to je užitečné například pro export korelačních tabulek) Balík `stargazer` obsahuje jedinou funkci `stargazer()`, která má velké množství parametrů -- srozumitelný tutoriál najdete například zde: http://jakeruss.com/cheatsheets/stargazer.html Vstupem do funkce může být: - jeden nebo více odhadnutých modelů různých tříd (viz seznam v helpu) - data.framy, vektory nebo matice pro tvorbu popisných statistik nebo přímý export Vstupy mohou být vloženy jednotlivě nebo jako `list`, popřípadě `list` v `listu`. ### Tabulka s popisnými statistikami Tvorbu tabulky s popisnými statistikami je možné ilustrovat na minimalistickém příkladu. Vstupem je v tomto případě `data.frame`. Klíčovým parametrem je `summary = TRUE`, který říká `stargazeru`, že má vytvořit popisné statistiky. V případě `FALSE` by exportoval vstupní tabulku tak jak je. Parametr `type` udává formát, do kterého `stargazer` exportuje. Možná nastavení jsou `html`, `text` a `latex` (default). Zde nevyužitým parametrem je `out` -- do něj se vkládá jméno výstupního souboru. ```r stargazer( VGAMdata::oly12, summary = TRUE, type = "latex" ) ```
StatisticNMeanSt. Dev.MinPctl(25)Pctl(75)Max
Age10,38426.0695.44113222971
Height9,8231.7690.1131.3201.6901.8502.210
Weight9,10472.85316.06736.00061.00081.000218.000
Gold10,3840.0170.1360002
Silver10,3840.0170.1330002
Bronze10,3840.0180.1360002
Total10,3840.0520.2500005
V základním nastavení `stargazer` vrací počet pozorování, průměr, směrodatnou odchylku, minimum a maximum. Výčet statistik, stejná jako formátování tabulky, lze měnit pomocí parametrů funkce `stargazer()`. **Pozor, stargazer neumí zpracovat tibble!** Tabulka `oly12` je tibble: ```r stargazer( oly12, summary = TRUE, type = "latex" ) ``` Výstup bude vypadat následovně:
StatisticNMeanSt. Dev.MinPctl(25)Pctl(75)Max
Pro správnou funkci je potřeba tabulku odtibblovat -- například pomocí `as.data.frame(oly12)`. ### Tabulka s regresními modely Pokud jsou vstupem stagazeru odhadnuté modely nebo jejich list, potom `stargazer` automaticky vytvoří regresní tabulku. Pro ilustraci můžeme použít list `Spec_models`. Pro vytvoření regresní tabulky ze `Spec_models` je vhodné použít další z mnoha parametrů funkce `stargazer()`. Jeden z modelů ve `Spec_models` obsahuje velké množství fixních efektů. Jejich odhady nás v podstatě nezajímají, ale jejich výpis by tabulku zvětšil do obrovských rozměrů. Použijeme proto parametr `omit`, který umožňuje potlačit výpis parametrů, jejichž jména odpovídají zadanému regulárnímu výrazu. V parametru `omit.labels` je možné nastavit jméno, které se má použít pro signalizaci přítomnosti nevytisknutých proměnných v regresi. ```r stargazer( Spec_models, type = pandoc.output.format(), omit = "factor\\(Country\\)", omit.labels = "Country FE" ) ```
Dependent variable:
Weight
(1)(2)(3)
Height94.635***78.902***94.887***
(1.141)(1.814)(1.234)
Age0.410***0.392***0.491***
(0.110)(0.109)(0.110)
I(Age2)-0.004**-0.004**-0.005***
(0.002)(0.002)(0.002)
SexM5.592***-39.474***5.502***
(0.256)(4.065)(0.266)
Height:SexM25.695***
(2.313)
Constant-106.013***-78.891***-112.691***
(2.385)(3.402)(6.240)
Country FENoNoNo
Observations9,0389,0389,038
R20.6030.6080.628
Adjusted R20.6030.6080.620
Residual Std. Error10.126 (df = 9033)10.058 (df = 9032)9.909 (df = 8837)
F Statistic3,430.928*** (df = 4; 9033)2,806.601*** (df = 5; 9032)74.631*** (df = 200; 8837)
Note:*p<0.1; **p<0.05; ***p<0.01
Za povšimnutí stojí, že `stargazer` defaultně používá jinou hvězdičkovou konvenci, než je tomu ve zbytku R. Ve v tabulce vytvořené stargazerem vidíte stejnou konvenci, na kterou jste zvyklí z Gretlu. R defaultně hvězdičkami více šetří -- viz `summary` výše. > Stargazer poskytuje extrémně užitečnou funkcionalitu, nicméně celkově se jedná o dost nemoderní balík s velmi složitým kódem, který je náchylný k chybám a divnému chování. Výše ukázané tabulky neobsahují robustní chyby. Protože z testů víme, že rezidua jsou heteroskedastická, musíme je do tabulek dostat. V defaultním nastavení `stargazer` získává odhad standardních chyb vnitřním voláním `summary`. Uživatel však může vložit vlastní odhady robustních chyb do argumentu `se`. Ten očekává `list` numerických vektorů. (V nápovědě k parametru `se` je v aktuální verzi chyba.) Postup je následující: 1. Nejprve vytvoříme funkci `get.se`, která odhadne robustní chyby a vrátí je jako vektor. 2. Pomocí `lapply` aplikujeme funkci `get.se` na všechny modely v listu `Spec_models`. 3. Výsledný list vektorů vložíme do parametru `se` `stargazeru`. Zároveň totéž provedeme pro p-hodnoty. ```r get.se <- function(x) coeftest(x, vcov. = vcovHC) %>% tidy %>% pull(std.error) get.pval <- function(x) coeftest(x, vcov. = vcovHC) %>% tidy %>% pull(p.value) se.list <- lapply(Spec_models, get.se) pval.list <- lapply(Spec_models, get.pval) stargazer( Spec_models, type = pandoc.output.format(), omit = "Country", omit.labels = "Country FE", se = se.list, p = pval.list ) ```
Dependent variable:
Weight
(1)(2)(3)
Height94.635***78.902***94.887
(1.212)(1.682)
Age0.410***0.392***0.491
(0.092)(0.091)
I(Age2)-0.004***-0.004**-0.005
(0.002)(0.002)
SexM5.592***-39.474***5.502
(0.244)(4.125)
Height:SexM25.695***
(2.352)
Constant-106.013***-78.891***-112.691
(2.453)(3.149)
Country FENoNoNo
Observations9,0389,0389,038
R20.6030.6080.628
Adjusted R20.6030.6080.620
Residual Std. Error10.126 (df = 9033)10.058 (df = 9032)9.909 (df = 8837)
F Statistic3,430.928*** (df = 4; 9033)2,806.601*** (df = 5; 9032)74.631*** (df = 200; 8837)
Note:*p<0.1; **p<0.05; ***p<0.01
Pokud bychom chtěli vytvořit regresní tabulku z modelů odhadnutých pro jednotlivé sporty, můžeme využít toho, že `data.frame` je ve své podstatě `list`: ```r stargazer( Sport_models$est_model, type = "html" ) ```