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.
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
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.
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()
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.
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()
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.
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é.
acf(residuals(ols_time), main = "ACF rezíduí časového modelu")
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
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
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.