Intuice (zkušenost, teorie) nám říká, že existuje vztah mezi:

• velikostí mzdy a pohlavím
• výškou a váhou člověka
• počtem hvězd na obloze a velikostí střešního okna
• …

Ekonometrická analýza nám umožňuje otestovat naše hypotézy na datech.

• BPE_ZAEK Základy ekonometrie
• MPE_EKON Ekonometrie
• BPE_CARA Časové řady
• MPE_BAAN Bayesiánská analýza
• BPE_AVED Analýza a vizualizace ekonomických dat (pouze potřebné technické nástroje)

Fungování ekonometrických nástrojů můžeme ilustrovat na příkladu datasetu s výškou, váhou, věkem a pohlavím sportovců, kteří se zůčastnili olympijských her v Londýně v roce 2012 (dostupný v VGAMdata::oly12):

print(oly12, n=5)

## # A tibble: 9,038 × 4
##   Height Weight    Sex   Age
##    <dbl>  <int> <fctr> <int>
## 1   1.70     60      M    23
## 2   1.93    125      M    33
## 3   1.87     76      M    30
## 4   1.78     85      F    26
## 5   1.82     80      M    27
## # ... 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í:

Našim cílem je najít takové parametry \(\alpha\) a \(\beta\) (tj. úrovňovou konstantu a sklon) přímky, které zajistí "nejlepší proložení".

Data již máme v tabulce oly12. Máme i představu o rovnici, kterou chceme odhadnout:

\[Weight = \alpha + \beta Height + \varepsilon\] Tuto rovnici však musíme zapsat tak, aby byla srozumitelná i pro R.

K tomu slouží objekty třídy formula. Jejich formulace následuje tuto syntax:

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:

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:

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:

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:

model <- Weight ~ log(Height)
model <- Weight ~ Height >= 1.8

První varianta například odpovídá specifikaci:

\[Weight = \alpha + \beta \log(Height) + \varepsilon\]

V rovnici může být přimo přítomna celá řada funkcí – s výjimkou těch, jejichž interpratace 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:

model <- Weight ~ I(Height + Age)

Třída forrmula umožňuje jednoduše jednoduše definovat i interakční členy pomocí operátorů : a *:

model <- Weight ~ Height:Sex

Chápe R jako

\[Weight = \alpha + \beta Height\times Sex_{male} + \gamma Height\times Sex_{female} + \varepsilon\]

Další možností jak popsat interakci proměnných je *:

model <- Weight ~ Height*Sex

Tato formula odpovídá jiné rovnici:

\[Weight = \alpha + \beta Height + \epsilon Sex_{male} + \gamma Height\times Sex_{male} + \varepsilon\]

Pro třídu formula jsou naimplementovány zajímavé metody. Velmi užitečná je například funkce update():

update(old,new,...)

Kde old je původní rovnice a new pravidla, která stanovují, jak se původní rovnice má upravit.

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:

model %>% 
  update(. ~ . -Age)

## Weight ~ Height

Nebo je možné nahradit celou stranu rovnice:

model %>% 
  update(. ~ Age)

## Weight ~ Age

Pro odhad použijeme model v následující specifikaci:

model <- Weight ~ Height + Age + Sex

Proměnné Weight, Height a Age jsou numerické spojité proměnné. Sex je zcela jiný případ. Je to kategoriální proměnná, která nemůže vstoupit do regrese přímo.

Před samotným odhadem je jí nutné transformovat na vektor (matici v případě více než dvou úrovní) nul a jedniček.

Estimační funkce si všechny vstupní faktory transformuje sama pomocí interního volání funkce model.matrix():

model.matrix(model, data = oly12) %>% head

##   (Intercept) Height Age SexM
## 1           1   1.70  23    1
## 2           1   1.93  33    1
## 3           1   1.87  30    1
## 4           1   1.78  26    0
## 5           1   1.82  27    1
## 6           1   1.82  30    0

Pokud je vstupní proměnná character, potom si ji model.matrix() vnitřně konvertuje na faktor a provede transformaci na vektor/matici dummy proměnných:

oly12 %>% 
    mutate(
        Sex = as.character(Sex)
    ) %>% 
    model.matrix(model, data = .) %>% 
    head

##   (Intercept) Height Age SexM
## 1           1   1.70  23    1
## 2           1   1.93  33    1
## 3           1   1.87  30    1
## 4           1   1.78  26    0
## 5           1   1.82  27    1
## 6           1   1.82  30    0

Přesto existuje dobrý důvod, proč může uživatel chtít provést konverzi do faktoru sám.

Počet sloupců, které vzniknou po transformaci faktoru je (v případě, že modelová specifikace obsahuje konstantu) o jeden nižší, než je počet úrovní faktoru (přítomných v datech). Jedna úroveň faktoru je totiž použita jako referenční a v modelu není explicitně obsažena. Odhadnuté koeficienty pro zachované úrovně následně ukazují rozdíl v průměrné hodnotě vůči referenční úrovni.

V předchozích příkladech se proměnná Sex transformovala do vektoru s 1 pro muže a 0 pro ženy. Toto chování se dá změnit nastavením referenční úrovně faktoru.

S pomocí base funkce:

oly12 %>% 
    mutate(
        Sex = relevel(Sex, ref="M")
    ) %>% 
    model.matrix(model, data = .)

Nebo funkce fct_relevel() z balíku forcats (součást tidyverse):

oly12 %>% 
    mutate(
        Sex = forcats::fct_relevel(Sex, "M")
    ) %>% 
    model.matrix(model, data = .) %>% 
    head

##   (Intercept) Height Age SexF
## 1           1   1.70  23    0
## 2           1   1.93  33    0
## 3           1   1.87  30    0
## 4           1   1.78  26    1
## 5           1   1.82  27    0
## 6           1   1.82  30    1

Pokud je proměnná numerická (typicky integer), ale chceme, aby její hodnoty byly chápány jako úrovně faktoru, je nutné provést transformaci. Buď v datech, nebo přímo v rovnici, např.:

Weight ~ Height + factor(Age) + Sex

Konverze numerických hodnot do faktoru je poměrně častá. Zvláště u dotazníkových dat čísla mohou pouze kódovat hodnoty kategoriální proměnné. Dalším připadem je například vytvoření fixních efektů pro roky.

R umožňuje tvorbu tzv. ordered faktorů. Ty doporučujeme nepoužívat – pokud si nejste jisti, co děláte.

S připraveným modelem a daty je možné volat funkci lm():

lm(model, data = oly12) -> est_model

print(est_model)

## 
## Call:
## lm(formula = model, data = oly12)
## 
## Coefficients:
## (Intercept)       Height          Age         SexM  
##    -103.040       94.898        0.167        5.573

Produkt funkce lm() je teď přiřazen do proměnné est_model. To je S3 objekt zvláštní třídy lm:

class(est_model)

## [1] "lm"

Tento objekt si můžeme představit jako list se standardizovanou strukturou:

names(est_model)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"

Její základní položky jsou:

• coefficients – pojmenovaný vektor odhadnutých koeficientů

est_model$coefficients

## (Intercept)      Height         Age        SexM 
## -103.040347   94.898203    0.166981    5.573154

• residuals – vektor reziduí – tj. odhadů náhodné složky \(\varepsilon\)

est_model$residuals %>% head

##         1         2         3         4         5         6 
## -7.700316 33.803288 -9.001877 14.780039  0.243976 -1.683813

• fitted.values – vyrovnané hodnoty

est_model$fitted.values %>% head

##        1        2        3        4        5        6 
## 67.70032 91.19671 85.00188 70.21996 79.75602 74.68381

Třídu lm dědí i výsledky jiných estimačních funkcí. Například glm() vrací objekty s třídami glm a lm:

glm(model, data = oly12) %>% class

## [1] "glm" "lm"

Uživatelé Gretlu a jiného statistického software jsou zvyklí, že odhad nevrací jen velmi střídmý přehled odhadnutých koeficientů, ale že ukazuje velmi bohatou tabulku s výsledky různých testů atp.

Takový výstup můžeme získat pomocí funkce summary:

summary(est_model)

## 
## Call:
## lm(formula = model, data = oly12)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.277  -5.770  -1.352   3.754 135.731 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -103.0404     1.9867 -51.864   <2e-16 ***
## Height        94.8982     1.1355  83.572   <2e-16 ***
## Age            0.1670     0.0196   8.518   <2e-16 ***
## SexM           5.5732     0.2560  21.774   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.13 on 9034 degrees of freedom
## Multiple R-squared:  0.6028, Adjusted R-squared:  0.6027 
## F-statistic:  4571 on 3 and 9034 DF,  p-value: < 2.2e-16

Funkce tidy() z balíku broom vytvoří z objektu třídy lm data.frame s odhady parametrů, jejich směrodatnými chybami, t-statistikami a p-hodnotami:

est_model %>% tidy

##          term    estimate  std.error  statistic       p.value
## 1 (Intercept) -103.040347 1.98673555 -51.864149  0.000000e+00
## 2      Height   94.898203 1.13552521  83.572080  0.000000e+00
## 3         Age    0.166981 0.01960294   8.518163  1.881305e-17
## 4        SexM    5.573154 0.25595405  21.774042 1.711081e-102

Funkce broom::glance() umožňuje získat jednořátkovou summarizaci "statistik" pro odhad modelu:

est_model %>% glance

##   r.squared adj.r.squared    sigma statistic p.value df    logLik      AIC
## 1 0.6028396     0.6027077 10.12833  4570.826       0  4 -33748.38 67506.76
##        BIC deviance df.residual
## 1 67542.31 926736.1        9034

Funkce tidy() a glance() jsou zvláště užitečné při strojovém zpracování velkého množstvé modelů.

Abychom mohli srovnat pozorované a vyrovnané hodnoty potřebujeme data.frame ve kterém budou obě dvě. Momentálně máme pozorované hodnoty v tabulce oly12, popřípadě v est_model$model, a predikované hodnoty v est_model$fitted.values. Spojit je dohromady nám pomůže funkce augment z balíku broom:

augment(est_model) %>% as_tibble() %>% print(n=3)

## # A tibble: 9,038 × 11
##   Weight Height   Age    Sex  .fitted   .se.fit    .resid         .hat
##    <int>  <dbl> <int> <fctr>    <dbl>     <dbl>     <dbl>        <dbl>
## 1     60   1.70    23      M 67.70032 0.2155292 -7.700316 0.0004528311
## 2    125   1.93    33      M 91.19671 0.2172001 33.803288 0.0004598795
## 3     76   1.87    30      M 85.00188 0.1622649 -9.001877 0.0002566690
## # ... with 9,035 more rows, and 3 more variables: .sigma <dbl>,
## #   .cooksd <dbl>, .std.resid <dbl>

Ta vrací původní datovou tabulku rozšířenou o dodatečné sloupce (rezidua, predikované hodnoty a další). Momentálně je podstatné i to, že vrací tidy data.frame. Její výstup můžeme použít pro vytvoření obrázku:

augment(est_model) %>% 
  ggplot(
    aes(
      x = Weight,
      y = .fitted
    )
  ) +
  geom_point(
    alpha = 0.1
  ) +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  geom_smooth(method = "lm") +
  theme_linedraw() -> p

Do prostého bodového grafu je na obrázku pomocí geom_abline() vložena i červená přímka zobrazující osu kvadrantu. V ideálním případě by všechny hodnoty měly ležet právě na ní.

Modrou barvou je vyznačeno proložení vztahu mezi pozorovanou a predikovanou hodnotou. V idálním případě by měly obě přímky splývat do jedné. Tady se to ovšem neděje. Zejména u vyšších vah je proložení vychylováno shlukem vlivných pozorování.

Pro popis kvality proložení (předpovědi) se používají i statistiky jako je například Mean squared error (MSE) – tedy průměr čtverců reziduí. Tato statistika se dá vypočítat ze součástí objektu lm. Buď je možné přímo adresovat vektor reziduí (est_model$residuals), nebo použít speciální funkci residuals(), která je extrahuje:

residuals(est_model)^2 %>% mean

## [1] 102.5377

Popřípadě lze použít funkci deviance(). Výsledek musí být identický:

deviance(est_model)/nobs(est_model)

## [1] 102.5377

To znamená, že průměrná vzdálenost pozorované a vyrovnané hodnoty je 10.13 Kg.

Doporučeníhodné je použití funkci residuals() a to zejména proto, že umí extrahovat různé typy reziduí. To je užitečné zejména u jiných estimačních funkcí – typicky glm() (logit, probit, poisson, PPML,…).

Kromě vykreslení vyrovnaných hodnot (fitted values) je možné vypočítat pomocí funkce predict() i predikované hodnoty. Při zavolání funkce predict() bez nastavení dalších parametrů vrací vyrovnané hodnoty.

all.equal(predict(est_model),fitted(est_model))

## [1] TRUE

To znamená, že predikované hodnoty jsou vypočítány pro pozorování z tabulka, která byla využita pro odhad modelu – jedná se tzv. o in-sample predikce.

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ž ktrerý 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:

# Náhodně vybereme polovinu řádků
idx <- sample(nrow(oly12), round(nrow(oly12)/2))

oly12[idx,]  -> oly12_a

oly12[-idx,]  -> oly12_b

Pro odhad parametrů použijeme tabulku oly12_a:

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:

predict(est_model_a, newdata = oly12_b) -> prediction_oly12_b

Prvním způsobem, jak ověřit normalitu reziduí je pohledem na jejich rozdělení:

augment(est_model) %>% 
    ggplot(
        aes(x = .resid)
    ) +
    geom_histogram(
        binwidth = 1
    ) +
    geom_density(
        aes(y = 1 * ..count..),
        color = "blue"
    ) +
    theme_linedraw() -> p

Prostou okonometrií se zdá, že rezidua mohou být od normálního rozdělení odchýlena.

Jiný pohled poskytuje Q-Q plot, který porovnává kvantily normálního rozdělení a standardizovaných reziduí. V případě, že je rozdělení reizudí normální, by pozorování měla ležet na ose kvadrantu:

augment(est_model) %>% 
    ggplot(
        aes(sample = .std.resid)
    ) +
    geom_qq() +
    geom_abline(
        slope = 1, 
        intercept = 0,
        color = "red"
        ) +
    theme_linedraw() -> p

Výstupem jb.norm.test() je objekt typu htest. Tato třída je oblíbená pro výsledky statistických testů. Pro tuto ptřídu existují metody v balíku broom, např.:

tidy(jb)

##   statistic p.value                         method
## 1  82320.78       0 Jarque-Bera test for normality

Nulovou hypotézou Jarque-Bera testu je normalita, kterou tento test zamítá.

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 vzroku.

V R je Rainbow test dostupný v balíku lmtest ve funkci raintest:

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.

library(lmtest)
raintest(model, order.by = ~ Height, data = oly12)

## 
##  Rainbow test
## 
## data:  model
## Rain = 1.1299, df1 = 4519, df2 = 4515, p-value = 2.047e-05

Rain test v tomto případě zamítá nulovou hypotézu. To znamená, že fit na neomezeném vzroku je významně horší, než na podmožině "prostředních" pozorování. Tento výsledek indikuje možný problém s linearitou.

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ů.

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 šloužit přibuzné testy bptest() z balíku lmtest a ncvTest z balíku car:

library(lmtest)
bptest(est_model)

## 
##  studentized Breusch-Pagan test
## 
## data:  est_model
## BP = 70.186, df = 3, p-value = 3.895e-15

library(car)
ncvTest(est_model)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 535.4902    Df = 1     p = 1.805286e-118

Oba testy testují nulovou hypotézu o homoskedasticitě oba ji suveréně zamítají – p-hodnoty obou testů jsou v podsatatě nerozlišitelné od nuly.

To je problém. Nyní už víme, že nelze věřit testům významnosti parametrů, které vrací volání summary() – jsou totiž vychýlené.

Jejich odhad jde naštestí opravit použitím tzv. HC (heteroskedasticity consistent) estimátorů kovarianční matice. Ty jsou naimplementovány například v balících car (car::hccm) nebo sandwich (sandwich::vcovHC):

vcovHC(x,
  type = c("HC3", "const", "HC", "HC0", "HC1", "HC2", "HC4", "HC4m", "HC5"),
  omega = NULL, sandwich = TRUE, ...)

Pokud obsahuje definice funkce více možných nastavení parametru (zde například vektor v parametru type), potom se první z nich chápe jako defaultní nastavení.

Vstupem funkce je odhadnutý model třídy lm. Funkce také umožňuje zvolit si v parametru type způsob korekce heteroskedasticity. Výběr je velký (pro bližší popis viz help). Balík sandwich obsahuje i HAC (heteroskedasticity and autocorrelation consistent) estimátory kovarianční matice ve vcovHAC.

Výstupem funkce je kovarianční matice, kterou můžeme použít pro proveení testů významnosti odhadnutých parametrů. Ten provedeme pomocí funkce coeftest() z balíku lmtest.

coeftest(x, vcov. = NULL, df = NULL, ...)

coeftest() očekává jako první argument odhadnutý model třídy lm. V, pro nás monetálně klíčovém, parametru vcov. se vkládá robustní kovarianční matice, nebo funkce, která se použije pro její odhad.

library(sandwich)
coeftest(est_model, vcov. = vcovHC)

## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) -103.040347    2.133328 -48.300 < 2.2e-16 ***
## Height        94.898203    1.213206  78.221 < 2.2e-16 ***
## Age            0.166981    0.018707   8.926 < 2.2e-16 ***
## SexM           5.573154    0.243771  22.862 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

V této konfigurace je použita defaultní hodnota type u funkce vcovHC().

Č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šchny specifikace pomocí lapply()

Najprve vytvoříme různé modelové specifikace s pomocí funkce update():

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:

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 neumí pracovat s listem odhadnutých modelů, který byl vytvořen prvním callem.

Výsledné objekty jsou nyní uložené jako položky listu Spec_models. K nim můžeme přistupovat pomocí subsetování:

Spec_models[[2]] %>% summary

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.

library(stargazer)

Stargazer dokaž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é nožství parametrů – srozumitelný tutorial 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.

V defaultní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:

Výstup bude vypdat následovně:

Pro správnou funkci je potřeba tabulku odtibblovat – například pomocí as.data.frame(oly12).

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.

stargazer(
    Spec_models,
    type = "html",
    omit = "factor\\(Country\\)",
    omit.labels = "Country FE"
)

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.

Postup je náseldující:

1. Nejprve vytvoříme funkci get.se, která odhadne robustní chyby a vrátí je jako vektor.

get.se <- function(x) coeftest(x, vcov. = vcovHC(x, type = "HC0")) %>% 
    tidy %>% extract2("std.error")

1. Pomocí lapply aplikujeme funci get.se na všechny modely v listu Spec_models.

se.list <- lapply(Spec_models, get.se)

1. Výsledný list vektorů vložíme do parametru se stargazeru.

stargazer(
    Spec_models,
    type = "html",
    omit = "factor\\(Country\\)",
    omit.labels = "Country FE",
    se = se.list
)

Pro odhad koeficientů \(\beta\) použijte OLS estimátor a data ze souboru HW_wages_select.Rdata. Ty obsahují následující proměnné:

• surveyy – rok sběru dat
• gender – pohlaví (1 \(=\) žena) – kategoriální proměnná
• wagegrhr – hrubá hodinová mzda
• tenuexpe – délka praxe v letech
• educalmh – úroveň vzdělání (1 = nízká, 2 = střední, 3 = vysoká) – kategoriální proměnná
• age – věk respondenta

V proměnné gender nastavte hodnotu žena jako referenční. Pozor, některé proměnné jsou kategoriální – ve skutečnosti nejde o čísla. Ta jen kódují různé úrovně.

Pamatujte:

• Na pořadí řádků a sloupců v tabulce emodel nezáleží.
• Úkol bude hodnocen s pomocí jiných dat. Ruční vložení adjR2 a cor_OF bude mít za výsledek nula bodů.