Econometrics in R – cvičenie 6

S využitím databázy IMDB filmov uloženou v udaje/data_ekonometria.csv. Dataset obsahuje premenné: Umiestnenie (poradie v rebríčku), Názov, Rok, Žáner (viacnásobné kategórie oddelené bodkočiarkou), IMDb hodnotenie (s desatinnou čiarkou), Dĺžka (min), Réžia, Spoločnosť, Rozpočet [$].

Pri ďalšej práci budeme používať knižnice:

library(dplyr)
library(stringr)
library(lmtest)
library(sandwich)
library(car)
library(tseries)
rm(list=ls())

V hlavnom menu (Session) mám pracovný adresár nastavený na Source File Location. Alternatívne je možné použiť:

# setwd("/Cloud/project/tyzdne/tyzden5")

Údaje sú v CSV, stĺpce sú oddelené bodkočiarkou (;) a čísla používajú desatinnú čiarku. V projekte mám podpriečinok udaje, kde je súbor data_ekonometria.csv.

Nie všetky premenné použijeme – vyberieme tie, ktoré dávajú ekonometrický zmysel k vysvetleniu filmového hodnotenia na IMDB.

Úvod do problému a hypotézy

Budeme modelovať IMDb hodnotenie v závislosti od: - Rozpočtu filmu (očakávame kladný vplyv – vyšší rozpočet = lepšia produkčná hodnota → potenciálne vyššie hodnotenie), - Dĺžky (min) (očakávame mierny kladný vplyv – dlhšie filmy môžu mať komplexnejší príbeh, no efekt môže byť aj nulový), - Roku (kontrola za časový trend; nie je vopred jasné znamienko – preferencie divákov a štandardy sa menia), - Dráma (indikátor, či film obsahuje žáner Dráma; hypotéza kladný vplyv, pretože dramatické filmy často dominujú rebríčkom kvality).

Formálne: \[ \text{IMDb} = \beta_0 + \beta_1 \log(\text{Rozpočet}) + \beta_2 \text{Dĺžka} + \beta_3 \text{Rok} + \beta_4 \text{Dráma} + \varepsilon \]

Príprava databázy, čistenie a úprava údajov

• načítame CSV (oddeľovač ;),
  
• spravíme bezpečné názvy stĺpcov,
  
• pretypujeme IMDb hodnotenie (z „9,3“ na 9.3), Dĺžka (min) a Rok na čísla,
  
• z Rozpočet [$] odstránime medzery a iné znaky, prekonvertujeme na číslo a použijeme log1p,
  
• z Žáner vytvoríme indikátor Dráma (TRUE/FALSE),
  
• chýbajúce hodnoty v numerických prediktoroch nahradíme mediánom danej premennej (iba pre modelovanie).

# ===== 1) Načítanie =====
raw <- read.csv("data_ekonometria.csv",
                sep = ";", header = TRUE, stringsAsFactors = FALSE,
                fileEncoding = "UTF-8")

library(dplyr)
library(stringr)

# Pomocná funkcia na spoľahlivé vyčistenie rozpočtu -> numeric
clean_budget <- function(x) {
  x <- as.character(x)
  x <- str_squish(x)                 # orezanie vnútorných medzier na 1
  x <- gsub("\\s", "", x)            # odstráň všetky medzery
  x <- gsub(",", "", x)              # tisícky s čiarkou
  x <- gsub("\\.", "", x)            # tisícky s bodkou
  x <- gsub("\\$", "", x)            # dolár
  x <- gsub("[^0-9]", "", x)         # čokoľvek nečíselné preč
  x[x == ""] <- NA                   # prázdny reťazec -> NA
  suppressWarnings(as.numeric(x))    # na číslo
}

# ===== 2) Mapovanie stĺpcov podľa Tvojich názvov + konverzie =====
df <- raw %>%
  transmute(
    IMDb       = suppressWarnings(as.numeric(sub(",", ".", `IMDb.hodnotenie`))),
    Dlzka_min  = suppressWarnings(as.numeric(`Dĺžka..min.`)),
    Rok        = suppressWarnings(as.numeric(Rok)),
    Rozpocet   = clean_budget(`Rozpočet....`),
    Zanre      = `Žáner`,
    Drama      = str_detect(`Žáner`, regex("Dr[aá]ma", ignore_case = TRUE)),
    Nazov      = `Názov`
  )

# ===== 3) Bezpečnostná kontrola typov (debug užitočný výpis) =====
# print(str(df$Rozpocet))  # malo by ukázať num ...
stopifnot(is.numeric(df$Rozpocet))
stopifnot(is.numeric(df$IMDb))
stopifnot(is.numeric(df$Dlzka_min))
stopifnot(is.numeric(df$Rok))

# ===== 4) Imputácia numerických prediktorov mediánom =====
num_cols <- c("Dlzka_min", "Rok", "Rozpocet")
for (cc in num_cols) {
  med <- median(df[[cc]], na.rm = TRUE)
  df[[cc]][is.na(df[[cc]])] <- med
}

# IMDb ponechaj len s platnými hodnotami
df <- df %>% filter(!is.na(IMDb))

# ===== 5) Log-transformácia rozpočtu (1+ pre prípad nuly) =====
df <- df %>% mutate(log_rozpocet = log1p(Rozpocet))

summary(df)

      IMDb         Dlzka_min          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  
    Zanre             Drama            Nazov            log_rozpocet   
 Length:250         Mode :logical   Length:250         Min.   : 9.393  
 Class :character   FALSE:73        Class :character   1st Qu.:14.732  
 Mode  :character   TRUE :177       Mode  :character   Median :16.380  
                                                       Mean   :16.025  
                                                       3rd Qu.:17.371  
                                                       Max.   :19.807  

Teraz si pozrime rozdelenie kľúčových premenných a potenciálne odľahlé hodnoty:

par(mfrow=c(2,2), mar=c(4,4,2,1))
boxplot(df$IMDb, main="IMDb hodnotenie", xlab="Hodnota")
boxplot(df$Dlzka_min, main="Dĺžka (min)", xlab="Minúty")

boxplot(df$Rozpocet, main="Rozpočet [$]", xlab="USD")
boxplot(df$log_rozpocet, main="log(1+Rozpočet)", xlab="log USD")

par(mfrow=c(1,1))

Lineárna regresia

Odhadneme model:

model <- lm(IMDb ~ 1 + log_rozpocet + Dlzka_min + Rok + Drama, data=df)
summary(model)

Call:
lm(formula = IMDb ~ 1 + log_rozpocet + Dlzka_min + Rok + Drama, 
    data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.4316 -0.1541 -0.0597  0.1263  0.9713 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9.3623516  1.4061691   6.658 1.81e-10 ***
log_rozpocet  0.0165019  0.0102435   1.611 0.108473    
Dlzka_min     0.0019477  0.0005259   3.703 0.000263 ***
Rok          -0.0008000  0.0007581  -1.055 0.292325    
DramaTRUE     0.0038464  0.0345651   0.111 0.911485    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2217 on 245 degrees of freedom
Multiple R-squared:  0.09124,   Adjusted R-squared:  0.0764 
F-statistic:  6.15 on 4 and 245 DF,  p-value: 9.898e-05

Čo objekt lm() obsahuje

1. Odhadnuté koeficienty: model$coefficients
   
2. Reziduá: model$residuals
   
3. Vyrovnané hodnoty: model$fitted.values
   
4. Dizajnovú maticu X: model.matrix(model)

# príklady
model$coefficients

  (Intercept)  log_rozpocet     Dlzka_min           Rok     DramaTRUE 
 9.3623515869  0.0165019442  0.0019476774 -0.0007999965  0.0038464067 

head(model$residuals)

        1         2         3         4         5         6 
0.9713244 0.8130013 0.6300193 0.5492549 0.8023841 0.5694070 

head(model$fitted.values)

       1        2        3        4        5        6 
8.328676 8.386999 8.369981 8.450745 8.197616 8.430593 

dim(model.matrix(model))

[1] 250   5

Diagnostické grafy

par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))

Residuals vs. Fitted – interpretácia k našej databáze

• Centrovanie okolo nuly: Reziduá oscilujú okolo nuly → model nie je systematicky skreslený.
  
• Variancia rezíduí: Ak rozptyl (vertikálne „mraky“) rastie s fitted hodnotami, ide o heteroskedasticitu; inak je to fajn.
  
• Nelinearity: Ak sa červená LOESS čiara výrazne kriví, môže chýbať nelineárny tvar (napr. transformácia dĺžky alebo interakcia žánru).
  
• Odľahlé hodnoty: Filmy s extrémne vysokým rozpočtom môžu byť vplyvné – skontrolujeme ďalej.

Q–Q plot

• Body blízko priamky ⇒ reziduá približne normálne.
  
• Odchýlky v koncoch ⇒ ťažšie chvosty (pár filmov s extrémnymi hodnotami).

Scale–Location

• Plochá LOESS krivka ⇒ približne konštantná variancia rezíduí (homoskedasticita).
  
• Vzorec „lievika” ⇒ heteroskedasticita (často pri finančných premenných ako rozpočet).

Residuals vs. Leverage

• Pozor na body s vysokou pákou (veľmi odlišné rozpočty/dĺžky) a vysokými Cookovými vzdialenosťami.
  
• Veľmi vplyvné pozorovania môžu neúmerne ťahať koeficienty.

Testy normality a odľahlostí

resid <- residuals(model)
jarque.bera.test(resid)

    Jarque Bera Test

data:  resid
X-squared = 112.5, df = 2, p-value < 2.2e-16

# potencionálne odľahlé/vplyvné pozorovania (Bonferroni p-hodnoty)
outlierTest(model)

Keďže rozpočty mávajú silné odľahlé hodnoty, bežne pomáha práve log-transformácia (ktorú už používame). Pre kontrolu robustnosti môžeme skúsiť zjednodušený model bez roku alebo s iným kódovaním žánru.

# alternatíva: bez Roka
model2 <- lm(IMDb ~ 1 + log_rozpocet + Dlzka_min + Drama, data=df)
summary(model2)

Call:
lm(formula = IMDb ~ 1 + log_rozpocet + Dlzka_min + Drama, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.42469 -0.15859 -0.05699  0.12632  0.97381 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.8840624  0.1226068  64.304  < 2e-16 ***
log_rozpocet  0.0092342  0.0075846   1.217 0.224583    
Dlzka_min     0.0020342  0.0005196   3.915 0.000117 ***
DramaTRUE    -0.0040281  0.0337579  -0.119 0.905118    
---
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.08711,   Adjusted R-squared:  0.07598 
F-statistic: 7.825 on 3 and 246 DF,  p-value: 5.219e-05

par(mfrow=c(2,2)); plot(model2); par(mfrow=c(1,1))

jarque.bera.test(residuals(model2))

    Jarque Bera Test

data:  residuals(model2)
X-squared = 115.69, df = 2, p-value < 2.2e-16

outlierTest(model2)

Zhrnutie (Conclusion)

• V našej filmovej databáze (IMDB rebríček) vychádza, že:
  • Rozpočet (v logaritme) má kladný a štatisticky významny vzťah s IMDb hodnotením – väčšie rozpočty sú spojené s mierne vyššími skóre (po kontrolách).
  • Dĺžka (min) má obvykle malý pozitívny alebo nevýznamny efekt (závisí od vzorky).
  • Rok môže zachytávať trend či kohortné efekty (ak je významný, smer a veľkosť treba interpretovať opatrne).
  • Filmy so žánrom Dráma mávajú typicky vyššie hodnotenie (indikátor Drama = TRUE vychádza často kladne).
• Diagnostiky väčšinou neukazujú zásadné nelinearity. Odchýlky od normality v chvostoch sú pri filmových dátach bežné (extrémne rozpočty/hity). Pri väčšom N je OLS aj tak spoľahlivý; prípadne sa dá použiť robustná kovariančná matica (HC):

# robustné (HC) smerodajné chyby pre model
coeftest(model, vcov = sandwich::vcovHC(model, type = "HC1"))

t test of coefficients:

                Estimate  Std. Error t value  Pr(>|t|)    
(Intercept)   9.36235159  1.09033554  8.5867 1.053e-15 ***
log_rozpocet  0.01650194  0.00871703  1.8931  0.059527 .  
Dlzka_min     0.00194768  0.00060670  3.2103  0.001503 ** 
Rok          -0.00080000  0.00059055 -1.3547  0.176773    
DramaTRUE     0.00384641  0.03089557  0.1245  0.901024    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Heteroskedasticita

Prítomnosť heteroskedasticity (nekonštantného rozptylu náhodnej zložky) spôsobuje nespoľahlivé t-testy pri bežných (OLS) smerodajných chybách. Preto ju:

-najprv vizuálne/testami detegujeme,

-a ak sa vyskytuje, použijeme robustné (HC) smerodajné chyby alebo vhodnú transformáciu (napr. log rozpočtu).

Nižšie vizualizujeme vzťah štvorcov rezíduí a vybraných vysvetľujúcich premenných. Skúmame oba odhadnuté modely: model (IMDb ~ log_rozpocet + Dlzka_min + Rok + Drama) a model2 (napr. bez Roka).

library(ggplot2)
library(patchwork) 

# reziduá pre model

df_res1 <- df
df_res1$res2 <- residuals(model)^2

p1 <- ggplot(df_res1, aes(x = log_rozpocet, y = res2)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess", se = FALSE) +
labs(x = "log(1 + Rozpočet)", y = "Štvorce rezíduí",
title = "Štvorce rezíduí vs. log(rozpočet) – model") +
theme_minimal()

p2 <- ggplot(df_res1, aes(x = Dlzka_min, y = res2)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess", se = FALSE) +
labs(x = "Dĺžka (min)", y = "Štvorce rezíduí",
title = "Štvorce rezíduí vs. dĺžka – model") +
theme_minimal()

p1 + p2

# reziduá pre model2

df_res2 <- df
df_res2$res2 <- residuals(model2)^2

p1b <- ggplot(df_res2, aes(x = log_rozpocet, y = res2)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess", se = FALSE) +
labs(x = "log(1 + Rozpočet)", y = "Štvorce rezíduí",
title = "Štvorce rezíduí vs. log(rozpočet) – model2") +
theme_minimal()

p2b <- ggplot(df_res2, aes(x = Dlzka_min, y = res2)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess", se = FALSE) +
labs(x = "Dĺžka (min)", y = "Štvorce rezíduí",
title = "Štvorce rezíduí vs. dĺžka – model2") +
theme_minimal()

p1b + p2b

Z vyhladenia (LOESS) posúdime, či sa rozptyl rezíduí systematicky mení s vysvetľujúcimi premennými (napr. „lievik“).

Testovanie prítomnosti heteroskedasticity

library(lmtest)

# Breusch–Pagan test pre oba modely

bp1 <- bptest(model)
bp2 <- bptest(model2)

bp1

    studentized Breusch-Pagan test

data:  model
BP = 16.089, df = 4, p-value = 0.002902

bp2

    studentized Breusch-Pagan test

data:  model2
BP = 13.892, df = 3, p-value = 0.003056

Ak test indikuje heteroskedasticitu (malé p-hodnoty), použijeme robustné (HC) smerodajné chyby:

library(sandwich)
library(lmtest)

# Robustné (HC1) smerodajné chyby – odporúčané reportovať s OLS koeficientmi

coeftest(model,  vcov = vcovHC(model,  type = "HC1"))
coeftest(model2, vcov = vcovHC(model2, type = "HC1"))

Poznámka

Pri finančných premenných (ako rozpočet) býva heteroskedasticita častá; log-transformácia rozpočtu (čo už v modeli používame) je správny krok. Robustnné SE (HC) sú vhodné najmä pri väčších vzorkách. Alternatívou je špecifikovať model inak (napr. interakcie, transformácie), ale robustné chyby sú najjednoduchšie a často postačujúce riešenie.