AVED: Ekonometrie v R

Zdroje

Kleiber, C., & Zeileis, A. (2008). Applied econometrics with R. Springer Science & Business Media.

Úvod

R obsahuje mnoho různých nástrojů, které téměř kompletně pokrývají potřeby empirického výzkumu. Základní estimátory (lm() a glm()) jsou implementovány v základním R, ale pro mnoho funkcí je nutné sahat do dalších balíků, které nejsou bohužel příliš standardizované.

Účelem tohoto handoutu není předvést všechny ekonometrické nástroje, které jsou v R dostupné – to by díky jejich počtu ani nebylo možné – ale vysvětlit logiku fungování základní sady nástrojů. Praktický příklad ukazuje v mnoha ohledech nejjednoduší druh regresní analýzy – odhad modelu pomocí OLS na průřezových datech. Její pochopení Vám umožní lehce pochopit fungování i složitějších nástrojů.

---

Pro každou ekonometrickou analýzu potřebujete několik věcí:

1. Modelovou specifikaci (tj. rovnicí nebo rovnicemi popsanou představu o tom, jak svět funguje)
2. Data
3. Estimátor

Když dáte tyto věci dohromady získáte výsledky, které můžete dále zpracovávat, analyzovat nebo rovnou zahodit…

Používání ekonometrických nátrojů v R se řídí stejnou logikou:

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:

oly12 %>% 
  ggplot(
    aes(
      x = Height,
      y = Weight
    )
  ) +
  geom_point(
    alpha = 0.1
  ) +
  theme_linedraw()

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í”.

Pro jejich nalezení (odhad modelu) je potřeba najít vhodný estimátor. Pro odhad parametrů použijeme metodu nejmenších čtverců (OLS), která je implementována ve funkci lm():

lm(formula, data, subset, weights, na.action,
   method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
   singular.ok = TRUE, contrasts = NULL, offset, ...)

Funkce lm() má mnoho parametrů:

• formula je objekt třídy “formula”, který obsahuje symbolický popis rovnice, která se má odhadnout
• data je data.frame který obsahuje data, která se mají použít. Jména sloupců musí odpovídat proměnným použitým v parametru formula
• 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.
• a další…

Model: Formula

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\]

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 (lm,…) transformována na vektor s hondnotami 1 a 0. Viz příklad takovéto transformace v tabulce:

oly12 %>% 
  mutate(
    cond_l = Height >= 1.8
  ) %>% 
  mutate(
    cond_n = as.numeric(cond_l)
  ) %>% print(n=5)

## # A tibble: 9,038 × 6
##   Height Weight    Sex   Age cond_l cond_n
##    <dbl>  <int> <fctr> <int>  <lgl>  <dbl>
## 1   1.70     60      M    23  FALSE      0
## 2   1.93    125      M    33   TRUE      1
## 3   1.87     76      M    30   TRUE      1
## 4   1.78     85      F    26  FALSE      0
## 5   1.82     80      M    27   TRUE      1
## # ... with 9,033 more rows

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

Výše uvedený příklad by se tak R tlumočil jako

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\] 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 *:

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 \(^\). 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:

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 naimplementovány i některé zajímavé metody. Velmi užitečná je například funkce update() s velmi prostou syntaxí:

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

Odhad modelu

Při volání funkce lm() nebo jiné estimační funkce se dohromady spojí všechny tři komponenty: modelová specifikace, data a estimátor.

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

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.

Volání lm()

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"

Výsledky a predikce

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

Všimněte si, že summary() používá méně obvyklé – a přísnější – hvězdičkové značení!

V R najdete řadu dalších funkcí, které pracují s objekty třídy lm a extrahují (nebo dopočítávají) z něj užitečné informace:

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í

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

Často nás zajímá, jak dobře model predikuje. Kromě různých statistik je zajímá se podívat na bodový graf, kde na jedné osa jsou pozorované hodnoty vysvětlované proměnné (Weight) a na druhé její očekávané hodnoty – tj. predikce modelu.

Abychom mohli takový obrázek vytvořit, potřebujeme data.frame ve kterém budou pozorované i predikované hodnoty. 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()

## # 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.70031552 0.0004528311
## 2     125   1.93    33      M 91.19671 0.2172001 33.80328752 0.0004598795
## 3      76   1.87    30      M 85.00188 0.1622649 -9.00187722 0.0002566690
## 4      85   1.78    26      F 70.21996 0.1802440 14.78003929 0.0003166983
## 5      80   1.82    27      M 79.75602 0.1434618  0.24397601 0.0002006302
## 6      73   1.82    30      F 74.68381 0.2186361 -1.68381294 0.0004659807
## 7      75   1.87    23      M 83.83301 0.1684282 -8.83301004 0.0002765372
## 8      80   1.90    27      M 87.34788 0.1639915 -7.34788024 0.0002621600
## 9      60   1.73    28      F 65.80901 0.1670191 -5.80901261 0.0002719296
## 10     64   1.71    28      F 63.91105 0.1651660  0.08895145 0.0002659287
## # ... with 9,028 more rows, and 3 more variables: .sigma <dbl>,
## #   .cooksd <dbl>, .std.resid <dbl>

Pro podobnou funkcionalitu jakou poskytuje balík broom můžete použít i balík modelr. broom však obsahuje nástroje pro mnohem širší škálu tříd.

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()

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 klastrem 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

Výsledkem je vektor predikovaných hodnot. Podívejme se, jak si model vedl při predikci out-of-sample:

oly12_b %<>%
    cbind(Weight_predicted = prediction_oly12_b) %>% 
    mutate(
        prediction_error = Weight - Weight_predicted
    )

oly12_b %>% 
    ggplot(
        aes(
            x = Weight,
            y = Weight_predicted
        )
    ) +
    geom_point(
        alpha = 0.1
    ) +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  geom_smooth(method = "lm") +
  theme_linedraw()

Můžeme opět vypočítat MSE – v tomto případě tedy MSPE (mean squared prediction error):

oly12_b$prediction_error^2 %>% mean

## [1] 103.8298

To znamená, že v out-of-sample predikci je průměrná vzdálenost mezi pozorovanou a očekávanou hodnotou 10.19 Kg. Model si tak vede stejně dobře (špatně) v in-sample a out-of-sample predikcích. (To je obecně dobré znamení.)

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 naimplementovány v influence.measures():

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:

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)") +
    theme_bw()

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

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

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()

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:

library(normtest)

residuals(est_model) %>% 
    jb.norm.test() -> jb

class(jb)

## [1] "htest"

print(jb)

## 
##  Jarque-Bera test for normality
## 
## data:  .
## JB = 82321, 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 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á.

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

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

est_model %>% 
    augment %>% 
    arrange(Height) %>% 
    ggplot(
        aes(
            x = Height,
            y = .resid
            )
    ) +
    geom_point(
        alpha = 0.1
    ) +
    ylab("Residuals") +
    theme_bw()

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é:

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

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(). Pokud chceme zvolit jiný type musíme do parametru vcov. vložit již odhadnutou kovarianční matici:

coeftest(est_model, vcov. = vcovHC(est_model, type="HC0"))

## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept) -103.040347    2.131765 -48.3357 < 2.2e-16 ***
## Height        94.898203    1.212318  78.2783 < 2.2e-16 ***
## Age            0.166981    0.018689   8.9348 < 2.2e-16 ***
## SexM           5.573154    0.243641  22.8744 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

V tomto případě jsou výsledky při použití defaultní HC3 a HC0 v podstatě identické. Odhadnuté koeficienty zůstavají statisticky významné i při použití robustních chyb.

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š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

## 
## Call:
## lm(formula = x, data = VGAMdata::oly12)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -59.075  -5.639  -1.235   3.645 135.519 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -76.03684    3.12692 -24.317   <2e-16 ***
## Height       79.11129    1.81191  43.662   <2e-16 ***
## Age           0.16446    0.01947   8.446   <2e-16 ***
## SexM        -39.60004    4.06565  -9.740   <2e-16 ***
## Height:SexM  25.75678    2.31361  11.133   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.06 on 9033 degrees of freedom
##   (1346 observations deleted due to missingness)
## Multiple R-squared:  0.6082, Adjusted R-squared:  0.608 
## F-statistic:  3506 on 4 and 9033 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).

VGAMdata::oly12 %>%
    # Vyřadíme pozorování s chybějícími pozorováními na LHS 
    filter(!is.na(Weight)) %>% 
    # Pokud atlet nastoupil ve více sportech tak se zkopíruje
    # na více řádků -- každý pro jiný sport (tidyr)
    separate_rows(Sport,sep=",") %>% 
    # Na některých řádcích jsou chyby -- mezery na koncích jmen 
    # sportů. Ty je pro správné grupování potřeba odstranit:
    rowwise() %>% 
    mutate(
        Sport = str_trim(Sport)
    ) %>% 
    # Některé sporty (moderní gymnastika a synchronizované plavání)
    # jsou otevřeny pouze ženám. Ty musíme kvůli specifikaci modelu
    # vyloučit:
    group_by(Sport) %>% 
    mutate(
        Males_total = sum(Sex=="M")
    ) %>% 
    filter(Males_total != 0) %>%
    # S pomocí funkce dplyr::do se model odhadne pro každou 
    # skupinu pozorování
    do(
        est_model = lm(Weight ~ Height + Age + Sex, data = .)
        ) -> Sport_models

Výsledná proměnná je tibble se dvěma sloupci: groupující proměnnou Sport a est_model, který obsahuje odhadnuté modely:

Sport_models %>% class

## [1] "rowwise_df" "tbl_df"     "tbl"        "data.frame"

Sport_models %>% print

## Source: local data frame [30 x 2]
## Groups: <by row>
## 
## # A tibble: 30 × 2
##                      Sport est_model
## *                    <chr>    <list>
## 1                  Archery  <S3: lm>
## 2                Athletics  <S3: lm>
## 3                Badminton  <S3: lm>
## 4               Basketball  <S3: lm>
## 5         Beach Volleyball  <S3: lm>
## 6             Canoe Slalom  <S3: lm>
## 7             Canoe Sprint  <S3: lm>
## 8            Cycling - BMX  <S3: lm>
## 9  Cycling - Mountain Bike  <S3: lm>
## 10          Cycling - Road  <S3: lm>
## # ... with 20 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() vytvoří tabulku se jménem sportu a adjustovaným koeficientem determinace:

Sport_models %>% 
    do(glance(.$est_model)) %>% 
    bind_cols(Sport_models,.) %>% 
    select(Sport,adj.r.squared) %>% 
    arrange(desc(adj.r.squared)) %>% 
    print

## # A tibble: 30 × 2
##                Sport adj.r.squared
##                <chr>         <dbl>
## 1          Triathlon     0.8741160
## 2  Modern Pentathlon     0.8707596
## 3   Beach Volleyball     0.8694573
## 4       Canoe Sprint     0.8492032
## 5           Handball     0.8484207
## 6             Hockey     0.8355383
## 7         Basketball     0.8300793
## 8             Tennis     0.8142285
## 9             Rowing     0.8135938
## 10          Swimming     0.8044436
## # ... with 20 more rows

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.

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.

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žítým parametrem je out – do něj se vkládá jméno výstupního souboru.

stargazer(
    VGAMdata::oly12,
    summary = TRUE,
    type = "html"
)

Statistic	N	Mean	St. Dev.	Min	Max
Age	10,384	26.069	5.441	13	71
Height	9,823	1.769	0.113	1.320	2.210
Weight	9,104	72.853	16.067	36	218
Gold	10,384	0.017	0.136	0	2
Silver	10,384	0.017	0.133	0	2
Bronze	10,384	0.018	0.136	0	2
Total	10,384	0.052	0.250	0	5

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:

stargazer(
    oly12,
    summary = TRUE,
    type = "html"
)

Výstup bude vypdat následovně:

Statistic	N	Mean	St. Dev.	Min	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íxh 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"
)

	Dependent variable:
	Weight
	(1)	(2)	(3)
Height	94.898***	79.111***	95.116***
	(1.136)	(1.812)	(1.232)
Age	0.167***	0.164***	0.178***
	(0.020)	(0.019)	(0.020)
SexM	5.573***	-39.600***	5.495***
	(0.256)	(4.066)	(0.266)
Height:SexM		25.757***
		(2.314)
Constant	-103.040***	-76.037***	-108.614***
	(1.987)	(3.127)	(6.081)
Country FE	No	No	No
Observations	9,038	9,038	9,038
R2	0.603	0.608	0.628
Adjusted R2	0.603	0.608	0.619
Residual Std. Error	10.128 (df = 9034)	10.060 (df = 9033)	9.913 (df = 8838)
F Statistic	4,570.826*** (df = 3; 9034)	3,505.755*** (df = 4; 9033)	74.902*** (df = 199; 8838)
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áseldující:

1. Nejprve vytvoříme funkci get.se, která odhadne robustní chyby a vrátí je jako vektor.
2. Pomocí lapply aplikujeme funci get.se na všechny modely v listu Spec_models.
3. Výsledný list vektorů vložíme do parametru se stargazeru.

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

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

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

	Dependent variable:
	Weight
	(1)	(2)	(3)
Height	94.898***	79.111***	95.116***
	(1.212)	(1.680)	(1.359)
Age	0.167***	0.164***	0.178***
	(0.019)	(0.019)	(0.019)
SexM	5.573***	-39.600***	5.495***
	(0.244)	(4.118)	(0.264)
Height:SexM		25.757***
		(2.349)
Constant	-103.040***	-76.037***	-108.614***
	(2.132)	(2.889)	(3.091)
Country FE	No	No	No
Observations	9,038	9,038	9,038
R2	0.603	0.608	0.628
Adjusted R2	0.603	0.608	0.619
Residual Std. Error	10.128 (df = 9034)	10.060 (df = 9033)	9.913 (df = 8838)
F Statistic	4,570.826*** (df = 3; 9034)	3,505.755*** (df = 4; 9033)	74.902*** (df = 199; 8838)
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:

stargazer(
    Sport_models$est_model,
    type = "html"
)

Další estimační funkce

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	sampleSelection::heckit()

…a mnoho, mnoho dalších.

Hint: Funkce pro regresní analýzu jsou roztříštěny do mnoha balíků. O jejich unfikaci se snažil balík AER, který s sebou natáhne většinu užitečných balíků.

Homework

Odhadněte průměrný rozdíl mezi mzdami mužů a žen v ČR.

Velikost hodinové hrubé mzdy \(w_{it}\) respondenta \(i\) je vysvětlena modelem:

\[\log(w_{i}) = \beta_0 + \beta_1log(age_{i}) + \beta_2\log(experience_{i} + 0.1) + \beta_3educ_i + \beta_4gender_i + \varepsilon_i\]

kde \(age_i\) je věk respondenta; \(experience_i\) délka jeho praxe; \(educ_i\) úroveň zdělání a \(gender_i\) je pohlaví respondenta. K hodnotě proměnné \(experience\) je přičteno 0,1, aby bylo možné tuto proměnnou logaritmovat.

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

Máte za úkol:

1. Odhadněte model a výsledný objekt třídy lm uložte do proměnné emodel.

2. Vytvořte tabulku etab (třídy data.frame), která bude ve sloupcích obsahovat:

• jméno proměnné ve sloupci term – pro identifikaci musí jméno proměnné obsahovat jméno sloupce z tabulky data a Intercept pro konstantu – např. SEXgender je v pořádku, pouze SEX nikoliv. Opravný skript by v druhém případě vyhodnotil řešení jako špatné.
• robustní směrodané chyby ve sloupci std.error – použijte White-corrected standardní chyby (tj. HC0)
• tabulka může obsahovat i další sloupce – jejich přítomnost nebude předmětem hodnocení

3. Do proměnné adjR2 uložte adjustovaný koeficient determinace odhadnutého modelu.

4. Do proměnné cor_OF uložte korelační koeficient pozorovaných a vyrovnaných (fitted) hodnot

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