# 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 World Bank World Development Indicators.
V nadväznosti na predchádzajúce cvičenie budeme modelovať očakávanú dĺžku života pri narodení (life_expectancy_at_birth) pre rok 2013. Budeme postupovať od jednoduchej otázky k regresnému modelu:
- Súvisí miera digitalizácie a prístupu k informáciám (podiel používateľov internetu) s očakávanou dĺžkou života?
- Grafické znázornenie tejto 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 world_bank_development_indicators.csv uložím 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)
library(dplyr)
Присоединяю пакет: 'dplyr'Следующий объект скрыт от 'package:kableExtra':
group_rowsСледующие объекты скрыты от 'package:stats':
filter, lagСледующие объекты скрыты от 'package:base':
intersect, setdiff, setequal, unionlibrary(janitor)
Присоединяю пакет: 'janitor'Следующие объекты скрыты от 'package:stats':
chisq.test, fisher.testlibrary(lubridate)
Присоединяю пакет: 'lubridate'Следующие объекты скрыты от 'package:base':
date, intersect, setdiff, unionPri 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 mojom súbore sú stĺpce oddelené čiarkou a desatinné čísla používajú bodku.
cesta <- "world_bank_development_indicators.csv"
if (!file.exists(cesta)) cesta <- file.path("udaje", "world_bank_development_indicators.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)
# Uprava nazvov pre R
names(udaje) <- make.names(trimws(names(udaje)), unique = TRUE)
# Filtrovanie pre rok 2013 a vyber stlpcov
rok <- substr(as.character(udaje$date), 1, 4)
udaje.2013 <- udaje[!is.na(rok) & rok == "2013",
c("country",
"life_expectancy_at_birth",
"individuals_using_internet.",
"government_health_expenditure.",
"GDP_current_US")]
# Premenovanie stlpcov na jednoduchsie nazvy pre model
names(udaje.2013) <- c("country", "life_expectancy", "internet", "health_exp", "gdp")
head(udaje.2013) |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Krajina", "Dĺžka života", "Internet (%)", "Zdravotníctvo (%)", "HDP (USD)"),
caption = "Prvých šesť krajín vo vybraných údajoch") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Krajina | Dĺžka života | Internet (%) | Zdravotníctvo (%) | HDP (USD) |
|---|---|---|---|---|
| Afghanistan | 62.42 | 5.90 | 0.44 | 20146416876 |
| Africa Eastern and Southern | 60.71 | 10.13 | 2.64 | 986342932619 |
| Africa Western and Central | 55.67 | 12.58 | 0.72 | 834097010135 |
| Albania | 78.12 | 51.80 | 2.58 | 12776224497 |
| Algeria | 74.61 | 22.50 | 4.29 | 209755003251 |
| American Samoa | NA | NA | NA | 638000000 |
Jeden riadok predstavuje jednu krajinu v roku 2013. 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 pri narodení (v rokoch); vysvetľovaná premenná. |
| internet | Podiel obyvateľov používajúcich internet (v %); kľúčová vysvetľujúca premenná. |
| health_exp | Vládne výdavky na zdravotníctvo (v % z celkových výdavkov). |
| gdp | HDP v bežných cenách (v USD). |
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.2013),
Pocet_NA = colSums(is.na(udaje.2013))
)
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 | 16 |
| internet | 25 |
| health_exp | 39 |
| gdp | 15 |
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.
# Doplnenie mediánom pre vysvetľujúce premenné
for (premenna in c("internet", "health_exp", "gdp")) {
median_premennej <- median(udaje.2013[[premenna]], na.rm = TRUE)
if (!is.finite(median_premennej)) stop("Premennú nemožno doplniť mediánom.")
chyba <- is.na(udaje.2013[[premenna]])
udaje.2013[[premenna]][chyba] <- median_premennej
}
premenne <- c("life_expectancy", "internet", "health_exp", "gdp")
opis <- data.frame(
Premenna = c("Dĺžka života (roky)", "Internet (%)", "Zdravotníctvo (%)", "HDP (USD)"),
Minimum = sapply(udaje.2013[premenne], min),
Median = sapply(udaje.2013[premenne], median),
Priemer = sapply(udaje.2013[premenne], mean),
Maximum = sapply(udaje.2013[premenne], max)
)
opis |>
kable(format = "html", digits = 2, row.names = FALSE,
caption = "Opisná štatistika po doplnení chýbajúcich hodnôt") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Premenna | Minimum | Median | Priemer | Maximum |
|---|---|---|---|---|
| Dĺžka života (roky) | NA | NA | NA | NA |
| Internet (%) | 0.00 | 3.826000e+01 | 4.03800e+01 | 9.655000e+01 |
| Zdravotníctvo (%) | 0.34 | 2.670000e+00 | 3.18000e+00 | 9.890000e+00 |
| HDP (USD) | 38615890.63 | 5.094967e+10 | 2.34943e+12 | 7.771472e+13 |
3 Graf závislosti
Formulujeme pracovnú hypotézu: krajiny s vyššou mierou penetrácie internetu (ktorá reflektuje celkový rozvoj infraštruktúry a dostupnosť informácií) vykazujú v priemere vyššiu očakávanú dĺžku života.
graf_internet <- ggplot(udaje.2013, aes(x = internet, y = life_expectancy)) +
geom_point(color = "steelblue", size = 2, alpha = 0.75) +
labs(x = "Používatelia internetu (% populácie)",
y = "Očakávaná dĺžka života (roky)") +
theme_minimal(base_size = 13)
graf_internetWarning: Removed 16 rows containing missing values or values outside the scale range
(`geom_point()`).
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 ~ internet, data = udaje.2013)
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 miera používania internetu") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 61.510 |
| internet | 0.242 |
Po odhade získavame rovnicu:
\[ \widehat{LifeExpectancy}_i = 61.51 + 0.242 \cdot Internet_i \].
Smernica (\(\beta_1\)): Nárast podielu používateľov internetu o 1 percentuálny bod sa spája s priemerným nárastom predikovanej dĺžky života o približne r round(coef(model1)[2], 2) roka.
Konštanta (\(\beta_0\)): V hypotetickej krajine s 0 % používateľov internetu by model predikoval dĺžku života približne r round(coef(model1)[1], 2) roka.
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_internet +
geom_abline(intercept = coef(model1)["(Intercept)"],
slope = coef(model1)["internet"],
color = "firebrick", linewidth = 1)Warning: Removed 16 rows containing missing values or values outside the scale range
(`geom_point()`).
5 Viacnásobná lineárna regresia
Dĺžka života môže súvisieť s viacerými charakteristikami súčasne. Rozšírime model o vládne výdavky na zdravotníctvo (health_exp) a HDP (gdp):
\[LifeExpectancy_i = \beta_0 + \beta_1 Internet_i + \beta_2 HealthExp_i + \beta_3 GDP_i + u_i.\]
Keďže celkový HDP v dolároch dosahuje vysoké hodnoty, pre lepšiu čitateľnosť koeficientu ho pred odhadom vyjadríme v miliardách USD (vydelíme \(10^9\)).
# Prevod HDP na miliardy USD pre zrozumitelnejsie koeficienty
udaje.2013$gdp_mld <- udaje.2013$gdp / 1e9
model <- lm(life_expectancy ~ internet + health_exp + gdp_mld,
data = udaje.2013)
vysledky <- tidy(model)
vysledky[c("term", "estimate")] |>
kable(format = "html", digits = 6, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Viacnásobná regresia: internet, zdravotníctvo a HDP (v mld. USD)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 61.044126 |
| internet | 0.223663 |
| health_exp | 0.356779 |
| gdp_mld | 0.000026 |
V tabuľke ponechávame šesť desatinných miest, aby sa malý koeficient pri HDP nezaokrúhlil na nulu. Pri použitom postupe dostaneme približne:
\[ \widehat{LifeExpectancy}_i = 61.044 + 0.224\,Internet_i + 0.357\,HealthExp_i + 0.000026\,GDP\_mld_i. \]
Koeficient internet sa teraz interpretuje pri rovnakých hodnotách HDP a výdavkov na zdravotníctvo. Rozdiel jedného percentuálneho bodu v podiele používateľov internetu sa spája s rozdielom približne 0.22 roka v predikovanej dĺžke života.
Koeficient gdp_mld je číselne malý, pretože sa vzťahuje na jednu miliardu USD celkového HDP krajiny. Rozdiel 1 000 miliárd (teda 1 bilióna) USD v HDP sa pri nezmenených ostatných premenných spája s rozdielom približne:
\[ 1000 \times 0.000026 = 0.0263 \]
roka v predikovanej dĺžke života.
Všimnite si: Koeficient pri premennej internet sa po pridaní ďalších premenných mierne zmenil. Jednoduchý a viacnásobný model odpovedajú na odlišné otázky: druhý zohľadňuje aj ďalšie zahrnuté charakteristiky krajín.
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) |
|---|---|---|
| 258 | 0.693 | 4.459 |
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.693. Model teda zachytáva približne 69.3 % výberovej variability dĺžky života.
Ú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 vypočítame rezíduum:
\[ e_i = y_i - \widehat y_i. \]
pouzite_riadky <- names(fitted(model))
udaje_modelu <- udaje.2013[pouzite_riadky, ]
porovnanie <- data.frame(
Krajina = udaje_modelu$country,
Skutocnost = udaje_modelu$life_expectancy,
Predikcia = as.numeric(fitted(model)),
Reziduum = as.numeric(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 | 62.42 | 62.52 | -0.11 |
| Africa Eastern and Southern | 60.71 | 64.28 | -3.57 |
| Africa Western and Central | 55.67 | 64.14 | -8.46 |
| Albania | 78.12 | 73.55 | 4.57 |
| Algeria | 74.61 | 67.61 | 7.00 |
| Angola | 59.31 | 64.56 | -5.25 |
| Antigua and Barbuda | 77.58 | 76.33 | 1.26 |
| Arab World | 70.20 | 69.60 | 0.60 |
Overenie pre prvú krajinu v tabuľke (Afghanistan):
\[ 62.42 - 62.52 = -0.11 \quad (\text{Reziduum: } -0.11) \]
8 Predikcia pre zadanú hodnotu
Využijeme jednoduchý model model1 a predikujeme očakávanú dĺžku života pre krajiny s hypotetickou penetráciou internetu na úrovni 20 %, 50 % a 80 %:
nove_udaje <- data.frame(internet = c(20, 50, 80))
nove_udaje$Predikcia <- predict(model1, newdata = nove_udaje)
nove_udaje |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Používatelia internetu (%)", "Predikovaná dĺžka života (roky)"),
caption = "Bodové predikcie jednoduchého regresného modelu") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Používatelia internetu (%) | Predikovaná dĺžka života (roky) |
|---|---|
| 20 | 66.35 |
| 50 | 73.61 |
| 80 | 80.87 |
9 Malé cvičenia
Ktorá premenná je v našom modeli vysvetľovaná a ktoré sú vysvetľujúce? • Vysvetľovaná premenná (\(Y\), závislá): life_expectancy (očakávaná dĺžka života pri narodení v rokoch). • Vysvetľujúce premenné (\(X\), nezávislé): internet (% používateľov internetu), health_exp (% vládnych výdavkov na zdravotníctvo), gdp (HDP v USD).
Odhadnite jednoduchý model
Life.expectancy ~ GDP. Pomocoutidy()získajte koeficienty, zobrazte ich pomocoukable()akable_styling()a pomocouggplot2nakreslite bodový graf s regresnou priamkou.
model_gdp <- lm(life_expectancy ~ gdp, data = udaje.2013)
tidy(model_gdp)[c("term", "estimate")] |>
kable(format = "html", digits = 10, row.names = FALSE,
col.names = c("Člen", "Koeficient")) |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen | Koeficient |
|---|---|
| (Intercept) | 70.84932 |
| gdp | 0.00000 |
ggplot(udaje.2013, aes(x = gdp, y = life_expectancy)) +
geom_point(color = "darkgreen", alpha = 0.6) +
geom_abline(intercept = coef(model_gdp)[1], slope = coef(model_gdp)[2], color = "red") +
labs(x = "HDP (USD)", y = "Dĺžka života (roky)") +
theme_minimal()Warning: Removed 16 rows containing missing values or values outside the scale range
(`geom_point()`).
Vysvetlite koeficient Schooling v jednoduchom aj viacnásobnom modeli. Prečo pri viacnásobnom modeli dopĺňame „pri rovnakých ostatných premenných“? V jednoduchom modeli zachytáva koeficient celkovú asociáciu medzi internetom a dĺžkou života (vrátane vplyvu vynechaných premenných). Vo viacnásobnom modeli klauzula ceteris paribus („pri nezmenených ostatných premenných“) izoluje vplyv internetu od vplyvu HDP a výdavkov na zdravotníctvo.
Pre Internet = 60 % vypočítajte predikciu pomocou
model1a overte ju ručným dosadením do odhadnutej rovnice:
# 1. Vypocet pomocou funkcie predict()
pred_funkcia <- predict(model1, newdata = data.frame(internet = 60))
# 2. Rucne dosadenie do rovnice
b0 <- coef(model1)[1]
b1 <- coef(model1)[2]
pred_rovnica <- b0 + b1 * 60- Pomocou funkcie
predict(model1): 76.03 roka. - Ručným dosadením do rovnice: \[ \widehat{LifeExpectancy} = 61.51 + 0.242 \times 60 = 76.03 \text{ roka.} \]
Oba postupy dávajú rovnaký výsledok.
Nájdite v tabuľke
porovnaniejednu krajinu s kladným a jednu so záporným rezíduom. Vysvetlite význam znamienka. • Kladné rezíduum (\(e_i > 0\)): Skutočná dĺžka života v danej krajine je vyššia, než model očakával na základe jej parametrov. • Záporné rezíduum (\(e_i < 0\)): Skutočná dĺžka života zaostáva za predikciou modelu.Vysvetlite, prečo z kladného koeficienta nemožno automaticky vyvodiť kauzálny účinok. Kladný koeficient nedokazuje, že surfovanie na internete priamo predlžuje život. Penetrácia internetu koreluje s celkovou infraštruktúrou, hygienou, vzdelaním a prístupom k zdravotnej starostlivosti. Ide o štatistickú asociáciu, nie o priamu kauzalitu.