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

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

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

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

18.5.2 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ášť.

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:

Sport_models %>% print
## # A tibble: 26 x 2
## # Rowwise: 
##    Sport      est_model
##    <chr>      <list>   
##  1 Archery    <lm>     
##  2 Athletics  <lm>     
##  3 Badminton  <lm>     
##  4 Basketball <lm>     
##  5 Beach      <lm>     
##  6 Canoe      <lm>     
##  7 Cycling    <lm>     
##  8 Diving     <lm>     
##  9 Equestrian <lm>     
## 10 Fencing    <lm>     
## # … 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:

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
##    <chr> <list>        <dbl>         <dbl> <dbl>     <dbl>     <dbl> <int>
##  1 Tria… <lm>          0.879         0.874  2.98      179. 4.32e- 44     5
##  2 Mode… <lm>          0.878         0.870  3.81      110. 3.90e- 27     5
##  3 Beach <lm>          0.874         0.869  4.56      153. 9.52e- 39     5
##  4 Hand… <lm>          0.850         0.848  6.01      423. 9.22e-122     5
##  5 Hock… <lm>          0.837         0.835  4.29      485. 2.43e-147     5
##  6 Bask… <lm>          0.833         0.830  6.50      315. 6.52e- 97     5
##  7 Tenn… <lm>          0.819         0.815  4.49      171. 4.96e- 55     5
##  8 Rowi… <lm>          0.815         0.813  5.84      557. 5.13e-184     5
##  9 Swim… <lm>          0.806         0.805  5.04      887. 2.58e-302     5
## 10 Canoe <lm>          0.797         0.794  5.08      293. 3.55e-102     5
## # … with 16 more rows, and 5 more variables: logLik <dbl>, AIC <dbl>,
## #   BIC <dbl>, deviance <dbl>, df.residual <int>

Je vidět, že nejvyšší \(\bar{R}^2\) má model odhadnutý na vzorku triatlonistů.