Načítanie balíkov

library(lmtest)
library(car)
library(tseries)

##Predošlý dataset - divné grafy

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
udaje <- read.csv2("udaje/productivity.csv",header=TRUE,sep=",",dec=".")
udaje <- udaje %>%
  rename_with(~ abbreviate(.x, strict = FALSE))
model_stary<- lm(Pr_S ~ Mn_S + E_S_+ Sc_D, data = udaje)
summary(model_stary)
## 
## Call:
## lm(formula = Pr_S ~ Mn_S + E_S_ + Sc_D, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.10432 -0.94079 -0.04562  0.97725  2.61627 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.777e-01  2.180e-02 -17.323   <2e-16 ***
## Mn_S         5.257e-04  2.804e-06 187.503   <2e-16 ***
## E_S_         1.408e-03  3.344e-03   0.421    0.674    
## Sc_D         3.769e-04  8.884e-04   0.424    0.671    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.217 on 99996 degrees of freedom
## Multiple R-squared:  0.2601, Adjusted R-squared:  0.2601 
## F-statistic: 1.172e+04 on 3 and 99996 DF,  p-value: < 2.2e-16

Kedže som s predošlým datasetom mala problém, že mi tam vychdádzal iba mesačný plat ako štatisticky významný a robilo mi to problémy v grafoch, rozhodla som sa zmeniť dataset, ak by s tým bol problém, rada to zmením, ale nevedela som si poradiť s grafmi, ktoré mi vznikli pri predošlom modeli.

plot(model_stary, which = 1)
Residuals vs Fitted- povodny

Residuals vs Fitted- povodny

Nadčítanie údajov a model

library(dplyr)
udaje_nove <- read.csv2("udaje/company_data.csv",header=TRUE,sep=",",dec=".")
model <- lm(Profit ~ Price + Order.Qty + Unit.Cost, data = udaje_nove)
summary(model)
## 
## Call:
## lm(formula = Profit ~ Price + Order.Qty + Unit.Cost, data = udaje_nove)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -17243   -365   -100    -23  42237 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -28.49945   20.64798   -1.38    0.168    
## Price        11.62553    0.06444  180.40   <2e-16 ***
## Order.Qty     9.48104    0.33894   27.97   <2e-16 ***
## Unit.Cost   -11.82672    0.13620  -86.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1799 on 14996 degrees of freedom
## Multiple R-squared:  0.7179, Adjusted R-squared:  0.7178 
## F-statistic: 1.272e+04 on 3 and 14996 DF,  p-value: < 2.2e-16

Uvažujeme lineárny regresný model, v ktorom vysvetľujeme hodnotu Profit (dosiahnutý zisk z objednávky) pomocou jednotkovej ceny produktu, objednaného množstva a jednotkových nákladov. Model má nasledovný tvar:\(Profit_i = \beta_0 + \beta_1 Price_i + \beta_2 Order.Qty_i + \beta_3 Unit.Cost_i + \varepsilon_i\) Pri jednotkovej cene predpokladáme pozitívny vplyv na zisk \((\beta_1 > 0)\). Vyššia predajná cena pri konštantných nákladoch priamo zvyšuje maržu, čo vedie k vyššiemu zisku z každej predanej jednotky. Pri objednanom množstve očakávame pozitívny vplyv na celkový zisk \((\beta_2 > 0)\). S narastajúcim počtom predaných kusov v rámci jednej objednávky rastie celkový objem vytvoreného zisku. Pri premennej jednotkové náklady na jednotku predpokladáme negatívny vplyv \((\beta_3 < 0)\). Vyššie náklady na výrobu produktu priamo znižujú rozdiel medzi predajnou cenou a nákladmi, čím znižujú výsledný zisk.

Diagnostické grafy

Residuals vs Fitted

plot(model, which = 1)
Residuals vs Fitted

Residuals vs Fitted

Rezíduá sú rovnomerne rozložené okolo nulovej hodnoty, čo potvrdzuje, že model nemá výrazné systematické skreslenie v predikciách zisku. Väčšina bodov sa nachádza v relatívne úzkom pásme okolo nuly, čo je priaznivý znak.

Červená vyhladzovacia (LOESS) krivka je takmer dokonale lineárna a kopíruje nulovú líniu bez výrazných výkyvov. Nie je tu prítomná nelinearita, čo znamená, že lineárny vzťah medzi cenou, množstvom, nákladmi a ziskom je zvolený správne.

V grafe sa nachádza niekoľko označených extrémnych hodnôt (napr. pozorovanie č. 13812, 11352, 16508), ktoré môžu predstavovať potenciálne odľahlé alebo vplyvné pozorovania. Tieto by bolo vhodné ďalej preskúmať, napríklad pomocou Cookovej vzdialenosti alebo testu na odľahlé hodnoty.

Celkovo je však model špecifikovaný veľmi dobre a popisuje štruktúru firemných dát.

Q-Q plot

plot(model, which = 2)
Normal Q-Q plot

Normal Q-Q plot

Q-Q graf porovnáva empirické rozdelenie štandardizovaných rezíduí s teoretickým normálnym rozdelením. Väčšina bodov v strednej časti grafu takmer dokonale kopírujú referenčnú priamku, čo naznačuje, že pre väčšinu bežných objednávok model spĺňa predpoklad normality chýb.

Na pravej strane body stúpajú nad priamku a na ľavej klesajú pod ňu, čo indikuje prítomnosť extrémnych hodnôt v datasete.

Scale-Location plot

plot(model, which = 3)
Scale-Location plot

Scale-Location plot

Scale-Location graf slúži na posúdenie homoskedasticity, teda konštantnosti rozptylu rezíduí v závislosti od prispôsobených hodnôt. Koncentrácia väčšiny pozorovaní v ľavej časti grafu odzrkadľuje realitu predajných dát organizácie, kde dominujú objednávky s nižšou jednotkovou cenou a ziskom. Červená krivka má mierne rastúci trend. To naznačuje, že s rastúcimi fitted hodnotami sa mierne zvyšuje aj variabilita rezíduí. Ide teda o náznak slabej heteroskedasticity, avšak nie o výrazné porušenie predpokladov lineárneho regresného modelu.

Residuals vs Leverage

plot(model, which = 5)
Residuals vs Leverage

Residuals vs Leverage

Graf Residuals vs Leverage slúži na identifikáciu vplyvných pozorovaní. Väčšina pozorovaní má veľmi nízky leverage a je sústredená blízko ľavej časti grafu, čo naznačuje, že väčšina dát je dobre vyvážená.

Poznámka k LOESS krivke

Červená LOESS krivka predstavuje lokálne vyhladený trend medzi premennými bez predpokladu konkrétneho funkčného tvaru. Slúži na vizuálnu identifikáciu systematických odchýlok od náhodného rozloženia rezíduí. Ak je približne horizontálna, model je špecifikovaný vhodne; ak sa zakrivuje, môže to naznačovať nelinearitu alebo iný problém v modeli.

Test normality rezíduí

Shapiro-Wilkov test

shapiro.test(sample(residuals(model), 5000))
## 
##  Shapiro-Wilk normality test
## 
## data:  sample(residuals(model), 5000)
## W = 0.41804, p-value < 2.2e-16

Shapiro-Wilkov test overuje nulovú hypotézu, že reziduály pochádzajú z normálneho rozdelenia. Ak je p-hodnota menšia než zvolená hladina významnosti (napríklad 0.05), nulovú hypotézu zamietame a usudzujeme, že reziduály nie sú normálne rozdelené.

Jarque-Bera test

jarque.bera.test(residuals(model))
## 
##  Jarque Bera Test
## 
## data:  residuals(model)
## X-squared = 5036890, df = 2, p-value < 2.2e-16

Jarque-Bera test je ďalší test normality, ktorý je založený na šikmosti a špicatosti rezíduí. Malá p-hodnota signalizuje porušenie predpokladu normálneho rozdelenia.

Test heteroskedasticity

Breusch-Pagan test

bptest(model)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 1106.1, df = 3, p-value < 2.2e-16

Breusch-Pagan test skúma nulovú hypotézu homoskedasticity, teda konštantného rozptylu rezíduí. Ak je p-hodnota malá, nulovú hypotézu zamietame a usudzujeme, že v modeli je prítomná heteroskedasticita.

Voliteľne: Whiteov test

bptest(model, ~ fitted(model) + I(fitted(model)^2))
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 1051.2, df = 2, p-value < 2.2e-16

Táto verzia testu rozširuje Breusch-Paganov prístup a umožňuje zachytiť aj zložitejšie formy heteroskedasticity.

Test autokorelácie

Durbin-Watson test

dwtest(model)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 2.02, p-value = 0.8898
## alternative hypothesis: true autocorrelation is greater than 0

Durbin-Watson test skúma nulovú hypotézu, že reziduály nie sú autokorelované prvého rádu. Hodnota testovej štatistiky blízka 2 naznačuje neprítomnosť autokorelácie. Hodnoty výrazne nižšie než 2 naznačujú kladnú autokoreláciu a hodnoty výrazne vyššie než 2 zápornú autokoreláciu.

Breusch-Godfrey test

bgtest(model, order = 1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model
## LM test = 1.505, df = 1, p-value = 0.2199

Breusch-Godfrey test skúma prítomnosť autokorelácie rezíduí. Je všeobecnejší než Durbin-Watson test a umožňuje testovať aj vyššie rády autokorelácie. Ak je p-hodnota malá, nulovú hypotézu o neprítomnosti autokorelácie zamietame.

Odľahlé a vplyvné pozorovania

Test na odľahlé hodnoty

outlierTest(model)
##       rstudent unadjusted p-value Bonferroni p
## 9854  23.93236        2.9789e-124  4.4683e-120
## 2048  22.61852        2.0653e-111  3.0980e-107
## 13762 17.94021         3.1753e-71   4.7630e-67
## 11181 15.65730         8.0466e-55   1.2070e-50
## 6837  15.02581         1.1627e-50   1.7441e-46
## 13488 13.85692         2.1386e-43   3.2080e-39
## 1916  13.39387         1.1248e-40   1.6873e-36
## 10302 12.58128         4.0739e-36   6.1109e-32
## 2541  12.42508         2.8492e-35   4.2738e-31
## 8962  11.77404         7.3376e-32   1.1006e-27

Tento test pomáha identifikovať pozorovania s mimoriadne veľkými štandardizovanými rezíduami.

Cookova vzdialenosť

cd <- cooks.distance(model)
head(sort(cd, decreasing = TRUE), 10)
##      5958     10635      3199     13762       758      5093      9854      9073 
## 1.5937009 0.5363113 0.2255723 0.2201817 0.1621086 0.1596306 0.1581160 0.1468108 
##     11921      5515 
## 0.1394411 0.1209062

Cookova vzdialenosť meria, nakoľko jednotlivé pozorovanie ovplyvňuje odhad regresných koeficientov. Väčšie hodnoty si zaslúžia dodatočnú pozornosť.

Krátke teoretické poznámky

Štandardizované rezíduá sú rezíduá vydelené ich odhadovanou smerodajnou odchýlkou, vďaka čomu sú medzi sebou porovnateľné a pri splnení predpokladov modelu sa správajú približne ako hodnoty zo štandardného normálneho rozdelenia.

Leverage vyjadruje, ako veľmi sa pozorovanie líši v hodnotách vysvetľujúcich premenných od ostatných a aký potenciál má ovplyvniť odhad regresného modelu. Formálne ide o diagonálne prvky tzv. hat matice, pričom hat matica je definovaná ako \(H = X(X'X)^{-1}X'\) a platí \(\hat{y} = Hy\).

Pre i-te pozorovanie je leverage daný vzťahom \(h_i = x_i'(X'X)^{-1}x_i\), kde \(x_i\) je stĺpcový vektor vysvetľujúcich premenných pre i-te pozorovanie.

Záverečná interpretačná poznámka

Pri interpretácii diagnostických testov je vhodné kombinovať grafické a formálne prístupy. Diagnostické grafy často odhalia problém intuitívne, zatiaľ čo testy poskytujú formálnejšie štatistické rozhodnutie. Pri väčších vzorkách môžu byť testy veľmi citlivé, preto je rozumné posudzovať aj ekonomický a praktický význam zistených porušení predpokladov.