| Premenná | Význam | Jednotka |
|---|---|---|
| diabetes | Výskyt cukrovky v populácii; vysvetľovaná premenná. | % |
| spotreba_na_obyv | Ročná spotreba cukru na obyvateľa. | kg / osoba / rok |
| obezita | Miera obezity. | % |
| hdp_tis | HDP na obyvateľa v tisícoch USD. | tis. USD |
| urbanizacia | Miera urbanizácie. | % |
Spotreba cukru a výskyt cukrovky
Lineárna regresia
1 Predmet výskumu
Pri nasledujúcom zadaní vychádzam z otázky, ktorú som si položila už pri opise údajov:
- Súvisí vyššia spotreba cukru na obyvateľa s vyšším výskytom cukrovky?
2 Údaje
Pracujem so súborom sugar_consumption_dataset.csv, ktorý som vyčistila v minulom zadaní. Jeden riadok je jeden záznam (krajina + rok). Na jednu dvojicu krajina + rok pripadá viac záznamov s odlišnými hodnotami, preto nevyberám jeden rok ako vo vzore. Všetkých 10000 záznamov beriem ako samostatné pozorovania, teda ako združené prierezové údaje (pooled cross-section). `r if (pocet_na == 0)
| Premenná | Minimum | Medián | Priemer | Maximum |
|---|---|---|---|---|
| diabetes | 2.14 | 9.33 | 9.35 | 16.55 |
| spotreba_na_obyv | 5.02 | 37.88 | 37.63 | 69.99 |
| obezita | 5.00 | 22.26 | 22.36 | 40.00 |
| hdp_tis | 0.50 | 25.16 | 25.24 | 50.00 |
| urbanizacia | 10.02 | 50.20 | 49.97 | 90.00 |
3 Graf
Hypotéza: Očakávam, že záznamy s vyššou spotrebou cukru majú spravidla aj vyšší
graf_cukor <- ggplot(udaje.model,
aes(x = spotreba_na_obyv, y = diabetes)) +
geom_point(color = farba, size = 1, alpha = 0.15) +
labs(x = "Spotreba cukru (kg na osobu za rok)",
y = "Výskyt cukrovky (%)") +
theme_minimal(base_size = 13)
graf_cukor
Odpovede na otázky zo vzoru:
- Je závislosť rastúca, alebo klesajúca? Korelačný koeficient je 0,122, takže závislosť je rastúca a slabá. Body pokrývajú takmer celú plochu grafu, takže sklon voľným okom takmer nevidno.
- Majú všetky záznamy s podobnou spotrebou aj rovnaký výskyt cukrovky? Nie. Napríklad pri spotrebe 35 až 40 kg ročne (763 záznamov) sa výskyt cukrovky pohybuje od 2,9 % do 15,9 %.
4 Jednoduchá lineárna regresia
4.1 Zápis modelu
Najprv použijem iba jednu vysvetľujúcu premennú:
\[ Diabetes_i = \beta_0 + \beta_1 Spotreba_i + u_i. \]
- \(i\) označuje záznam (krajina + rok).
- \(\beta_0\) je konštanta, teda priesečník regresnej priamky s osou Y.
- \(\beta_1\) je smernica priamky: vyjadruje rozdiel predikovaného výskytu cukrovky pri rozdiele spotreby o 1 kg ročne.
- \(u_i\) zachytáva tú časť výskytu cukrovky, 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(diabetes ~ spotreba_na_obyv, data = udaje.model)
vysledky1 <- tidy(model1)
vysledky1 |>
mutate(odhad = sapply(estimate, koef)) |> # desatinná čiarka
select(term, odhad) |>
kable(format = "html", row.names = FALSE, align = c("l", "r"),
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Jednoduchá regresia: výskyt cukrovky a spotreba cukru") |>
styl()| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 8,419 |
| spotreba_na_obyv | 0,0247 |
b0 <- coef(model1)["(Intercept)"]
b1 <- coef(model1)["spotreba_na_obyv"]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. Zobrazujem iba tieto dva stĺpce.
Znak ~ čítam ako „vysvetľujem pomocou“. Vľavo je vysvetľovaná premenná a vpravo vysvetľujúca premenná. Konštantu R zahrnie automaticky.
Odhadnutá regresná priamka:
\[ \widehat{Diabetes}_i = 8{,}419 + 0{,}0247\,Spotreba_i. \]
Interpretácia smernice: Výskyt cukrovky je v percentách, preto rozdiely vyjadrujem v percentuálnych bodoch (p. b.), nie v percentách. Rozdiel 1 kg v ročnej spotrebe cukru na obyvateľa sa v tomto modeli spája s približne o 0,0247 p. b. vyšším predikovaným výskytom cukrovky. Pri rozdiele 10 kg (asi 27 g cukru denne) je to približne 0,25 p. b.
Interpretácia konštanty: Pri nulovej spotrebe cukru by priamka predikovala výskyt cukrovky približne 8,42 %. Najnižšia spotreba v údajoch je však 5,0 kg, takže nula leží mimo pozorovaného rozsahu. Konštantu preto chápem predovšetkým ako parameter polohy priamky.
Poznámka
Súvis medzi záznamami sám osebe nedokazuje, že ak niekto zje o 10 kg cukru ročne viac, jeho riziko cukrovky stúpne o 0,25 p. b. Krajiny a roky sa líšia aj v mnohých ďalších charakteristikách a údaje sú navyše pravdepodobne simulované.
4.3 Regresná priamka v grafe
K bodovému grafu pridám priamku pomocou geom_abline(). Jej konštantu (intercept) a smernicu (slope) prevezmem priamo z odhadnutého modelu model1.
graf_cukor +
geom_abline(intercept = coef(model1)["(Intercept)"],
slope = coef(model1)["spotreba_na_obyv"],
color = "firebrick", linewidth = 1)
5 Viacnásobná lineárna regresia
Výskyt cukrovky môže súvisieť s viacerými charakteristikami súčasne. Model rozšírim o tri premenné:
- obezita je jeden z hlavných rizikových faktorov cukrovky 2. typu,
- HDP na obyvateľa zachytáva rozdiely v životnej úrovni, strave aj v tom, ako často sa cukrovka diagnostikuje,
- urbanizácia súvisí s menším pohybom a vyššou spotrebou spracovaných potravín.
\[ Diabetes_i = \beta_0 + \beta_1 Spotreba_i + \beta_2 Obezita_i + \beta_3 HDP_i + \beta_4 Urbanizacia_i + u_i. \]
Do modelu nezaraďujem denný príjem cukru (denny_prijem). V minulom zadaní som overila, že je len prepočtom spotreby (g za deň = kg × 1000 / 365). Obe premenné teda nesú tú istú informáciu a R by nevedel oddeliť ich vplyv (dokonalá multikolinearita).
model <- lm(diabetes ~ spotreba_na_obyv + obezita + hdp_tis + urbanizacia,
data = udaje.model)
vysledky <- tidy(model)
vysledky |>
mutate(odhad = sapply(estimate, koef)) |>
select(term, odhad) |>
kable(format = "html", row.names = FALSE, align = c("l", "r"),
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Viacnásobná regresia: spotreba cukru, obezita, HDP a urbanizácia") |>
styl()| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 8,459 |
| spotreba_na_obyv | 0,0247 |
| obezita | -0,00165 |
| hdp_tis | -0,00230 |
| urbanizacia | 0,00111 |
b <- coef(model)Vo vzore bolo potrebných osem desatinných miest, aby sa malý koeficient pri GDP nezaokrúhlil na nulu. Tu tento problém nemám, lebo HDP som už v minulom zadaní prepočítala na tisíce USD (hdp_tis). Počet desatinných miest sa navyše nastavuje automaticky podľa veľkosti koeficienta. Odhadnutá rovnica:
\[ \widehat{Diabetes}_i = 8{,}459 + 0{,}0247\,Spotreba_i - 0{,}00165\,Obezita_i - 0{,}00230\,HDP_i + 0{,}00111\,Urbanizacia_i. \]
Každý koeficient sa teraz interpretuje pri rovnakých hodnotách ostatných premenných v modeli:
- Spotreba cukru: rozdiel 1 kg ročne sa pri rovnakej obezite, HDP a urbanizácii spája s približne o 0,0247 p. b. vyšším predikovaným výskytom cukrovky.
- Obezita: rozdiel 1 p. b. v miere obezity sa spája s približne o 0,00165 p. b. nižším predikovaným výskytom cukrovky.
- HDP: rozdiel 1 000 USD v HDP na obyvateľa sa spája s približne o 0,00230 p. b. nižším predikovaným výskytom cukrovky. Vďaka jednotke „tisíc USD“ netreba koeficient násobiť tisícom ako vo vzore.
- Urbanizácia: rozdiel 1 p. b. v miere urbanizácie sa spája s približne o 0,00111 p. b. vyšším predikovaným výskytom cukrovky.
5.1 Zmenil sa koeficient spotreby cukru?
Vo vzore sa koeficient Schooling po pridaní ďalších premenných zmenšil (z 2,23 na 1,92). Vzdelanie totiž súvisí s HDP, a tak časť súvisu, ktorú jednoduchý model pripisoval vzdelaniu, prevzalo HDP. Pozriem sa preto, ako spolu súvisia moje vysvetľujúce premenné:
vysvetlujuce <- c("spotreba_na_obyv", "obezita", "hdp_tis", "urbanizacia")
korelacie <- cor(udaje.model[vysvetlujuce])
korelacie |>
kable(format = "html", digits = 3,
caption = "Korelácie medzi vysvetľujúcimi premennými") |>
styl()| spotreba_na_obyv | obezita | hdp_tis | urbanizacia | |
|---|---|---|---|---|
| spotreba_na_obyv | 1.000 | 0.002 | -0.002 | 0.012 |
| obezita | 0.002 | 1.000 | 0.008 | 0.008 |
| hdp_tis | -0.002 | 0.008 | 1.000 | -0.008 |
| urbanizacia | 0.012 | 0.008 | -0.008 | 1.000 |
max_kor <- max(abs(korelacie[upper.tri(korelacie)]))
zmena <- (b["spotreba_na_obyv"] - b1) / b1Koeficient spotreby cukru je v jednoduchom modeli 0,0247 a vo viacnásobnom 0,0247, takže sa takmer nezmenil (zmena o -0,1 %). Dôvod vidno v tabuľke korelácií: najväčšia korelácia medzi vysvetľujúcimi premennými (v absolútnej hodnote) je iba 0,012. Keď pridaná premenná so spotrebou cukru nesúvisí, nemôže od nej „prevziať“ časť súvisu a odhad sa takmer nezmení. Aj toto zodpovedá tomu, že údaje sú simulované a premenné boli generované nezávisle.
Jednoduchý a viacnásobný model aj tak odpovedajú na odlišné otázky: druhý porovnáva záznamy s rovnakou obezitou, HDP a urbanizáciou.
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. Vyberám počet pozorovaní, \(R^2\) a odhad štandardnej odchýlky chyby. Pre porovnanie uvádzam oba modely.
suhrn1 <- glance(model1)
suhrn <- glance(model)
suhrny <- data.frame(
Model = c("jednoduchý", "viacnásobný"),
rbind(suhrn1[c("nobs", "r.squared", "sigma")],
suhrn[c("nobs", "r.squared", "sigma")])
)
suhrny |>
kable(format = "html", digits = c(0, 0, 4, 3), row.names = FALSE,
col.names = c("Model", "Počet záznamov", "R²", "Odhad SD chyby (p. b.)"),
caption = "Základné charakteristiky modelov") |>
styl()| Model | Počet záznamov | R² | Odhad SD chyby (p. b.) |
|---|---|---|---|
| jednoduchý | 10000 | 0.0149 | 3.767 |
| viacnásobný | 10000 | 0.0150 | 3.767 |
r2 <- suhrn$r.squarednobsje počet záznamov použitých pri odhade.r.squaredje podiel výberovej variability výskytu cukrovky zachytený modelom.sigmaje odhad štandardnej odchýlky modelovej chyby, tu vyjadrený v percentuálnych bodoch.
Hodnota \(R^2\) je vo viacnásobnom modeli 0,0150. Model teda zachytáva približne 1,5 % výberovej variability výskytu cukrovky. To je veľmi málo: takmer všetky rozdiely vo výskyte cukrovky medzi záznamami zostávajú nevysvetlené. Zodpovedá to zisteniu z minulého zadania, že premenné v týchto (pravdepodobne simulovaných) údajoch spolu takmer nesúvisia. Pridanie troch premenných zvýšilo \(R^2\) z 0,0149 na 0,0150.
Nízke \(R^2\) neznamená, že odhady koeficientov sú „zlé“. Znamená len, že výskyt cukrovky ovplyvňuje najmä niečo, čo v modeli nie je (alebo náhoda). Hodnota \(R^2\) tiež neznamená podiel správnych predikcií a nedokazuje kauzalitu.
Úplný textový súhrn je dostupný príkazom summary(model). tidy(model) okrem odhadov obsahuje aj std.error, statistic a p.value. Ich vysvetlenie patrí na hodinu venovanú štatistickému testovaniu, preto ich tu nezobrazujem.
7 Skutočná hodnota, predikcia a rezíduum
Pre každý záznam môžem porovnať skutočný výskyt cukrovky s hodnotou na regresnej ploche:
\[ e_i = y_i - \widehat y_i. \]
Tento rozdiel sa nazýva rezíduum. Kladné rezíduum znamená, že skutočná hodnota je vyššia než predikcia, záporné znamená opak. Rezíduum viem vypočítať, zatiaľ čo skutočnú modelovú chybu \(u_i\) nepoznám.
porovnanie <- data.frame(
Krajina = udaje.model$country,
Rok = udaje.model$year,
Skutocnost = udaje.model$diabetes,
Predikcia = fitted(model),
Reziduum = residuals(model)
)
head(porovnanie, 8) |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Krajina", "Rok", "Skutočnosť", "Predikcia", "Rezíduum"),
caption = "Skutočný a predikovaný výskyt cukrovky (%)") |>
styl()| Krajina | Rok | Skutočnosť | Predikcia | Rezíduum |
|---|---|---|---|---|
| France | 1972 | 11.50 | 8.74 | 2.76 |
| Australia | 2003 | 6.52 | 8.98 | -2.46 |
| Germany | 1963 | 8.25 | 9.23 | -0.98 |
| France | 1965 | 11.00 | 9.65 | 1.34 |
| Germany | 2010 | 14.89 | 8.93 | 5.96 |
| South Africa | 1997 | 14.89 | 10.07 | 4.82 |
| India | 1965 | 11.99 | 8.98 | 3.01 |
| Mexico | 1970 | 8.88 | 9.84 | -0.96 |
prvy <- porovnanie[1, ]Úloha zo vzoru (ručné overenie): V prvom riadku (France, 1972) je skutočný výskyt cukrovky 11,4951 % a predikcia 8,7378 %. Ich rozdiel 11,4951 − 8,7378 = 2,7573 sa rovná rezíduu 2,7573.
Funkcia coef(model) vracia koeficienty, fitted(model) vyrovnané hodnoty a residuals(model) rezíduá. Maticu vysvetľujúcich premenných by som získala príkazom model.matrix(model).
8 Predikcia pre zadanú hodnotu
Na ukážku predikcie sa vrátim k jednoduchému modelu. Vypočítam, aký výskyt cukrovky priamka predikuje pre tri hodnoty spotreby. Prvá zodpovedá hornej hranici odporúčania WHO (50 g denne × 365 dní = 18,25 kg ročne), ktorú som používala v minulom zadaní:
nove_udaje <- data.frame(spotreba_na_obyv = c(18.25, 40, 60))
nove_udaje$Predikcia <- predict(model1, newdata = nove_udaje)
nove_udaje |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Spotreba cukru (kg/rok)", "Predikcia výskytu cukrovky (%)"),
caption = "Predikcie jednoduchého regresného modelu") |>
styl()| Spotreba cukru (kg/rok) | Predikcia výskytu cukrovky (%) |
|---|---|
| 18.25 | 8.87 |
| 40.00 | 9.41 |
| 60.00 | 9.90 |
Rozdiel medzi spotrebou 18,25 kg a 60 kg ročne sa v modeli spája s rozdielom 1,03 p. b. vo výskyte cukrovky. Ide o bodové predikcie, nie o isté výsledky. Skutočný záznam môže mať iný výskyt cukrovky. Neistotu predikcie (predikčné intervaly) preberieme neskôr.
9 Malé cvičenia
9.1 Vysvetľovaná a vysvetľujúce premenné
Vysvetľovaná premenná je výskyt cukrovky (diabetes). V jednoduchom modeli je vysvetľujúcou premennou spotreba cukru na obyvateľa (spotreba_na_obyv). Vo viacnásobnom modeli k nej pribudli obezita, HDP na obyvateľa a urbanizácia.
9.2 Jednoduchý model s obezitou
Namiesto modelu s GDP zo vzoru odhadnem model s obezitou, lebo pri cukrovke dáva väčší zmysel.
model_ob <- lm(diabetes ~ obezita, data = udaje.model)
tidy(model_ob) |>
mutate(odhad = sapply(estimate, koef)) |>
select(term, odhad) |>
kable(format = "html", row.names = FALSE, align = c("l", "r"),
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Jednoduchá regresia: výskyt cukrovky a obezita") |>
styl()| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 9,384 |
| obezita | -0,00157 |
ggplot(udaje.model, aes(x = obezita, y = diabetes)) +
geom_point(color = farba, size = 1, alpha = 0.15) +
geom_abline(intercept = coef(model_ob)["(Intercept)"],
slope = coef(model_ob)["obezita"],
color = "firebrick", linewidth = 1) +
labs(x = "Miera obezity (%)", y = "Výskyt cukrovky (%)") +
theme_minimal(base_size = 13)
Rozdiel 1 p. b. v miere obezity sa spája s približne o 0,00157 p. b. nižším predikovaným výskytom cukrovky.
9.3 Koeficient spotreby v oboch modeloch
V jednoduchom modeli koeficient 0,0247 porovnáva záznamy, ktoré sa líšia spotrebou cukru, bez ohľadu na ostatné charakteristiky. Vo viacnásobnom modeli koeficient 0,0247 porovnáva záznamy, ktoré sa líšia spotrebou, ale majú rovnakú obezitu, HDP a urbanizáciu. Preto dopĺňam „pri rovnakých ostatných premenných“: bez toho by sa súvis spotreby mohol miešať so súvisom premenných, ktoré sa so spotrebou menia spolu.
9.4 Predikcia pre spotrebu 30 kg
predikcia_30 <- predict(model1, newdata = data.frame(spotreba_na_obyv = 30))
rucne_30 <- coef(model1)["(Intercept)"] + coef(model1)["spotreba_na_obyv"] * 30
c(predict = unname(predikcia_30), rucne = unname(rucne_30))predict rucne
9.16056 9.16056 Ručne: 8,419 + 0,0247 × 30 = 9,161 %, čo sa zhoduje s výsledkom funkcie predict() (9,161 %).
9.5 Kladné a záporné rezíduum
kladne <- porovnanie[porovnanie$Reziduum > 0, ][1, ]
zaporne <- porovnanie[porovnanie$Reziduum < 0, ][1, ]
rbind(kladne, zaporne) |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Krajina", "Rok", "Skutočnosť", "Predikcia", "Rezíduum"),
caption = "Príklad kladného a záporného rezídua") |>
styl()| Krajina | Rok | Skutočnosť | Predikcia | Rezíduum |
|---|---|---|---|---|
| France | 1972 | 11.50 | 8.74 | 2.76 |
| Australia | 2003 | 6.52 | 8.98 | -2.46 |
- France (1972) má kladné rezíduum 2,76: skutočný výskyt cukrovky je vyšší, než predikuje model.
- Australia (2003) má záporné rezíduum -2,46: skutočný výskyt cukrovky je nižší, než predikuje model.
9.6 Prečo koeficient neznamená kauzálny účinok
Koeficient opisuje, ako sa v údajoch spolu menia spotreba cukru a výskyt cukrovky. Nehovorí, čo by sa stalo, keby sme spotrebu cukru zmenili. Súvis môžu spôsobovať iné faktory, ktoré v modeli nie sú (strava, vek populácie, pohyb, zdravotníctvo), alebo môže ísť o náhodu. Na kauzálny záver by bol potrebný experiment alebo špeciálne metódy. Pri simulovaných údajoch navyše nejde o skutočné správanie ľudí.