1. Úvod

Po autokorelácii a heteroskedasticite rezíduí je multikolinearita tretím závažným porušením predpokladov použitia metódy najmenších štvorcov.

Pri klasickej lineárnej regresii predpokladáme okrem iného, že matica vysvetľujúcich premenných X má lineárne nezávislé stĺpce. To znamená, že žiadny stĺpec (žiadna vysvetľujúca premenná) sa nedá presne vyjadriť ako lineárna kombinácia ostatných. Tento predpoklad zabezpečuje, že matica X’X je regulárna (má inverziu), a teda vieme vypočítať odhady regresných koeficientov beta = (X’X)^(-1) X’y.

V praxi sa však často nestane, že by bola kombinácia premenných presne lineárne závislá – skôr ide o situáciu, keď sú niektoré premenné „takmer“ lineárne závislé. Vtedy hovoríme o multikolinearite. V takom prípade je matica (X’X)^(-1) veľmi nestabilná, jej prvky na hlavnej diagonále sú veľké a to sa prejaví na veľkých štandardných chybách odhadovaných koeficientov.

Multikolinearita teda neskreslí priemer odhadov, ale zhorší ich presnosť a stabilitu.

2. Dôsledky multikolinearity

Hlavné dôsledky multikolinearity:

Nezavádza bias – odhady koeficientov sú v priemere správne.
Zvyšuje štandardné chyby koeficientov – intervaly sú širokejšie, p-hodnoty väčšie.
Odhady sú nestabilné – malé zmeny v dátach vedú k veľkým zmenám koeficientov.
Zhoršuje interpretáciu modelu – ťažko sa vysvetľuje samostatný vplyv jednej premennej.

3. Východiskový model a údaje

Budeme pracovať s databázou 250 najlepšie hodnotených filmov z IMDb (súbor data_ekonometria.csv). Zaujíma nás, ako sa hodnotenie IMDb (0–10) dá vysvetliť pomocou týchto premenných:

rok uvedenia filmu,
dĺžka filmu v minútach,
rozpočet filmu v dolároch.

Model:

IMDb = beta0 + beta1 * rok + beta2 * dlzka + beta3 * rozpocet + chyba

Najprv načítame a pripravíme údaje:

udaje <- read.csv("data_ekonometria.csv",
                  dec = ",", sep = ";", header = TRUE)

# premenovanie stĺpcov na jednoduchšie názvy
names(udaje) <- c("umiestnenie", "nazov", "rok", "zaner",
                  "imdb", "dlzka", "rezia", "spolocnost", "rozpocet_text")

# prevod rozpočtu z textu s medzerami na číselnú premennú v dolároch
udaje$rozpocet <- as.numeric(gsub(" ", "", udaje$rozpocet_text))

# vyber premenných, s ktorými budeme pracovať v regresii
udaje <- udaje[, c("imdb", "rok", "dlzka", "rozpocet")]

# jednoduchá imputácia – prípadné chýbajúce hodnoty nahradíme mediánom
column_medians <- sapply(udaje, median, na.rm = TRUE)
for (col in names(udaje)) {
  udaje[[col]][is.na(udaje[[col]])] <- column_medians[col]
}

nrow(udaje)

## [1] 250

head(udaje)

##   imdb  rok dlzka rozpocet
## 1  9.3 1994   142 2.50e+07
## 2  9.2 1972   175 6.00e+06
## 3  9.0 2008   152 1.85e+08
## 4  9.0 1974   202 1.30e+07
## 5  9.0 1957    96 3.37e+05
## 6  9.0 1993   195 2.20e+07

Komentár:

v databáze je 250 filmov – ide o prierezové údaje,
imdb je hodnotenie, typicky medzi 8 a 9.5,
rok sa pohybuje približne od polovice 20. storočia po súčasnosť,
dlzka je v minútach, väčšinou 90–200,
rozpocet je v dolároch, rádovo milióny až stovky miliónov.

4. Odhad základného regresného modelu

model <- lm(imdb ~ rok + dlzka + rozpocet,
            data = udaje)
summary(model)

## 
## Call:
## lm(formula = imdb ~ rok + dlzka + rozpocet, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42001 -0.15545 -0.05865  0.13091  0.98723 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.678e+00  1.238e+00   7.010 2.28e-11 ***
## rok         -3.391e-04  6.267e-04  -0.541    0.589    
## dlzka        2.126e-03  4.689e-04   4.535 9.00e-06 ***
## rozpocet     3.424e-10  2.624e-10   1.305    0.193    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2218 on 246 degrees of freedom
## Multiple R-squared:  0.08716,    Adjusted R-squared:  0.07603 
## F-statistic: 7.829 on 3 and 246 DF,  p-value: 5.186e-05

Interpretácia základného modelu:

R-squared je približne 0.087 – model vysvetľuje asi 8.7 % rozdielov v hodnotení IMDb. Nie je to veľa, ale vzhľadom na to, že máme len tri jednoduché vysvetľujúce premenné, je to realistické.
F-test pre celý model má veľmi malú p-hodnotu – aspoň jedna z premenných (rok, dlzka, rozpocet) má štatisticky významný vplyv na IMDb.
Koeficient pri dlzka je kladný, štatisticky významný: dlhšie filmy majú v priemere o niečo vyššie hodnotenie. Napríklad nárast dĺžky o 60 minút zvýši očakávané hodnotenie približne o 0.12 bodu.
Koeficient pri rok je blízko nule a štatisticky nevýznamny – po zohľadnení dĺžky a rozpočtu rok uvedenia filmu už nehrá veľkú rolu.
Koeficient pri rozpocet je extrémne malý (lebo pracujeme v dolároch) a tiež nie je štatisticky významny – vo vzorke top 250 filmov nevidíme jasný vzťah, že vyšší rozpočet znamená lepšie hodnotenie.
Intercept (konštanta) má hodnotu okolo 8.7 – berieme ho ako technický parameter, nie ako reálne „hodnotenie filmu z roku 0, s dĺžkou 0 a rozpočtom 0“.

Zhrnutie: najdôležitejšou premennou v tomto modeli je dĺžka filmu, rok a rozpočet už veľa navyše nevysvetľujú.

5. Korelačná matica a párové vzťahy

xvars <- udaje[, c("rok", "dlzka", "rozpocet")]
round(cor(xvars), 3)

##            rok dlzka rozpocet
## rok      1.000 0.111    0.452
## dlzka    0.111 1.000    0.110
## rozpocet 0.452 0.110    1.000

Interpretácia korelačnej matice:

rok a rozpocet: mierne kladná korelácia – novšie filmy majú často vyšší rozpočet,
rok a dlzka: korelácia je nízka – novšie filmy nie sú zásadne dlhšie či kratšie,
dlzka a rozpocet: korelácia je skôr slabá – drahšie filmy nie sú automaticky dlhšie.

Nevidíme extrémne korelácie typu 0.9, takže multikolinearita nie je na prvý pohľad zrejmá len z korelačnej matice.

pairs(xvars,
      main = "Scatterplotová matica – rok, dlzka, rozpocet")

Zo scatterplotov vidíme:

dĺžka filmu má mierne rastúci vzťah k IMDb (dlhšie filmy majú mierne vyššie hodnotenie),
vzťah rozpočtu k IMDb je veľmi rozptýlený – existujú lacnejšie, ale veľmi dobre hodnotené filmy.

6. VIF – Variance Inflation Factor

library(car)
vif(model)

##      rok    dlzka rozpocet 
## 1.262283 1.017010 1.262082

Interpretácia VIF:

VIF približne 1–2 znamená, že multikolinearita nie je silná z pohľadu danej premennej.
Pri našom modeli VIF vychádza okolo 1–1.3 – to je veľmi nízke.
To naznačuje, že žiadna premenná nie je extrémne lineárne závislá od ostatných.

Záver: VIF nám neindikuje vážny problém multikolinearity. To však ešte neznamená, že nemôže existovať problém s číselnou stabilitou kvôli rozdielnym rádom premenných – to zachytíme Condition Number.

7. Condition Number

X <- model.matrix(model)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)

condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number

## [1] 2345673

Pravidlo:

menej ako 10 – nízka multikolinearita,
10–30 – mierna,
30–100 – silná,
viac ako 100 – veľmi vážna.

V našom prípade je Condition Number veľmi vysoký (rádovo státisíce až milióny). To znamená:

matica X’X je zle podmienená,
malé zmeny údajov môžu spôsobiť výrazné zmeny v koeficientoch,
rozdielne rády premenných (rozpočet v desiatkach až stovkách miliónov, rok okolo 2000, dlzka okolo 100) spôsobujú numerickú nestabilitu.

VIF teda nie je problém, ale Condition Number ukazuje na číselnú nestabilitu modelu.

8. Riešenia multikolinearity a nestability

8.1 Vynechanie niektorej premennej

model_noRok <- lm(imdb ~ dlzka + rozpocet, data = udaje)
summary(model_noRok)

## 
## Call:
## lm(formula = imdb ~ dlzka + rozpocet, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41837 -0.16151 -0.06026  0.12843  0.98413 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.009e+00  6.216e-02 128.859  < 2e-16 ***
## dlzka       2.109e-03  4.671e-04   4.515 9.81e-06 ***
## rozpocet    2.793e-10  2.347e-10   1.190    0.235    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2214 on 247 degrees of freedom
## Multiple R-squared:  0.08607,    Adjusted R-squared:  0.07867 
## F-statistic: 11.63 on 2 and 247 DF,  p-value: 1.488e-05

model_noDlzka <- lm(imdb ~ rok + rozpocet, data = udaje)
summary(model_noDlzka)

## 
## Call:
## lm(formula = imdb ~ rok + rozpocet, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29853 -0.18305 -0.08277  0.11290  1.01132 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.565e+00  1.286e+00   6.661 1.75e-10 ***
## rok         -1.438e-04  6.495e-04  -0.221    0.825    
## rozpocet     4.228e-10  2.720e-10   1.554    0.121    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2304 on 247 degrees of freedom
## Multiple R-squared:  0.01084,    Adjusted R-squared:  0.00283 
## F-statistic: 1.353 on 2 and 247 DF,  p-value: 0.2603

model_noRozpocet <- lm(imdb ~ rok + dlzka, data = udaje)
summary(model_noRozpocet)

## 
## Call:
## lm(formula = imdb ~ rok + dlzka, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42572 -0.16175 -0.05905  0.12769  0.98021 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.963e+00  1.111e+00   7.165 8.94e-12 ***
## rok         2.476e-05  5.620e-04   0.044    0.965    
## dlzka       2.168e-03  4.685e-04   4.627 5.99e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2221 on 247 degrees of freedom
## Multiple R-squared:  0.08084,    Adjusted R-squared:  0.0734 
## F-statistic: 10.86 on 2 and 247 DF,  p-value: 3.011e-05

Typický obraz:

model bez dĺžky má výrazne nižší R-squared – dĺžka je kľúčová premenná,
model bez roku má podobný R-squared ako pôvodný model – rok veľa nepridáva,
model bez rozpočtu má tiež podobný R-squared – rozpočet tiež nepridáva veľa informácie navyše k dĺžke a roku.

Ak chceme jednoduchší model, dáva zmysel ponechať dĺžku (najdôležitejšia premenná) a pokojne vynechať rok alebo rozpočet, pretože ich prínos je malý.

8.2 Škálovanie premenných (standardizácia)

Teraz odstránime problém rozdielnych rádov premenných tým, že ich štandardizujeme (odčítame priemer a delíme smerodajnou odchýlkou).

udaje$rok_c       <- scale(udaje$rok,      center = TRUE, scale = TRUE)
udaje$dlzka_c     <- scale(udaje$dlzka,    center = TRUE, scale = TRUE)
udaje$rozpocet_c  <- scale(udaje$rozpocet, center = TRUE, scale = TRUE)

model_centered <- lm(imdb ~ rok_c + dlzka_c + rozpocet_c,
                     data = udaje)
summary(model_centered)

## 
## Call:
## lm(formula = imdb ~ rok_c + dlzka_c + rozpocet_c, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42001 -0.15545 -0.05865  0.13091  0.98723 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.294400   0.014025 591.400   <2e-16 ***
## rok_c       -0.008543   0.015789  -0.541    0.589    
## dlzka_c      0.064272   0.014172   4.535    9e-06 ***
## rozpocet_c   0.020600   0.015788   1.305    0.193    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2218 on 246 degrees of freedom
## Multiple R-squared:  0.08716,    Adjusted R-squared:  0.07603 
## F-statistic: 7.829 on 3 and 246 DF,  p-value: 5.186e-05

vif(model_centered)

##      rok_c    dlzka_c rozpocet_c 
##   1.262283   1.017010   1.262082

Interpretácia:

R-squared ostáva prakticky rovnaké – vysvetlená variabilita sa nezmenila,
koeficienty teraz hovoria o zmene IMDb pri zmene premennej o 1 smerodajnú odchýlku,
VIF sú stále nízke (okolo 1–1.3).

Skontrolujeme Condition Number po škálovaní:

X <- model.matrix(model_centered)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)

condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number

## [1] 1.654121

Teraz je Condition Number blízko hodnoty 1–2 → matica je dobre podmienená, číselná stabilita je výrazne lepšia.

Nevýhoda: koeficienty sú menej intuitívne (pracujeme v jednotkách „1 smerodajná odchýlka“).

8.3 Transformácia rozpočtu na milióny dolárov

Ak chceme zachovať interpretáciu v reálnych jednotkách, môžeme upraviť len rozpočet: namiesto dolárov použijeme milióny dolárov.

udaje$rozpocet_mil <- udaje$rozpocet / 1e6
head(udaje)

##   imdb  rok dlzka rozpocet      rok_c    dlzka_c rozpocet_c rozpocet_mil
## 1  9.3 1994   142 2.50e+07  0.3048280  0.3845886 -0.1815834       25.000
## 2  9.2 1972   175 6.00e+06 -0.5683772  1.4764179 -0.4974294        6.000
## 3  9.0 2008   152 1.85e+08  0.8605040  0.7154460  2.4781724      185.000
## 4  9.0 1974   202 1.30e+07 -0.4889949  2.3697328 -0.3810651       13.000
## 5  9.0 1957    96 3.37e+05 -1.1637444 -1.1373553 -0.5915681        0.337
## 6  9.0 1993   195 2.20e+07  0.2651369  2.1381326 -0.2314538       22.000

model_rozp_mil <- lm(imdb ~ rok + dlzka + rozpocet_mil,
                     data = udaje)
summary(model_rozp_mil)

## 
## Call:
## lm(formula = imdb ~ rok + dlzka + rozpocet_mil, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42001 -0.15545 -0.05865  0.13091  0.98723 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.6783808  1.2379574   7.010 2.28e-11 ***
## rok          -0.0003391  0.0006267  -0.541    0.589    
## dlzka         0.0021265  0.0004689   4.535 9.00e-06 ***
## rozpocet_mil  0.0003424  0.0002624   1.305    0.193    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2218 on 246 degrees of freedom
## Multiple R-squared:  0.08716,    Adjusted R-squared:  0.07603 
## F-statistic: 7.829 on 3 and 246 DF,  p-value: 5.186e-05

vif(model_rozp_mil)

##          rok        dlzka rozpocet_mil 
##     1.262283     1.017010     1.262082

Interpretácia:

koeficient pri rozpocet_mil hovorí, o koľko bodov sa zmení IMDb pri zvýšení rozpočtu o 1 milión dolárov,
VIF sú podobné ako predtým (multikolinearita medzi rokom a dĺžkou môže mierne pretrvávať),
Condition Number je nižší ako pri pôvodnom rozpočte v dolároch, ale stále vyšší ako po úplnom škálovaní všetkých premenných.

Je to kompromis: stále lepšie čísla ako pri pôvodnom rozpočte a koeficienty sú dobre interpretovateľné.

library(dplyr)
udaje <- udaje %>%
  dplyr::select(imdb, rok, dlzka, rozpocet, rozpocet_mil,
                rok_c, dlzka_c, rozpocet_c)

8.4 PCA (voliteľné)

PCA (Principal Component Analysis) vytvorí nové premenné – hlavné komponenty – ako lineárne kombinácie pôvodných premenných. Môžeme napríklad spojiť rok a rozpocet do jednej komponenty (PC1).

X_pca <- scale(udaje[, c("rok", "rozpocet")],
               center = TRUE, scale = TRUE)
pca_res <- prcomp(X_pca)

summary(pca_res)

## Importance of components:
##                           PC1    PC2
## Standard deviation     1.2049 0.7405
## Proportion of Variance 0.7258 0.2742
## Cumulative Proportion  0.7258 1.0000

udaje$PC1 <- pca_res$x[,1]

model_pca <- lm(imdb ~ dlzka + PC1, data = udaje)
summary(model_pca)

## 
## Call:
## lm(formula = imdb ~ dlzka + PC1, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42087 -0.16565 -0.05762  0.12685  0.98014 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.0172014  0.0627437 127.777  < 2e-16 ***
## dlzka        0.0021261  0.0004691   4.533 9.08e-06 ***
## PC1         -0.0085266  0.0117667  -0.725    0.469    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2218 on 247 degrees of freedom
## Multiple R-squared:  0.08278,    Adjusted R-squared:  0.07536 
## F-statistic: 11.15 on 2 and 247 DF,  p-value: 2.319e-05

vif(model_pca)

##    dlzka      PC1 
## 1.017009 1.017009

Interpretácia:

PC1 zachytáva väčšinu variability dvojice (rok, rozpocet) – je to akýsi index „modernosti a finančnej náročnosti“ filmu,
dĺžka ostáva významným vysvetľujúcim faktorom,
PC1 nemusí byť významný – kombinácia roku a rozpočtu stále nehrá veľkú rolu pri vysvetľovaní IMDb,
VIF sú nízke, pretože PC1 a dlzka sú nekorelované (ortogonálne).

Nevýhoda: PC1 je ťažšie interpretovateľná ako rok alebo rozpocet.

8.5 Hrebeňová regresia (ridge regression – voliteľné)

Hrebeňová regresia je modifikácia OLS, pri ktorej pridávame do matice X’X tzv. penalizáciu lambda * I. Odhady koeficientov sú mierne skreslené, ale stabilnejšie pri multikolinearite.

library(MASS)

X <- as.matrix(udaje[, c("rok", "dlzka", "rozpocet")])
y <- udaje$imdb

ridge_fit <- lm.ridge(y ~ X, lambda = seq(0, 10, 1))
ridge_fit

##                      Xrok      Xdlzka    Xrozpocet
##  0 8.678381 -0.0003390826 0.002126497 3.424405e-10
##  1 8.670991 -0.0003347735 0.002118001 3.407312e-10
##  2 8.663724 -0.0003305316 0.002109575 3.390442e-10
##  3 8.656578 -0.0003263555 0.002101216 3.373792e-10
##  4 8.649551 -0.0003222438 0.002092924 3.357357e-10
##  5 8.642640 -0.0003181951 0.002084698 3.341132e-10
##  6 8.635843 -0.0003142083 0.002076538 3.325112e-10
##  7 8.629156 -0.0003102819 0.002068443 3.309294e-10
##  8 8.622578 -0.0003064147 0.002060411 3.293672e-10
##  9 8.616106 -0.0003026056 0.002052443 3.278244e-10
## 10 8.609738 -0.0002988534 0.002044537 3.263005e-10

Interpretácia:

pri lambda = 0 dostaneme klasickú OLS regresiu,
ako lambda rastie, koeficienty sa zmenšujú smerom k nule a model je stabilnejší,
dôležitejšie premenné (ako dlzka) si držia nenulové koeficienty,
menej dôležité (rok, rozpocet) sa pri väčších lambda blížia k nule.

Ridge regresia je už pokročilejší nástroj, ale pekne ilustruje, ako sa dá technicky bojovať s multikolinearitou a nestabilitou odhadov.

9. Zhrnutie

Na databáze 250 najlepšie hodnotených filmov z IMDb sme skúmali vplyv roku, dĺžky a rozpočtu na hodnotenie IMDb.
Základný model ukázal, že najdôležitejšou premennou je dĺžka filmu – dlhšie filmy majú mierne vyššie hodnotenie.
Rok a rozpočet nie sú po zohľadnení dĺžky štatisticky významné.
VIF neukazuje vážnu multikolinearitu, ale Condition Number poukazuje na číselnú nestabilitu modelu kvôli rozdielnym rádom premenných.
Riešenia:
- vynechanie málo dôležitých premenných (rok, rozpocet),
- škálovanie premenných (výrazne zlepší Condition Number),
- transformácia rozpočtu na milióny dolárov,
- PCA a hrebeňová regresia (pokročilejšie metódy).

Celkovo vidíme, že multikolinearita a číselná nestabilita nie sú len teoretický problém, ale reálne ovplyvňujú interpretáciu regresného modelu aj pri praktických dátach, ako sú filmy z IMDb.

Cvičenie 10 – Multikolinearita v regresných modeloch

Diana Ďuriančiková & ChatGPT

1. Úvod

2. Dôsledky multikolinearity

3. Východiskový model a údaje

4. Odhad základného regresného modelu

5. Korelačná matica a párové vzťahy

6. VIF – Variance Inflation Factor

7. Condition Number

8. Riešenia multikolinearity a nestability

8.1 Vynechanie niektorej premennej

8.2 Škálovanie premenných (standardizácia)

8.3 Transformácia rozpočtu na milióny dolárov

8.4 PCA (voliteľné)

8.5 Hrebeňová regresia (ridge regression – voliteľné)

9. Zhrnutie