Špecifikácia modelu – hodnotenie filmov

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

V tejto úlohe budeme pomocou regresie skúmať, nakoľko je hodnotenie filmu podľa IMDb ovplyvnené:

• dĺžkou filmu v minútach,
  
• rokom uvedenia (staršie vs. novšie filmy),
  
• rozpočtom filmu.

Na to použijeme databázu najlepšie hodnotených filmov z portálu IMDb, ktorá obsahuje napríklad umiestnenie v rebríčku, názov filmu, rok, žáner, hodnotenie IMDb, dĺžku, režiséra, štúdio a rozpočet.

#######################################################################
# PRÍPRAVA ÚDAJOV
#######################################################################

# načítanie dát (oddelené bodkočiarkou, desatinná čiarka)
udaje_raw <- read.csv("data_ekonometria.csv",
                      sep = ";",
                      dec = ",",
                      header = TRUE,
                      check.names = FALSE)

# pozrime si názvy stĺpcov
names(udaje_raw)

[1] "Umiestnenie"     "Názov"           "Rok"             "Žáner"          
[5] "IMDb hodnotenie" "Dĺžka (min)"     "Réžia"           "Spoločnosť"     
[9] "Rozpočet [$]"   

# vytvoríme si jednoduchšie (bez diakritiky) a numerické premenné
udaje <- udaje_raw

# IMDb hodnotenie
udaje$imdb <- udaje[["IMDb hodnotenie"]]

# dĺžka filmu v minútach
udaje$dlzka <- udaje[["Dĺžka (min)"]]

# rok uvedenia
udaje$rok <- udaje[["Rok"]]

# rozpočet – je v texte (medzery a znak $), preto ho musíme očistiť a prekonvertovať na číslo
udaje$rozpocet <- as.numeric(gsub(" |\\$|,", "", udaje[["Rozpočet [$]"]]))

# vyberieme len premenné, ktoré budeme používať v modeli
udaje <- udaje[, c("imdb", "dlzka", "rok", "rozpocet")]

summary(udaje)

      imdb           dlzka            rok          rozpocet        
 Min.   :8.000   Min.   : 45.0   Min.   :1921   Min.   :    12000  
 1st Qu.:8.100   1st Qu.:110.0   1st Qu.:1966   1st Qu.:  2500000  
 Median :8.200   Median :128.0   Median :1994   Median : 13000000  
 Mean   :8.294   Mean   :130.4   Mean   :1986   Mean   : 35923313  
 3rd Qu.:8.400   3rd Qu.:148.0   3rd Qu.:2006   3rd Qu.: 35000000  
 Max.   :9.300   Max.   :238.0   Max.   :2022   Max.   :400000000  

Teraz máme pripravenú menšiu databázu, kde každý riadok predstavuje film a premenné:

• imdb – hodnotenie filmu na IMDb (závislá premenná),
  
• dlzka – dĺžka filmu v minútach,
  
• rok – rok uvedenia filmu,
  
• rozpocet – rozpočet filmu v USD.

Môže sa stať, že niektoré filmy nebudú mať uvedený rozpočet; tieto chýbajúce hodnoty nahradíme mediánmi.

# IMPUTÁCIA CHÝBAJÚCICH HODNÔT

# vypočítame mediány pre jednotlivé stĺpce (ignorujeme NA)
column_medians <- sapply(udaje, median, na.rm = TRUE)

# nahradíme chýbajúce hodnoty mediánom danej premennej
udaje_imputed <- udaje
for (col in names(udaje_imputed)) {
  udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}

udaje <- udaje_imputed
summary(udaje)

      imdb           dlzka            rok          rozpocet        
 Min.   :8.000   Min.   : 45.0   Min.   :1921   Min.   :    12000  
 1st Qu.:8.100   1st Qu.:110.0   1st Qu.:1966   1st Qu.:  2500000  
 Median :8.200   Median :128.0   Median :1994   Median : 13000000  
 Mean   :8.294   Mean   :130.4   Mean   :1986   Mean   : 35923313  
 3rd Qu.:8.400   3rd Qu.:148.0   3rd Qu.:2006   3rd Qu.: 35000000  
 Max.   :9.300   Max.   :238.0   Max.   :2022   Max.   :400000000  

1. Základná regresia

Základný model bude mať tvar:

• vysvetľovaná premenná: IMDb hodnotenie filmu (imdb),
• vysvetľujúce premenné: dĺžka filmu (dlzka), rok uvedenia (rok), rozpočet (rozpocet).

Formálne:

\[ \text{imdb}_i = \beta_0 + \beta_1 \text{dlzka}_i + \beta_2 \text{rok}_i + \beta_3 \text{rozpocet}_i + u_i \]

################################################################################
# ZÁKLADNÁ REGRESIA
################################################################################
attach(udaje)

model <- lm(imdb ~ 1 + dlzka + rok + rozpocet, data = udaje)
summary(model)

Call:
lm(formula = imdb ~ 1 + dlzka + rok + rozpocet, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42001 -0.15545 -0.05865  0.13091  0.98723 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.678e+00  1.238e+00   7.010 2.28e-11 ***
dlzka        2.126e-03  4.689e-04   4.535 9.00e-06 ***
rok         -3.391e-04  6.267e-04  -0.541    0.589    
rozpocet     3.424e-10  2.624e-10   1.305    0.193    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2218 on 246 degrees of freedom
Multiple R-squared:  0.08716,   Adjusted R-squared:  0.07603 
F-statistic: 7.829 on 3 and 246 DF,  p-value: 5.186e-05

Interpretácia tohto modelu: snažíme sa zistiť, či sú filmy s vyšším rozpočtom, dlhšou stopážou a novším rokom uvedenia na IMDb hodnotené systematicky inak (lepšie alebo horšie).

2. Test správnej funkčnej formy – Ramsey RESET test

Chceme otestovať, či lineárna špecifikácia je pre naše dáta (hodnotenie filmov) vhodná, alebo či by sme mali napríklad zaviesť kvadratické členy, logaritmy a podobne.

Myšlienka: nech pôvodný model má tvar

\[ \text{imdb}_i = \beta_0 + \beta_1 \text{dlzka}_i + \beta_2 \text{rok}_i + \beta_3 \text{rozpočet}_i + u_i \]

Ak je model správne špecifikovaný, potom pridaním mocnín vyrovnaných hodnôt (\(\hat y^2, \hat y^3\)) sa model podstatne nezlepší.

Budeme testovať:

• \(H_0:\) model je správne špecifikovaný,
  
• \(H_1:\) model je špecifikovaný nesprávne (chýbajú dôležité nelineárne alebo iné premenné).

{r # RESET test resettest(model)

Interpretácia:

• Ak p-hodnota < 0.05 → odmietame \(H_0\) → model filmového hodnotenia je pravdepodobne zle špecifikovaný (chýbajú nelineárne členy alebo ďalšie premenné).
• Ak p-hodnota ≥ 0.05 → nemáme dôkaz o nesprávnej špecifikácii.

3. Grafická analýza

3.1 Graf Residuals vs. Fitted

Pozrieme sa na vzťah medzi vyrovnanými hodnotami IMDb hodnotenia a rezíduami.

plot(model, which = 1)

Ak rezíduá vykazujú systematický (nezhlukovaný) vzor – napríklad zakrivenie –, môže to znamenať, že model závislosti IMDb hodnotenia na dĺžke, roku a rozpočte nie je úplne lineárny. Vtedy môže pomôcť transformácia niektorej premennej (napr. zavedenie kvadratického člena) alebo doplnenie novej premennej (napr. dummy pre „staršie“/„novšie“ filmy).

3.2 C+R grafy (Component + Residual plots)

C+R grafy nám pomôžu zisťovať, či je vzťah medzi IMDb hodnotením a jednotlivými vysvetľujúcimi premennými lineárny, alebo zakrivený.

car::crPlots(model)

Ak sa krivka výrazne odchyľuje od priamky, zvážime transformáciu danej premennej (napr. kvadratický člen, logaritmus, odmocnina). Predstavme si, že najväčší odklon od linearity pozorujeme najmä pri premennej dĺžka filmu (dlzka) a možno aj pri rozpočte (rozpocet).

4. Nelineárna špecifikácia – kvadratické členy

Zložitejšie nelineárne vzťahy môžeme približovať polynómom, napr. zavedením kvadrátu vybraných vysvetľujúcich premenných.

Predpokladajme, že sa rozhodneme zaviesť kvadrát premennej dĺžka a rozpočet, lebo C+R grafy naznačujú zakrivenie ich vzťahu k IMDb hodnoteniu.

model <- lm(imdb ~ 1 + dlzka + rok + rozpocet, data = udaje)

model_kvadr <- lm(imdb ~ 1 + dlzka + rok + rozpocet +
                    I(dlzka^2) + I(rozpocet^2),
                  data = udaje)

summary(model_kvadr)

Call:
lm(formula = imdb ~ 1 + dlzka + rok + rozpocet + I(dlzka^2) + 
    I(rozpocet^2), data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.41110 -0.15480 -0.05242  0.11408  0.98026 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    9.245e+00  1.308e+00   7.066 1.66e-11 ***
dlzka          4.049e-03  3.010e-03   1.345   0.1799    
rok           -6.942e-04  6.844e-04  -1.014   0.3114    
rozpocet       9.925e-10  5.942e-10   1.670   0.0961 .  
I(dlzka^2)    -6.907e-06  1.062e-05  -0.650   0.5160    
I(rozpocet^2) -2.352e-18  1.967e-18  -1.196   0.2330    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2218 on 244 degrees of freedom
Multiple R-squared:  0.09387,   Adjusted R-squared:  0.07531 
F-statistic: 5.056 on 5 and 244 DF,  p-value: 0.0001988

# porovnanie lineárneho a kvadratického modelu
anova(model, model_kvadr)

Analysis of Variance Table

Model 1: imdb ~ 1 + dlzka + rok + rozpocet
Model 2: imdb ~ 1 + dlzka + rok + rozpocet + I(dlzka^2) + I(rozpocet^2)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    246 12.097                           
2    244 12.008  2  0.088998 0.9042 0.4062

# RESET test pre kvadratický model
resettest(model_kvadr)

    RESET test

data:  model_kvadr
RESET = 2.5842, df1 = 2, df2 = 242, p-value = 0.07754

Ak má transformovaný model vyšší upravený koeficient determinácie \(R^2_{adj}\) a zároveň sa zlepší aj výsledok RESET testu, môžeme konštatovať, že kvadratické členy pomohli lepšie vystihnúť vzťah medzi hodnotením filmu a jeho dĺžkou/rozpočtom.

Ak sa ukáže, že niektoré kvadratické členy sú štatisticky nevýznamné, môžeme ich postupne vylučovať a zvoliť „úspornejší“ model (bez zbytočných členov).

5. Rozšírený model s nelineárnymi členmi

Skúsme model, kde pre všetky tri vysvetľujúce premenné (dĺžka, rok, rozpočet) zavedie­me aj kvadratické členy a následne vyhodnotíme ich štatistickú významnosť:

model_rozsireny <- lm(imdb ~ dlzka + rok + rozpocet +
                        I(dlzka^2) + I(rok^2) + I(rozpocet^2),
                      data = udaje)
summary(model_rozsireny)

Call:
lm(formula = imdb ~ dlzka + rok + rozpocet + I(dlzka^2) + I(rok^2) + 
    I(rozpocet^2), data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.43993 -0.14709 -0.05205  0.12665  0.96420 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
(Intercept)   -2.064e+02  9.073e+01  -2.275   0.0238 *
dlzka          2.905e-03  3.020e-03   0.962   0.3370  
rok            2.176e-01  9.183e-02   2.369   0.0186 *
rozpocet       1.108e-09  5.906e-10   1.876   0.0618 .
I(dlzka^2)    -3.559e-06  1.061e-05  -0.335   0.7376  
I(rok^2)      -5.520e-05  2.322e-05  -2.377   0.0182 *
I(rozpocet^2) -2.186e-18  1.950e-18  -1.121   0.2632  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2198 on 243 degrees of freedom
Multiple R-squared:  0.1145,    Adjusted R-squared:  0.0926 
F-statistic: 5.235 on 6 and 243 DF,  p-value: 4.33e-05

anova(model, model_rozsireny)

Analysis of Variance Table

Model 1: imdb ~ 1 + dlzka + rok + rozpocet
Model 2: imdb ~ dlzka + rok + rozpocet + I(dlzka^2) + I(rok^2) + I(rozpocet^2)
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    246 12.097                              
2    243 11.735  3   0.36183 2.4974 0.06034 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

resettest(model_rozsireny)

    RESET test

data:  model_rozsireny
RESET = 4.0892, df1 = 2, df2 = 241, p-value = 0.01793

Interpretácia:

• štatisticky významne kvadratické koeficienty naznačujú zakrivený vzťah (napr. veľmi dlhé filmy už nemusia byť hodnotené lepšie),
• porovnaním modelov pomocou ANOVA sledujeme, či rozšírený model prináša štatisticky významné zlepšenie oproti pôvodnému lineárnemu modelu.

6. Dummy premenná a lineárna lomená funkcia

Predpokladajme teraz, že C+R graf pre premennú dĺžka filmu (dlzka) naznačuje, že približne pri 150 minútach nastáva zmena trendu – kratšie filmy sa správajú inak ako veľmi dlhé filmy.

Zavedieme preto dummy premennú:

\[ DUM_i = \begin{cases} 0, & \text{ak } \text{dlzka}_i < 150 \\ 1, & \text{ak } \text{dlzka}_i \ge 150 \end{cases} \]

Túto dummy možno využiť:

1. na zlom v autonómnom (konštantnom) člene,
   
2. na zlom v sklone voči premenným (napr. voči dĺžke).

6.1 Zlom v autonómnom člene

udaje$DUM <- ifelse(udaje$dlzka < 150, 0, 1)

modelD_auto <- lm(imdb ~ 1 + DUM + dlzka + rok + rozpocet,
                  data = udaje)
summary(modelD_auto)

Call:
lm(formula = imdb ~ 1 + DUM + dlzka + rok + rozpocet, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42320 -0.15874 -0.06193  0.12451  0.99780 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.670e+00  1.240e+00   6.993 2.54e-11 ***
DUM          3.008e-02  5.289e-02   0.569   0.5701    
dlzka        1.799e-03  7.436e-04   2.419   0.0163 *  
rok         -3.168e-04  6.288e-04  -0.504   0.6148    
rozpocet     3.382e-10  2.629e-10   1.287   0.1995    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2221 on 245 degrees of freedom
Multiple R-squared:  0.08836,   Adjusted R-squared:  0.07348 
F-statistic: 5.937 on 4 and 245 DF,  p-value: 0.0001416

Tu testujeme, či sa pri prechode z kratších na dlhé filmy (nad 150 minút) posunie priemerné IMDb hodnotenie o nejakú konštantnú hodnotu (bez zmeny sklonu).

6.2 Zlom v sklone (lineárna lomená funkcia)

Teraz vytvoríme model so zlomeným sklonom voči dĺžke – čiže iný vplyv dĺžky na hodnotenie pri kratších a pri dlhších filmoch:

modelD_sklon <- lm(imdb ~ 1 + dlzka + I(DUM * dlzka) + rok + rozpocet,
                   data = udaje)
summary(modelD_sklon)

Call:
lm(formula = imdb ~ 1 + dlzka + I(DUM * dlzka) + rok + rozpocet, 
    data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42400 -0.15910 -0.06235  0.12705  0.99543 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     8.654e+00  1.242e+00   6.970 2.92e-11 ***
dlzka           1.861e-03  8.020e-04   2.320   0.0212 *  
I(DUM * dlzka)  1.337e-04  3.269e-04   0.409   0.6829    
rok            -3.118e-04  6.313e-04  -0.494   0.6218    
rozpocet        3.388e-10  2.630e-10   1.288   0.1989    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2221 on 245 degrees of freedom
Multiple R-squared:  0.08778,   Adjusted R-squared:  0.07289 
F-statistic: 5.894 on 4 and 245 DF,  p-value: 0.0001522

Interpretácia:

• koeficient pri dlzka vyjadruje vplyv dĺžky na IMDb hodnotenie pre kratšie filmy (kde DUM = 0),
  
• súčet koeficientov dlzka + I(DUM * dlzka) vyjadruje vplyv dĺžky na hodnotenie pre dlhšie filmy (DUM = 1).

Porovnanie pôvodného lineárneho modelu a modelu so zlomeným sklonom urobíme pomocou ANOVA a RESET testu:

anova(model, modelD_sklon)

Analysis of Variance Table

Model 1: imdb ~ 1 + dlzka + rok + rozpocet
Model 2: imdb ~ 1 + dlzka + I(DUM * dlzka) + rok + rozpocet
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    246 12.097                           
2    245 12.089  1 0.0082517 0.1672 0.6829

resettest(modelD_sklon)

    RESET test

data:  modelD_sklon
RESET = 0.60852, df1 = 2, df2 = 243, p-value = 0.545

Ak ANOVA test ukáže, že model so zlomom v sklone je štatisticky lepší, môžeme konštatovať, že dĺžka filmu ovplyvňuje hodnotenie IMDb inak pri kratších filmoch a inak pri veľmi dlhých.

7. Box–Coxov transformačný test (voliteľné)

Box–Coxova transformácia sa používa na systematické hľadanie vhodnej transformácie závislej premenné, v našom prípade IMDb hodnotenia filmu.

\[ Y_i^{(\lambda)} = \begin{cases} \dfrac{Y_i^{\lambda} - 1}{\lambda}, & \text{ak } \lambda \neq 0, \\[10pt] \ln(Y_i), & \text{ak } \lambda = 0. \end{cases} \]

library(MASS)
boxcox(model)

Z grafu Box–Cox vyhodnotíme, ktorá hodnota \(\lambda\) by mohla byť vhodná. Podľa toho:

• ak \(\lambda \approx 1\) → netreba transformovať,
  
• ak \(\lambda \approx 0\) → vhodná logaritmická transformácia \(\log(\text{imdb})\),
  
• ak \(\lambda \approx 0{,}5\) → vhodná odmocninová transformácia,
  
• ak \(\lambda \approx -1\) → transformácia \(1/\text{imdb}\).

Pre ilustráciu si vyskúšajme nejakú vybranú hodnotu \(\lambda\) (napríklad \(\lambda = 1.8\)) a znova odhadnime model:

lambda <- 1.8

model_lambda <- lm(I((imdb^lambda - 1) / lambda) ~ 1 + dlzka + rok + rozpocet,
                   data = udaje)
summary(model_lambda)

Call:
lm(formula = I((imdb^lambda - 1)/lambda) ~ 1 + dlzka + rok + 
    rozpocet, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2858 -0.8483 -0.3328  0.7056  5.6102 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.658e+01  6.836e+00   3.889  0.00013 ***
dlzka        1.177e-02  2.589e-03   4.544 8.66e-06 ***
rok         -1.859e-03  3.461e-03  -0.537  0.59160    
rozpocet     1.869e-09  1.449e-09   1.290  0.19843    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.225 on 246 degrees of freedom
Multiple R-squared:  0.08726,   Adjusted R-squared:  0.07613 
F-statistic: 7.839 on 3 and 246 DF,  p-value: 5.118e-05

resettest(model_lambda)

    RESET test

data:  model_lambda
RESET = 0.7661, df1 = 2, df2 = 244, p-value = 0.4659

# anova porovnanie s pôvodným modelom nemá zmysel, pretože sme zmenili vysvetľovanú premennú

Ak by transformácia výraznejšie zlepšila vlastnosti modelu (vyšší \(R^2_{adj}\), lepší RESET test), bolo by možné ju ďalej používať. Zároveň však strácame priamočiaru interpretáciu „o koľko bodov IMDb sa zmení hodnotenie“, preto je pri praktickej interpretácii často výhodnejšie zostať pri pôvodnej mierke IMDb a nelinearity riešiť napríklad kvadratickými členmi alebo lomenými funkciami.