Špecifikácia modelu – cestovné náklady

knitr::opts_chunk$set(echo = TRUE)

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

#Popis dát a cieľ ##V tejto časti budeme modelovať, ako sú celkové náklady na zájazd (TotalCost) ovplyvnené vybranými charakteristikami cesty – najmä dĺžkou pobytu (Duration) a vekom cestovateľa (Age). Vytvoríme si teda premennú TotalCost a potom budeme krok za krokom testovať, či je lineárna špecifikácia modelu vhodná, podobne ako pri príklade s očakávanou dĺžkou života.

#PRÍPRAVA ÚDAJOV

# PRÍPRAVA ÚDAJOV
# načítanie pôvodných dát (uprav cestu k súboru podľa seba)
udaje_raw <- read.csv("Travel_data.csv",
dec = ".",
sep = ";",
header = TRUE)

# pozri si štruktúru, aby si videla presné názvy stĺpcov
str(udaje_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 ...

#Výber relevantných stĺpcov

udaje <- udaje_raw[, c("Duration..days.",
"Traveler.age",
"Accommodation.cost",
"Transportation.cost")]

# premenovanie na jednoduchšie názvy
names(udaje) <- c("Duration", "Age", "Accommodation", "Transport")

# vytvorenie celkových nákladov na zájazd
udaje$TotalCost <- udaje$Accommodation + udaje$Transport

# rýchla kontrola
summary(udaje)

    Duration           Age        Accommodation 
 Min.   : 5.000   Min.   :20.00   Min.   : 100  
 1st Qu.: 7.000   1st Qu.:28.00   1st Qu.: 600  
 Median : 7.000   Median :31.00   Median : 900  
 Mean   : 7.606   Mean   :33.18   Mean   :1245  
 3rd Qu.: 8.000   3rd Qu.:38.00   3rd Qu.:1200  
 Max.   :14.000   Max.   :60.00   Max.   :8000  
                                                
   Transport        TotalCost    
 Min.   :  20.0   Min.   :  200  
 1st Qu.: 200.0   1st Qu.: 1000  
 Median : 550.0   Median : 1400  
 Mean   : 645.2   Mean   : 1899  
 3rd Qu.: 800.0   3rd Qu.: 1900  
 Max.   :3000.0   Max.   :10500  
 NA's   :1        NA's   :1      

#Identifikácia číselných premenných

numeric_cols <- sapply(udaje, is.numeric)

# výpočet mediánov pre číselné stĺpce
column_medians <- sapply(udaje[, numeric_cols], median, na.rm = TRUE)

# imputácia chýbajúcich hodnôt
udaje_imputed <- udaje
for (col in names(column_medians)) {
udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}

udaje <- udaje_imputed
summary(udaje)

    Duration           Age        Accommodation 
 Min.   : 5.000   Min.   :20.00   Min.   : 100  
 1st Qu.: 7.000   1st Qu.:28.00   1st Qu.: 600  
 Median : 7.000   Median :31.00   Median : 900  
 Mean   : 7.606   Mean   :33.18   Mean   :1245  
 3rd Qu.: 8.000   3rd Qu.:38.00   3rd Qu.:1200  
 Max.   :14.000   Max.   :60.00   Max.   :8000  
   Transport        TotalCost    
 Min.   :  20.0   Min.   :  200  
 1st Qu.: 200.0   1st Qu.: 1000  
 Median : 550.0   Median : 1400  
 Mean   : 644.5   Mean   : 1895  
 3rd Qu.: 800.0   3rd Qu.: 1900  
 Max.   :3000.0   Max.   :10500  

#Základný lineárny model ##Ako základný model zvolíme: ##vysvetľovaná premenná: TotalCost (celkové náklady na zájazd), ##vysvetľujúce premenné: Duration (počet dní) a Age (vek cestovateľa).

# ZÁKLADNÁ REGRESIA
attach(udaje)
model <- lm(TotalCost ~ +1 + Duration + Age, data = udaje)
summary(model)

Call:
lm(formula = TotalCost ~ +1 + Duration + Age, data = udaje)

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

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 2401.686   1123.981   2.137   0.0344 *
Duration    -103.226     98.862  -1.044   0.2983  
Age            8.395     22.155   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

##Koeficient pri Duration ukazuje, o koľko sa v priemere zmenia celkové náklady na zájazd pri predĺžení pobytu o 1 deň (pri fixnom veku). Koeficient pri Age zachytáva, ako sa menia priemerné náklady s vekom cestovateľa (pri rovnakej dĺžke pobytu). R^2 a upravené R^2 hovoria, aká časť variability nákladov je vysvetlená týmito dvomi premennými.

#Test špecifikácie – Ramsey RESET test ##Najprv otestujeme, či lineárna špecifikácia modelu vyzerá byť vhodná alebo naznačuje, že chýbajú nejaké nelineárne členy (napr. štvorce, mocniny) alebo premenné. ##Formálne testujeme: ##H₀: model je správne špecifikovaný, ##H₁: model je nesprávne špecifikovaný.

resettest(model)

    RESET test

data:  model
RESET = 2.0638, df1 = 2, df2 = 132, p-value = 0.131

##Ak je p-hodnota < 0.05, zamietame H₀ → lineárna špecifikácia pravdepodobne nestačí,bolo by vhodné uvažovať o nelineárnej transformácii alebo doplnení premenných. ##Ak je p-hodnota ≥ 0.05, nemáme dôkaz o zlej funkčnej forme (aspoň podľa tohto testu).

#Grafická analýza špecifikácie ##Graf Residuals vs Fitted

plot(model, which = 1)

##Ak sú rezíduá rozptýlené náhodne okolo horizontálnej osi bez výrazného vzoru, lineárny model môže byť vhodný. ##Ak vidíš zakrivenie, lievikovitý tvar alebo systematický vzor, naznačuje to, že by mohla pomôcť transformácia (napr. log(TotalCost)) alebo pridanie nelineárnych členov (napr. Duration^2).

#C+R (Component + Residual) grafy

car::crPlots(model)

##Pre každú vysvetľujúcu premennú vidíš na osi X jej hodnoty a na osi Y kombináciu jej príspevku v modeli + rezíduá. ##Ak sa hladká krivka výrazne odchyľuje od priamky, pri danej premennej môže byť vhodné zvážiť napr. kvadratický člen (I(Duration^2)), logaritmus alebo inú transformáciu.

#Nelineárna špecifikácia – kvadratické členy ##Predpokladajme, že C+R graf naznačuje nelinearitu najmä pri Duration (napr. krátke vs. dlhé pobyty).

model_quad <- lm(TotalCost ~ +1 + Duration + I(Duration^2) + Age, data = udaje)
summary(model_quad)

Call:
lm(formula = TotalCost ~ +1 + Duration + I(Duration^2) + Age, 
    data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-1771.2  -883.7  -468.7    59.5  8608.2 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)    3161.74    2751.78   1.149    0.253
Duration       -292.97     634.42  -0.462    0.645
I(Duration^2)    11.32      37.38   0.303    0.763
Age               8.38      22.23   0.377    0.707

Residual standard error: 1839 on 133 degrees of freedom
Multiple R-squared:  0.01065,   Adjusted R-squared:  -0.01167 
F-statistic: 0.4771 on 3 and 133 DF,  p-value: 0.6987

anova(model, model_quad)

Analysis of Variance Table

Model 1: TotalCost ~ +1 + Duration + Age
Model 2: TotalCost ~ +1 + Duration + I(Duration^2) + Age
  Res.Df       RSS Df Sum of Sq      F Pr(>F)
1    134 450070864                           
2    133 449760797  1    310067 0.0917 0.7625

resettest(model_quad)

    RESET test

data:  model_quad
RESET = 2.6961, df1 = 2, df2 = 131, p-value = 0.07121

##Ak sa uprav. R² zvýši a zároveň Anova test ukáže, že nový model je štatisticky lepší (p-hodnota < 0.05), kvadratický člen Duration^2 pravdepodobne zlepšuje model. ##RESET test na kvadratickom modeli nám ukáže, či sme problém nesprávnej špecifikácie aspoň čiastočne odstránili.

#Rozšírený model (viac kvadratických členov)

model_rozsireny <- lm(TotalCost ~ Duration + I(Duration^2) + Age + I(Age^2),
data = udaje)
summary(model_rozsireny)

Call:
lm(formula = TotalCost ~ Duration + I(Duration^2) + Age + I(Age^2), 
    data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-1786.9  -894.7  -452.8   115.9  8627.0 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)   2242.4507  4269.9393   0.525    0.600
Duration      -273.1854   640.4787  -0.427    0.670
I(Duration^2)   10.2043    37.7185   0.271    0.787
Age             56.2824   171.1706   0.329    0.743
I(Age^2)        -0.6541     2.3173  -0.282    0.778

Residual standard error: 1845 on 132 degrees of freedom
Multiple R-squared:  0.01124,   Adjusted R-squared:  -0.01872 
F-statistic: 0.3753 on 4 and 132 DF,  p-value: 0.826

anova(model, model_rozsireny)

Analysis of Variance Table

Model 1: TotalCost ~ +1 + Duration + Age
Model 2: TotalCost ~ Duration + I(Duration^2) + Age + I(Age^2)
  Res.Df       RSS Df Sum of Sq      F Pr(>F)
1    134 450070864                           
2    132 449489502  2    581362 0.0854 0.9182

resettest(model_rozsireny)

    RESET test

data:  model_rozsireny
RESET = 2.1746, df1 = 2, df2 = 130, p-value = 0.1178

#Zlom pomocou dummy premennej – lineárna lomená funkcia ##Predpokladajme, že dlhšie pobyty (napr. od 9 dní vyššie) sa správajú inak ako kratšie pobyty. ##Zavedieme si dummy premennú DUM_long: ##DUM_long = 0, ak Duration < 9, ##DUM_long = 1, ak Duration ≥ 9.

udaje$DUM_long <- ifelse(udaje$Duration < 9, 0, 1)
table(udaje$DUM_long)

  0   1 
104  33 

#Zlom v autonómnom člene
modelD_auto <- lm(TotalCost ~ +1 + DUM_long + Duration + Age, data = udaje)
summary(modelD_auto)

Call:
lm(formula = TotalCost ~ +1 + DUM_long + Duration + Age, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-1743.1  -827.2  -471.0   107.0  8576.0 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 2052.678   1486.734   1.381    0.170
DUM_long    -223.505    620.489  -0.360    0.719
Duration     -54.615    167.483  -0.326    0.745
Age            9.393     22.399   0.419    0.676

Residual standard error: 1839 on 133 degrees of freedom
Multiple R-squared:  0.01093,   Adjusted R-squared:  -0.01138 
F-statistic: 0.4899 on 3 and 133 DF,  p-value: 0.6899

##Koeficient pri DUM_long ukazuje, o koľko sú v priemere celkové náklady iné pri dlhších pobytoch (Duration ≥ 9) oproti kratším (Duration < 9), pri rovnakej dĺžke a veku.

#Zlom v sklone – interakcia DUM_long * Duration

modelD_sklon <- lm(TotalCost ~ +1 + Duration + I(DUM_long * Duration) + Age,
data = udaje)
summary(modelD_sklon)

Call:
lm(formula = TotalCost ~ +1 + Duration + I(DUM_long * Duration) + 
    Age, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-1731.3  -839.6  -462.0    99.8  8581.7 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)
(Intercept)            2131.267   1593.423   1.338    0.183
Duration                -65.350    186.261  -0.351    0.726
I(DUM_long * Duration)  -16.636     69.238  -0.240    0.810
Age                       9.056     22.403   0.404    0.687

Residual standard error: 1839 on 133 degrees of freedom
Multiple R-squared:  0.01039,   Adjusted R-squared:  -0.01193 
F-statistic: 0.4657 on 3 and 133 DF,  p-value: 0.7067

#Porovnanie modelov
anova(model, modelD_sklon)

Analysis of Variance Table

Model 1: TotalCost ~ +1 + Duration + Age
Model 2: TotalCost ~ +1 + Duration + I(DUM_long * Duration) + Age
  Res.Df       RSS Df Sum of Sq      F Pr(>F)
1    134 450070864                           
2    133 449875593  1    195271 0.0577 0.8105

resettest(modelD_sklon)

    RESET test

data:  modelD_sklon
RESET = 2.0851, df1 = 2, df2 = 131, p-value = 0.1284

##Sklon pri Duration pre krátke pobyty (DUM_long = 0) je daný koeficientom pri Duration. ##Pre dlhé pobyty (DUM_long = 1) je sklon = beta_Duration + beta_interakcia.

#Box-Coxova transformácia závislej premennej

boxcox(model)

##Z grafu Box-Coxovej transformácie zistíme optimálnu hodnotu λ (lambda), pri ktorej je log-vierohodnosť najvyššia. ##Približne: ##λ ≈ 1 → netreba transformovať, ##λ ≈ 0 → logaritmus: log(TotalCost), ##λ ≈ 0.5 → odmocnina: sqrt(TotalCost), ##λ ≈ -1 → recipročná transformácia: 1/TotalCost.

lambda_opt <- 0.5   
model_lambda <- lm(I((TotalCost^lambda_opt - 1) / lambda_opt) ~ Duration + Age,
data = udaje)
summary(model_lambda)

Call:
lm(formula = I((TotalCost^lambda_opt - 1)/lambda_opt) ~ Duration + 
    Age, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-50.878 -16.598  -6.179   7.122 125.395 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  76.6500    19.5941   3.912 0.000145 ***
Duration     -1.3080     1.7234  -0.759 0.449205    
Age           0.3722     0.3862   0.964 0.336870    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 31.95 on 134 degrees of freedom
Multiple R-squared:  0.01257,   Adjusted R-squared:  -0.002166 
F-statistic: 0.853 on 2 and 134 DF,  p-value: 0.4284

resettest(model_lambda)

    RESET test

data:  model_lambda
RESET = 0.59218, df1 = 2, df2 = 132, p-value = 0.5546

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

plot(cars)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.