S využitím databázy hotel_booking.csv database. V databáze sa nachádzajú údaje o hotelových rezerváciach, počte hostí, pošte dní pobytu, časový interval medzi rezerváciou a pobytom…

Úvod do problému, stanovenie hypotéz

Rozdhodli sme sa zamerať len na najzaujímavejšie číselné údaje spomedzi dostupných dát, a to is_canceled (či bola rezervácia zrušená), lead_time (počet dní medzi rezerváciou a pobytom), arrival_date_week_number (číslo týždňa pobytu v roku), stays_in_weekend_nights (počet víkondových dní pobytu), stays_in_week_nights (počet pracovných dní pobytu), previous_cancellations (počet predošlých zrušených rezervácií) a previous_bookings_not_canceled (počet nezrušených predošlých rezervácií).

udaje <- read.csv("Ekonometria/hotel_bookings.csv", header = TRUE, sep = ",", dec = ".")

udaje.new <- udaje[udaje$is_canceled == 1,
  c("lead_time","arrival_date_week_number","stays_in_weekend_nights",
    "stays_in_week_nights","previous_cancellations","previous_bookings_not_canceled")]

column_medians <- sapply(udaje.new, median, na.rm = TRUE)

udaje_imputed <- udaje.new
for (col in names(udaje.new)) {
  udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}
udaje.new <- udaje_imputed

sapply(udaje.new, is.numeric)  # rýchla kontrola typov

##                      lead_time       arrival_date_week_number 
##                           TRUE                           TRUE 
##        stays_in_weekend_nights           stays_in_week_nights 
##                           TRUE                           TRUE 
##         previous_cancellations previous_bookings_not_canceled 
##                           TRUE                           TRUE

Teraz chceme vidieť tvar údajov.

head(udaje.new)

# Determine number of plots
num_plots <- length(names(udaje.new))
num_plots

## [1] 6

# Set the layout: 2 rows × 2 columns
par(mfrow = c(2, 3))
par(mar = c(4, 4, 2, 1))  # Adjust margins (optional)

# Loop through columns and plot each boxplot
for (col in names(udaje.new)) {
  boxplot(udaje.new[[col]], main = col, xlab = "Value", col = "lightblue")
}

# Add a global caption / title
mtext("Boxploty jednotlivých premenných", outer = TRUE, cex = 1.4, font = 2)

# Reset layout to default (1 plot per figure)
#par(mfrow = c(1, 1))

Boxploty krásne vykreslujú stav jednotlivých premenných. Najviac výrazné sú outliery, ľudia, ktorí príliš často rušia rezervácie, alebo využívajú hotel namiesto prenájmu apartmánu na dlhú dobu.

colnames(udaje.new)

## [1] "lead_time"                      "arrival_date_week_number"      
## [3] "stays_in_weekend_nights"        "stays_in_week_nights"          
## [5] "previous_cancellations"         "previous_bookings_not_canceled"

model <- lm(lead_time ~ +1 + arrival_date_week_number + stays_in_weekend_nights + stays_in_week_nights + previous_cancellations + previous_bookings_not_canceled, data=udaje.new)
model

## 
## Call:
## lm(formula = lead_time ~ +1 + arrival_date_week_number + stays_in_weekend_nights + 
##     stays_in_week_nights + previous_cancellations + previous_bookings_not_canceled, 
##     data = udaje.new)
## 
## Coefficients:
##                    (Intercept)        arrival_date_week_number  
##                        101.342                           1.498  
##        stays_in_weekend_nights            stays_in_week_nights  
##                         -8.162                           3.395  
##         previous_cancellations  previous_bookings_not_canceled  
##                          7.649                          -5.576

Objekt triedy lm() nám poskytuje niekoľko výsledkov:

Vector odhadnutých koeficientov model$coefficients
Vektor rezíduí model$ residuals
Vektor vyrovnaných hodnôt vysvetľovanej veličiny model$fitted.values
Maticu X model$x

summary(model)

## 
## Call:
## lm(formula = lead_time ~ +1 + arrival_date_week_number + stays_in_weekend_nights + 
##     stays_in_week_nights + previous_cancellations + previous_bookings_not_canceled, 
##     data = udaje.new)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -265.77  -89.10  -24.55   68.37  504.79 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    101.34182    1.50034  67.546  < 2e-16 ***
## arrival_date_week_number         1.49792    0.04234  35.375  < 2e-16 ***
## stays_in_weekend_nights         -8.16160    0.62500 -13.059  < 2e-16 ***
## stays_in_week_nights             3.39482    0.33535  10.123  < 2e-16 ***
## previous_cancellations           7.64900    0.41774  18.310  < 2e-16 ***
## previous_bookings_not_canceled  -5.57589    0.81776  -6.818 9.32e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 116.2 on 44218 degrees of freedom
## Multiple R-squared:  0.04101,    Adjusted R-squared:  0.0409 
## F-statistic: 378.2 on 5 and 44218 DF,  p-value: < 2.2e-16

Podľa výsledku lineárnej regresie vidíme, že všetky nami vybrané údaje štatisticky významné pre sledovanie počtu dní medzi rezerváciou a pobytom.

# Nastaviť rozloženie 2 x 2
par(mfrow = c(2, 2))

# Vykresliť všetky 4 diagnostické grafy modelu
plot(model)

Diagnostické grafy regresného modelu

# (Voliteľné) pridať spoločný nadpis
#mtext("Diagnostické grafy regresného modelu", outer = TRUE, cex = 1.2, font = 2)

# Resetovať layout
par(mfrow = c(1, 1))

1. Residuals vs Fitted

Účel: overenie linearity a konštantného rozptylu reziduí
Body by mali byť náhodne rozptýlené okolo červenej čiary
V tomto prípade vidno zhlukovanie a mierny trend, čo naznačuje porušenie linearity a možnú nehomoskedasticitu (rozptyl reziduí nie je konštantný).

2. Normal Q-Q

Účel: kontrola, či majú reziduá približne normálne rozdelenie
Ak sú body na grafe v blízkosti priamky, predpoklad normality je splnený
Tu pozorujeme odchýlky v chvostoch grafu, čo poukazuje na narušenie normality reziduí

3. Scale-Location (Spread-Location)

Účel: overenie, či je rozptyl reziduí rovnaký pre všetky predikované hodnoty
Červená čiara by mala byť približne vodorovná
V našom prípade má červená čiara stúpajúci trend, čo signalizuje heteroskedasticitu.

4. Residuals vs Leverage

Účel: identifikovať vplyvné pozorovania, ktoré môžu výrazne ovplyvniť odhad parametrov
Body s vysokým „leverage“ a veľkou Cookovou vzdialenosťou (nad 0.5 alebo 1) sú potenciálne problematické
V našom prípade sa v grafe nachádza niekoľko takýchto bodov, čo znamená, že model je citlivý na niektoré pozorovania

Zhrnutie

Na základe diagnostických grafov možno konštatovať, že: - predpoklady lineárneho modelu nie sú úplne splnené - reziduá nevykazujú konštantný rozptyl - rozdelenie reziduí sa odchyľuje od normálneho
- niektoré pozorovania majú výrazný vplyv na výsledky

Čo ďalej?

Podľa výsledku lm modelu sú všetky naše dáta relevantné, ale grafy nám nevychádzajú až tak pekne, ako by sme si želali. Preto si upravíme outliery.

# normality tests
residuals <- residuals(model)
jb_test <- jarque.bera.test(residuals)
jb_test

## 
##  Jarque Bera Test
## 
## data:  residuals
## X-squared = 9237.4, df = 2, p-value < 2.2e-16

# outlier test (see p-value for Bonferroni correction)
outlier_test <- outlierTest(model)
outlier_test

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##       rstudent unadjusted p-value Bonferroni p
## 65232 4.346172         1.3884e-05      0.61402

Aby sme zistili, či sú reziduá z nášho modelu normálne rozdelené, použili sme Jarque–Bera test. Tento test v podstate skúma, či má rozdelenie reziduí „normálny tvar“. Nulová hypotéza hovorí, že reziduá sú normálne. Keďže p-hodnota vyšla extrémne malá (prakticky nulová), túto hypotézu musíme zamietnuť. Inými slovami, reziduá sa spravajú inak než pri normálnom rozdelení, čo vidno aj na Q–Q grafe, kde sa body odchyľujú od referenčnej priamky.

Potom sme pomocou funkcie outlierTest() skontrolovali, či v dátach nemáme nejaké extrémne pozorovania, ktoré by mohli mať výrazný vplyv na výsledky regresie. Tento test sleduje študentizované reziduá a kontroluje, či niektoré nie sú „príliš veľké“. V našom prípade sa však žiadne pozorovanie neukázalo ako štatisticky významný outlier. Aj keď najväčšie študentizované reziduum bolo pomerne vysoké (okolo 4.35), po Bonferroniho korekcii vyšla p-hodnota približne 0.61, čo je vysoko nad hranicou významnosti 0.05. Znamená to, že podľa tohto testu nemáme v dátach žiadny problémový outlier.

Keďže sa nepreukázala normalita rezíduí, pokúsme sa eliminovať odľahlé hodnoty. Použijeme transformáciu log1p(x), ktorá znamená log(1 + x), keďže naše premenné často obsahujú hodnotu 0, ktorá je vo všetkých prípadoch relevantným údajom.

model2 <- lm(
  log1p(lead_time) ~ 
    arrival_date_week_number +
    log1p(stays_in_weekend_nights) +
    log1p(stays_in_week_nights) +
    log1p(previous_cancellations) +
    log1p(previous_bookings_not_canceled),
  data = udaje.new
)
summary(model2)

## 
## Call:
## lm(formula = log1p(lead_time) ~ arrival_date_week_number + log1p(stays_in_weekend_nights) + 
##     log1p(stays_in_week_nights) + log1p(previous_cancellations) + 
##     log1p(previous_bookings_not_canceled), data = udaje.new)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4872 -0.5620  0.2424  0.7885  4.6538 
## 
## Coefficients:
##                                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)                            3.7290929  0.0186363 200.099   <2e-16
## arrival_date_week_number               0.0088597  0.0004159  21.302   <2e-16
## log1p(stays_in_weekend_nights)         0.0200643  0.0108917   1.842   0.0655
## log1p(stays_in_week_nights)            0.4025148  0.0118566  33.948   <2e-16
## log1p(previous_cancellations)          0.6887648  0.0186501  36.931   <2e-16
## log1p(previous_bookings_not_canceled) -1.1953641  0.0489749 -24.408   <2e-16
##                                          
## (Intercept)                           ***
## arrival_date_week_number              ***
## log1p(stays_in_weekend_nights)        .  
## log1p(stays_in_week_nights)           ***
## log1p(previous_cancellations)         ***
## log1p(previous_bookings_not_canceled) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.132 on 44218 degrees of freedom
## Multiple R-squared:  0.07479,    Adjusted R-squared:  0.07469 
## F-statistic: 714.9 on 5 and 44218 DF,  p-value: < 2.2e-16

Po logaritmickej úprave premenných sa viaceré koeficienty v modeli stále ukázali ako štatisticky významné. Najsilnejší vplyv má premenná stays_in_week_nights, ktorej koeficient je pozitívny a veľmi významný, čo naznačuje, že dlhší pobyt počas pracovných dní súvisí s dlhším časovým predstihom rezervácie. Podobne aj previous_cancellations a previous_bookings_not_canceled majú výrazný a štatisticky významný efekt, pričom prvá premenná pôsobí pozitívne a druhá negatívne. To znamená, že ľudia, ktorí v minulosti veľa rušili, rezervujú s väčším predstihom, kým tí, ktorí rezervácie nerušia, rezervujú bližšie k dátumu pobytu.

Premenná arrival_date_week_number je tiež vysoko významná, čo naznačuje, že týždeň v roku hrá určitú rolu pri dĺžke predstihu rezervácie. Jediná premenná, ktorá nevyšla ako štatisticky významná na bežnej úrovni 5 %, je log1p(stays_in_weekend_nights), hoci jej koeficient je kladný.

# Nastaviť rozloženie 2 x 2
par(mfrow = c(2, 2))

# Vykresliť všetky 4 diagnostické grafy modelu
plot(model2)

Diagnostické grafy regresného modelu 2

Po logaritmickej transformácii sa správanie reziduí vo viacerých ohľadoch zlepšilo, ale nie sú splnené úplne všetky predpoklady lineárnej regresie. V grafe Residuals vs Fitted vidíme, že rozptyl reziduí je síce vyrovnanejší než v pôvodnom modeli, no stále sa objavuje určitý náznak zakrivenia, čo naznačuje mierne porušenie linearity. Q–Q graf ukazuje, že reziduá sú bližšie k teoretickej priamke než predtým, ale odchýlky v chvostoch stále zostávajú. Graf Scale–Location naznačuje, že heteroskedasticita sa čiastočne zmiernila, hoci rozptyl reziduí nie je úplne rovnomerný pre všetky hodnoty. V grafe Residuals vs Leverage sa ukazuje niekoľko bodov s vyšším vplyvom, ktoré však nie sú extrémne problematické.

Celkovo logaritmická transformácia pomohla modelu stabilizovať variabilitu a mierne zlepšila normalitu reziduí, no model stále nie je ideálny a určité nedostatky v diagnostike pretrvávajú.

Heteroskedasticita

Prítomnosť heteroskedasticity (t. j. nekonštantného rozptylu náhodnej zložky) môže spôsobovať skreslené vyhodnocovanie t-testov významnosti jednotlivých regresných koeficientov. Preto je dôležité heteroskedasticitu odhaliť — či už vizuálne alebo pomocou formálnych testov — a v prípade jej výskytu model vhodne upraviť.

Aj v našom prípade sa snažíme vizuálne posúdiť, či sa heteroskedasticita vyskytuje. Slúžia nám na to diagnostické grafy, najmä graf Scale–Location, ktorý ukazuje, či sa rozptyl reziduí mení spolu s veľkosťou predikovaných hodnôt.

V ďalšom kroku sa preto pozrieme na to, ako sa heteroskedasticita správa v dvoch rôznych modeloch — v pôvodnom modeli (model) a v modeli so zlogaritmovanými premennými (model2). Model2 využíva logaritmickú transformáciu vybraných premenných, aby sa znížil vplyv extrémnych hodnôt a stabilizoval rozptyl. Táto úprava bola zvolená práve kvôli tomu, že pôvodné dáta obsahovali veľa veľmi nízkych a zároveň veľmi vysokých hodnôt, čo často vedie k heteroskedasticite. Model je teda pôvodný lineárny model, zatiaľ čo model2 je jeho transformovaná verzia.

library(ggplot2)

# riadky, ktoré model reálne použil
mf  <- model.frame(model)
idx <- as.integer(rownames(mf))

# dáta na kreslenie (len použité riadky) + štvorce reziduí
df_plot <- data.frame(
  previous_cancellations          = udaje.new$previous_cancellations[idx],
  previous_bookings_not_canceled  = udaje.new$previous_bookings_not_canceled[idx],
  sqres = residuals(model)^2
)

p1 <- ggplot(df_plot, aes(previous_cancellations, sqres)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "previous_cancellations", y = "Squared Residuals",
       title = "Squared Residuals vs previous_cancellations") +
  theme_minimal()

p2 <- ggplot(df_plot, aes(previous_bookings_not_canceled, sqres)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "previous_bookings_not_canceled", y = "Squared Residuals",
       title = "Squared Residuals vs previous_bookings_not_canceled") +
  theme_minimal()

# bez patchworku: vytlač po sebe, alebo vedľa seba s gridExtra
print(p1); print(p2)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 31739 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at -0.13

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 0.0169

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at -0.13

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.13

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 1

## Warning: Failed to fit group -1.
## Caused by error in `predLoess()`:
## ! NA/NaN/Inf in foreign function call (arg 5)

## Warning: Removed 31739 rows containing missing values or values outside the scale range
## (`geom_point()`).

Skúmanie heteroskedasticity

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 31739 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at -0.22

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 0.0484

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at -0.22

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.22

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger

## Warning: Failed to fit group -1.
## Caused by error in `predLoess()`:
## ! NA/NaN/Inf in foreign function call (arg 5)

## Warning: Removed 31739 rows containing missing values or values outside the scale range
## (`geom_point()`).