model <- lm(diabetes ~ 1 + spotreba_na_obyv + obezita + hdp_tis + urbanizacia,
data = udaje.model)Spotreba cukru a výskyt cukrovky
Cvičenie 5 – základná diagnostika regresného modelu, hat matica a vplyvné pozorovania
TipMožnosti dokumentu
Tlačidlo Skryť kód vpravo dole skryje všetok R kód aj technické výpisy konzoly. Zostanú len opisy, vysvetlenia, tabuľky a grafy. Opätovným kliknutím sa kód zobrazí. Rovnakú funkciu má aj ponuka Kód vpravo hore.
1 Predmet a údaje
Nadväzujem na minulé zadanie, v ktorom som odhadla viacnásobný regresný model pre výskyt cukrovky:
\[ Diabetes_i=\beta_0+\beta_1\,Spotreba_i+\beta_2\,Obezita_i+\beta_3\,HDP_i+\beta_4\,Urbanizacia_i+u_i . \]
Tentoraz model neinterpretujem len bodovo, ale testujem jeho koeficienty, posudzujem kvalitu vyrovnania a robím diagnostiku rezíduí: normalitu, odľahlé hodnoty, grafickú analýzu a vplyv jednotlivých pozorovaní pomocou hat matice.
1.1 Príprava údajov a pracovného prostredia
Pracujem s rovnakými údajmi ako v minulom zadaní (sugar_consumption_dataset.csv). Na začiatku som načítala uložený workspace pamat.Rdata, takže som nemusela údaje znovu importovať a čistiť; ak súbor nie je k dispozícii, údaje sa pripravia priamo z CSV súboru. Pre model som vybrala výskyt cukrovky, spotrebu cukru na obyvateľa, obezitu, HDP na obyvateľa (v tis. USD) a urbanizáciu – spolu 10 000 záznamov bez chýbajúcich hodnôt. Pracovné prostredie som uložila príkazom save.image(), nepotrebné objekty zmazala príkazom rm() a uložený obsah opäť načítala príkazom load(). Rovnako ho na konci dokumentu uložím pre ďalšie cvičenie.
2 Lineárna regresia
Člen 1 + vo vzorci znamená absolútny člen. R ho zahrnie aj bez neho, tu ho píšem kvôli prehľadnosti.
tidy(model) |>
mutate(
p.value = p_hod(p.value),
Vyznamnost = case_when(
tidy(model)$p.value < 0.001 ~ "***",
tidy(model)$p.value < 0.01 ~ "**",
tidy(model)$p.value < 0.05 ~ "*",
tidy(model)$p.value < 0.10 ~ ".",
TRUE ~ ""
)
) |>
kable(digits = 5, align = c("l", "r", "r", "r", "r", "c"),
col.names = c("Premenná", "Odhad", "Smer. chyba", "t hodnota",
"p-hodnota", "Významnosť"),
caption = "Odhady regresných koeficientov a t-testy") |>
styl()| Premenná | Odhad | Smer. chyba | t hodnota | p-hodnota | Významnosť |
|---|---|---|---|---|---|
| (Intercept) | 8.45873 | 0.15789 | 53.57387 | < 0,001 | *** |
| spotreba_na_obyv | 0.02471 | 0.00201 | 12.28516 | < 0,001 | *** |
| obezita | -0.00165 | 0.00373 | -0.44306 | 0,658 | |
| hdp_tis | -0.00230 | 0.00265 | -0.86551 | 0,387 | |
| urbanizacia | 0.00111 | 0.00163 | 0.68592 | 0,493 |
koef_tab <- tidy(model)
b <- coef(model)
n <- nobs(model)
k <- length(b) # počet koeficientov vrátane absolútneho člena2.1 Vysvetlenie stĺpcov tabuľky
| Stĺpec | Vysvetlenie |
|---|---|
| Premenná | Názov premennej. (Intercept) je absolútny člen \(\hat\beta_0\). |
| Odhad | Bodový odhad \(\hat\beta_j\): zmena predikovanej hodnoty \(y\) pri jednotkovej zmene \(x_j\), ak ostatné premenné ostanú nezmenené (ceteris paribus). |
| Smer. chyba | \(SE(\hat\beta_j)\) – neistota odhadu. Čím menšia, tým presnejší odhad. |
| t hodnota | Testovacia štatistika \(t_j=\hat\beta_j/SE(\hat\beta_j)\). |
| p-hodnota | Pravdepodobnosť, že pri platnosti \(H_0:\beta_j=0\) dostaneme \(|t|\) aspoň taký veľký ako pozorovaný. Malá p-hodnota (spravidla \(<0{,}05\)) je dôvodom zamietnuť \(H_0\). |
| Významnosť | Hviezdičky podľa p-hodnoty (pozri tabuľku nižšie). |
2.2 t-test štatistickej významnosti koeficientov
Pre každý koeficient testujem
\[ H_0:\beta_j=0 \qquad\text{proti}\qquad H_1:\beta_j\neq 0, \qquad t_j=\frac{\hat\beta_j}{SE(\hat\beta_j)} \overset{H_0}{\sim} t(n-k). \]
Tu \(n=10000\) je počet pozorovaní a \(k=5\) počet koeficientov vrátane absolútneho člena, takže stupne voľnosti sú \(n-k=9995\). Smerodajné chyby sú odmocninami diagonály matice \(\widehat{\operatorname{Var}}(\hat{\boldsymbol\beta})=\hat\sigma^2(\mathbf X^\mathsf T\mathbf X)^{-1}\).
Rozhodovacie pravidlo: ak \(p<\alpha\) (zvyčajne \(\alpha=0{,}05\)), zamietam \(H_0\) a koeficient považujem za štatisticky významný.
| Symbol | p-hodnota | Interpretácia |
|---|---|---|
*** |
\(p<0{,}001\) | významný na hladine 0,1 % |
** |
\(p<0{,}01\) | významný na hladine 1 % |
* |
\(p<0{,}05\) | významný na hladine 5 % |
. |
\(p<0{,}10\) | významný na hladine 10 % |
| bez symbolu | \(p\geq 0{,}10\) | nevýznamný |
2.2.1 Ručné overenie t-testu
Smerodajné chyby, t-štatistiky aj p-hodnoty si prepočítam maticovo a porovnám s výsledkom lm():
X <- model.matrix(model)
sigma2 <- sum(residuals(model)^2) / (n - k)
V <- sigma2 * solve(t(X) %*% X) # kovariančná matica odhadov
se <- sqrt(diag(V))
t_stat <- b / se
p_val <- 2 * pt(abs(t_stat), df = n - k, lower.tail = FALSE)
t_krit <- qt(0.975, df = n - k)
stopifnot(isTRUE(all.equal(unname(se), koef_tab$std.error)),
isTRUE(all.equal(unname(t_stat), koef_tab$statistic)),
isTRUE(all.equal(unname(p_val), koef_tab$p.value)))
data.frame(Premenna = names(b), SE = se, t = t_stat, p = p_hod(p_val),
Zamietam_H0 = ifelse(abs(t_stat) > t_krit, "áno", "nie")) |>
kable(digits = 4, row.names = FALSE,
col.names = c("Premenná", "SE (ručne)", "t (ručne)", "p (ručne)",
"Zamietam H0 pri 5 %?"),
caption = "Ručný výpočet t-testov") |>
styl()| Premenná | SE (ručne) | t (ručne) | p (ručne) | Zamietam H0 pri 5 %? |
|---|---|---|---|---|
| (Intercept) | 0.1579 | 53.5739 | < 0,001 | áno |
| spotreba_na_obyv | 0.0020 | 12.2852 | < 0,001 | áno |
| obezita | 0.0037 | -0.4431 | 0,658 | nie |
| hdp_tis | 0.0027 | -0.8655 | 0,387 | nie |
| urbanizacia | 0.0016 | 0.6859 | 0,493 | nie |
Ručný výpočet sa zhoduje s lm(). Kritická hodnota na hladine 5 % je \(t_{0{,}975}(9995)=1,960\). Pri takmer desiatich tisícoch stupňov voľnosti je t-rozdelenie prakticky totožné s normovaným normálnym (\(z_{0{,}975}=1{,}960\)).
2.3 Interpretácia koeficientov
Spotreba cukru (spotreba_na_obyv): rozdiel 1 kg v ročnej spotrebe na obyvateľa sa pri rovnakej obezite, HDP a urbanizácii spája s predikovaným výskytom cukrovky vyšším o 0,0247 percentuálneho bodu (p. b.). Pri rozdiele 10 kg je to 0,25 p. b. Koeficient je štatisticky významný (\(t=12,29\), \(p < 0,001\)), hoci jeho praktická veľkosť je malá.
Obezita, HDP a urbanizácia: ich p-hodnoty sú 0,658, 0,387 a 0,493. Ani na hladine 10 % nezamietam nulovú hypotézu. Odhady (-0,00165, -0,00230, 0,00111) sú teda od nuly štatisticky nerozlíšiteľné. Záporné znamienko pri obezite by bolo v rozpore s medicínskymi poznatkami, no pri nevýznamnom koeficiente ho neinterpretujem ako skutočný súvis.
Absolútny člen: \(\hat\beta_0=8,459\) je predikcia výskytu cukrovky pri nulových hodnotách všetkých vysvetľujúcich premenných. Nulová spotreba cukru ani nulová urbanizácia sa v údajoch nevyskytujú, preto má len technický význam (poloha regresnej nadroviny).
3 Kvalita vyrovnania modelu
suhrn <- glance(model)
suhrn |>
select(-deviance) |>
mutate(p.value = p_hod(p.value)) |>
kable(digits = 4,
col.names = c("R²", "Upravené R²", "Sigma", "F štat.", "p-hodnota",
"df", "logLik", "AIC", "BIC", "df.residual", "N"),
caption = "Súhrnné charakteristiky modelu") |>
styl() |>
scroll_box(width = "100%")| R² | Upravené R² | Sigma | F štat. | p-hodnota | df | logLik | AIC | BIC | df.residual | N |
|---|---|---|---|---|---|---|---|---|---|---|
| 0.015 | 0.0146 | 3.7673 | 38.1489 | < 0,001 | 4 | -27450.46 | 54912.93 | 54956.19 | 9995 | 10000 |
3.1 Vysvetlenie stĺpcov
| Stĺpec | Vysvetlenie |
|---|---|
| R² | Koeficient determinácie – podiel variability \(y\) vysvetlený modelom. |
| Upravené R² | \(\bar R^2\) – penalizuje počet premenných; vhodné na porovnanie modelov s rôznym počtom premenných. |
| Sigma | Reziduálna štandardná chyba \(\hat\sigma\) – odhad smerodajnej odchýlky náhodnej zložky. |
| F štat. | Štatistika celkového F-testu spoločnej významnosti vysvetľujúcich premenných. |
| p-hodnota | p-hodnota F-testu. |
| df | Počet vysvetľujúcich premenných bez absolútneho člena (\(k-1\)), stupne voľnosti čitateľa F-testu. |
| logLik | Maximálna hodnota logaritmu vierohodnostnej funkcie. |
| AIC, BIC | Informačné kritériá – nižšia hodnota znamená lepší kompromis medzi presnosťou a zložitosťou. |
| df.residual | Reziduálne stupne voľnosti \(n-k\). |
| N | Počet pozorovaní použitých pri odhade. |
| deviance | RSS – súčet štvorcov rezíduí (z tabuľky vynechaný, použijem ho pri ručných výpočtoch). |
3.2 Interpretácia pre náš model
Všetky charakteristiky si najprv overím ručne zo vzorcov:
e <- residuals(model)
y <- udaje.model$diabetes
RSS <- sum(e^2)
TSS <- sum((y - mean(y))^2)
R2 <- 1 - RSS / TSS
R2_adj <- 1 - (1 - R2) * (n - 1) / (n - k)
sig <- sqrt(RSS / (n - k))
F_stat <- (R2 / (k - 1)) / ((1 - R2) / (n - k))
p_F <- pf(F_stat, k - 1, n - k, lower.tail = FALSE)
F_krit <- qf(0.95, k - 1, n - k)
logL <- -n / 2 * (log(2 * pi) + log(RSS / n) + 1)
K_par <- k + 1 # koeficienty + rozptyl náhodnej zložky
AIC_r <- -2 * logL + 2 * K_par
BIC_r <- -2 * logL + K_par * log(n)
rucne <- c(R2, R2_adj, sig, F_stat, p_F, logL, AIC_r, BIC_r, RSS)
z_R <- with(suhrn, c(r.squared, adj.r.squared, sigma, statistic, p.value,
logLik, AIC, BIC, deviance))
stopifnot(isTRUE(all.equal(unname(rucne), unname(z_R))))
data.frame(
Ukazovatel = c("R²", "Upravené R²", "Sigma", "F", "p-hodnota F",
"logLik", "AIC", "BIC", "RSS (deviance)"),
Rucne = rucne, glance = z_R
) |>
kable(digits = 5, row.names = FALSE,
col.names = c("Ukazovateľ", "Ručný výpočet", "glance()"),
caption = "Ručné overenie súhrnných charakteristík") |>
styl()| Ukazovateľ | Ručný výpočet | glance() |
|---|---|---|
| R² | 0.01504 | 0.01504 |
| Upravené R² | 0.01464 | 0.01464 |
| Sigma | 3.76730 | 3.76730 |
| F | 38.14894 | 38.14894 |
| p-hodnota F | 0.00000 | 0.00000 |
| logLik | -27450.46471 | -27450.46471 |
| AIC | 54912.92941 | 54912.92941 |
| BIC | 54956.19145 | 54956.19145 |
| RSS (deviance) | 141854.36998 | 141854.36998 |
1. Koeficient determinácie. \(R^2=1-RSS/TSS=0,0150\). Model vysvetľuje len 1,50 % variability výskytu cukrovky, zvyšných 98,50 % zostáva nevysvetlených. To je veľmi málo. Zodpovedá to zisteniu z minulých zadaní, že premenné v týchto údajoch spolu takmer nesúvisia. Nízke \(R^2\) však neznamená, že odhady sú skreslené – len to, že výskyt cukrovky určuje hlavne niečo mimo modelu.
2. Upravené \(R^2\). \(\bar R^2=1-(1-R^2)\frac{n-1}{n-k}=0,0146\). Je o niečo nižšie ako \(R^2\) – penalizácia za tri premenné, ktoré k vysvetleniu prakticky neprispievajú. Pri veľkom \(n\) je rozdiel malý.
3. Reziduálna štandardná chyba. \(\hat\sigma=\sqrt{RSS/(n-k)}=3,767\) p. b. Typická odchýlka skutočného výskytu cukrovky od predikcie je teda asi 3,8 percentuálneho bodu. Pre porovnanie: smerodajná odchýlka samotnej premennej diabetes je 3,795 p. b. Model teda neistotu takmer nezmenšil.
4.–5. F-test. Testujem
\[ H_0:\beta_1=\beta_2=\beta_3=\beta_4=0 \qquad\text{proti}\qquad H_1:\text{aspoň jeden }\beta_j\neq 0,\ j=1,\dots,4, \]
\[ F=\frac{R^2/(k-1)}{(1-R^2)/(n-k)}=38,15 \overset{H_0}{\sim}F(4;\,9995). \]
Kritická hodnota na hladine 5 % je \(F_{0{,}95}=2,373\) a p-hodnota je < 0,001. Zamietam \(H_0\) – model ako celok je štatisticky významný. Je to vďaka spotrebe cukru. Paradox model je významný, ale vysvetľuje len 1,5 % je pri veľkých vzorkách bežný: F-test hovorí, či je \(R^2\) rozlíšiteľné od nuly, nie či je veľké.
6. df. \(k-1=4\) vysvetľujúce premenné, čo sú stupne voľnosti čitateľa F-štatistiky.
7. logLik. Pri normálnom lineárnom modeli platí \(\ln L=-\frac n2\left[\ln(2\pi)+\ln\frac{RSS}{n}+1\right]=-27 450,46\). Samotná hodnota nemá interpretáciu; používa sa pri porovnaní modelov na rovnakých údajoch.
8.–9. AIC a BIC. \(AIC=-2\ln L+2K=54 912,93\) a \(BIC=-2\ln L+K\ln n=54 956,19\), kde \(K=6\) (5 koeficientov + rozptyl \(\sigma^2\)). Samotné čísla nič nehovoria – porovnám ich s upraveným modelom v kapitole Úprava modelu.
10.–11. df.residual a N. \(n-k=10000-5=9995\); použitých bolo všetkých \(N=10000\) záznamov, žiadny nevypadol pre chýbajúce hodnoty.
12. deviance. \(RSS=141 854,4\). Práve túto veličinu metóda najmenších štvorcov minimalizuje.
DôležitéVýsledok
Model je štatisticky významný (F-test), významná je len spotreba cukru (t-testy), ale vysvetľovacia schopnosť je veľmi nízka (\(R^2\approx 0,015\)).
4 Grafická analýza rezíduí – plot(model)
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))Červená čiara v grafoch je vyhladená krivka LOWESS (locally weighted scatterplot smoothing). Pre každý bod ju R počíta váženou metódou najmenších štvorcov z približne 2/3 najbližších pozorovaní, pričom bližšie pozorovania majú vyššiu váhu. Čísla pri bodoch označujú tri najextrémnejšie pozorovania v danom grafe.
fit <- fitted(model)
r_std <- rstandard(model)
rozsah_fit <- range(fit)
rozsah_e <- range(e)4.1 Residuals vs Fitted
- Os X: vyrovnané hodnoty \(\hat y_i\), os Y: rezíduá \(e_i\).
- Vyrovnané hodnoty ležia len v úzkom pásme 8,45 až 10,25 %, kým rezíduá od -6,62 do 6,64 p. b. Graf preto vyzerá ako rovnomerne vyplnený obdĺžnik – to je dôsledok nízkeho \(R^2\).
- Červená čiara je takmer vodorovná pri nule. To naznačuje, že model nemá systematické skreslenie predikcií a nevidím nelinearitu.
- Šírka pásma sa so zmenou \(\hat y\) nemení, nevzniká lievik – nič nenasvedčuje heteroskedasticite.
- Pásmo má ostré horné aj dolné okraje, nie postupne rednúce body ako pri normálnom rozdelení. To je prvý náznak ľahkých chvostov (pozri Q–Q graf).
4.2 Q–Q Residuals (normálny Q–Q graf)
- Os X: teoretické kvantily normálneho rozdelenia, os Y: štandardizované rezíduá.
- V strede (približne −1 až +1) body ležia blízko priamky.
- Na oboch koncoch sa body ohýbajú smerom k vodorovnej osi: vľavo ležia nad priamkou, vpravo pod ňou. Vzniká tvar písmena S.
- Je to typický obraz symetrického rozdelenia s ľahkými chvostmi (pozri prílohu textu o hat maticiach – rovnomerné rozdelenie). Extrémne rezíduá sú menšie, ako by sme pri normalite čakali. Najväčšie štandardizované rezíduum má absolútnu hodnotu len 1,76, hoci pri 10000 normálnych pozorovaniach by sme čakali hodnoty okolo 3,5 – 4.
4.3 Scale–Location
- Os X: vyrovnané hodnoty, os Y: \(\sqrt{|r_i|}\) (odmocnina absolútnych štandardizovaných rezíduí).
- Červená čiara je takmer vodorovná, rozptyl bodov sa po osi X nemení – rozptyl rezíduí je približne konštantný (homoskedasticita).
- Žiadny bod nepresahuje hodnotu \(\sqrt{2}\approx1{,}41\), teda žiadne \(|r_i|>2\). Neexistujú teda ani odľahlé hodnoty z pohľadu rozptylu.
4.4 Residuals vs Leverage
- Os X: leverage \(h_{ii}\) (pákový účinok, pozri ďalšiu kapitolu), os Y: štandardizované rezíduá.
- Prerušované červené kontúry Cookovej vzdialenosti (0,5 a 1) sa pri tak malých leverage do grafu ani nezmestia. Žiadne pozorovanie nemá podstatný vplyv na odhad.
- Všetky leverage sú menšie ako 0,0012. Aj najvzdialenejší záznam teda tvorí len nepatrný zlomok celkovej informácie.
PoznámkaVýsledky grafickej analýzy
Linearita a homoskedasticita sú v poriadku, vplyvné pozorovania nevidím. Problémom je tvar rozdelenia rezíduí – ľahké chvosty. Overím to formálnym testom normality.
5 Hat matica, leverage a vplyvné pozorovania
Diagnostické grafy 3 a 4 aj test odľahlých hodnôt pracujú s hat maticou. Táto kapitola ukazuje, čo sa za nimi skrýva.
5.1 Hat matica a jej vlastnosti
Vyrovnané hodnoty dostaneme ako
\[ \widehat{\mathbf y}=\mathbf X\widehat{\boldsymbol\beta} =\underbrace{\mathbf X(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf X^\mathsf T}_{\mathbf H}\mathbf y =\mathbf H\mathbf y . \]
Matica \(\mathbf H\) nasadí \(\mathbf y\) klobúk (hat), odtiaľ jej názov. Je symetrická (\(\mathbf H^\mathsf T=\mathbf H\)) a idempotentná (\(\mathbf H^2=\mathbf H\)). Jej diagonálne prvky \(h_{ii}\) sú leverage (pákový účinok) pozorovaní:
\[ \frac1n\le h_{ii}\le 1,\qquad \sum_{i=1}^n h_{ii}=\operatorname{tr}(\mathbf H)=k,\qquad \bar h=\frac kn . \]
Dolná hranica \(1/n\) platí pre model s absolútnym členom. Leverage závisí len od vysvetľujúcich premenných \(\mathbf X\), nie od \(y\). Meria, ako ďaleko leží pozorovanie od ťažiska ostatných v priestore vysvetľujúcich premenných.
5.1.1 Ukážka na malej vzorke
Celá hat matica pre 10000 záznamov by mala \(10000\times10000\) prvkov (asi 800 MB). Vlastnosti preto najprv ukážem na prvých desiatich záznamoch a jednoduchšom modeli so spotrebou a obezitou:
mala <- udaje.model[1:10, ]
X_m <- model.matrix(~ spotreba_na_obyv + obezita, data = mala)
H_m <- X_m %*% solve(t(X_m) %*% X_m) %*% t(X_m)
m_mala <- lm(diabetes ~ spotreba_na_obyv + obezita, data = mala)
round(H_m, 3) 1 2 3 4 5 6 7 8 9 10
1 0.353 0.138 0.060 0.188 0.311 -0.145 0.103 -0.030 0.155 -0.132
2 0.138 0.242 0.192 -0.110 0.057 -0.018 0.298 0.003 0.140 0.059
3 0.060 0.192 0.176 -0.074 0.013 0.081 0.241 0.065 0.114 0.132
4 0.188 -0.110 -0.074 0.496 0.295 0.149 -0.222 0.182 0.066 0.031
5 0.311 0.057 0.013 0.295 0.316 -0.059 -0.005 0.034 0.128 -0.090
6 -0.145 -0.018 0.081 0.149 -0.059 0.386 -0.022 0.272 0.027 0.329
7 0.103 0.298 0.241 -0.222 -0.005 -0.022 0.386 -0.014 0.147 0.087
8 -0.030 0.003 0.065 0.182 0.034 0.272 -0.014 0.211 0.053 0.223
9 0.155 0.140 0.114 0.066 0.128 0.027 0.147 0.053 0.120 0.048
10 -0.132 0.059 0.132 0.031 -0.090 0.329 0.087 0.223 0.048 0.314data.frame(
Vlastnost = c("Symetria H = Hᵀ", "Idempotentnosť H² = H",
"Stopa tr(H) = k", "Zhoda diag(H) s hatvalues()",
"Zhoda Hy s fitted()"),
Vysledok = c(isTRUE(all.equal(H_m, t(H_m))),
isTRUE(all.equal(H_m %*% H_m, H_m)),
paste0(cislo(sum(diag(H_m)), 4), " (k = ", ncol(X_m), ")"),
isTRUE(all.equal(unname(diag(H_m)), unname(hatvalues(m_mala)))),
isTRUE(all.equal(as.numeric(H_m %*% mala$diabetes),
unname(fitted(m_mala)))))
) |>
kable(col.names = c("Vlastnosť", "Výsledok"),
caption = "Overenie vlastností hat matice na 10 záznamoch") |>
styl()| Vlastnosť | Výsledok |
|---|---|
| Symetria H = Hᵀ | TRUE |
| Idempotentnosť H² = H | TRUE |
| Stopa tr(H) = k | 3,0000 (k = 3) |
| Zhoda diag(H) s hatvalues() | TRUE |
| Zhoda Hy s fitted() | TRUE |
data.frame(zaznam = 1:10,
krajina = mala$country, rok = mala$year,
spotreba = mala$spotreba_na_obyv, obezita = mala$obezita,
leverage = diag(H_m)) |>
kable(digits = 3, col.names = c("Záznam", "Krajina", "Rok", "Spotreba",
"Obezita", "Leverage h_ii"),
caption = "Leverage v malej vzorke (priemer k/n = 3/10 = 0,3)") |>
styl() |>
row_spec(which.max(diag(H_m)), bold = TRUE, background = "#fde8d8")| Záznam | Krajina | Rok | Spotreba | Obezita | Leverage h_ii |
|---|---|---|---|---|---|
| 1 | France | 1972 | 12.828 | 34.094 | 0.353 |
| 2 | Australia | 2003 | 21.363 | 16.040 | 0.242 |
| 3 | Germany | 1963 | 32.077 | 14.406 | 0.176 |
| 4 | France | 1965 | 47.649 | 38.499 | 0.496 |
| 5 | Germany | 2010 | 23.214 | 36.191 | 0.316 |
| 6 | South Africa | 1997 | 66.799 | 14.448 | 0.386 |
| 7 | India | 1965 | 19.314 | 10.976 | 0.386 |
| 8 | Mexico | 1970 | 55.629 | 19.811 | 0.211 |
| 9 | Mexico | 1986 | 28.699 | 22.690 | 0.120 |
| 10 | Russia | 2013 | 59.171 | 10.891 | 0.314 |
Najvyššiu leverage (0,496) má záznam 4. Jeho kombinácia spotreby a obezity je od ostatných deviatich najvzdialenejšia. Pri desiatich pozorovaniach môže jediný bod výrazne potiahnuť regresnú rovinu k sebe.
5.1.2 Leverage v celom modeli
Pre celý model netreba počítať celú maticu \(\mathbf H\). Stačí jej diagonála: \(h_{ii}=\mathbf x_i^\mathsf T(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf x_i\), čo je riadkový súčet súčinu \(\mathbf X(\mathbf X^\mathsf T\mathbf X)^{-1}\) a \(\mathbf X\) po prvkoch.
h <- rowSums((X %*% solve(t(X) %*% X)) * X)
stopifnot(isTRUE(all.equal(unname(h), unname(hatvalues(model)))))
data.frame(
Ukazovatel = c("Súčet h_ii (= k)", "Priemer h_ii (= k/n)", "Dolná hranica 1/n",
"Minimum", "Medián", "Maximum",
"Orientačný prah 2k/n", "Počet h_ii > 2k/n"),
Hodnota = c(cislo(sum(h), 4), cislo(mean(h), 6), cislo(1/n, 6),
cislo(min(h), 6), cislo(median(h), 6), cislo(max(h), 6),
cislo(2*k/n, 6), as.character(sum(h > 2*k/n)))
) |>
kable(col.names = c("Ukazovateľ", "Hodnota"), align = c("l", "r"),
caption = "Leverage v modeli so všetkými záznamami") |>
styl()| Ukazovateľ | Hodnota |
|---|---|
| Súčet h_ii (= k) | 5,0000 |
| Priemer h_ii (= k/n) | 0,000500 |
| Dolná hranica 1/n | 0,000100 |
| Minimum | 0,000104 |
| Medián | 0,000488 |
| Maximum | 0,001218 |
| Orientačný prah 2k/n | 0,001000 |
| Počet h_ii > 2k/n | 45 |
Súčet leverage je presne \(k=5\) a priemer \(k/n=0,0005\). Orientačný prah \(2k/n\) prekračuje 45 záznamov. Aj najvyššia leverage 0,00122 je však v absolútnom vyjadrení zanedbateľná: tento záznam tvorí asi 0,12 % vlastnej predikcie, zvyšok predikcie určujú ostatné pozorovania. Prahy \(2k/n\) a \(3k/n\) sú len orientačné.
i_h <- which.max(h)
udaje.model[order(h, decreasing = TRUE)[1:5], ] |>
mutate(leverage = sort(h, decreasing = TRUE)[1:5]) |>
kable(digits = c(0, 0, 2, 2, 2, 2, 2, 6), row.names = TRUE,
col.names = c("Krajina", "Rok", "Cukrovka", "Spotreba", "Obezita",
"HDP (tis.)", "Urbanizácia", "Leverage"),
caption = "Päť záznamov s najvyššou leverage") |>
styl()| Krajina | Rok | Cukrovka | Spotreba | Obezita | HDP (tis.) | Urbanizácia | Leverage | |
|---|---|---|---|---|---|---|---|---|
| 4740 | Germany | 2013 | 9.51 | 6.88 | 39.94 | 48.68 | 88.11 | 0.001218 |
| 5427 | Japan | 1970 | 9.25 | 5.13 | 5.99 | 49.29 | 84.90 | 0.001195 |
| 6510 | China | 1984 | 9.70 | 7.02 | 5.60 | 4.09 | 89.54 | 0.001159 |
| 8992 | Brazil | 2021 | 8.63 | 61.57 | 38.98 | 0.74 | 11.17 | 0.001130 |
| 1328 | Japan | 2000 | 6.32 | 66.83 | 6.31 | 1.71 | 14.29 | 0.001109 |
priemery <- colMeans(udaje.model[c("spotreba_na_obyv", "obezita", "hdp_tis", "urbanizacia")])Najvyššiu leverage má záznam Germany (2013). Má takmer najnižšiu spotrebu (6,9 kg oproti priemeru 37,6 kg) a súčasne takmer maximálnu obezitu (39,9 %), HDP (48,7 tis. USD) aj urbanizáciu (88,1 %). Leží teda v rohu štvorrozmerného priestoru vysvetľujúcich premenných.
5.2 Štandardizované a studentizované rezíduá
Rezíduá \(\mathbf e=(\mathbf I-\mathbf H)\mathbf y\) nemajú rovnaký rozptyl:
\[ \operatorname{Var}(e_i\mid\mathbf X)=\sigma^2(1-h_{ii}). \]
Pozorovania s vysokou leverage priťahujú regresnú rovinu, a preto majú menšie rezíduá. Aby sa dali porovnávať, štandardizujem ich:
\[ \boxed{r_i=\frac{e_i}{\hat\sigma\sqrt{1-h_{ii}}}}\quad(\texttt{rstandard()}), \qquad t_i=\frac{e_i}{\hat\sigma_{(i)}\sqrt{1-h_{ii}}}\quad(\texttt{rstudent()}), \]
kde \(\hat\sigma_{(i)}\) je odhad z modelu bez \(i\)-teho pozorovania. Externe studentizované rezíduá \(t_i\) používa test odľahlých hodnôt.
r_rucne <- e / (sig * sqrt(1 - h))
sig_i <- sqrt((RSS - e^2 / (1 - h)) / (n - k - 1)) # sigma bez i-teho pozorovania
t_rucne <- e / (sig_i * sqrt(1 - h))
stopifnot(isTRUE(all.equal(unname(r_rucne), unname(rstandard(model)))),
isTRUE(all.equal(unname(t_rucne), unname(rstudent(model)))))
data.frame(
Ukazovatel = c("Počet |r_i| > 2", "Počet |r_i| > 3", "max |r_i|",
"Očakávaný počet |r_i| > 2 pri normalite"),
Hodnota = c(sum(abs(r_rucne) > 2), sum(abs(r_rucne) > 3),
cislo(max(abs(r_rucne)), 3), cislo(n * 2 * pnorm(-2), 0))
) |>
kable(col.names = c("Ukazovateľ", "Hodnota"), align = c("l", "r"),
caption = "Štandardizované rezíduá") |>
styl()| Ukazovateľ | Hodnota |
|---|---|
| Počet |r_i| > 2 | 0 |
| Počet |r_i| > 3 | 0 |
| max |r_i| | 1,764 |
| Očakávaný počet |r_i| > 2 pri normalite | 455 |
Ručný výpočet sa zhoduje s rstandard() aj rstudent(). Pri normálnom rozdelení by približne 4,6 % pozorovaní, teda asi 455 záznamov, malo \(|r_i|>2\). V našom modeli takých záznamov nie je ani jeden, čo opäť potvrdzuje ľahké chvosty. Faktor \(\sqrt{1-h_{ii}}\) je tu takmer 1, preto sa \(r_i\) od \(e_i/\hat\sigma\) prakticky nelíšia.
5.3 Cookova vzdialenosť
Cookova vzdialenosť spája veľkosť rezídua a leverage. Meria, ako by sa zmenili všetky vyrovnané hodnoty, keby sme \(i\)-te pozorovanie vynechali:
\[ \boxed{D_i=\frac{\sum_j(\hat y_j-\hat y_{j(i)})^2}{k\,\hat\sigma^2} =\frac{e_i^2}{k\,\hat\sigma^2}\,\frac{h_{ii}}{(1-h_{ii})^2} =\frac{r_i^2}{k}\,\frac{h_{ii}}{1-h_{ii}}.} \]
D <- r_rucne^2 / k * h / (1 - h)
stopifnot(isTRUE(all.equal(unname(D), unname(cooks.distance(model)))))
i_D <- which.max(D)par(mfrow = c(1, 2))
plot(h, type = "h", col = farba, xlab = "Záznam", ylab = "Leverage h_ii",
main = "Pákový účinok")
abline(h = c(k/n, 2*k/n), lty = c(1, 2), col = c("grey40", "red"))
legend("topright", c("k/n", "2k/n"), lty = c(1, 2), col = c("grey40", "red"),
bty = "n", cex = .8)
plot(D, type = "h", col = farba, xlab = "Záznam", ylab = "Cookova vzdialenosť D_i",
main = "Vplyv pozorovaní")
abline(h = 4/n, lty = 2, col = "red")
legend("topright", "4/n", lty = 2, col = "red", bty = "n", cex = .8)
par(mfrow = c(1, 1))Najväčšia Cookova vzdialenosť je \(D_{max}=0,000589\) (záznam India, 2014). Orientačný prah \(4/n=0,0004\) prekračuje 113 záznamov. Prahy 0,5 a 1, ktoré kreslí plot(model), sú však vzdialené o tri rády. Žiadne pozorovanie teda nemá podstatný vplyv. Prah \(4/n\) pri veľkých vzorkách klesá k nule a označí vždy niekoľko percent záznamov, preto ho beriem len ako upozornenie.
5.4 Je vysoká leverage to isté ako vplyvné pozorovanie?
Nie. Záznam s najvyššou leverage (Germany, 2013) má štandardizované rezíduum len 0,25, preto jeho Cookova vzdialenosť 0,000016 nie je najväčšia. Naopak, najvplyvnejší záznam (India, 2014) kombinuje nadpriemernú leverage (0,00103, t. j. 2,1-násobok priemeru) s veľkým rezíduom (\(r_i=1,69\)). Vplyv vzniká až kombináciou neobvyklého \(\mathbf x_i\) a neobvyklého \(y_i\).
data.frame(h = h, r = r_rucne, D = D) |>
ggplot(aes(x = h, y = r, size = D)) +
geom_point(alpha = 0.25, color = farba) +
geom_point(data = \(d) d[c(i_h, i_D), ], color = "firebrick", alpha = 1) +
annotate("text", x = h[c(i_h, i_D)], y = r_rucne[c(i_h, i_D)],
label = c("najvyššia leverage", "najväčšie D"),
hjust = 1.1, vjust = -0.6, color = "firebrick", size = 3.5) +
geom_vline(xintercept = 2*k/n, linetype = 2) +
scale_size_area(max_size = 4, guide = "none") +
labs(x = "Leverage h_ii", y = "Štandardizované rezíduum r_i") +
theme_minimal(base_size = 13)
5.4.1 Citlivosť odhadov na vynechanie pozorovaní
Ako najprísnejšiu skúšku vynechám naraz všetkých 113 záznamov s \(D_i>4/n\) a porovnám koeficienty:
model_bez <- lm(formula(model), data = udaje.model[D <= 4/n, ])
data.frame(
Parameter = names(b),
Vsetky = b, Bez = coef(model_bez),
Zmena = 100 * (coef(model_bez) - b) / abs(b),
SE = se
) |>
kable(digits = c(0, 5, 5, 1, 5), row.names = FALSE,
col.names = c("Parameter", "Všetky záznamy", "Bez záznamov D > 4/n",
"Zmena (%)", "SE pôvodného odhadu"),
caption = "Koeficienty s vplyvnými záznamami a bez nich") |>
styl()| Parameter | Všetky záznamy | Bez záznamov D > 4/n | Zmena (%) | SE pôvodného odhadu |
|---|---|---|---|---|
| (Intercept) | 8.45873 | 8.43855 | -0.2 | 0.15789 |
| spotreba_na_obyv | 0.02471 | 0.02455 | -0.6 | 0.00201 |
| obezita | -0.00165 | -0.00072 | 56.5 | 0.00373 |
| hdp_tis | -0.00230 | -0.00318 | -38.4 | 0.00265 |
| urbanizacia | 0.00111 | 0.00170 | 52.8 | 0.00163 |
Koeficient spotreby cukru sa zmenil z 0,02471 na 0,02455. Rozdiel je menší ako jeho smerodajná chyba (0,00201), takže záver o spotrebe cukru nezávisí od niekoľkých záznamov. Percentuálne zmeny pri nevýznamných koeficientoch sú väčšie len preto, že ich odhady ležia blízko nuly.
6 Testy normality rezíduí
Náhodné zložky \(\mathbf u\) nepozorujeme. Rezíduá sú však ich lineárnou transformáciou:
\[ \mathbf e=\mathbf y-\mathbf H\mathbf y=(\mathbf I-\mathbf H)(\mathbf X\boldsymbol\beta+\mathbf u) =(\mathbf I-\mathbf H)\mathbf u , \]
pretože \((\mathbf I-\mathbf H)\mathbf X=\mathbf X-\mathbf X=\mathbf 0\). Lineárna kombinácia normálne rozdelených veličín je opäť normálne rozdelená. Ak sú \(\mathbf u\) normálne, sú normálne aj \(\mathbf e\). Normalitu nepozorovateľných chýb preto testujem sprostredkovane cez rezíduá.
6.1 Jarqueov–Berov test
Test porovnáva výberovú šikmosť \(S\) a špicatosť \(K\) rezíduí s hodnotami normálneho rozdelenia (\(S=0\), \(K=3\)):
\[ H_0:\text{rezíduá sú normálne rozdelené},\qquad JB=\frac n6\left(S^2+\frac{(K-3)^2}{4}\right)\overset{a}{\sim}\chi^2(2). \]
# Jarqueov–Berov test priamo podľa vzorca
# (dáva rovnaký výsledok ako tseries::jarque.bera.test())
jarque_bera <- function(x) {
xc <- x - mean(x)
S <- mean(xc^3) / mean(xc^2)^1.5 # výberová šikmosť
K <- mean(xc^4) / mean(xc^2)^2 # výberová špicatosť
JB <- length(x) / 6 * (S^2 + (K - 3)^2 / 4)
list(S = S, K = K, statistic = JB,
p.value = pchisq(JB, df = 2, lower.tail = FALSE))
}
jb_test <- jarque_bera(e)
S <- jb_test$S
K <- jb_test$K
JB <- jb_test$statistic
data.frame(S = S, K = K, JB = JB, krit = qchisq(0.95, 2),
p = p_hod(jb_test$p.value)) |>
kable(digits = 4,
col.names = c("Šikmosť S", "Špicatosť K", "JB štatistika",
"Kritická hodnota χ²(2)", "p-hodnota"),
caption = "Jarqueov–Berov test normality rezíduí") |>
styl()| Šikmosť S | Špicatosť K | JB štatistika | Kritická hodnota χ²(2) | p-hodnota |
|---|---|---|---|---|
| 0.0049 | 1.7817 | 618.4602 | 5.9915 | < 0,001 |
Interpretácia: \(JB=618,5\) výrazne prekračuje kritickú hodnotu 5,991 a \(p < 0,001\). Zamietam hypotézu o normalite rezíduí. Rozklad štatistiky ukazuje príčinu:
- šikmosť \(S=0,0049\) je prakticky nulová – rozdelenie je symetrické, člen \(S^2\) prispieva k JB len 0,04,
- špicatosť \(K=1,782\) je výrazne nižšia ako 3. Rovnomerné rozdelenie má \(K=9/5=1{,}8\), takže rezíduá sú svojím tvarom takmer rovnomerné. Celú hodnotu JB tvorí druhý člen.
Výsledok testu sa zhoduje s obrazom Q–Q grafu (tvar S – ľahké chvosty).
set.seed(123)
rovn <- runif(n, -sqrt(3) * sig, sqrt(3) * sig) # rovnomerné s rovnakým rozptylom
par(mfrow = c(1, 2))
hist(e, breaks = 40, freq = FALSE, col = "lightblue", border = "white",
main = "Histogram rezíduí", xlab = "Rezíduum (p. b.)")
curve(dnorm(x, 0, sig), add = TRUE, col = "red", lwd = 2)
qqnorm(e, pch = ".", col = farba, main = "Q–Q graf: rezíduá vs. rovnomerné")
qqline(e, col = "red", lwd = 2)
q <- qqnorm(sort(rovn), plot.it = FALSE)
lines(q$x, q$y, col = "orange", lwd = 2, lty = 2)
legend("topleft", c("rezíduá modelu", "rovnomerné rozdelenie"),
col = c(farba, "orange"), lty = c(1, 2), lwd = 2, bty = "n", cex = .8)
par(mfrow = c(1, 1))Histogram je plochý – namiesto zvonu tvorí takmer obdĺžnik. Rezíduá v Q–Q grafe sa takmer kryjú s krivkou rovnomerného rozdelenia s rovnakým rozptylom. Vysvetlenie je jednoduché: aj samotná premenná diabetes je v údajoch rozložená takmer rovnomerne (od 2,1 do 16,5 %, špicatosť 1,82). Model z nej vysvetľuje len 1,5 %, preto rezíduá zdedia jej tvar.
PoznámkaJe porušenie normality vážny problém?
Normalita chýb je potrebná na presné rozdelenie t- a F-štatistík v malých vzorkách. Pri \(n=10000\) zaručuje centrálna limitná veta, že odhady \(\hat\beta_j\) sú približne normálne aj bez normálnych chýb. Ľahké chvosty sú navyše neškodná odchýlka – neprodukujú extrémne hodnoty. Závery t-testov a F-testu preto považujem za spoľahlivé. Jarqueov–Berov test je pri veľkých vzorkách veľmi citlivý a zamietne normalitu aj pri malých odchýlkach. Preto ho vždy dopĺňam Q–Q grafom. Ďalšie testy normality sú napríklad Shapirov–Wilkov, Andersonov–Darlingov alebo Kolmogorovov–Smirnovov.
7 Test odľahlých pozorovaní – outlierTest()
Funkcia car::outlierTest() (alebo rovnaký výpočet priamo podľa vzorca, ktorý používam nižšie) testuje najväčšie externe studentizované rezíduum
\[ t_i=\frac{e_i}{\hat\sigma_{(i)}\sqrt{1-h_{ii}}}\overset{H_0}{\sim}t(n-k-1), \]
kde \(H_0\): pozorovanie nie je odľahlé vzhľadom na model. Keďže tým vlastne testujem všetkých \(n\) pozorovaní naraz, p-hodnoty sa upravia Bonferroniho korekciou \(p_{i,\mathrm{Bonf}}=\min(1,\,n\cdot p_i)\). Ak \(p_{i,\mathrm{Bonf}}<0{,}05\), pozorovanie považujem za štatisticky významne odľahlé.
# rovnaký výpočet ako car::outlierTest(), bez potreby balíka car
tabulka_outlier <- function(m) {
t_i <- rstudent(m) # externe studentizované rezíduá
i <- which.max(abs(t_i)) # najextrémnejšie pozorovanie
df <- df.residual(m) - 1 # n - k - 1
p <- 2 * pt(abs(t_i[i]), df = df, lower.tail = FALSE)
data.frame(
Zaznam = names(t_i)[i],
Krajina = paste0(udaje.model$country[i], " (", udaje.model$year[i], ")"),
rstudent = unname(t_i[i]),
p = unname(p),
bonf = min(1, nobs(m) * unname(p)) # Bonferroniho korekcia
)
}
ot1 <- tabulka_outlier(model)
ot1 |>
kable(digits = 4, row.names = FALSE,
col.names = c("Záznam", "Krajina", "Studentizované rezíduum",
"Neupravená p-hodnota", "Bonferroniho p-hodnota"),
caption = "Test odľahlých pozorovaní s Bonferroniho korekciou") |>
styl()| Záznam | Krajina | Studentizované rezíduum | Neupravená p-hodnota | Bonferroniho p-hodnota |
|---|---|---|---|---|
| 7458 | Indonesia (2006) | 1.7642 | 0.0777 | 1 |
| Stĺpec | Význam |
|---|---|
| Záznam | Poradové číslo (riadok) záznamu s najväčším \(|t_i|\). |
| Studentizované rezíduum | \(t_i\) – rozptyl odhadnutý bez tohto pozorovania. |
| Neupravená p-hodnota | p-hodnota jedného t-testu, bez ohľadu na to, že testujeme \(n\) pozorovaní. |
| Bonferroniho p-hodnota | \(\min(1, n\cdot p_i)\) – podľa nej rozhodujem. |
Interpretácia: Tabuľka ukazuje záznam s najväčším \(|t_i|=1,764\). Jeho neupravená p-hodnota 0,0777 by síce bola hraničná, po korekcii na 10000 testov je však \(n\cdot p_i=777>1\), teda Bonferroniho p-hodnota je 1. (Funkcia car::outlierTest() by v takom prípade vypísala NA.) V modeli nie sú žiadne odľahlé pozorovania. Pri ľahkých chvostoch rezíduí to ani nemožno čakať: aby bol záznam pri 10000 pozorovaniach významne odľahlý, potreboval by \(|t_i|\) aspoň okolo 4,57.
Poznámka
Odľahlé pozorovanie nemusí byť chybné a nemali by sme ho automaticky odstraňovať. Test hodnotí len neobvyklosť \(y_i\) vzhľadom na model; na posúdenie vplyvu na koeficienty slúži Cookova vzdialenosť.
8 Úprava modelu
Diagnostika ukázala dva poznatky:
- tri vysvetľujúce premenné (obezita, HDP, urbanizácia) sú nevýznamné,
- problém s normalitou nespôsobujú odľahlé hodnoty ani nelinearita, ale tvar rozdelenia vysvetľovanej premennej.
Logaritmická transformácia, ktorá pomáha pri nelinearite a pri odľahlých hodnotách v pravom chvoste, tu preto nie je opodstatnená. Rozumnou úpravou je vynechať nevýznamné premenné:
\[ Diabetes_i=\beta_0+\beta_1\,Spotreba_i+u_i . \]
8.1 F-test spoločnej nevýznamnosti vynechaných premenných
Predtým, než premenné vynechám, otestujem ich spoločne (jednotlivé t-testy nestačia):
\[ H_0:\beta_2=\beta_3=\beta_4=0,\qquad F=\frac{(RSS_R-RSS_U)/q}{RSS_U/(n-k)}\overset{H_0}{\sim}F(q,\,n-k), \]
kde \(RSS_R\) je súčet štvorcov rezíduí obmedzeného (menšieho) modelu, \(RSS_U\) neobmedzeného modelu a \(q=3\) počet obmedzení.
model2 <- lm(diabetes ~ spotreba_na_obyv, data = udaje.model)
test_obm <- anova(model2, model)
test_obm |>
as.data.frame() |>
mutate(`Pr(>F)` = ifelse(is.na(`Pr(>F)`), NA, p_hod(`Pr(>F)`))) |>
kable(digits = 4, row.names = FALSE,
col.names = c("Rez. df", "RSS", "q", "Rozdiel RSS", "F", "p-hodnota"),
caption = "F-test: obmedzený model (1. riadok) vs. pôvodný model (2. riadok)") |>
styl()| Rez. df | RSS | q | Rozdiel RSS | F | p-hodnota |
|---|---|---|---|---|---|
| 9998 | 141874.6 | ||||
| 9995 | 141854.4 | 3 | 20.2475 | 0.4755 | 0,699 |
\(F=0,476\) s p-hodnotou 0,699. Nulovú hypotézu nezamietam – obezita, HDP a urbanizácia ani spoločne neprispievajú k vysvetleniu výskytu cukrovky. Ich vynechanie je štatisticky opodstatnené.
tidy(model2) |>
mutate(p.value = p_hod(p.value)) |>
kable(digits = 5,
col.names = c("Premenná", "Odhad", "Smer. chyba", "t hodnota", "p-hodnota"),
caption = "Upravený model: výskyt cukrovky a spotreba cukru") |>
styl()| Premenná | Odhad | Smer. chyba | t hodnota | p-hodnota |
|---|---|---|---|---|
| (Intercept) | 8.41877 | 0.08452 | 99.60484 | < 0,001 |
| spotreba_na_obyv | 0.02473 | 0.00201 | 12.29605 | < 0,001 |
Koeficient spotreby (0,02473) sa oproti pôvodnému modelu (0,02471) prakticky nezmenil. Vynechané premenné totiž so spotrebou nekorelujú, ako som zistila minule.
8.2 Grafická diagnostika upraveného modelu
par(mfrow = c(2, 2))
plot(model2)
par(mfrow = c(1, 1))Grafy sú takmer totožné s pôvodnými. Pribudla len jedna zaujímavosť v grafe Residuals vs Leverage: pri jedinej vysvetľujúcej premennej je leverage priamo funkciou vzdialenosti spotreby od priemeru, \(h_{ii}=\frac1n+\frac{(x_i-\bar x)^2}{\sum_j(x_j-\bar x)^2}\). Body preto začínajú presne na hodnote \(1/n=0,0001\) (záznamy so spotrebou blízko priemeru) a najvyššiu leverage majú záznamy s najnižšou a najvyššou spotrebou. Keďže je spotreba rozložená takmer rovnomerne, maximum je len 0,00040, teda približne \(4/n\).
8.3 Normalita a odľahlé hodnoty upraveného modelu
jb2 <- jarque_bera(residuals(model2))
ot2 <- tabulka_outlier(model2)
data.frame(
Test = c("Jarque–Bera: JB štatistika", "Jarque–Bera: p-hodnota",
"outlierTest: max |t_i|", "outlierTest: Bonferroniho p-hodnota"),
Povodny = c(cislo(jb_test$statistic, 2), p_hod(jb_test$p.value),
cislo(abs(ot1$rstudent[1]), 3), cislo(ot1$bonf[1], 3)),
Upraveny = c(cislo(jb2$statistic, 2), p_hod(jb2$p.value),
cislo(abs(ot2$rstudent[1]), 3), cislo(ot2$bonf[1], 3))
) |>
kable(col.names = c("Test", "Pôvodný model", "Upravený model"),
align = c("l", "r", "r"),
caption = "Testy normality a odľahlých hodnôt") |>
styl()| Test | Pôvodný model | Upravený model |
|---|---|---|
| Jarque–Bera: JB štatistika | 618,46 | 618,72 |
| Jarque–Bera: p-hodnota | < 0,001 | < 0,001 |
| outlierTest: max |t_i| | 1,764 | 1,742 |
| outlierTest: Bonferroniho p-hodnota | 1,000 | 1,000 |
Normalitu rezíduí zamietam aj v upravenom modeli a odľahlé pozorovania opäť nie sú. Vynechanie premenných diagnostiku nezhoršilo, model je však jednoduchší pri rovnakej vysvetľovacej schopnosti.
9 Malé úlohy
9.1 Porovnanie modelov pomocou AIC
Okrem pôvodného a upraveného modelu som vyskúšala ďalšie dve modifikácie, aby som overila, či nepomôže iný funkčný tvar spotreby:
- logaritmický model \(Diabetes_i=\beta_0+\beta_1\ln(Spotreba_i)+u_i\), kde \(\beta_1/100\) je zmena výskytu cukrovky v p. b. pri raste spotreby o 1 %,
- kvadratický model \(Diabetes_i=\beta_0+\beta_1 Spotreba_i+\beta_2 Spotreba_i^2+u_i\), ktorý testuje nelinearitu.
model3 <- lm(diabetes ~ log(spotreba_na_obyv), data = udaje.model)
model4 <- lm(diabetes ~ spotreba_na_obyv + I(spotreba_na_obyv^2), data = udaje.model)
modely <- list(
"1: pôvodný (4 premenné)" = model,
"2: len spotreba" = model2,
"3: ln(spotreba)" = model3,
"4: spotreba + spotreba²" = model4
)
porovnanie_modelov <- lapply(modely, glance) |>
bind_rows(.id = "Model") |>
select(Model, df, r.squared, adj.r.squared, sigma, AIC, BIC) |>
mutate(dAIC = AIC - min(AIC), dBIC = BIC - min(BIC))
porovnanie_modelov |>
kable(digits = c(0, 0, 5, 5, 4, 2, 2, 2, 2),
col.names = c("Model", "Premenné", "R²", "Upravené R²", "Sigma",
"AIC", "BIC", "ΔAIC", "ΔBIC"),
caption = "Porovnanie modelov (ΔAIC, ΔBIC = rozdiel oproti najlepšiemu)") |>
styl() |>
row_spec(which.min(porovnanie_modelov$AIC), bold = TRUE, background = "#e3f1e6")| Model | Premenné | R² | Upravené R² | Sigma | AIC | BIC | ΔAIC | ΔBIC |
|---|---|---|---|---|---|---|---|---|
| 1: pôvodný (4 premenné) | 4 | 0.01504 | 0.01464 | 3.7673 | 54912.93 | 54956.19 | 4.57 | 26.20 |
| 2: len spotreba | 1 | 0.01490 | 0.01480 | 3.7670 | 54908.36 | 54929.99 | 0.00 | 0.00 |
| 3: ln(spotreba) | 1 | 0.01343 | 0.01333 | 3.7698 | 54923.24 | 54944.87 | 14.89 | 14.89 |
| 4: spotreba + spotreba² | 2 | 0.01490 | 0.01470 | 3.7672 | 54910.34 | 54939.19 | 1.99 | 9.20 |
naj <- porovnanie_modelov$Model[which.min(porovnanie_modelov$AIC)]
p_kv <- tidy(model4)$p.value[3]Odpoveď: Podľa Akaikeho informačného kritéria je najlepší model 2: len spotreba s najnižším AIC. Pôvodný model má AIC vyššie o 4,57. Tri pridané premenné teda zvýšili \(R^2\) len nepatrne, menej, ako je penalizácia \(2\times3=6\) za tri parametre navyše. BIC, ktoré penalizuje zložitosť ešte prísnejšie (\(\ln n\approx9,21\) za parameter), dáva rovnaké poradie s ešte väčším rozdielom (26,20). Pravidlo palca: rozdiel AIC väčší ako 2 už znamená citeľnú prednosť modelu s nižším AIC.
9.2 Iná modifikácia modelu
- Logaritmus spotreby model zhoršil (ΔAIC = 14,89). Údaje teda nepodporujú predstavu, že s rastúcou spotrebou sa prírastok výskytu cukrovky zmenšuje (to by predpokladal logaritmus).
- Kvadratický člen je nevýznamný (\(p=0,912\)) a AIC je vyššie o 1,99. Potvrdzuje to, čo ukázal graf Residuals vs Fitted: vzťah je lineárny.
Ani jedna modifikácia teda nezlepšila technické charakteristiky. Podobne ako vo vzore, kde logaritmická transformácia HDP nepomohla, aj tu sa ukazuje, že transformácia má zmysel len vtedy, keď ju podporuje diagnostika.
10 Najdôležitejšie zistenia
Odhadla som viacnásobný regresný model výskytu cukrovky a podrobila ho diagnostike:
- t-testy: štatisticky významná je len spotreba cukru. Obezita, HDP a urbanizácia sú nevýznamné jednotlivo aj spoločne (F-test obmedzení).
- Kvalita vyrovnania: model ako celok je významný (F-test), vysvetľuje však len asi 1,5 % variability (\(R^2\)).
- Grafická analýza: linearita a homoskedasticita sú v poriadku. Q–Q graf ukazuje ľahké chvosty.
- Hat matica: všetky leverage sú veľmi malé (max. 0,00122), Cookove vzdialenosti sú o tri rády pod kritickými hodnotami. Ani vynechanie všetkých záznamov s \(D_i>4/n\) nezmení závery – vplyvné pozorovania v modeli nie sú.
- Normalita: Jarqueov–Berov test normalitu zamieta. Príčinou je nízka špicatosť (\(K\approx1,78\), takmer rovnomerné rozdelenie), nie šikmosť ani extrémne hodnoty. Pri 10000 pozorovaniach to spoľahlivosť testov výrazne neohrozuje.
- Odľahlé hodnoty: žiadne (Bonferroniho p-hodnota = 1).
- Úprava modelu: najlepší podľa AIC aj BIC je jednoduchý model so spotrebou cukru. Logaritmická ani kvadratická transformácia ho nezlepšili.
Diagnostika teda neodhalila chybu v údajoch ani v špecifikácii. Hlavným obmedzením modelu je jeho nízka vysvetľovacia schopnosť: výskyt cukrovky v týchto údajoch takmer nesúvisí s premennými, ktoré mám k dispozícii.