# Spustite raz v konzole, iba ak tieto balíky chýbajú.
install.packages(c("knitr", "rmarkdown", "broom", "kableExtra", "ggplot2"))Econometrics in R – cvičenie 4
Prvé zoznámenie s lineárnou regresiou
1 Predmet výskumu
V tomto cvičení si prvýkrát ukážeme, ako pomocou R odhadnúť a interpretovať lineárnu regresiu. Použijeme údaje z databázy Life Expectancy (WHO).
Budeme postupovať od jednoduchej otázky k regresnému modelu:
- Súvisí vyššia hodnota ukazovateľa vzdelania s vyššou očakávanou dĺžkou života?
- Grafické znázornenie tejt závislosti.
- Odhad regresnej priamky a jej interpretácia.
- Pridanie ďalších vysvetľujúcich premenných.
- Predikcia.
V tomto cvičení sa nevenujeme diagnostickým testom. Cieľom je pochopiť základnú myšlienku regresie a prácu s funkciou
lm().
2 Príprava údajov
2.1 Umiestnenie súboru
Súbor Life_Expectancy_Data.csv uložte do podpriečinka udaje. Kód umožňuje aj uloženie CSV priamo k dokumentu. Ak existujú obe kópie, použije sa kópia v udaje.
V RStudio alebo Posit Cloud otvorte .qmd súbor a stlačte Render. Modely odhadneme funkciou lm(). Výsledky prevedieme na tabuľky pomocou balíka broom a zobrazíme pomocou známych funkcií kable() a kable_styling(). Na renderovanie dokumentu je potrebné prostredie Quarto s balíkmi knitr a rmarkdown, ktoré môžu byť už nainštalované.
library(broom)
library(knitr)
library(kableExtra)
library(ggplot2)Pri formátovaní použijeme zápis |>, ktorý odovzdá výsledok predchádzajúceho príkazu ako prvý argument nasledujúcej funkcie. Napríklad tabulka |> kable() zobrazí objekt tabulka pomocou kable().
Nepotrebujeme v dokumente nastavovať pevnú absolútnu cestu pomocou setwd(). Pri renderovaní samostatného dokumentu sa relatívne cesty štandardne odvodzujú od jeho umiestnenia; projektové nastavenia môžu toto správanie zmeniť.
2.2 Načítanie CSV
V našom súbore sú stĺpce oddelené čiarkou a desatinné čísla používajú bodku.
cesta <- "udaje/Life_Expectancy_Data.csv"
if (!file.exists(cesta)) cesta <- "Life_Expectancy_Data.csv"
if (!file.exists(cesta)) {
stop("Vložte CSV do podpriečinka udaje alebo vedľa QMD dokumentu.")
}
udaje <- read.csv(cesta, sep = ",", dec = ".",
na.strings = c("", "NA"),
check.names = FALSE, stringsAsFactors = FALSE)
# Odstránime medzery na začiatku a konci názvov stĺpcov.
# Názvy upravíme do podoby vhodnej pre R.
names(udaje) <- make.names(trimws(names(udaje)), unique = TRUE)
# Vyberieme rok 2015 a potrebné stĺpce.
udaje.2015 <- udaje[!is.na(udaje$Year) & udaje$Year == 2015,
c("Country", "Life.expectancy", "BMI", "GDP", "Schooling")]
head(udaje.2015) |>
kable(format = "html", digits = 2, row.names = FALSE,
caption = "Prvých šesť krajín vo vybraných údajoch") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Country | Life.expectancy | BMI | GDP | Schooling |
|---|---|---|---|---|
| Afghanistan | 65.0 | 19.1 | 584.26 | 10.1 |
| Albania | 77.8 | 58.0 | 3954.23 | 14.2 |
| Algeria | 75.6 | 59.5 | 4132.76 | 14.4 |
| Angola | 52.4 | 23.3 | 3695.79 | 11.4 |
| Antigua and Barbuda | 76.4 | 47.7 | 13566.95 | 13.9 |
| Argentina | 76.3 | 62.8 | 13467.12 | 17.3 |
Jeden riadok predstavuje jednu krajinu v roku 2015. Nemodelujeme vývoj v čase, takejto organizácii údajov hovoríme prierezové údaje (cross-sectional data)
| Premenná | Význam |
|---|---|
| Country | Názov krajiny. |
| Life.expectancy | Očakávaná dĺžka života v rokoch; vysvetľovaná premenná. |
| Schooling | Ukazovateľ vzdelania vyjadrený v rokoch. |
| GDP | Ukazovateľ HDP na obyvateľa podľa označenia databázy. |
| BMI | Ukazovateľ označený ako BMI; definíciu a kvalitu treba preveriť. |
2.3 Chýbajúce hodnoty
Chýbajúcu hodnotu R označuje ako NA. Najprv zistíme ich počty:
chybajuce <- data.frame(
Premenna = names(udaje.2015),
Pocet_NA = colSums(is.na(udaje.2015))
)
chybajuce |>
kable(format = "html", col.names = c("Premenná", "Počet chýbajúcich hodnôt"),
row.names = FALSE, caption = "Chýbajúce hodnoty pred doplnením") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Premenná | Počet chýbajúcich hodnôt |
|---|---|
| Country | 0 |
| Life.expectancy | 0 |
| BMI | 2 |
| GDP | 29 |
| Schooling | 10 |
Pri priložených údajoch máme 183 krajín. Chýbajú 2 hodnoty BMI, 29 hodnôt GDP a 10 hodnôt Schooling. Dĺžka života v tomto roku nechýba.
Aby sme sa mohli sústrediť na regresiu, ponecháme jednoduchý postup z pôvodného textu: chýbajúce hodnoty vysvetľujúcich premenných doplníme mediánom. Je to didaktické zjednodušenie, ktoré môže ovplyvniť výsledky. Chýbajúcu vysvetľovanú premennú takto nenahrádzame.
# Zachováme pôvodné údaje pre prípadnú neskoršiu prácu.
udaje.2015.povodne <- udaje.2015
# Ak by chýbala vysvetľovaná premenná, taký riadok nepoužijeme.
udaje.2015 <- udaje.2015[!is.na(udaje.2015$Life.expectancy), ]
for (premenna in c("BMI", "GDP", "Schooling")) {
median_premennej <- median(udaje.2015[[premenna]], na.rm = TRUE)
if (!is.finite(median_premennej)) stop("Premennú nemožno doplniť mediánom.")
chyba <- is.na(udaje.2015[[premenna]])
udaje.2015[[premenna]][chyba] <- median_premennej
}
premenne <- c("Life.expectancy", "BMI", "GDP", "Schooling")
opis <- data.frame(
Premenna = premenne,
Minimum = sapply(udaje.2015[premenne], min),
Median = sapply(udaje.2015[premenne], median),
Priemer = sapply(udaje.2015[premenne], mean),
Maximum = sapply(udaje.2015[premenne], max)
)
opis |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Premenná", "Minimum", "Medián", "Priemer", "Maximum"),
caption = "Opisná štatistika po doplnení chýbajúcich hodnôt") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Premenná | Minimum | Medián | Priemer | Maximum |
|---|---|---|---|---|
| Life.expectancy | 51.00 | 73.90 | 71.62 | 88.00 |
| BMI | 2.50 | 48.60 | 42.82 | 77.60 |
| GDP | 33.68 | 2916.23 | 6508.81 | 66346.52 |
| Schooling | 4.90 | 13.10 | 12.94 | 20.40 |
3 Graf
Očakávame, že krajiny s vyššou hodnotou ukazovateľa Schooling majú spravidla aj vyššiu očakávanú dĺžku života. Ide zatiaľ o pracovnú hypotézu.
graf_vzdelanie <- ggplot(udaje.2015,
aes(x = Schooling, y = Life.expectancy)) +
geom_point(color = "steelblue", size = 2, alpha = 0.75) +
labs(x = "Schooling (roky)",
y = "Očakávaná dĺžka života (roky)") +
theme_minimal(base_size = 13)
graf_vzdelanie
Otázky: Vidíte rastúcu alebo klesajúcu závislosť? Majú všetky krajiny s podobným Schooling aj rovnakú dĺžku života?
4 Jednoduchá lineárna regresia
4.1 Zápis modelu
Najprv použijeme iba jednu vysvetľujúcu premennú:
\[ LifeExpectancy_i = \beta_0 + \beta_1 Schooling_i + u_i. \]
- \(i\) označuje krajinu.
- \(\beta_0\) je konštanta, teda priesečník regresnej priamky s osou Y.
- \(\beta_1\) je smernica priamky: vyjadruje rozdiel predikovanej dĺžky života pri jednotkovom rozdiele Schooling.
- \(u_i\) zachytáva tú časť dĺžky života, ktorú tento jednoduchý model nevysvetľuje.
Regresná priamka nemusí prechádzať všetkými bodmi. Metóda najmenších štvorcov vyberie takú priamku, ktorá minimalizuje súčet štvorcov vertikálnych rozdielov medzi bodmi a priamkou.
4.2 Odhad pomocou lm()
model1 <- lm(Life.expectancy ~ Schooling, data = udaje.2015)
vysledky1 <- tidy(model1)
vysledky1[c("term", "estimate")] |>
kable(format = "html", digits = 3, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Jednoduchá regresia: dĺžka života a Schooling") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 42.816 |
| Schooling | 2.226 |
Funkcia tidy(model1) z balíka broom prevedie výsledky modelu na tabuľku: jeden riadok zodpovedá jednému koeficientu. Stĺpec term obsahuje názov člena modelu a estimate jeho odhad. Označenie (Intercept) znamená konštantu. Na úvod zobrazujeme iba tieto dva stĺpce.
tidy() pripravuje údaje o výsledkoch; kable() vytvára tabuľku a kable_styling() upravuje jej vzhľad. Odhad modelu sa formátovaním nemení.
Znak ~ čítame ako „vysvetľujeme pomocou“. Vľavo je vysvetľovaná premenná a vpravo vysvetľujúca premenná. Konštantu R zahrnie automaticky.
Pri priložených údajoch po uvedenom doplnení hodnôt získame približne:
\[ \widehat{LifeExpectancy}_i = 42{,}816 + 2{,}226\,Schooling_i. \]
Interpretácia smernice: Rozdiel jedného roka v ukazovateli Schooling sa v tomto jednoduchom modeli spája s rozdielom približne 2,23 roka v predikovanej dĺžke života.
Interpretácia konštanty: Pri Schooling = 0 by priamka predikovala približne 42,82 roka. Nula však leží mimo pozorovaného rozsahu Schooling v tejto vzorke, preto konštantu chápeme predovšetkým ako parameter polohy priamky.
Poznámka
Kladná závislosť medzi krajinami sama osebe nedokazuje, že ďalší rok vzdelania konkrétnemu človeku predĺži život o 2,23 roka. Krajiny sa líšia aj v mnohých ďalších charakteristikách.
4.3 Regresná priamka v grafe
K bodovému grafu pridáme priamku pomocou geom_abline(). Jej konštantu (intercept) a smernicu (slope) prevezmeme priamo z odhadnutého modelu model1.
graf_vzdelanie +
geom_abline(intercept = coef(model1)["(Intercept)"],
slope = coef(model1)["Schooling"],
color = "firebrick", linewidth = 1)
5 Viacnásobná lineárna regresia
Dĺžka života môže súvisieť s viacerými charakteristikami súčasne. Rozšírime model o GDP a BMI, ako v pôvodnom zadaní:
\[ LifeExpectancy_i = \beta_0 + \beta_1 BMI_i + \beta_2 GDP_i + \beta_3 Schooling_i + u_i. \]
model <- lm(Life.expectancy ~ BMI + GDP + Schooling,
data = udaje.2015)
vysledky <- tidy(model)
vysledky[c("term", "estimate")] |>
kable(format = "html", digits = 8, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Viacnásobná regresia: BMI, GDP a Schooling") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 44.19364375 |
| BMI | 0.04947912 |
| GDP | 0.00006972 |
| Schooling | 1.92097373 |
V tabuľke ponechávame osem desatinných miest, aby sa malý koeficient pri GDP nezaokrúhlil na nulu. Pri použitom postupe dostaneme približne:
\[ \widehat{LifeExpectancy}_i = 44{,}194 + 0{,}0495\,BMI_i + 0{,}00006972\,GDP_i + 1{,}921\,Schooling_i. \]
Koeficient Schooling sa teraz interpretuje pri rovnakých hodnotách GDP a BMI. Rozdiel jedného roka Schooling sa spája s rozdielom približne 1,92 roka v predikovanej dĺžke života.
Koeficient GDP je číselne malý aj preto, že sa vzťahuje na jednu jednotku GDP. Rozdiel 1 000 jednotiek GDP sa pri nezmenených ostatných premenných spája s rozdielom približne
\[ 1000\times 0{,}00006972 = 0{,}06972 \]
roka v predikovanej dĺžke života. Pred ekonomickou interpretáciou treba overiť presné jednotky GDP v zdroji.
Všimnite si: Koeficient Schooling sa po pridaní ďalších premenných zmenil. Jednoduchý a viacnásobný model odpovedajú na odlišné otázky: druhý zohľadňuje aj ďalšie zahrnuté charakteristiky krajín.
BMI zatiaľ ponechávame na ukážku syntaxe viacnásobnej regresie. V tomto súbore obsahuje hodnoty, ktorých význam treba overiť. Jeho kladné znamienko preto neinterpretujeme ako dôkaz zdravotného prínosu vyššieho BMI a premennú nevyraďujeme iba preto, že znamienko nezodpovedá očakávaniu.
6 Stručný súhrn modelu pomocou glance()
Zatiaľ čo tidy() poskytuje riadok pre každý koeficient, glance() poskytuje jeden riadok so súhrnnými charakteristikami celého modelu. Na prvom cvičení vyberieme iba počet pozorovaní, \(R^2\) a odhad štandardnej odchýlky chyby.
suhrn <- glance(model)
suhrn[c("nobs", "r.squared", "sigma")] |>
kable(format = "html", digits = c(0, 3, 3), row.names = FALSE,
col.names = c("Počet krajín", "R²", "Odhad SD chyby (roky)"),
caption = "Základné charakteristiky viacnásobného modelu") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Počet krajín | R² | Odhad SD chyby (roky) |
|---|---|---|
| 183 | 0.622 | 5.033 |
nobsje počet krajín použitých pri odhade.r.squaredje podiel výberovej variability dĺžky života zachytený modelom.sigmaje odhad štandardnej odchýlky modelovej chyby, tu vyjadrený v rokoch.
Hodnota \(R^2\) je v našom viacnásobnom modeli približne 0,622. Model teda zachytáva približne 62,2 % výberovej variability dĺžky života. Neznamená to, že „62,2 % predikcií je správnych“, ani to nedokazuje kauzalitu.
Úplný textový súhrn je stále dostupný príkazom summary(model). tidy(model) okrem názvov a odhadov obsahuje aj std.error, statistic a p.value. Ich vysvetlenie ponecháme na hodinu venovanú štatistickému testovaniu. Preto ich v úvodných tabuľkách nezobrazujeme.
7 Skutočná hodnota, predikcia a rezíduum
Pre každú krajinu môžeme porovnať skutočnú dĺžku života s hodnotou na regresnej ploche:
\[ e_i = y_i - \widehat y_i. \]
Tento rozdiel nazývame rezíduum. Kladné rezíduum znamená, že skutočná hodnota je vyššia než predikcia; záporné znamená opak. Rezíduum vieme vypočítať, zatiaľ čo skutočnú modelovú chybu \(u_i\) nepoznáme.
porovnanie <- data.frame(
Krajina = udaje.2015$Country,
Skutocnost = udaje.2015$Life.expectancy,
Predikcia = fitted(model),
Reziduum = residuals(model)
)
head(porovnanie, 8) |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Krajina", "Skutočnosť", "Predikcia", "Rezíduum"),
caption = "Skutočné a predikované hodnoty dĺžky života v rokoch") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Krajina | Skutočnosť | Predikcia | Rezíduum |
|---|---|---|---|
| Afghanistan | 65.0 | 64.58 | 0.42 |
| Albania | 77.8 | 74.62 | 3.18 |
| Algeria | 75.6 | 75.09 | 0.51 |
| Angola | 52.4 | 67.50 | -15.10 |
| Antigua and Barbuda | 76.4 | 74.20 | 2.20 |
| Argentina | 76.3 | 81.47 | -5.17 |
| Armenia | 74.8 | 71.33 | 3.47 |
| Australia | 82.8 | 90.62 | -7.82 |
Úloha: Vyberte jeden riadok a ručne overte, že Skutocnost - Predikcia = Reziduum.
Funkcia coef(model) vracia koeficienty, fitted(model) vyrovnané hodnoty a residuals(model) rezíduá. Ak budeme neskôr potrebovať maticu vysvetľujúcich premenných, získame ju príkazom model.matrix(model).
8 Predikcia pre zadanú hodnotu
Na prvú ukážku predikcie sa vrátime k jednoduchému modelu so Schooling. Vypočítame, akú dĺžku života priamka predikuje pre hodnoty 10, 12 a 14 rokov:
nove_udaje <- data.frame(Schooling = c(10, 12, 14))
nove_udaje$Predikcia <- predict(model1, newdata = nove_udaje)
nove_udaje |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Schooling (roky)", "Predikcia dĺžky života (roky)"),
caption = "Predikcie jednoduchého regresného modelu") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Schooling (roky) | Predikcia dĺžky života (roky) |
|---|---|
| 10 | 65.08 |
| 12 | 69.53 |
| 14 | 73.98 |
Ide o bodové predikcie, nie o isté výsledky. Skutočná krajina môže mať inú dĺžku života. O neistote predikcie a predikčných intervaloch sa budeme učiť neskôr.
9 Malé cvičenia
- Ktorá premenná je v našom modeli vysvetľovaná a ktoré sú vysvetľujúce?
- Odhadnite jednoduchý model
Life.expectancy ~ GDP. Pomocoutidy()získajte koeficienty, zobrazte ich pomocoukable()akable_styling()a pomocouggplot2nakreslite bodový graf s regresnou priamkou. - Vysvetlite koeficient Schooling v jednoduchom aj viacnásobnom modeli. Prečo pri viacnásobnom modeli dopĺňame „pri rovnakých ostatných premenných“?
- Pomocou
model1vypočítajte predikciu pre Schooling = 13. Overte ju ručným dosadením do odhadnutej rovnice. - Nájdite v tabuľke
porovnaniejednu krajinu s kladným a jednu so záporným rezíduom. Vysvetlite význam znamienka. - Vysvetlite, prečo z kladného koeficienta nemožno automaticky vyvodiť kauzálny účinok.