R Notebook

# 0) Balíčky (len ak chýbajú, doinštalujú sa)
needs <- c("dplyr","ggplot2","tseries","car")
to_install <- needs[!needs %in% rownames(installed.packages())]
if (length(to_install)) install.packages(to_install, quiet = TRUE)
lapply(needs, library, character.only = TRUE)

[[1]]
[1] "dplyr"     "stats"     "graphics"  "grDevices" "utils"     "datasets" 
[7] "methods"   "base"     

[[2]]
[1] "ggplot2"   "dplyr"     "stats"     "graphics"  "grDevices" "utils"    
[7] "datasets"  "methods"   "base"     

[[3]]
 [1] "tseries"   "ggplot2"   "dplyr"     "stats"     "graphics"  "grDevices"
 [7] "utils"     "datasets"  "methods"   "base"     

[[4]]
 [1] "car"       "carData"   "tseries"   "ggplot2"   "dplyr"     "stats"    
 [7] "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"     

# 1) Robustné načítanie dát
safe_read <- function(path){
  df <- tryCatch(read.csv(path, header=TRUE, sep=";", dec=",",
                          check.names=FALSE, stringsAsFactors=FALSE), error=function(e) NULL)
  if (is.null(df) || ncol(df) == 1)
    df <- read.csv(path, header=TRUE, sep=",", dec=".",
                   check.names=FALSE, stringsAsFactors=FALSE)
  df
}
udaje <- safe_read("dataEKONOMETRIA.csv")
if (is.null(udaje)) stop("Súbor 'dataEKONOMETRIA.csv' sa nepodarilo načítať. Cesta: ", getwd())

cat("Nájdené stĺpce:\n"); print(colnames(udaje)); cat("\n")

Nájdené stĺpce:
[1] "Nazov"       "Kategoria"   "Forma"       "ROE"         "ROA"        
[6] "EBIT"        "EBITDAmarza" "M"           "Z"          

head(data)

                                                                            
1 function (..., list = character(), package = NULL, lib.loc = NULL,        
2     verbose = getOption("verbose"), envir = .GlobalEnv, overwrite = TRUE) 
3 {                                                                         
4     fileExt <- function(x) {                                              
5         db <- grepl("\\\\.[^.]+\\\\.(gz|bz2|xz)$", x)                     
6         ans <- sub(".*\\\\.", "", x)                                      

install.packages(c("readr", "dplyr", "ggplot2", "lmtest", "sandwich", "car", "tseries"))

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/readr_2.1.5.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/dplyr_1.1.4.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/ggplot2_4.0.0.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/lmtest_0.9-40.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/sandwich_3.1-1.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/car_3.1-3.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/tseries_0.10-58.tar.gz'

The downloaded source packages are in
    ‘/tmp/RtmpUc0CYs/downloaded_packages’

# Načítanie potrebných balíčkov

library(readr)
library(dplyr)
library(ggplot2)
library(lmtest)
library(sandwich)
library(car)
library(tseries)

# Načítanie dát a základný opis

firmy <- read_csv("dataEKONOMETRIA.csv")

glimpse(firmy)

Rows: 53
Columns: 9
$ Nazov       <chr> "Accenture, s.r.o.", "Beiersdorf Slovakia, s.r.o.", "Be…
$ Kategoria   <dbl> 3, 1, 4, 5, 3, 4, 4, 3, 4, 4, 2, 4, 3, 4, 4, 3, 3, 3, 3…
$ Forma       <chr> "sro", "sro", "sro", "sro", "as", "sro", "as", "sro", "…
$ ROE         <dbl> 0.270192860, 0.293773272, 0.664287656, 0.095621306, 0.1…
$ ROA         <dbl> 0.078790333, 0.080062193, 0.055303738, 0.036313948, 0.0…
$ EBIT        <dbl> 1250080, 1072061, 12586000, 19677000, 8126000, 9582736,…
$ EBITDAmarza <dbl> 0.13899269, 0.03386505, 0.05668215, 0.08028712, 0.00000…
$ M           <dbl> 1, 2, 6, 5, 6, 4, 7, 0, 10, 1, 6, 4, 8, 5, 3, 11, 2, 1,…
$ Z           <dbl> 0, 1, 0, 3, 2, 0, 4, 2, 7, 0, 0, 0, 1, 2, 0, 2, 0, 1, 1…

summary(firmy)

    Nazov             Kategoria        Forma                ROE          
 Length:53          Min.   :1.000   Length:53          Min.   :-1.84304  
 Class :character   1st Qu.:2.000   Class :character   1st Qu.: 0.08053  
 Mode  :character   Median :3.000   Mode  :character   Median : 0.12639  
                    Mean   :3.057                      Mean   : 0.20327  
                    3rd Qu.:4.000                      3rd Qu.: 0.31271  
                    Max.   :5.000                      Max.   : 0.99781  
      ROA                EBIT            EBITDAmarza             M     
 Min.   :-0.16527   Min.   :  -790635   Min.   :-0.04082   Min.   : 0  
 1st Qu.: 0.02724   1st Qu.:  2008948   1st Qu.: 0.03528   1st Qu.: 2  
 Median : 0.06000   Median :  6643237   Median : 0.05559   Median : 4  
 Mean   : 0.09227   Mean   : 77168704   Mean   : 0.11727   Mean   : 6  
 3rd Qu.: 0.11555   3rd Qu.: 30844065   3rd Qu.: 0.13260   3rd Qu.: 9  
 Max.   : 0.70448   Max.   :932605000   Max.   : 0.83243   Max.   :23  
       Z        
 Min.   :0.000  
 1st Qu.:0.000  
 Median :1.000  
 Mean   :1.717  
 3rd Qu.:3.000  
 Max.   :7.000  

# Úprava premenných a vytvorenie nových ukazovateľov

firmy <- firmy %>%
mutate(
# Prevod EBIT na číselnú premennú odstránením oddeľovača tisícok
EBIT_num = as.numeric(gsub(",", "", EBIT)),
# Celkový počet členov vedenia
pocet_clenov = M + Z,
# Podiel žien vo vedení
podiel_zien = Z / pocet_clenov,
# Podiel mužov vo vedení
podiel_muzov = M / pocet_clenov,
# Kategória a forma ako faktorové premenné
Kategoria = factor(Kategoria),
Forma = factor(Forma)
)

# Kontrola upravených dát

firmy %>%
select(Nazov, Kategoria, Forma, ROE, ROA, EBIT_num, EBITDAmarza, M, Z, podiel_zien, podiel_muzov) %>%
head()

NA

# Základné vizualizácie vybraných premenných

ggplot(firmy, aes(x = ROE)) +
geom_histogram(bins = 15) +
labs(title = "Rozdelenie ukazovateľa ROE")

ggplot(firmy, aes(x = podiel_zien)) +
geom_histogram(bins = 10) +
labs(title = "Rozdelenie podielu žien vo vedení")

ggplot(firmy, aes(x = podiel_zien, y = ROE)) +
geom_point() +
labs(title = "Vzťah medzi ROE a podielom žien vo vedení")

# Odhad lineárneho regresného modelu s ROE ako závislou premennou

model1 <- lm(ROE ~ podiel_zien + ROA + EBITDAmarza, data = firmy)

# Výstup modelu obsahuje odhady parametrov, ich smerodajné odchýlky a štatistickú významnosť

summary(model1)

Call:
lm(formula = ROE ~ podiel_zien + ROA + EBITDAmarza, data = firmy)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.25674 -0.07073 -0.00477  0.05133  0.82942 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.05891    0.08947   0.658   0.5133   
podiel_zien  0.21318    0.28974   0.736   0.4654   
ROA          1.55811    0.50071   3.112   0.0031 **
EBITDAmarza -0.36215    0.41051  -0.882   0.3820   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4214 on 49 degrees of freedom
Multiple R-squared:  0.2117,    Adjusted R-squared:  0.1634 
F-statistic: 4.385 on 3 and 49 DF,  p-value: 0.008219

# Rozšírený model s kategóriou a formou ako dodatočnými vysvetľujúcimi premennými

model2 <- lm(ROE ~ podiel_zien + ROA + EBITDAmarza + Kategoria + Forma, data = firmy)

summary(model2)

Call:
lm(formula = ROE ~ podiel_zien + ROA + EBITDAmarza + Kategoria + 
    Forma, data = firmy)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.15800 -0.10117 -0.03845  0.10002  0.89147 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.233991   0.275110  -0.851  0.39963   
podiel_zien  0.255261   0.301150   0.848  0.40124   
ROA          1.812357   0.536228   3.380  0.00153 **
EBITDAmarza -0.550542   0.445212  -1.237  0.22280   
Kategoria2   0.158131   0.249659   0.633  0.52976   
Kategoria3   0.299437   0.245102   1.222  0.22833   
Kategoria4   0.397898   0.254715   1.562  0.12542   
Kategoria5   0.226616   0.380624   0.595  0.55464   
Formasro     0.008296   0.142221   0.058  0.95375   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4272 on 44 degrees of freedom
Multiple R-squared:  0.2726,    Adjusted R-squared:  0.1403 
F-statistic: 2.061 on 8 and 44 DF,  p-value: 0.06081

# Diagnostické grafy pre model1

par(mfrow = c(2, 2))
plot(model1)
par(mfrow = c(1, 1))

# Test normality rezíduí (Shapiro-Wilk)

shapiro.test(residuals(model1))

    Shapiro-Wilk normality test

data:  residuals(model1)
W = 0.65623, p-value = 7.067e-10

# Jarque-Bera test normality rezíduí

jarque.bera.test(residuals(model1))

    Jarque Bera Test

data:  residuals(model1)
X-squared = 661.83, df = 2, p-value < 2.2e-16

# Test heteroskedasticity rezíduí (Breusch-Paganov test)

bptest(model1)

    studentized Breusch-Pagan test

data:  model1
BP = 4.9401, df = 3, p-value = 0.1762

# Odhad koeficientov s robustnými smerodajnými chybami (HC1)

coeftest(model1, vcov = vcovHC(model1, type = "HC1"))

t test of coefficients:

             Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.058906   0.079550  0.7405  0.46254  
podiel_zien  0.213180   0.385813  0.5525  0.58308  
ROA          1.558106   0.738731  2.1092  0.04006 *
EBITDAmarza -0.362146   0.227874 -1.5892  0.11844  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Test prítomnosti autokorelácie rezíduí (Durbin-Watsonov test)

durbinWatsonTest(model1)

 lag Autocorrelation D-W Statistic p-value
   1       0.3252467      1.347163    0.03
 Alternative hypothesis: rho != 0

# Identifikácia pozorovaní s vysokou vplyvnosťou a odľahlosťou

influencePlot(model1, main = "Vplyvné a odľahlé pozorovania")

NA

# Prehľad vplyvných pozorovaní podľa Cookovej vzdialenosti

cooks <- cooks.distance(model1)

# Zobrazenie pozorovaní s vysokou hodnotou Cookovej vzdialenosti

vplyvne <- which(cooks > 4 * mean(cooks, na.rm = TRUE))
firmy[vplyvne, ]

NA

summary(model1)

Call:
lm(formula = ROE ~ podiel_zien + ROA + EBITDAmarza, data = firmy)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.25674 -0.07073 -0.00477  0.05133  0.82942 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.05891    0.08947   0.658   0.5133   
podiel_zien  0.21318    0.28974   0.736   0.4654   
ROA          1.55811    0.50071   3.112   0.0031 **
EBITDAmarza -0.36215    0.41051  -0.882   0.3820   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4214 on 49 degrees of freedom
Multiple R-squared:  0.2117,    Adjusted R-squared:  0.1634 
F-statistic: 4.385 on 3 and 49 DF,  p-value: 0.008219

Zhrnutie výsledkov

Na základe výsledkov regresnej analýzy možno konštatovať, že:

• Premenná podiel žien vo vedení má síce pozitívny koeficient, no jej vplyv na ROE nie je štatisticky významný (p = 0.4654).
  
• Premenná ROA má pozitívny a štatisticky významný vplyv na ROE (p = 0.0031). Vyššia rentabilita aktív teda vedie k vyššej rentabilite vlastného kapitálu.
  
• Premenná EBITDA marža má negatívny, ale štatisticky nevýznamný vplyv na ROE (p = 0.3820).
  
• Model ako celok je štatisticky významný (p = 0.0082) a vysvetľuje približne 16 % variability ROE (Adjusted R² = 0.1634).

Celkovo možno povedať, že rentabilita aktív (ROA) má preukázateľne najsilnejší vplyv na rentabilitu vlastného kapitálu, zatiaľ čo podiel žien vo vedení a EBITDA marža nevykazujú štatisticky významný efekt.

shapiro.test(residuals(model1))

    Shapiro-Wilk normality test

data:  residuals(model1)
W = 0.65623, p-value = 7.067e-10

jarque.bera.test(residuals(model1))

    Jarque Bera Test

data:  residuals(model1)
X-squared = 661.83, df = 2, p-value < 2.2e-16

bptest(model1)

    studentized Breusch-Pagan test

data:  model1
BP = 4.9401, df = 3, p-value = 0.1762

durbinWatsonTest(model1)

 lag Autocorrelation D-W Statistic p-value
   1       0.3252467      1.347163   0.036
 Alternative hypothesis: rho != 0

Diagnostika modelu

Na základe vykonaných diagnostických testov možno zhrnúť, že:

• Podľa Shapiro-Wilk a Jarque-Bera testu rezíduá modelu nie sú normálne rozdelené, čo naznačuje porušenie predpokladu normality.
• Breusch-Paganov test nepotvrdil prítomnosť heteroskedasticity (p = 0.1762), teda rozptyl rezíduí možno považovať za približne konštantný.
• Durbin-Watsonov test ukázal hodnotu 1.35 s p-hodnotou 0.036, čo naznačuje autokoreláciu prvého rádu medzi rezíduami.

Celkovo možno konštatovať, že model ako celok je štatisticky významný, no niektoré predpoklady klasického lineárneho modelu (najmä normalita a nezávislosť rezíduí) nie sú úplne splnené. Tieto zistenia je potrebné zohľadniť pri interpretácii výsledkov.

Heteroskedasticita

Prítomnosť heteroskedasticity (nekonštantného rozptylu rezíduí) môže viesť k nespoľahlivým odhadom smerodajných chýb a následne k nesprávnym záverom o štatistickej významnosti regresných koeficientov.
V tejto časti overujeme, či je v modeli prítomná heteroskedasticita, a ak áno, ukážeme spôsoby jej riešenia.

library(lmtest)
library(sandwich)
library(ggplot2)
library(patchwork)

# vizuálne posúdenie - závislosť štvorcov rezíduí od vysvetľujúcich premenných
p1 <- ggplot(firmy, aes(x = ROA, y = resid(model1)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(title = "Squared residuals vs ROA", x = "ROA", y = "Squared residuals") +
  theme_minimal()

p2 <- ggplot(firmy, aes(x = podiel_zien, y = resid(model1)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(title = "Squared residuals vs podiel žien", x = "Podiel žien vo vedení", y = "Squared residuals") +
  theme_minimal()

p1 + p2

Z grafov je možné posúdiť, či existuje systematický trend medzi rozptylom rezíduí a vysvetľujúcimi premennými. Ak by červená krivka mala zreteľný rastúci alebo klesajúci tvar, mohlo by to naznačovať prítomnosť heteroskedasticity.

# Breusch–Paganov test heteroskedasticity
bptest(model1)

    studentized Breusch-Pagan test

data:  model1
BP = 4.9401, df = 3, p-value = 0.1762

Na základe výsledku testu (p = 0.1762) možno konštatovať, že heteroskedasticita v modeli nebola preukázaná. Pre účely demonštrácie však môžeme ukázať, ako by sa heteroskedasticita riešila, ak by bola prítomná.
Najčastejším postupom je použitie tzv. Whiteovho (robustného) odhadu kovariančnej matice, ktorý upravuje smerodajné chyby koeficientov.

# Odhad koeficientov s robustnými smerodajnými chybami
coeftest(model1, vcov = vcovHC(model1, type = "HC1"))

t test of coefficients:

             Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.058906   0.079550  0.7405  0.46254  
podiel_zien  0.213180   0.385813  0.5525  0.58308  
ROA          1.558106   0.738731  2.1092  0.04006 *
EBITDAmarza -0.362146   0.227874 -1.5892  0.11844  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Graf štvorcov rezíduí voči predikovaným hodnotám ROE
plot(
  fitted(model1),
  residuals(model1)^2,
  xlab = "Odhadované hodnoty ROE",
  ylab = "Štvorce rezíduí",
  main = "Štvorce rezíduí vs. odhadované hodnoty"
)

library(sandwich)
library(lmtest)

# Odhad koeficientov s robustnými smerodajnými chybami (Whiteova korekcia)
coeftest(model1, vcov = vcovHC(model1, type = "HC1"))

t test of coefficients:

             Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.058906   0.079550  0.7405  0.46254  
podiel_zien  0.213180   0.385813  0.5525  0.58308  
ROA          1.558106   0.738731  2.1092  0.04006 *
EBITDAmarza -0.362146   0.227874 -1.5892  0.11844  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

``` Na základe výsledku Breusch–Paganovho testu (p-hodnota 0.1762) možno konštatovať, že heteroskedasticita rezíduí v modeli nebola preukázaná. Z tohto dôvodu predpoklad konštantného rozptylu náhodnej zložky považujeme za splnený. Pre demonštráciu postupu však boli vypočítané aj robustné smerodajné chyby, ktoré by bolo vhodné používať v prípade, ak by heteroskedasticita bola prítomná.