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:
- Vytvořit si
listse všemi rovnicemi (specifikacemi) - Odhadnout model pro všechny specifikace pomocí funkce
map()z balíku purrr
Nejprve vytvoříme různé modelové specifikace s pomocí funkce update():
Models <- 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))
)Nyní pomocí purrr::map() zavoláme lm() na všechny prvky listu Models – tedy na všechny modelové specifikace:
Spec_models <- purrr::map(Models, lm, data = VGAMdata::oly12)
Spec_models <- purrr::map(Models, function(x) lm(x, data = VGAMdata::oly12))
Spec_models <- purrr::map(Models, ~lm(., data = VGAMdata::oly12))Všechny tři výše uvedené způsoby použití map() 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 se seznamem 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 = ., 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-1618.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 podobnou operaci existuje více přístupů. Historicky byla pro takové operace určena funkce dplyr::do(). Ta však bude v budoucnosti v balíku dplyr nahrazena a není ještě zcela jasné jak.
Dobrou variantou je postup analogický k odhadu více modelů na stejných datech. Nyní chceme odhadnout stejný model na různých datech. V předchozím případě jsme přes list rovnic iterovali pomocí funkce map(). Nyní bychom chtěli iterovat přes list tibblů. Ten vytvoříme pomocí funkce split() z base R.
Vstupy split() bude úplný tibble a pravostranná formule se jménem vektoru, podle kterého se má tento tibble rozdělit na menší tibblíky (bylo by možné použít přímo i tento vektor ve tvaru .$Sport). Typicky je tento vektor přímo součástí vstupního tibble. To není problém – můžeme využí placeholder .. Dejme tomu, že chceme odhadnout parametry pro každý sport zvlášť:
Sport_models <- oly12 %>%
split(~ Sport) %>%
map(
# Zde si pomáháme anonymní funkcí, která nám pomůže doručit data
# do správného parametru funkce lm()
function(x) lm(Models[[1]], data = x)
)Výsledná proměnná je list, který obsahuje odhadnuté modely:
Sport_models[1:3] %>% print()## $Archery
##
## Call:
## lm(formula = Models[[1]], data = x)
##
## Coefficients:
## (Intercept) Height Age I(Age^2) SexM
## -45.146045 56.714558 0.459166 0.002153 10.191336
##
##
## $Athletics
##
## Call:
## lm(formula = Models[[1]], data = x)
##
## Coefficients:
## (Intercept) Height Age I(Age^2) SexM
## -145.12558 123.04350 -0.40467 0.01024 2.15605
##
##
## $Badminton
##
## Call:
## lm(formula = Models[[1]], data = x)
##
## Coefficients:
## (Intercept) Height Age I(Age^2) SexM
## -70.285717 75.702484 0.093450 0.001788 3.787492S modely uloženými v seznamu. Následující kód například opět za využívá map(). Pomocí broom::glance() vytvoří tabulku modelů seřazenou sestupně podle adjustovaného koeficientu determinace:
Sport_models %>%
# Varinata map_dfr() je ekvivalentní k map() %>% bind_rows()
map_dfr(glance) %>%
arrange(desc(adj.r.squared)) %>%
print()## # A tibble: 26 × 12
## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.879 0.874 2.98 179. 4.32e- 44 4 -256. 524. 540.
## 2 0.878 0.870 3.81 110. 3.90e- 27 4 -179. 371. 384.
## 3 0.874 0.869 4.56 153. 9.52e- 39 4 -271. 553. 568.
## 4 0.850 0.848 6.01 423. 9.22e-122 4 -974. 1960. 1983.
## 5 0.837 0.835 4.29 485. 2.43e-147 4 -1099. 2209. 2233.
## 6 0.833 0.830 6.50 315. 6.52e- 97 4 -846. 1705. 1726.
## 7 0.819 0.815 4.49 171. 4.96e- 55 4 -453. 918. 936.
## 8 0.815 0.813 5.84 557. 5.13e-184 4 -1627. 3267. 3292.
## 9 0.806 0.805 5.04 887. 2.58e-302 4 -2605. 5223. 5251.
## 10 0.797 0.794 5.08 293. 3.55e-102 4 -923. 1858. 1880.
## # … with 16 more rows, and 3 more variables: deviance <dbl>, df.residual <int>,
## # nobs <int>(Pozn. takové řazení je užitečné jenom ve velmi omezeném počtu případů.)