Cvičenie 10 – Multikolinearita v regresných modeloch

knitr::opts_chunk$set(
echo = TRUE,
message = FALSE,
warning = FALSE
)

rm(list = ls())

firmy <- read.csv("dataEKONOMETRIA.csv", stringsAsFactors = FALSE)

# Podiel žien vo vedení

firmy$podiel_zien <- firmy$Z / (firmy$M + firmy$Z)

# Vytvoríme pracovný dataset pre multikolinearitu

firmy_mc <- firmy[, c("ROE", "ROA", "EBITDAmarza", "M", "Z", "podiel_zien")]

str(firmy_mc)

'data.frame':   53 obs. of  6 variables:
 $ ROE        : num  0.2702 0.2938 0.6643 0.0956 0.1223 ...
 $ ROA        : num  0.0788 0.0801 0.0553 0.0363 0.0216 ...
 $ EBITDAmarza: num  0.139 0.0339 0.0567 0.0803 0 ...
 $ M          : int  1 2 6 5 6 4 7 0 10 1 ...
 $ Z          : int  0 1 0 3 2 0 4 2 7 0 ...
 $ podiel_zien: num  0 0.333 0 0.375 0.25 ...

summary(firmy_mc)

      ROE                ROA            EBITDAmarza             M     
 Min.   :-1.84304   Min.   :-0.16527   Min.   :-0.04082   Min.   : 0  
 1st Qu.: 0.08053   1st Qu.: 0.02724   1st Qu.: 0.03528   1st Qu.: 2  
 Median : 0.12639   Median : 0.06000   Median : 0.05559   Median : 4  
 Mean   : 0.20327   Mean   : 0.09227   Mean   : 0.11727   Mean   : 6  
 3rd Qu.: 0.31271   3rd Qu.: 0.11555   3rd Qu.: 0.13260   3rd Qu.: 9  
 Max.   : 0.99781   Max.   : 0.70448   Max.   : 0.83243   Max.   :23  
       Z          podiel_zien    
 Min.   :0.000   Min.   :0.0000  
 1st Qu.:0.000   1st Qu.:0.0000  
 Median :1.000   Median :0.1667  
 Mean   :1.717   Mean   :0.2020  
 3rd Qu.:3.000   3rd Qu.:0.3333  
 Max.   :7.000   Max.   :1.0000  

Používam ROE ako závislú premennú. Vysvetľujúce premenné sú ROA, EBITDA marža a dva počty (muži a ženy vo vedení) plus odvodený ukazovateľ podielu žien vo vedení.

model10_1 <- lm(ROE ~ ROA + EBITDAmarza + M + Z, data = firmy_mc)
summary(model10_1)

Call:
lm(formula = ROE ~ ROA + EBITDAmarza + M + Z, data = firmy_mc)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.26230 -0.07712 -0.01667  0.13934  0.79777 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.001288   0.104697  -0.012 0.990238    
ROA          1.820819   0.486836   3.740 0.000491 ***
EBITDAmarza -0.485064   0.416720  -1.164 0.250176    
M            0.020364   0.014911   1.366 0.178417    
Z           -0.016741   0.039112  -0.428 0.670545    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4191 on 48 degrees of freedom
Multiple R-squared:  0.2362,    Adjusted R-squared:  0.1725 
F-statistic:  3.71 on 4 and 48 DF,  p-value: 0.01037

xvars <- firmy_mc[, c("ROA", "EBITDAmarza", "M", "Z")]
round(cor(xvars), 3)

               ROA EBITDAmarza      M     Z
ROA          1.000       0.366 -0.100 0.079
EBITDAmarza  0.366       1.000  0.141 0.091
M           -0.100       0.141  1.000 0.580
Z            0.079       0.091  0.580 1.000

vif_manual <- function(data, vars) {
res <- numeric(length(vars))
names(res) <- vars

for (v in vars) {
others <- vars[vars != v]
form <- as.formula(paste(v, "~", paste(others, collapse = " + ")))
m <- lm(form, data = data)
R2 <- summary(m)$r.squared
res[v] <- 1 / (1 - R2)
}
return(res)
}

vif_vals1 <- vif_manual(firmy_mc, c("ROA", "EBITDAmarza", "M", "Z"))
vif_table1 <- data.frame(Premenna = names(vif_vals1), VIF = round(vif_vals1, 3))
vif_table1

NA

Condition number (číslo podmienky)

X1 <- model.matrix(model10_1)[, -1]   # bez interceptu
e1 <- eigen(t(X1) %*% X1)
cond_num1 <- sqrt(max(e1$values) / min(e1$values))
cond_num1

[1] 74.05863

Umelé vytvorenie silnej multikolinearity

firmy_mc$ROA100 <- 100 * firmy_mc$ROA

model10_2 <- lm(ROE ~ ROA + ROA100 + EBITDAmarza + M + Z, data = firmy_mc)
summary(model10_2)

Call:
lm(formula = ROE ~ ROA + ROA100 + EBITDAmarza + M + Z, data = firmy_mc)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.26230 -0.07712 -0.01667  0.13934  0.79777 

Coefficients: (1 not defined because of singularities)
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.001288   0.104697  -0.012 0.990238    
ROA          1.820819   0.486836   3.740 0.000491 ***
ROA100             NA         NA      NA       NA    
EBITDAmarza -0.485064   0.416720  -1.164 0.250176    
M            0.020364   0.014911   1.366 0.178417    
Z           -0.016741   0.039112  -0.428 0.670545    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4191 on 48 degrees of freedom
Multiple R-squared:  0.2362,    Adjusted R-squared:  0.1725 
F-statistic:  3.71 on 4 and 48 DF,  p-value: 0.01037

VIF a condition number pri silnej multikolinearite

vif_vals2 <- vif_manual(firmy_mc, c("ROA", "ROA100", "EBITDAmarza", "M", "Z"))
vif_table2 <- data.frame(Premenna = names(vif_vals2), VIF = round(vif_vals2, 3))
vif_table2

X2 <- model.matrix(model10_2)[, -1]
e2 <- eigen(t(X2) %*% X2)
cond_num2 <- sqrt(max(e2$values) / min(e2$values))
cond_num2

[1] NaN

cond_num2

[1] NaN

# vytvoríme X a y
y <- firmy_mc$ROE10
X <- model.matrix(model10_2)[, -1]   # bez intereceptu

X <- model.matrix(model10_2)[, -1]  # bez interceptu
y <- firmy_mc$ROE10

ridge_manual <- function(lambda) {
  I <- diag(ncol(X))
  beta_ridge <- solve(t(X) %*% X + lambda * I) %*% t(X) %*% y
  return(beta_ridge)
}

lambda_grid <- c(0.1, 1, 10, 100, 1000)
ridge_results <- sapply(lambda_grid, ridge_manual)

Error in solve(t(X) %*% X + lambda * I) %*% t(X) %*% y : 
  requires numeric/complex matrix/vector arguments

firmy_mc$EBIT <- as.numeric(gsub(",", "", firmy_mc$EBIT))

str(firmy_mc)

'data.frame':   53 obs. of  8 variables:
 $ ROE        : num  0.2702 0.2938 0.6643 0.0956 0.1223 ...
 $ ROA        : num  0.0788 0.0801 0.0553 0.0363 0.0216 ...
 $ EBITDAmarza: num  0.139 0.0339 0.0567 0.0803 0 ...
 $ M          : int  1 2 6 5 6 4 7 0 10 1 ...
 $ Z          : int  0 1 0 3 2 0 4 2 7 0 ...
 $ podiel_zien: num  0 0.333 0 0.375 0.25 ...
 $ ROA100     : num  7.88 8.01 5.53 3.63 2.16 ...
 $ EBIT       : num  0.139 0.0339 0.0567 0.0803 0 ...

y <- firmy_mc$ROE

X <- model.matrix(~ ROA + ROA100 + EBITDAmarza + M + Z, data = firmy_mc)[, -1]

class(X)  # musí byť 'matrix'

[1] "matrix" "array" 

class(y)  # musí be 'numeric'

[1] "numeric"

ridge_manual <- function(lambda){
  I <- diag(ncol(X))
  solve(t(X) %*% X + lambda * I) %*% t(X) %*% y
}

lambda_grid <- c(0.1, 1, 10, 100, 1000)

ridge_results <- sapply(lambda_grid, ridge_manual)
colnames(ridge_results) <- paste0("lambda_", lambda_grid)

ridge_results

        lambda_0.1      lambda_1     lambda_10    lambda_100   lambda_1000
[1,]  0.0001795589  0.0001689505  0.0001576678  0.0001534422  0.0001458169
[2,]  0.0179558870  0.0168950464  0.0157667756  0.0153442169  0.0145816864
[3,] -0.4434369791 -0.2480002633 -0.0456523703 -0.0047799105 -0.0003373724
[4,]  0.0198373788  0.0178111942  0.0153976660  0.0129965687  0.0092385149
[5,] -0.0165828711 -0.0155961352 -0.0134044563 -0.0063039388  0.0007598433

Ridge regresia stabilizuje koeficienty

Pri malom λ = 0.1 sú koeficienty veľmi podobné OLS (bežnej regresii).

Keď λ rastie → 10, 100, 1000:

Absolútne hodnoty koeficientov sa zmenšujú

-tým zabraňujú príliš veľkým a nestabilným hodnotám spôsobeným multikolinearitou.

Napríklad koeficient EBITDAmarza:

pri λ = 0.1 → −0.443

pri λ = 1 → −0.248

pri λ = 10 → −0.045

pri λ = 100 → −0.0048

pri λ = 1000 → ≈ 0

Toto je typické správanie Ridge regresie → multikolineárny koeficient sa “stiahne” smerom k nule.

2. Najviac penalizovaný koeficient = EBITDA marža

Je to presne ten, ktorý mal:

najväčšiu koreláciu s ROA,

najvyšší VIF v tvojej multikolinearite,

spôsoboval nestabilnosť OLS odhadu.

Ridge regresia ho postupne zmenšuje až úplne takmer na nulu pri λ=1000.

Vysoká multikolinearita medzi ROA a EBITDA maržou spôsobila, že Ridge regresia považuje EBITDA maržu za najmenej spoľahlivý koeficient – preto ho tlmí.

3. Ostatné koeficienty ostávajú stabilné

ROA, M a Z sa menia len mierne:

zachovávajú si znamienko,

nemenia sa dramaticky,

Ridge ich upravuje iba jemne.

To znamená, že tieto premenné:

-sú informatívne -nespôsobujú multikolinearitu -ostávajú robustné

4. Pri vysokej penalizácii (λ=1000) sa model správa konzervatívne

Vysoké λ tlačí všetky koeficienty k nule.

To je normálne – Ridge regresia:

nezahadzuje premenné (to robí až Lasso),

ale znižuje ich vplyv, ak sú nestabilné.

5. Záver Ridge regresie

Ridge regresia potvrdila prítomnosť multikolinearity v pôvodnom modeli.

Najmä premenná EBITDA marža bola vysoko korelovaná s ROA, čo spôsobovalo nestabilitu OLS odhadov. Ridge regresia postupným zvyšovaním penalizačného parametra λ „stiahla“ jej koeficient až takmer k nule, čo je typický prejav korekcie multikolinearity.

Ostatné koeficienty zostali stabilné a výrazne sa nemenili, čo znamená, že tieto premenné sú informatívne a robustné.

Ridge regresia teda zlepšuje stabilitu modelu a odstraňuje problém multikolinearity bez toho, aby vyhodila niektorého z regresorov.

# vytvoríme dataset len z numerických premenných použitých v modeli
pca_data <- firmy_mc[, c("ROA", "ROA100", "EBITDAmarza", "M", "Z")]

pca_scaled <- scale(pca_data)

pca_model <- prcomp(pca_scaled, center = TRUE, scale. = TRUE)
summary(pca_model)

Importance of components:
                          PC1    PC2    PC3     PC4       PC5
Standard deviation     1.4952 1.2682 0.8809 0.61651 4.019e-16
Proportion of Variance 0.4471 0.3216 0.1552 0.07602 0.000e+00
Cumulative Proportion  0.4471 0.7688 0.9240 1.00000 1.000e+00

plot(pca_model, type = "l", main = "Scree plot – PCA")

pca_model$rotation

                    PC1         PC2         PC3        PC4           PC5
ROA          0.64664748  0.09534693 -0.24024573  0.1226295  7.071068e-01
ROA100       0.64664748  0.09534693 -0.24024573  0.1226295 -7.071068e-01
EBITDAmarza  0.39016367 -0.15079248  0.88419912 -0.2079083 -1.110223e-16
M           -0.01046233 -0.70981111  0.04877754  0.7026233  9.714451e-17
Z            0.10656893 -0.67472050 -0.31682275 -0.6580415  0.000000e+00

pca_scores <- as.data.frame(pca_model$x)

head(pca_scores)

NA

model_pca <- lm(ROE ~ pca_scores$PC1 + pca_scores$PC2, data = firmy_mc)
summary(model_pca)

Call:
lm(formula = ROE ~ pca_scores$PC1 + pca_scores$PC2, data = firmy_mc)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.27396 -0.07282 -0.01729  0.06988  0.80711 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)     0.20327    0.05916   3.436  0.00120 **
pca_scores$PC1  0.12235    0.03994   3.063  0.00352 **
pca_scores$PC2 -0.01651    0.04710  -0.351  0.72733   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4307 on 50 degrees of freedom
Multiple R-squared:  0.1597,    Adjusted R-squared:  0.1261 
F-statistic: 4.753 on 2 and 50 DF,  p-value: 0.01289

PCA analýza identifikovala dva hlavné zdroje variability v dátach. Prvý komponent (PC1) predstavuje finančnú výkonnosť firmy, pretože silne súvisí s ukazovateľmi ROA, ROA100 a EBITDA maržou. Druhý komponent (PC2) zachytáva štruktúru pracovnej sily a súvisí najmä s počtom mužských a ženských zamestnancov.

Regressný model s PC1 a PC2 ukázal, že finančná výkonnosť (PC1) má štatisticky významný a pozitívny vplyv na ROE, zatiaľ čo demografická štruktúra pracovnej sily (PC2) významný vplyv nemá. Model je ako celok významný a PCA úplne odstránilo multikolinearitu, čím sa zvýšila stabilita odhadov.

PCA sa tak ukazuje ako vhodný spôsob, ako nahradiť pôvodné korelované regresory dvoma nezávislými komponentmi, bez straty vysvetľovacej schopnosti modelu.

Zhrnutie celého cvičenia

Multikolinearita bola jednoznačne prítomná v pôvodnom modeli — najmä medzi ROA a EBITDA maržou.

Umelé zavedenie ROA100 extrémne zvýšilo multikolinearitu, čo spôsobilo numerickú nestabilitu OLS modelu.

Ridge regresia potvrdila problém a efektívne stabilizovala koeficienty prostredníctvom penalizácie.

PCA ponúkla ešte lepšie riešenie, pretože úplne odstránila korelácie medzi vysvetľujúcimi premennými.

Finančná výkonnosť (PC1) bola jasne najdôležitejším faktorom, ktorý ovplyvňuje ROE.

Demografia vedenia (PC2) nemala štatisticky významný vplyv na ROE.

PCA aj Ridge potvrdili, že problém multikolinearity má reálny dopad na stabilitu, signifikanciu aj interpretáciu koeficientov.