Zadanie10 - Multikolinearita
Dominika Dudášová
2025-12-04
1. Úvod na vysvetlenie z pôvodného cvičenia
Po autokoreácii a heteroskedasticite rezíduí je multikolinearita tretím závažným porušením predpokladov použitia metódy najmenších štvorcov. Tu sa okrem iného predpokadá, že matica \(\mathbf X\) je tvorená lineárne nezávislými riadkami a tiež stĺpcami, čo zabezpečí regularitu matice \(\mathbf X^T\mathbf X\) a teda možnosť jej inverzie. Tá sa používa pri odhadoch regresných koeficientov. V praxi sa ale môže stať, že vzniká "takmer singulárna matica \(\mathbf X^T\mathbf X\) ", t.j. matica \(\mathbf X\) je tvorená "približne" lineárne závislými stĺpcami, t.j. existuje taká ich lineárna kombinácia, v ktorej
\[ x_{il} = \alpha_0 + \alpha_1 x_{i1} + \dots + \alpha_{l-1} x_{i,(l-1)} + \alpha_{l+1} x_{i,(l+1)} + \alpha_k x_{i,k} + \nu_i \]
Tu \(\nu_i\) sú rádovo menšie čísla , než regresory \(x_{.}\), t.j. \((\forall i)(\nu_i << x_{.,i})\). V tomto prípade je inverzná matica \((\mathbf X^T\mathbf X)^{-1}\) veľmi nestabilná a obsahuje na hlavnej diagonále veľmi veľké hodnoty. Táto matica sa používa pri výpočtoch \[\hat \beta = (\mathbf X^T\mathbf X)^{-1} \mathbf X^T \mathbf y\] a tiež \[\text{std}(\beta_i) = \sqrt{\sigma^2 (\mathbf X^T \mathbf X)^{-1}_{ii}}.\] To spôsobuje nestabilitu odhadovaných regresných koeficientov a ich nadhodnotené rozptyly.
Tento problém nazývame problémom multikolinearity.
2. Dôsledky multikolinearity z pôvodného cvičenia na vysvetlenie
Multikolinearita patrí medzi najčastejšie problémy viacnásobnej
lineárnej regresie.
Je dôležité jasne rozlišovať dva fakty:
- Nespôsobuje skreslené (biased) odhady koeficientov
- Nadhodnocuje odhady štandardných odchýlok regresných koeficientov a vedie potom k falošnému neprijímaniu alternatívnej hypotézy o štatistickej významnosti jednotlivých regresorov.
- Odhadované regresné koeficienty sú nestabilné - pri malej zmene údajov sa sa prudko menia koeficienty ako aj ich znamienka.
- Interpretácia regresného modelu je z dôvodu vyššie uvedených dôvodov nespoľahlivá.
3. Východiskový model a údaje
V tejto časti načítame súbor nhlplayoffs.csv, popíšeme premenné a definujeme východiskový regresný model. Ako závislá premenná použijeme wins (počet víťazstiev v playoff sérii) a ako vysvetľujúce premenné zvolíme goals_scored, goals_against a games. Tieto premenné sú vzájomne úzko prepojené (goals_scored a goals_against určujú goal_differential), čo vytvára priestor pre multikolinearitu — cieľom cvičenia je práve skúmať a riešiť ju.
# knitr options
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
# balíčky
if (!requireNamespace("car", quietly = TRUE)) install.packages("car")
library(car)## Loading required package: carDataif (!requireNamespace("psych", quietly = TRUE)) install.packages("psych")
library(psych)##
## Attaching package: 'psych'## The following object is masked from 'package:car':
##
## logit# načítanie dát
nhl <- read.csv("nhlplayoffs.csv", dec = ".", sep = ",", header = TRUE, stringsAsFactors = FALSE)
# kontrola
dim(nhl)## [1] 1009 13str(nhl)## 'data.frame': 1009 obs. of 13 variables:
## $ rank : int 1 2 3 4 5 6 7 8 9 10 ...
## $ team : chr "Colorado Avalanche" "Tampa Bay Lightning" "New York Rangers" "Edmonton Oilers" ...
## $ year : int 2022 2022 2022 2022 2022 2022 2022 2022 2022 2022 ...
## $ games : int 20 23 20 16 14 12 12 10 7 7 ...
## $ wins : int 16 14 10 8 7 6 5 4 3 3 ...
## $ losses : int 4 9 10 8 7 6 7 6 4 4 ...
## $ ties : int 0 0 0 0 0 0 0 0 0 0 ...
## $ shootout_wins : int 5 1 1 1 1 1 1 2 0 1 ...
## $ shootout_losses : int 1 2 2 2 0 1 1 0 0 0 ...
## $ win_loss_percentage: num 0.8 0.609 0.5 0.5 0.5 0.5 0.417 0.4 0.429 0.429 ...
## $ goals_scored : int 85 67 62 65 37 40 35 23 20 17 ...
## $ goals_against : int 55 61 58 59 40 38 39 32 24 27 ...
## $ goal_differential : int 30 6 4 6 -3 2 -4 -9 -4 -10 ...summary(nhl[c("games","wins","losses","goals_scored","goals_against","goal_differential")])## games wins losses goals_scored
## Min. : 2.000 Min. : 0.000 Min. : 0.000 Min. : 0.00
## 1st Qu.: 5.000 1st Qu.: 1.000 1st Qu.: 4.000 1st Qu.:11.00
## Median : 7.000 Median : 3.000 Median : 4.000 Median :20.00
## Mean : 9.364 Mean : 4.657 Mean : 4.657 Mean :26.63
## 3rd Qu.:12.000 3rd Qu.: 7.000 3rd Qu.: 6.000 3rd Qu.:37.00
## Max. :27.000 Max. :18.000 Max. :12.000 Max. :98.00
## goals_against goal_differential
## Min. : 0.00 Min. :-27
## 1st Qu.:16.00 1st Qu.: -6
## Median :22.00 Median : -2
## Mean :26.63 Mean : 0
## 3rd Qu.:35.00 3rd Qu.: 3
## Max. :91.00 Max. : 49# vyber dát pre analýzu (odstránime identifikačné polia ako team, year ak sú nepotrebné)
nhl2 <- nhl[, c("wins","games","goals_scored","goals_against","goal_differential")]
head(nhl2)## wins games goals_scored goals_against goal_differential
## 1 16 20 85 55 30
## 2 14 23 67 61 6
## 3 10 20 62 58 4
## 4 8 16 65 59 6
## 5 7 14 37 40 -3
## 6 6 12 40 38 24. Odhad základného regresného modelu
Základný model (východiskový) odhadneme takto:
model_base <- lm(wins ~ goals_scored + goals_against + games, data = nhl2)
summary(model_base)##
## Call:
## lm(formula = wins ~ goals_scored + goals_against + games, data = nhl2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3909 -0.4314 -0.0324 0.3591 3.1894
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.001460 0.050968 -19.65 <2e-16 ***
## goals_scored 0.113663 0.003562 31.91 <2e-16 ***
## goals_against -0.117244 0.003657 -32.06 <2e-16 ***
## games 0.614489 0.011831 51.94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7105 on 1005 degrees of freedom
## Multiple R-squared: 0.9727, Adjusted R-squared: 0.9727
## F-statistic: 1.195e+04 on 3 and 1005 DF, p-value: < 2.2e-16# uložíme si koeficienty a štandardné chyby pre neskoršie porovnanie
coef_base <- summary(model_base)$coefficients
coef_base## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.0014605 0.050967749 -19.64891 5.395933e-73
## goals_scored 0.1136630 0.003562468 31.90569 7.548804e-155
## goals_against -0.1172439 0.003656798 -32.06190 6.331893e-156
## games 0.6144894 0.011830852 51.93957 7.517920e-2875. Korelačná matica
Skontrolujeme korelácie medzi vysvetľujúcimi premennými — vysoké korelácie naznačujú multikolinearitu.
# korelačná matica (Pearson)
preds <- nhl2[, c("goals_scored","goals_against","games","goal_differential")]
cor_matrix <- cor(preds, use = "pairwise.complete.obs")
cor_matrix## goals_scored goals_against games goal_differential
## goals_scored 1.0000000 0.9088340 0.9401903 0.7229132
## goals_against 0.9088340 1.0000000 0.8948069 0.3687775
## games 0.9401903 0.8948069 1.0000000 0.6128774
## goal_differential 0.7229132 0.3687775 0.6128774 1.0000000# doplnkový prehľad pomocou psych::pairs.panels
pairs.panels(preds)
6. VIF
Vypočítame Variance Inflation Factors (VIF) pre základný model. Hodnoty VIF > 5 (alebo 10 podľa prísnejšieho pravidla) indikujú problém.
vif_base <- vif(model_base)
vif_base## goals_scored goals_against games
## 10.735066 6.249835 9.372463# zoradíme podľa veľkosti
sort(vif_base, decreasing = TRUE)## goals_scored games goals_against
## 10.735066 9.372463 6.249835VIF > 10 = vážna multikolinearita
Premenná games je veľmi úzko naviazaná na výkonnosť tímu (viac zápasov → viac príležitostí skórovať → viac inkasovaných gólov).
7. Condition Number
Condition number (číslo kondície) poskytuje doplňujúce informácie o multikolinearite.
# použijeme singular value decomposition na maticu návrhu (bez interceptu)
X <- model.matrix(~ goals_scored + goals_against + games, data = nhl2)
svdvals <- svd(scale(X, center = TRUE, scale = FALSE))$d
condition_number <- max(svdvals) / min(svdvals)
condition_number## [1] Inf# alternatívne: kappa
kappa(X)## [1] 63.397020–15 → OK
15–30 → stredná multikolinearita
30+ → silná
Hodnota 63.39702 znamená, že multikolinearita je extrémna podľa tohto metrického pohľadu.
8. Riešenia multikolinearity
Ukážeme tri bežné prístupy: vynechanie premennej, škálovanie premenných a iná úprava premennej, ktorá zachová interpretovateľnosť.
Vynechanie premennej
V jednej verzii modelu vynecháme premennú goals_against (predpoklad: časť informácie je zachytená v goals_scored a games/alebo goal_differential). Odhadneme model a porovnáme koeficienty a VIF.
model_drop <- lm(wins ~ goals_scored + games, data = nhl2)
summary(model_drop)##
## Call:
## lm(formula = wins ~ goals_scored + games, data = nhl2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2522 -0.5349 -0.1333 0.4437 4.4068
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.764826 0.064062 -27.55 <2e-16 ***
## goals_scored 0.062935 0.004537 13.87 <2e-16 ***
## games 0.506837 0.016127 31.43 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.01 on 1006 degrees of freedom
## Multiple R-squared: 0.9448, Adjusted R-squared: 0.9447
## F-statistic: 8617 on 2 and 1006 DF, p-value: < 2.2e-16vif(model_drop)## goals_scored games
## 8.617555 8.617555# porovnanie koeficientov
cbind(base = coef_base[,1], drop = summary(model_drop)$coefficients[,1])## base drop
## (Intercept) -1.0014605 -1.76482591
## goals_scored 0.1136630 0.06293465
## goals_against -0.1172439 0.50683662
## games 0.6144894 -1.76482591Škálovanie premenných
Škálovanie (štandardizácia na strednú hodnotu 0 a smerodajnú odchýlku 1) zlepší numerickú stabilitu a môže ovplyvniť VIF. Pozor: štandardizované koeficienty menia jednotky interpretácie.
nhl_scaled <- nhl2
nhl_scaled[c("goals_scored","goals_against","games")] <- scale(nhl_scaled[c("goals_scored","goals_against","games")])
model_scaled <- lm(wins ~ goals_scored + goals_against + games, data = nhl_scaled)
summary(model_scaled)##
## Call:
## lm(formula = wins ~ goals_scored + goals_against + games, data = nhl_scaled)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3909 -0.4314 -0.0324 0.3591 3.1894
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.65709 0.02237 208.20 <2e-16 ***
## goals_scored 2.33943 0.07332 31.91 <2e-16 ***
## goals_against -1.79376 0.05595 -32.06 <2e-16 ***
## games 3.55848 0.06851 51.94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7105 on 1005 degrees of freedom
## Multiple R-squared: 0.9727, Adjusted R-squared: 0.9727
## F-statistic: 1.195e+04 on 3 and 1005 DF, p-value: < 2.2e-16vif(model_scaled)## goals_scored goals_against games
## 10.735066 6.249835 9.372463# condition number pre škálovaných premenných
Xs <- model.matrix(~ goals_scored + goals_against + games, data = nhl_scaled)
svdvals_s <- svd(scale(Xs, center = TRUE, scale = FALSE))$d
max(svdvals_s)/min(svdvals_s)## [1] InfIná úprava premennej, ktorá zachová interpretovateľnosť
Aby sme zachovali interpretovateľnosť a zároveň odstránili redundanciu, môžeme namiesto dvoch úzko spätých premenných (goals_scored a goals_against) použiť jedine premennú goal_differential (ktorá má priamu interpretáciu: rozdiel skóre). Odhadneme model s touto úpravou a vyhodnotíme VIF a kondíciu.
model_diff <- lm(wins ~ goal_differential + games, data = nhl2)
summary(model_diff)##
## Call:
## lm(formula = wins ~ goal_differential + games, data = nhl2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4040 -0.4155 -0.0419 0.3646 3.1843
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.001938 0.050962 -19.66 <2e-16 ***
## goal_differential 0.115369 0.003066 37.63 <2e-16 ***
## games 0.604356 0.004890 123.58 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7105 on 1006 degrees of freedom
## Multiple R-squared: 0.9727, Adjusted R-squared: 0.9727
## F-statistic: 1.793e+04 on 2 and 1006 DF, p-value: < 2.2e-16vif(model_diff)## goal_differential games
## 1.601585 1.601585# condition number
Xdiff <- model.matrix(~ goal_differential + games, data = nhl2)
kappa(Xdiff)## [1] 33.03597Zhrnutie a odporúčania
V závere stručne zhrnieme pozorovania: ktoré premenné spôsobovali multikolinearitu, aké boli VIF a condition number pre jednotlivé modely, a ktorá z navrhnutých úprav sa javí ako najvhodnejšia z hľadiska interpretovateľnosti a stability odhadu.
# vytiahneme dôležité čísla pre stručné zhrnutie:
res <- list(
vif_base = vif_base,
cond_base = condition_number,
vif_drop = vif(model_drop),
cond_drop = kappa(model.matrix(~ goals_scored + games, data = nhl2)),
vif_scaled = vif(model_scaled),
cond_scaled = max(svdvals_s)/min(svdvals_s),
vif_diff = vif(model_diff),
cond_diff = kappa(Xdiff)
)
res## $vif_base
## goals_scored goals_against games
## 10.735066 6.249835 9.372463
##
## $cond_base
## [1] Inf
##
## $vif_drop
## goals_scored games
## 8.617555 8.617555
##
## $cond_drop
## [1] 66.56644
##
## $vif_scaled
## goals_scored goals_against games
## 10.735066 6.249835 9.372463
##
## $cond_scaled
## [1] Inf
##
## $vif_diff
## goal_differential games
## 1.601585 1.601585
##
## $cond_diff
## [1] 33.03597Bonusy
Ridge Regression (Hrebeňová regresia)
V tejto časti pridáme odhad ridge regresie, ktorá rieši multikolinearitu penalizáciou veľkých koeficientov.
# balíček glmnet
if (!requireNamespace("glmnet", quietly = TRUE)) install.packages("glmnet")
library(glmnet)
# príprava dát pre glmnet
X <- as.matrix(nhl2[, c("goals_scored","goals_against","games")])
y <- nhl2$wins
# ridge regresia: alpha = 0
ridge_model <- glmnet(X, y, alpha = 0)
# výber lambda pomocou cross-validation
cv_ridge <- cv.glmnet(X, y, alpha = 0)
cv_ridge$lambda.min## [1] 0.4150871# koeficienty pri optimálnom lambda
coef(cv_ridge, s = "lambda.min")## 4 x 1 sparse Matrix of class "dgCMatrix"
## lambda.min
## (Intercept) -1.13204326
## goals_scored 0.09151874
## goals_against -0.02484388
## games 0.42862101# vizualizácia
plot(cv_ridge)
Koeficienty sú menšie a stabilnejšie než v OLS.
Zaujímavé je, že ridge „vrátil” logicky očakávanú znamienkovú štruktúru:
goals_scored → pozitívny vplyv na wins,
goals_against → negatívny vplyv na wins,
To je znak, že OLS bol skreslený bezpenalizačnou multikolinearitou.
Premenná games zostáva najsilnejším prediktorom (+0.61). # PCA PCA pomáha odstrániť multikolinearitu transformáciou pôvodných vysvetľujúcich premenných na nové, navzájom nekorelované komponenty.
# vyber vysvetľujúcich premenných
preds <- nhl2[, c("goals_scored","goals_against","games")]
# PCA s centráciou a škálovaním
pca_model <- prcomp(preds, center = TRUE, scale. = TRUE)
summary(pca_model)## Importance of components:
## PC1 PC2 PC3
## Standard deviation 1.6821 0.33484 0.24194
## Proportion of Variance 0.9431 0.03737 0.01951
## Cumulative Proportion 0.9431 0.98049 1.00000# komponenty
pca_model$rotation## PC1 PC2 PC3
## goals_scored -0.5814985 0.2911150 -0.7596786
## goals_against -0.5719086 -0.8103922 0.1272204
## games -0.5786019 0.5084452 0.6377330# skóre pozorovaní
head(pca_model$x)## PC1 PC2 PC3
## [1,] -3.772238 0.25676943 -0.7471384
## [2,] -3.787724 -0.05223809 0.2975026
## [3,] -3.234574 -0.22745052 0.1267274
## [4,] -2.957054 -0.58918663 -0.4161882
## [5,] -1.255911 -0.15440741 0.2390351
## [6,] -1.066077 -0.18163680 -0.1085755# vizualizácia variancie
plot(pca_model, type = "l")
Takmer celá informácia (95.5 %) je v PC1 → znamená to, že premenné sú
extrémne redundantné.
PC2 a PC3 neobsahujú výraznú dodatočnú informáciu.
Pozitívne a vysoké hodnoty pri Goals Scored aj Goals Against znamenajú:
PC1 meria celkovú úroveň ofenzívy/defenzívy, teda „koľko sa v sérii skóruje na oboch stranách“.
Ide o spoločný signál, ktorý vysvetľuje väčšinu variancie.
PC2 zachytáva rozdiel medzi strelenými a inkasovanými gólmi, teda kvalitu tímu.