Multikolinearita v regresných modeloch
Bartáková Terézia & ChatGPT
Východiskový model a údaje
udaje <- read.csv("world_population_data.csv", dec=".", sep=",", header = TRUE, stringsAsFactors = FALSE)
udaje2023 <- udaje[, c("X2023.population", "area..km..", "density..km..", "growth.rate")]
udaje2023$growth.rate <- as.numeric(gsub("%", "", udaje2023$growth.rate))
udaje_clean <- udaje2023[complete.cases(udaje2023), ]
column_medians <- sapply(udaje2023, median, na.rm = TRUE)
udaje_imputed <- udaje2023
for (col in names(udaje2023)) {udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]}
udaje2023 <- udaje_imputed
udaje <- udaje2023Odhad základného regresného modelu
model <- lm(X2023.population ~ area..km.. + density..km.. + growth.rate, data = udaje_clean)
summary(model)##
## Call:
## lm(formula = X2023.population ~ area..km.. + density..km.. +
## growth.rate, data = udaje_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -474570585 -15188922 -13460198 -4935689 1298489341
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.519e+07 1.092e+07 1.391 0.165
## area..km.. 3.530e+01 4.595e+00 7.683 4.46e-13 ***
## density..km.. 1.388e+02 4.101e+03 0.034 0.973
## growth.rate -1.438e+06 6.559e+06 -0.219 0.827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 123300000 on 230 degrees of freedom
## Multiple R-squared: 0.2051, Adjusted R-squared: 0.1947
## F-statistic: 19.78 on 3 and 230 DF, p-value: 1.931e-11Najvýznamnejšia je premenná area..km.., kde väčšia plocha krajiny súvisí s vyššou populáciou. Nevýznamné premenné ako density..km.. a growth.rate nemajú významný prínos k predikcii.
Korelačná matica
xvars <- udaje_clean[, c("area..km..", "density..km..", "growth.rate")]
round(cor(xvars), 3)## area..km.. density..km.. growth.rate
## area..km.. 1.000 -0.065 -0.007
## density..km.. -0.065 1.000 -0.074
## growth.rate -0.007 -0.074 1.000Môžeme vidieť, že tu nie je žiadna silná korelácia medzi premennými. Môžeme povedať, že multikolinearita tu pravdepodobne nie je problém. Na presné overenie multikolinearity pri regresii použijeme v nasledujúcom kroku VIF.
pairs(udaje_clean[, c("area..km..", "density..km..", "growth.rate")],
main = "Scatterplotová matica – area, density a growth.rate")
VIF
library(car)
vif(model)## area..km.. density..km.. growth.rate
## 1.004420 1.009959 1.005720Pri posudzovaní multikolinearity sa často používa ukazovateľ VIF (Variance Inflation Factor). Ak je jeho hodnota väčšia ako 5 (prísnejšie kritérium) alebo väčšia ako 10 (miernejšie kritérium), znamená to, že medzi premennými existuje vysoká závislosť a môže to ovplyvniť spoľahlivosť odhadov regresie. V našom prípade sú hodnoty VIF pre všetky vysvetľujúce premenné nižšie než tieto hranice.
Condition Number
Používame intuitívne pravidlo. Ak Conditional number je
- < 10 → nízka multikolinearita,
- 10–30 → mierna,
- 30–100 → silná,
- 100 → veľmi vážna
X <- model.matrix(model)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)
condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number## [1] 1204635Keďže v tomto prípade indikátor významne presahuje hranicu 100, signalizuje prítomnosť závažnej multikolinearity.
Riešenia multikolinearity
Vynechanie premennej
model_noArea <- lm(X2023.population ~ density..km.. + growth.rate, data = udaje_clean)
summary(model_noArea)##
## Call:
## lm(formula = X2023.population ~ density..km.. + growth.rate,
## data = udaje_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38444696 -34523949 -27941903 -9050031 1393974457
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37250343 11779538 3.162 0.00178 **
## density..km.. -1938 4577 -0.424 0.67231
## growth.rate -2055224 7335897 -0.280 0.77961
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 137900000 on 231 degrees of freedom
## Multiple R-squared: 0.001045, Adjusted R-squared: -0.007604
## F-statistic: 0.1208 on 2 and 231 DF, p-value: 0.8863model_noDensity <- lm(X2023.population ~ area..km.. + growth.rate, data = udaje_clean)
summary(model_noDensity)##
## Call:
## lm(formula = X2023.population ~ area..km.. + growth.rate, data = udaje_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -474482010 -15192317 -13398818 -4959566 1298518486
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.527e+07 1.060e+07 1.441 0.151
## area..km.. 3.529e+01 4.575e+00 7.714 3.63e-13 ***
## growth.rate -1.455e+06 6.526e+06 -0.223 0.824
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.23e+08 on 231 degrees of freedom
## Multiple R-squared: 0.2051, Adjusted R-squared: 0.1982
## F-statistic: 29.79 on 2 and 231 DF, p-value: 3.083e-12model_noGrowth <- lm(X2023.population ~ area..km.. + density..km.., data = udaje_clean)
summary(model_noGrowth)##
## Call:
## lm(formula = X2023.population ~ area..km.. + density..km.., data = udaje_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -473070872 -13923199 -13111941 -4878250 1298689351
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.375e+07 8.708e+06 1.579 0.116
## area..km.. 3.531e+01 4.585e+00 7.702 3.91e-13 ***
## density..km.. 2.063e+02 4.081e+03 0.051 0.960
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.23e+08 on 231 degrees of freedom
## Multiple R-squared: 0.2049, Adjusted R-squared: 0.198
## F-statistic: 29.76 on 2 and 231 DF, p-value: 3.157e-12Tu vidíme, že vynechanie premennej density..km.. alebo growth.rate prakticky nemení hodnotu upraveného koeficientu determinácie (Adjusted R²), zatiaľ čo vynechanie premennej area..km.. by úplne znížilo vysvetlenú variabilitu modelu.
Preto, ak by sme chceli zjednodušiť model, uprednostnili by sme vynechanie density..km.. alebo growth.rate, pretože ich odstránenie nijako neovplyvňuje presnosť predikcie, zatiaľ čo plocha krajiny (area..km..) je kľúčová premenná, ktorú treba zachovať.
Škálovanie premenných
udaje_clean$area_c <- scale(udaje_clean$area..km.., center = TRUE, scale = TRUE)
udaje_clean$density_c <- scale(udaje_clean$density..km.., center = TRUE, scale = TRUE)
udaje_clean$growth_c <- scale(udaje_clean$growth.rate, center = TRUE, scale = TRUE)
model_centered <- lm(X2023.population ~ area_c + density_c + growth_c,
data = udaje_clean)
summary(model_centered)##
## Call:
## lm(formula = X2023.population ~ area_c + density_c + growth_c,
## data = udaje_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -474570585 -15188922 -13460198 -4935689 1298489341
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34374425 8059678 4.265 2.92e-05 ***
## area_c 62191432 8094786 7.683 4.46e-13 ***
## density_c 274746 8117076 0.034 0.973
## growth_c -1776183 8100021 -0.219 0.827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 123300000 on 230 degrees of freedom
## Multiple R-squared: 0.2051, Adjusted R-squared: 0.1947
## F-statistic: 19.78 on 3 and 230 DF, p-value: 1.931e-11vif(model_centered)## area_c density_c growth_c
## 1.004420 1.009959 1.005720Conditional Number je
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.104844Z výsledkov vidíme, že ukazovateľ Condition number sa podstatne zlepšil a nesignalizuje už žiadnu multikolinearitu. Ukazovatele VIF aj naďalej ostali veľmi nízke, čo je dobrou vlastnosťou preškálovania premenných
Iná úprava premennej, ktorá zachová interpretovateľnosť
Premenná area..km.. v tisícoch km²
udaje_clean$area1000 <- udaje_clean$area..km.. / 1000
head(udaje_clean[, c("area..km..", "area1000")])## area..km.. area1000
## 1 3287590 3287.590
## 2 9706961 9706.961
## 3 9372610 9372.610
## 4 1904569 1904.569
## 5 881912 881.912
## 6 923768 923.768Potom lineárny model dosiahne výsledky:
model_scaled <- lm(X2023.population ~ area1000 + density..km.. + growth.rate,
data = udaje_clean)
summary(model_scaled)##
## Call:
## lm(formula = X2023.population ~ area1000 + density..km.. + growth.rate,
## data = udaje_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -474570585 -15188922 -13460198 -4935689 1298489341
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15187465.6 10915412.5 1.391 0.165
## area1000 35299.1 4594.5 7.683 4.46e-13 ***
## density..km.. 138.8 4100.8 0.034 0.973
## growth.rate -1438215.8 6558773.1 -0.219 0.827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 123300000 on 230 degrees of freedom
## Multiple R-squared: 0.2051, Adjusted R-squared: 0.1947
## F-statistic: 19.78 on 3 and 230 DF, p-value: 1.931e-11vif(model_scaled)## area1000 density..km.. growth.rate
## 1.004420 1.009959 1.005720Po škálovaní premennej area..km.. na tisíce km² (area1000) a zachovaní pôvodných jednotiek pre density..km.. a growth.rate sú regresné koeficienty porovnateľné z hľadiska rádu. Hodnoty VIF sú nízke vo všetkých prípadoch (area1000 ≈ 1.004, density..km.. ≈ 1.010, growth.rate ≈ 1.006), čo naznačuje, že multikolinearita nie je problém.
X <- model.matrix(model_scaled)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)
condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number## [1] 1318.354Conditional number modelu je 1318,35, čo je relatívne vysoké, ale nie také extrémne ako v pôvodnom prípade. Znamená to, že model je mierne citlivý na numerické chyby, ale stále je vhodný na interpretáciu a odhady koeficientov sú dostatočne spoľahlivé.