Econometrics in R – domáca úloha

Budeme odhadovať model celkových cestovných nákladov v závislosti od dĺžky pobytu a veku cestujúcej osoby:

\[ \textbf{TotalCost} = \text{Accommodation cost} + \text{Transportation cost} \]

\[ \text{TotalCost} = \beta_0 + \beta_1 \cdot \text{Duration (days)} + \beta_2 \cdot \text{Traveler age} + u \]

Knižnice a nastavenie prostredia

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
rm(list = ls())

Načítanie a príprava databázy

# načítanie dát 
data_raw <- read.csv("Travel_data.csv", sep = ";", header = TRUE, stringsAsFactors = FALSE)

# rýchla kontrola štruktúry
str(data_raw)

'data.frame':   137 obs. of  13 variables:
 $ Trip.ID             : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Destination         : chr  "London, UK" "Phuket, Thailand" "Bali, Indonesia" "New York, USA" ...
 $ Start.date          : chr  "01/05/2023" "15/06/2023" "01/07/2023" "15/08/2023" ...
 $ End.date            : chr  "08/05/2023" "20/06/2023" "08/07/2023" "29/08/2023" ...
 $ Duration..days.     : int  7 5 7 14 7 5 10 7 7 7 ...
 $ Traveler.name       : chr  "John Smith" "Jane Doe" "David Lee" "Sarah Johnson" ...
 $ Traveler.age        : int  35 28 45 29 26 42 33 25 31 39 ...
 $ Traveler.gender     : chr  "Male" "Female" "Male" "Female" ...
 $ Traveler.nationality: chr  "American" "Canadian" "Korean" "British" ...
 $ Accommodation.type  : chr  "Hotel" "Resort" "Villa" "Hotel" ...
 $ Accommodation.cost  : int  1200 800 1000 2000 700 1500 500 900 1200 2500 ...
 $ Transportation.type : chr  "Flight" "Flight" "Flight" "Flight" ...
 $ Transportation.cost : int  600 500 700 1000 200 800 1200 600 200 800 ...

head(data_raw, 5)

Tvorba analytickej tabuľky

• Vytvoríme TotalCost = Accommodation cost + Transportation cost
• Pre istotu premenujeme stĺpce na syntakticky bezpečné názvy
• Vyberieme len potrebné premenné do modelu
• Urobíme imputáciu mediánom pre chýbajúce hodnoty (ak sa vyskytnú)

# bezpečné názvy stĺpcov
names(data_raw) <- make.names(names(data_raw))

# konverzie číselných premenných (ak by pri importe ostali ako text)
num_cols <- c("Duration..days.", "Traveler.age", "Accommodation.cost", "Transportation.cost")
for (cl in num_cols) {
  if (cl %in% names(data_raw)) {
    data_raw[[cl]] <- suppressWarnings(as.numeric(data_raw[[cl]]))
  }
}

# výpočet TotalCost
data_raw$TotalCost <- with(data_raw, Accommodation.cost + Transportation.cost)

# pracovná tabuľka len s potrebnými premennými
dat <- subset(data_raw, select = c(TotalCost, Duration..days., Traveler.age))

# imputácia mediánom pre NA
col_meds <- sapply(dat, function(x) if (is.numeric(x)) median(x, na.rm = TRUE) else NA)
for (cl in names(dat)) {
  if (is.numeric(dat[[cl]])) {
    nas <- is.na(dat[[cl]])
    if (any(nas)) dat[[cl]][nas] <- col_meds[cl]
  }
}

summary(dat)

   TotalCost     Duration..days. 
 Min.   :  200   Min.   : 5.000  
 1st Qu.: 1000   1st Qu.: 7.000  
 Median : 1400   Median : 7.000  
 Mean   : 1895   Mean   : 7.606  
 3rd Qu.: 1900   3rd Qu.: 8.000  
 Max.   :10500   Max.   :14.000  
  Traveler.age  
 Min.   :20.00  
 1st Qu.:28.00  
 Median :31.00  
 Mean   :33.18  
 3rd Qu.:38.00  
 Max.   :60.00  

Rýchla vizualizácia – boxploty

par(mfrow = c(1, 3), mar = c(4, 4, 2, 1))
boxplot(dat$TotalCost, main = "TotalCost", xlab = "", col = "lightblue")
boxplot(dat$Duration..days., main = "Duration (days)", xlab = "", col = "lightgreen")

boxplot(dat$Traveler.age, main = "Traveler age", xlab = "", col = "lightpink")
par(mfrow = c(1, 1))

Lineárna regresia (Model 1)

model1 <- lm(TotalCost ~ 1 + Duration..days. + Traveler.age, data = dat)
summary(model1)

Call:
lm(formula = TotalCost ~ 1 + Duration..days. + Traveler.age, 
    data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-1742.6  -869.3  -483.8   111.0  8594.2 

Coefficients:
                Estimate Std. Error
(Intercept)     2401.686   1123.981
Duration..days. -103.226     98.862
Traveler.age       8.395     22.155
                t value Pr(>|t|)  
(Intercept)       2.137   0.0344 *
Duration..days.  -1.044   0.2983  
Traveler.age      0.379   0.7053  
---
Signif. codes:  
  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’
  0.1 ‘ ’ 1

Residual standard error: 1833 on 134 degrees of freedom
Multiple R-squared:  0.009965,  Adjusted R-squared:  -0.004811 
F-statistic: 0.6744 on 2 and 134 DF,  p-value: 0.5112

Diagnostické grafy

# Diagnostické grafy – vlastná verzia s lepšou mierkou
par(mfrow = c(2, 2), mar = c(4,4,2,1))

# 1️⃣ Residuals vs Fitted
plot(fitted(model1), residuals(model1), 
     pch = 19, col = "darkred",
     xlab = "Fitted values", ylab = "Residuals",
     main = "Residuals vs Fitted")
abline(h = 0, col = "gray40", lty = 2)

# 2️⃣ Q-Q plot
qqnorm(rstandard(model1), pch = 19, col = "darkblue",
       main = "Q-Q Plot of Standardized Residuals")
qqline(rstandard(model1), col = "gray40", lty = 2)

# 3️⃣ Scale-Location plot
plot(fitted(model1), sqrt(abs(rstandard(model1))),
     pch = 19, col = "darkgreen",
     xlab = "Fitted values", ylab = expression(sqrt("|Standardized residuals|")),
     main = "Scale-Location")
abline(h = 0, col = "gray40", lty = 2)

# 4️⃣ Residuals vs Leverage
plot(hatvalues(model1), rstandard(model1),
     pch = 19, col = "purple",
     xlab = "Leverage", ylab = "Standardized residuals",
     main = "Residuals vs Leverage")
abline(h = 0, col = "gray40", lty = 2)

par(mfrow = c(1, 1))

NA
NA

Testy normality a odľahlostí

res1 <- residuals(model1)

# Jarque–Bera test normality
jb1 <- jarque.bera.test(res1)
jb1

    Jarque Bera Test

data:  res1
X-squared = 628.91, df = 2, p-value <
2.2e-16

# Outlier test (Bonferroni)
out1 <- outlierTest(model1)
out1

Alternatívny model so zmiernením vplyvu odľahlých hodnôt (Model 2)

Ak sú náklady pozitívne a pravdepodobne prahovo posunuté a šikmé (pravostranná šikmosť), log-transformácia cieľa môže pomôcť.

# zabezpečenie, aby bolo všetko kladné (ak by boli nuly, pripočítame malú konštantu)
eps <- 1e-6
dat$TotalCost_pos <- ifelse(dat$TotalCost <= 0, dat$TotalCost + abs(min(dat$TotalCost)) + eps, dat$TotalCost)

model2 <- lm(log(TotalCost_pos) ~ 1 + Duration..days. + Traveler.age, data = dat)
summary(model2)

Call:
lm(formula = log(TotalCost_pos) ~ 1 + Duration..days. + Traveler.age, 
    data = dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.95511 -0.36180 -0.03388  0.27143  2.04096 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      7.024460   0.418866  16.770   <2e-16 ***
Duration..days. -0.017535   0.036842  -0.476    0.635    
Traveler.age     0.011721   0.008256   1.420    0.158    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.683 on 134 degrees of freedom
Multiple R-squared:  0.01788,   Adjusted R-squared:  0.003225 
F-statistic:  1.22 on 2 and 134 DF,  p-value: 0.2985

# Diagnostické grafy pre model2
par(mfrow = c(2, 2), mar = c(4,4,2,1))

# Residuals vs Fitted
plot(fitted(model2), residuals(model2), 
     pch = 19, col = "darkred",
     xlab = "Fitted values", ylab = "Residuals",
     main = "Residuals vs Fitted")
abline(h = 0, col = "gray40", lty = 2)

# Q-Q plot
qqnorm(rstandard(model2), pch = 19, col = "darkblue",
       main = "Q-Q Plot of Standardized Residuals")
qqline(rstandard(model2), col = "gray40", lty = 2)

# Scale-Location plot
plot(fitted(model2), sqrt(abs(rstandard(model2))),
     pch = 19, col = "darkgreen",
     xlab = "Fitted values", ylab = expression(sqrt("|Standardized residuals|")),
     main = "Scale-Location")
abline(h = 0, col = "gray40", lty = 2)

# Residuals vs Leverage
plot(hatvalues(model2), rstandard(model2),
     pch = 19, col = "purple",
     xlab = "Leverage", ylab = "Standardized residuals",
     main = "Residuals vs Leverage")
abline(h = 0, col = "gray40", lty = 2)

par(mfrow = c(1, 1))

# Testy pre model2
res2 <- residuals(model2)
jb2 <- jarque.bera.test(res2)
jb2

    Jarque Bera Test

data:  res2
X-squared = 18.827, df = 2, p-value = 8.16e-05

out2 <- outlierTest(model2)
out2

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:

Interpretácia

• Očakávame, že β1 pri Duration (days) bude kladný – dlhší pobyt zvyšuje celkové náklady.
• Koeficient pri Traveler age môže byť ± – ak s vekom skôr volíme komfortnejšie ubytovanie/dopravu, môže byť kladný; ak sú mladší cestujúci skôr „low-cost“, môže vychádzať záporný.
• Porovnajte Model 1 vs. Model 2 – ak log-transformácia zlepší diagnostiku (lepší QQ-plot, menšia heteroskedasticita), preferujte Model 2.

Záver

• Na našich dátach sme testovali, že dĺžka pobytu a vek cestujúcej osoby štatisticky významne formujú celkové náklady cesty.
• Diagnostické grafy a testy, najmä Jarque–Bera a outlierTest, pomohli posúdiť vhodnosť predpokladov a vplyv možných odľahlých hodnôt.