1 Cieľ práce

Cieľom práce je overiť kvalitu lineárneho regresného modelu na firemných dátach. Model vysvetľuje premennú Profit pomocou premenných Sales, Price, Order Qty a Unit Cost. Následne testujeme heteroskedasticitu a autokoreláciu a použijeme robustné štandardné chyby.

2 Načítanie dát

data <- read.csv("company_data.csv")
names(data) <- make.names(names(data))

# úprava dátumu
data$Order.Date <- as.Date(data$Order.Date, format = "%m/%d/%Y")

# ponecháme len potrebné premenné a odstránime chýbajúce hodnoty
data_model <- data %>%
  select(Order.Date, Profit, Sales, Price, Order.Qty, Unit.Cost, Cost.of.Sales) %>%
  na.omit()

head(data_model)
##   Order.Date     Profit     Sales  Price Order.Qty   Unit.Cost Cost.of.Sales
## 1 2017-09-13 2029.86529 2714.7200 304.00         9  76.0949677     684.85471
## 2 2016-08-20   20.17439   50.1414  12.99         4   7.4917528      29.96701
## 3 2016-07-08 1304.01176 1395.1128 159.99         9  10.1223377      91.10104
## 4 2018-08-11  452.04924  462.4200  25.69        18   0.5761533      10.37076
## 5 2017-07-15 1637.82101 2614.4000 304.00         9 108.5087768     976.57899
## 6 2018-10-02 5184.75931 7056.4000 299.00        24  77.9850288    1871.64069
summary(data_model[, c("Profit", "Sales", "Price", "Order.Qty", "Unit.Cost")])
##      Profit            Sales              Price           Order.Qty      
##  Min.   : -838.4   Min.   :    4.75   Min.   :   0.95   Min.   :   4.00  
##  1st Qu.:  242.2   1st Qu.:  645.00   1st Qu.:  59.00   1st Qu.:   9.00  
##  Median : 1012.7   Median : 2301.18   Median : 205.00   Median :  10.00  
##  Mean   : 2105.8   Mean   : 3692.78   Mean   : 296.51   Mean   :  16.74  
##  3rd Qu.: 2573.1   3rd Qu.: 4767.72   3rd Qu.: 366.00   3rd Qu.:  13.00  
##  Max.   :55692.6   Max.   :78312.00   Max.   :2899.99   Max.   :1560.00  
##    Unit.Cost        
##  Min.   :1.580e-03  
##  1st Qu.:1.685e+01  
##  Median :6.969e+01  
##  Mean   :1.244e+02  
##  3rd Qu.:1.643e+02  
##  Max.   :1.896e+03

3 Základný regresný model OLS

ols_model <- lm(Profit ~ Sales + Price + Order.Qty + Unit.Cost, data = data_model)
summary(ols_model)
## 
## Call:
## lm(formula = Profit ~ Sales + Price + Order.Qty + Unit.Cost, 
##     data = data_model)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12185.6    -70.5      4.2    103.2  10113.6 
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  11.408717   8.210604    1.390    0.165    
## Sales         0.591658   0.002094  282.614  < 2e-16 ***
## Price         4.846632   0.035098  138.090  < 2e-16 ***
## Order.Qty    -0.576843   0.139377   -4.139 3.51e-05 ***
## Unit.Cost   -12.199189   0.054168 -225.212  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 715.3 on 14995 degrees of freedom
## Multiple R-squared:  0.9554, Adjusted R-squared:  0.9554 
## F-statistic: 8.032e+04 on 4 and 14995 DF,  p-value: < 2.2e-16

Model OLS ukazuje, ako tržby, cena, počet objednaných kusov a jednotkové náklady súvisia so ziskom. Samotné koeficienty však nestačia — treba ešte overiť, či model neporušuje základné predpoklady.

4 Graf vzťahu medzi Sales a Profit

ggplot(data_model, aes(x = Sales, y = Profit)) +
  geom_point(alpha = 0.25) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    title = "Vzťah medzi tržbami a ziskom",
    x = "Sales",
    y = "Profit"
  ) +
  theme_minimal()

5 Test heteroskedasticity

Heteroskedasticita znamená, že rozptyl rezíduí nie je konštantný. Overíme ju Breusch-Paganovým testom.

bptest(ols_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  ols_model
## BP = 4682.2, df = 4, p-value < 2.2e-16

Ak je p-hodnota menšia ako 0,05, zamietame nulovú hypotézu o homoskedasticite a v modeli sa pravdepodobne nachádza heteroskedasticita.

5.1 Graf rezíduí

data_model$resid_ols <- residuals(ols_model)
data_model$fitted_ols <- fitted(ols_model)

ggplot(data_model, aes(x = fitted_ols, y = resid_ols)) +
  geom_point(alpha = 0.25) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_smooth(se = FALSE) +
  labs(
    title = "Rezíduá voči odhadovaným hodnotám",
    x = "Odhadované hodnoty",
    y = "Rezíduá"
  ) +
  theme_minimal()

6 Robustné štandardné chyby

Ak je prítomná heteroskedasticita, použijeme Whiteove robustné štandardné chyby.

coeftest(ols_model, vcov = vcovHC(ols_model, type = "HC1"))
## 
## t test of coefficients:
## 
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  11.408717  10.929724   1.0438   0.2966    
## Sales         0.591658   0.012394  47.7386   <2e-16 ***
## Price         4.846632   0.149748  32.3653   <2e-16 ***
## Order.Qty    -0.576843   0.442303  -1.3042   0.1922    
## Unit.Cost   -12.199189   0.153958 -79.2370   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robustné štandardné chyby nemenia samotné koeficienty modelu, ale upravujú štandardné chyby, t-hodnoty a p-hodnoty tak, aby boli spoľahlivejšie pri heteroskedasticite.

7 Autokorelácia rezíduí

Keďže dáta obsahujú dátum objednávky, vytvoríme denný dataset a otestujeme autokoreláciu rezíduí.

daily_data <- data_model %>%
  group_by(Order.Date) %>%
  summarise(
    Profit = sum(Profit, na.rm = TRUE),
    Sales = sum(Sales, na.rm = TRUE),
    Price = mean(Price, na.rm = TRUE),
    Order.Qty = sum(Order.Qty, na.rm = TRUE),
    Unit.Cost = mean(Unit.Cost, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(Order.Date) %>%
  na.omit()

ols_time <- lm(Profit ~ Sales + Price + Order.Qty + Unit.Cost, data = daily_data)
summary(ols_time)
## 
## Call:
## lm(formula = Profit ~ Sales + Price + Order.Qty + Unit.Cost, 
##     data = daily_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17099.1  -1464.7     27.7   1449.0  17094.9 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.920e+01  3.453e+02   0.056  0.95568    
## Sales        5.830e-01  5.145e-03 113.314  < 2e-16 ***
## Price        5.784e+01  1.539e+00  37.582  < 2e-16 ***
## Order.Qty   -1.842e+00  5.958e-01  -3.091  0.00204 ** 
## Unit.Cost   -1.404e+02  3.114e+00 -45.092  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3251 on 1091 degrees of freedom
## Multiple R-squared:  0.9506, Adjusted R-squared:  0.9504 
## F-statistic:  5244 on 4 and 1091 DF,  p-value: < 2.2e-16
bgtest(ols_time, order = 1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  ols_time
## LM test = 4.2701, df = 1, p-value = 0.03879

Ak je p-hodnota menšia ako 0,05, znamená to možný problém autokorelácie rezíduí. Vtedy môžu byť bežné štandardné chyby nespoľahlivé.

7.1 ACF graf rezíduí

acf(residuals(ols_time), main = "ACF rezíduí časového modelu")

8 Newey-West štandardné chyby

Pri možnej autokorelácii použijeme Newey-West korekciu.

coeftest(ols_time, vcov = NeweyWest(ols_time, lag = 1, prewhite = FALSE))
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value Pr(>|t|)    
## (Intercept)  1.9195e+01  3.4524e+02   0.0556  0.95567    
## Sales        5.8305e-01  8.6624e-03  67.3083  < 2e-16 ***
## Price        5.7839e+01  2.2664e+00  25.5199  < 2e-16 ***
## Order.Qty   -1.8417e+00  6.9923e-01  -2.6339  0.00856 ** 
## Unit.Cost   -1.4044e+02  4.8148e+00 -29.1680  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

9 Krátka ukážka GLS modelu

GLS model je alternatíva k OLS, keď sú porušené predpoklady o rezíduách. Tu použijeme model s AR(1) štruktúrou rezíduí.

gls_ar1 <- try(
  gls(Profit ~ Sales + Price + Order.Qty + Unit.Cost,
      data = daily_data,
      correlation = corAR1(form = ~ as.numeric(Order.Date))),
  silent = TRUE
)

if (inherits(gls_ar1, "try-error")) {
  print("GLS model sa nepodarilo odhadnúť na týchto dátach.")
} else {
  summary(gls_ar1)
}
## Generalized least squares fit by REML
##   Model: Profit ~ Sales + Price + Order.Qty + Unit.Cost 
##   Data: daily_data 
##        AIC      BIC    logLik
##   20832.41 20867.38 -10409.21
## 
## Correlation Structure: AR(1)
##  Formula: ~as.numeric(Order.Date) 
##  Parameter estimate(s):
##        Phi 
## 0.06362658 
## 
## Coefficients:
##                  Value Std.Error   t-value p-value
## (Intercept)   19.85637  348.4441   0.05699  0.9546
## Sales          0.58295    0.0052 112.25934  0.0000
## Price         57.76023    1.5441  37.40813  0.0000
## Order.Qty     -1.79216    0.5976  -2.99886  0.0028
## Unit.Cost   -140.30430    3.1059 -45.17410  0.0000
## 
##  Correlation: 
##           (Intr) Sales  Price  Ordr.Q
## Sales     -0.123                     
## Price     -0.348 -0.337              
## Order.Qty -0.402 -0.352  0.157       
## Unit.Cost -0.168 -0.033 -0.713  0.056
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -5.26192685 -0.45109790  0.01046962  0.44504211  5.26055941 
## 
## Residual standard error: 3251.535 
## Degrees of freedom: 1096 total; 1091 residual

10 Záver

V práci sme vytvorili lineárny regresný model, ktorý vysvetľuje zisk firmy pomocou tržieb, ceny, množstva objednávok a jednotkových nákladov. Následne sme overili dva dôležité problémy regresnej analýzy: heteroskedasticitu a autokoreláciu. Pri porušení predpokladov sme použili robustné štandardné chyby a Newey-West korekciu. Na záver sme ukázali aj GLS model ako alternatívny spôsob odhadu.