Obsah

Práca s databázou – import údajov, grafy, štatistiky
Lineárna regresia
Špecifikácia modelu – hodnotenie filmov
Zhluková analýza filmov
Autokorelácia rezíduí
Multikolinearita v regresných modeloch

Práca s databázou - import údajov, grafy, štatistiky

Úvod

Nastavenie

library(dplyr)
library(tidyr)
library(ggplot2)
library(readr)
library(stringr)
library(scales)
library(knitr)
library(broom)
library(kableExtra)

options(knitr.kable.NA = "")

Načítanie a čistenie dát

# Cesta k CSV
path <- "data_ekonometria.csv"

# Európske CSV s bodkočiarkou, UTF-8
df_raw <- read.csv2(path, sep = ";", fileEncoding = "UTF-8", encoding = "UTF-8", stringsAsFactors = FALSE)


# --- Očistenie názvov stĺpcov (rozšírené) ---
# zjednotenie medzier, trimming, fix znakov
names(df_raw) <- gsub("\\s+", " ", names(df_raw))
names(df_raw) <- trimws(names(df_raw))
names(df_raw) <- iconv(names(df_raw), from = "UTF-8", to = "UTF-8")

# Pomocná funkcia na premenovanie podľa vzoru (prvý jednoznačný match)
rename_if_present <- function(cols, pattern, new_name) {
  idx <- which(grepl(pattern, cols, ignore.case = TRUE, perl = TRUE))
  if (length(idx) == 1) cols[idx] <- new_name
  cols
}

nms <- names(df_raw)
# IMDb
nms <- rename_if_present(nms, "(?i)\\bimdb\\b", "IMDb hodnotenie")
# Dĺžka (min) – toleruj diakritiku a/alebo len (min)
nms <- rename_if_present(nms, "(?i)(d[lľĺ]žka)|(\\(\\s*min\\s*\\))", "Dĺžka (min)")
# Rozpočet
nms <- rename_if_present(nms, "(?i)rozpo[cč]et|budget", "Rozpočet [$]")
names(df_raw) <- nms

print(names(df_raw))

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

# Voliteľne: jasná hláška, ak niečo stále chýba
required_cols <- c("Umiestnenie","Názov","Rok","Žáner","IMDb hodnotenie","Dĺžka (min)","Réžia","Spoločnosť","Rozpočet [$]")
missing_cols <- setdiff(required_cols, names(df_raw))
if (length(missing_cols) > 0) {
  warning(paste0("Upozornenie: v údajoch chýbajú stĺpce: ", paste(missing_cols, collapse = ", "),
                 ". Dokument bude pokračovať so stĺpcami, ktoré sú k dispozícii."))
}


# --- Očistenie názvov stĺpcov ---
names(df_raw) <- gsub("\\s+", " ", names(df_raw))       # zjednotí viac medzier
names(df_raw) <- trimws(names(df_raw))                    # odstráni medzery na začiatku/konci
names(df_raw) <- iconv(names(df_raw), from = "UTF-8", to = "UTF-8")  # opraví skryté znaky
if (any(grepl("IMDb", names(df_raw)))) {
  names(df_raw)[grepl("IMDb", names(df_raw))] <- "IMDb hodnotenie"
}
print(names(df_raw))

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

# Čistenie a typy premenných
df <- df_raw %>%
  mutate(
    Rok = suppressWarnings(as.integer(Rok)),
    `IMDb hodnotenie` = str_replace_all(`IMDb hodnotenie`, ",", "."),
    `IMDb hodnotenie` = suppressWarnings(as.numeric(`IMDb hodnotenie`)),
    `Dĺžka (min)` = suppressWarnings(as.numeric(`Dĺžka (min)`)),
    `Rozpočet [$]` = str_remove_all(`Rozpočet [$]`, "[^0-9]"),
    `Rozpočet [$]` = suppressWarnings(as.numeric(`Rozpočet [$]`)),
    Dekáda = paste0(floor(Rok / 10) * 10, "s")
  )

glue <- function(...) paste0(...)

cat(glue("Počet riadkov: ", nrow(df), "\n"))

## Počet riadkov: 250

cat(glue("Stĺpce: ", paste(names(df), collapse = ", "), "\n"))

## Stĺpce: Umiestnenie, Názov, Rok, Žáner, IMDb hodnotenie, Dĺžka (min), Réžia, Spoločnosť, Rozpočet [$], Dekáda

Výber dvoch skupín

group_var <- "Spoločnosť"

# 2 najpočetnejšie kategórie
top2 <- df %>%
  filter(!is.na(.data[[group_var]]), .data[[group_var]] != "") %>%
  count(.data[[group_var]], sort = TRUE) %>%
  slice_head(n = 2) %>%
  pull(1)

df2 <- df %>% filter(.data[[group_var]] %in% top2)

# Premenujeme na "Skupina" kvôli konzistentnosti
names(df2)[names(df2) == group_var] <- "Skupina"

head(dplyr::select(df2, Skupina, Názov, Rok, Žáner, `IMDb hodnotenie`, `Dĺžka (min)`, `Rozpočet [$]`))

Grafy

Stĺpcový graf: priemerné IMDb podľa dekády a skupiny

agg_decade <- df2 %>%
  group_by(Skupina, Dekáda) %>%
  summarise(
    avg_imdb = mean(`IMDb hodnotenie`, na.rm = TRUE),
    n = n(), .groups = "drop"
  )

ggplot(agg_decade, aes(x = Dekáda, y = avg_imdb, fill = Skupina, group = Skupina)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.8) +
  geom_text(aes(label = round(avg_imdb, 2)),
            position = position_dodge(width = 0.8), vjust = -0.3, size = 3) +
  theme_minimal() +
  labs(title = "Priemerné IMDb hodnotenie podľa dekády",
       x = "Dekáda", y = "Priemer IMDb", fill = "Skupina") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Čiarový graf: medián rozpočtu podľa roku a skupiny

agg_year_budget <- df2 %>%
  filter(!is.na(Rok)) %>%
  group_by(Skupina, Rok) %>%
  summarise(median_budget = median(`Rozpočet [$]`, na.rm = TRUE), .groups = "drop")

ggplot(agg_year_budget, aes(x = Rok, y = median_budget, color = Skupina, group = Skupina)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2.5) +
  geom_text(aes(label = label_number(accuracy = 1, scale_cut = scales::cut_si(""))(median_budget)),
            vjust = -0.5, size = 2.7, show.legend = FALSE) +
  theme_minimal() +
  labs(title = "Medián rozpočtu podľa roku",
       x = "Rok", y = "Medián rozpočtu [$]", color = "Skupina") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Základné štatistiky

stats_tbl <- df2 %>%
  group_by(Skupina, Dekáda) %>%
  summarise(
    n = n(),
    mean_imdb = mean(`IMDb hodnotenie`, na.rm = TRUE),
    sd_imdb   = sd(`IMDb hodnotenie`, na.rm = TRUE),
    min_imdb  = min(`IMDb hodnotenie`, na.rm = TRUE),
    q25_imdb  = quantile(`IMDb hodnotenie`, 0.25, na.rm = TRUE),
    median_imdb = median(`IMDb hodnotenie`, na.rm = TRUE),
    q75_imdb  = quantile(`IMDb hodnotenie`, 0.75, na.rm = TRUE),
    max_imdb  = max(`IMDb hodnotenie`, na.rm = TRUE),
    mean_len  = mean(`Dĺžka (min)`, na.rm = TRUE),
    mean_budget = mean(`Rozpočet [$]`, na.rm = TRUE),
    .groups = "drop"
  )

kable(stats_tbl, digits = 2, caption = "Základné štatistiky podľa skupiny a dekády")

Základné štatistiky podľa skupiny a dekády
Skupina	Dekáda	n	mean_imdb	sd_imdb	min_imdb	q25_imdb	median_imdb	q75_imdb	max_imdb	mean_len	mean_budget
United Artists	1920s	2	8.10	0.00	8.1	8.10	8.10	8.10	8.1	81.00	836500
United Artists	1930s	2	8.50	0.00	8.5	8.50	8.50	8.50	8.5	87.00	1500000
United Artists	1940s	3	8.20	0.17	8.1	8.10	8.10	8.25	8.4	118.00	1496667
United Artists	1950s	4	8.50	0.35	8.2	8.35	8.40	8.55	9.0	105.25	1534250
United Artists	1960s	4	8.22	0.05	8.2	8.20	8.20	8.22	8.3	152.00	2600000
United Artists	1970s	4	8.35	0.30	8.1	8.10	8.30	8.55	8.7	130.25	9940000
United Artists	1980s	1	8.20		8.2	8.20	8.20	8.20	8.2	129.00	18000000
Warner Bros. Pictures	1940s	2	8.35	0.21	8.2	8.27	8.35	8.43	8.5	144.00	1750000
Warner Bros. Pictures	1950s	1	8.20		8.2	8.20	8.20	8.20	8.2	105.00	1400000
Warner Bros. Pictures	1960s	1	8.10		8.1	8.10	8.10	8.10	8.1	127.00	3200000
Warner Bros. Pictures	1970s	3	8.17	0.12	8.1	8.10	8.10	8.20	8.3	147.67	8433333
Warner Bros. Pictures	1980s	4	8.30	0.14	8.1	8.25	8.35	8.40	8.4	152.00	27250000
Warner Bros. Pictures	1990s	6	8.35	0.29	8.0	8.20	8.25	8.60	8.7	134.17	41233333
Warner Bros. Pictures	2000s	7	8.39	0.31	8.1	8.20	8.20	8.50	9.0	134.14	81000000
Warner Bros. Pictures	2010s	7	8.36	0.28	8.1	8.10	8.40	8.50	8.8	139.29	151871429

Testovanie hypotéz

t-test: rozdiel IMDb medzi 2 skupinami v najnovšej dekáde

latest_decade <- df2 %>%
  filter(!is.na(Dekáda)) %>%
  mutate(dec_num = as.integer(str_remove(Dekáda, "s"))) %>%
  slice_max(dec_num, n = 1, with_ties = FALSE) %>%
  pull(Dekáda) %>%
  unique()

df_latest <- df2 %>% filter(Dekáda == latest_decade)

ok_groups <- df_latest %>%
  group_by(Skupina) %>%
  filter(!is.na(`IMDb hodnotenie`)) %>%
  summarise(n = n(), .groups = "drop") %>%
  filter(n >= 2) %>% pull(Skupina)

if (length(ok_groups) == 2) {
  x <- df_latest %>% filter(Skupina == ok_groups[1]) %>% pull(`IMDb hodnotenie`)
  y <- df_latest %>% filter(Skupina == ok_groups[2]) %>% pull(`IMDb hodnotenie`)
  t_test_result <- t.test(x, y)
  t_test_result
} else {
  "t-test sa neuskutočnil: v najnovšej dekáde nie je dosť dát v oboch skupinách."
}

## [1] "t-test sa neuskutočnil: v najnovšej dekáde nie je dosť dát v oboch skupinách."

ANOVA: IMDb podľa žánru (top 5 žánrov)

top_genres <- df %>% count(Žáner, sort = TRUE) %>% slice_head(n = 5) %>% pull(Žáner)
anova_df <- df %>% filter(Žáner %in% top_genres, !is.na(`IMDb hodnotenie`))

anova_result <- aov(`IMDb hodnotenie` ~ Žáner, data = anova_df)
summary(anova_result)

##             Df Sum Sq Mean Sq F value Pr(>F)
## Žáner        4  0.594 0.14848   1.706  0.162
## Residuals   56  4.874 0.08703

Lineárna regresia

Model: IMDb ~ Rok + Skupina (dve vybrané skupiny z Spoločnosť alebo Žáner).

model <- lm(`IMDb hodnotenie` ~ Rok + Skupina, data = df2)
summary(model)

## 
## Call:
## lm(formula = `IMDb hodnotenie` ~ Rok + Skupina, data = df2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3341 -0.1684 -0.0516  0.1272  0.6818 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)
## (Intercept)                   4.729720   3.439589   1.375    0.175
## Rok                           0.001834   0.001759   1.042    0.302
## SkupinaWarner Bros. Pictures -0.061108   0.095570  -0.639    0.526
## 
## Residual standard error: 0.2414 on 48 degrees of freedom
## Multiple R-squared:  0.02238,    Adjusted R-squared:  -0.01836 
## F-statistic: 0.5494 on 2 and 48 DF,  p-value: 0.5809

coef.tbl <- tidy(model, conf.int = TRUE) %>%
  mutate(
    term = dplyr::case_when(
      term == "(Intercept)" ~ "Intercept",
      term == "Rok" ~ "Rok",
      TRUE ~ term
    ),
    stars = dplyr::case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      p.value < 0.1   ~ "·",
      TRUE            ~ ""
    )
  ) %>%
  transmute(
    Term = term,
    Estimate = estimate,
    `Std. Error` = std.error,
    `t value` = statistic,
    `p value` = p.value,
    `95% CI` = stringr::str_c("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"),
    Sig = stars
  )

coef.tbl %>%
  kable(
    digits = 2,
    caption = "OLS regresné koeficienty (IMDb ~ Rok + Skupina)"
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  column_spec(1, bold = TRUE) %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2") %>%
  footnote(
    general = "Signif. codes: *** p<0.001, ** p<0.01, * p<0.05, · p<0.1.",
    threeparttable = TRUE
  )

OLS regresné koeficienty (IMDb ~ Rok + Skupina)
Term	Estimate	Std. Error	t value	p value	95% CI
Intercept	4.73	3.44	1.38	0.18	[-2.19, 11.65]
Rok	0.00	0.00	1.04	0.30	[0, 0.01]
SkupinaWarner Bros. Pictures	-0.06	0.10	-0.64	0.53	[-0.25, 0.13]
Note:
Signif. codes: * p<0.001, p<0.01, * p<0.05, · p<0.1.

Lineárna regresia

Knižnice:

library(dplyr)
library(stringr)
library(lmtest)
library(sandwich)
library(car)
library(tseries)

# 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.

Ú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 \]

# ===== 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

Odhadnuté koeficienty: model$coefficients
Reziduá: model$residuals
Vyrovnané hodnoty: model$fitted.values
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.

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)

##   rstudent unadjusted p-value Bonferroni p
## 1 4.572282         7.6738e-06    0.0019184
## 2 3.794663         1.8653e-04    0.0466340

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)

##   rstudent unadjusted p-value Bonferroni p
## 1 4.582792         7.3138e-06    0.0018284
## 2 3.825446         1.6570e-04    0.0414260

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).

# 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:

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"))

## 
## 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

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

## 
## t test of coefficients:
## 
##                 Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)   7.88406236  0.11573503 68.1217 < 2.2e-16 ***
## log_rozpocet  0.00923415  0.00686918  1.3443 0.1800937    
## Dlzka_min     0.00203421  0.00060414  3.3671 0.0008814 ***
## DramaTRUE    -0.00402805  0.02982879 -0.1350 0.8926914    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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.

Špecifikácia modelu – hodnotenie filmov

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)

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.

#######################################################################
# 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.

# 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.

\[ \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)

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)

# RESET test pre kvadratický model
resettest(model_kvadr)

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

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) zavedieme 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)

resettest(model_rozsireny)

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

Interpretácia: - 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

Zavedieme 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ť:

na zlom v autonómnom (konštantnom) člene,
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

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

anova(model, modelD_sklon)

resettest(modelD_sklon)

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

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ú

Zhluková analýza filmov

Úvod

Keďže pôvodná databáza obsahuje pomerne veľa filmov a úplné názvy filmov by v grafoch a tabuľkách zhoršovali prehľadnosť, urobíme dve zjednodušenia:

Analyzujeme náhodne vybraných 20 filmov z rebríčka. Vďaka tomu zostanú výstupy prehľadné, ale zároveň bude možné dobre sledovať, ako algoritmus zhlukovania funguje.
Filmy budeme v analytických tabuľkách a grafoch označovať len ich umiestnením v rebríčku (premenná „Umiestnenie“). Plné názvy filmov ponecháme len v úvodnej tabuľke, ktorá opisuje použitú vzorku. V ďalšej analýze preto čitateľ pracuje najmä s číselnými identifikátormi filmov.

V ďalšom texte budeme pracovať s údajmi načítanými zo súboru data_ekonometria.csv. V Tab. 1 uvádzame náhodne vybranú vzorku 20 filmov, ktorá vstupuje do zhlukovej analýzy.

# Načítanie údajov o filmoch
filmy <- read.csv(
  "data_ekonometria.csv",
  sep = ";",
  stringsAsFactors = FALSE,
  check.names = FALSE   # ponechá pôvodné názvy stĺpcov s medzerami
)

# Pre didaktické účely vyberieme náhodnú vzorku 20 filmov
set.seed(123)  # aby bol výber reprodukovateľný
filmy_sample <- filmy %>%
  dplyr::slice_sample(n = 20)

# Prevod textových premenných na numerické -----------------------------

# IMDb hodnotenie je vo formáte s čiarkou ako desatinným oddeľovačom (napr. "9,3")
filmy_sample$IMDb_hodnotenie_num <- as.numeric(sub(",", ".", filmy_sample$`IMDb hodnotenie`))

# Dĺžka v minútach – premenujeme na jednoduchší názov
filmy_sample$Dlzka <- filmy_sample$`Dĺžka (min)`

# Rozpočet je zapísaný s medzerami ako oddeľovačmi tisícov (napr. "25 000 000")
filmy_sample$Rozpocet <- as.numeric(gsub(" ", "", filmy_sample$`Rozpočet [$]`))

# Výber premenných vhodných na zhlukovanie (A: IMDb + dĺžka + rozpočet)
filmy_num <- filmy_sample[, c("IMDb_hodnotenie_num", "Dlzka", "Rozpocet")]

# Ako identifikátor použijeme umiestnenie filmu v rebríčku IMDb
rownames(filmy_num) <- filmy_sample$Umiestnenie

# Odstránenie prípadných riadkov s chýbajúcimi hodnotami
filmy_num_complete <- na.omit(filmy_num)

# Úvodná prehľadová tabuľka (ukazujeme aj názov filmu)
tab1 <- filmy_sample[, c("Umiestnenie", "Názov", "Rok",
                         "Žáner", "IMDb hodnotenie", "Dĺžka (min)", "Rozpočet [$]")]

kable(tab1,
      caption = "Tab. 1: Náhodne vybraných 20 filmov z rebríčka IMDb") %>%
  kable_styling(full_width = FALSE)

Tab. 1: Náhodne vybraných 20 filmov z rebríčka IMDb
Umiestnenie	Názov	Rok	Žáner	IMDb hodnotenie	Dĺžka (min)	Rozpočet [$]
159	Gone with the Wind	1939	Dráma; Romantický; Vojnový	8,2	238	3 850 000
207	Tokyo Story	1953	Dráma	8,1	136	120 000
179	Klaus	2019	Animovaný; Dobrodružný; Komédia	8,1	96	40 000 000
14	Inception	2010	Akčný; Dobrodružný; Sci-Fi	8,8	148	160 000 000
195	Sherlock Jr.	1924	Akčný; Komédia; Romantický	8,1	45	200 000
170	Fargo	1996	Krimi; Triler	8,1	98	7 000 000
50	Cinema Paradiso	1988	Dráma; Romantický	8,5	155	5 000 000
118	Die Hard	1988	Akčný; Triler	8,2	132	33 000 000
43	Casablanca	1942	Dráma; Romantický; Vojnový	8,5	162	1 000 000
229	La haine	1995	Krimi; Dráma	8,1	98	2 600 000
247	The Help	2011	Dráma	8,0	146	25 000 000
243	Persona	1966	Dráma; Triler	8,1	85	78 000
153	The Thing	1982	Horor; Mysteriózny; Sci-Fi	8,2	109	15 000 000
90	Eternal Sunshine of the Spotless Mind	2004	Dráma; Romantický; Sci-Fi	8,3	108	20 000 000
91	2001: A Space Odyssey	1968	Dobrodružný; Sci-Fi	8,3	149	10 500 000
197	Mr. Smith Goes to Washington	1939	Komédia; Dráma	8,1	129	1 500 000
236	Amores perros	2000	Dráma; Triler	8,1	154	2 400 000
185	The Grand Budapest Hotel	2014	Dobrodružný; Komédia; Krimi	8,1	99	25 000 000
92	Reservoir Dogs	1992	Krimi; Triler	8,3	99	3 000 000
137	Pan’s Labyrinth	2006	Dráma; Fantazijný; Vojnový	8,2	118	19 000 000

V ďalšej analýze pracujeme už iba s číselným identifikátorom filmu – umiestnením v rebríčku IMDb. Všetky tabuľky a grafy preto používajú namiesto názvu filmu číslo Umiestnenie, aby zostali prehľadné. Plné názvy filmov je možné spätne dohľadať v Tab. 1.

Škálovanie premenných a boxploty

Hierarchická zhluková analýza pracuje s mierami vzdialenosti medzi pozorovaniami. Aby boli tieto vzdialenosti porovnateľné, je vhodné, aby všetky premenné boli na rovnakej škále. Preto použijeme tzv. z-škálovanie, pričom transformované $z$ hodnoty (skóre) vypočítame:

\[ z = \frac{x - \mu}{\sigma}, \]

kde $\mu$ je stredná hodnota a $\sigma$ je štandardná odchýlka pozorovaní $x$. Predpokladáme pritom, že súbor údajov už neobsahuje chýbajúce hodnoty, ktoré sme ošetrili v predchádzajúcom kroku.

Touto operáciou získame škálované pozorovania, ktorých rozloženie pre jednotlivé premenné znázorňujeme pomocou boxplotov.

# =======================================================
# 1) Príprava údajov a data.frame so škálovanými údajmi
# =======================================================

filmy_scaled <- scale(filmy_num_complete)

Obr. 1. Boxploty škálovaných numerických premenných (IMDb hodnotenie, dĺžka a rozpočet)

num_vars <- as.data.frame(filmy_scaled)
num_plots <- ncol(num_vars)

par(mfrow = c(ceiling(sqrt(num_plots)), ceiling(num_plots / ceiling(sqrt(num_plots)))))
par(mar = c(4, 4, 2, 1))

for (col in names(num_vars)) {
  boxplot(num_vars[[col]],
          main = col,
          col = "lightblue",
          horizontal = TRUE)
}

mtext("Boxploty numerických premenných (filmy)", outer = TRUE, cex = 1.3, font = 2)

Prípadné odľahlé hodnoty (napríklad extrémne vysoký rozpočet alebo dĺžka) nebudeme vylučovať, keďže predstavujú konkrétne filmy, ktoré sú z hľadiska interpretácie zaujímavé.

Korelačná matica premenných

Pri zhlukovej analýze je dôležitá korelačná matica premenných. Vysoká korelácia môže zvýhodňovať niektoré premenné pri tvorbe klastrov. Pri veľmi vysokej korelácii (napr. nad 0,8 alebo 0,9) by sme uvažovali o vylúčení jednej z dvojice premenných alebo o použití analýzy hlavných komponentov.

V Tab. 2 uvádzame korelačnú maticu troch použitých premenných: IMDb hodnotenia, dĺžky a rozpočtu.

cor_mat <- cor(filmy_scaled, use = "pairwise.complete.obs")
cor_mat <- round(cor_mat, 2)

kable(cor_mat,
      caption = "Tab. 2: Korelačná matica škálovaných premenných") %>%
  kable_styling(full_width = FALSE)

Tab. 2: Korelačná matica škálovaných premenných
	IMDb_hodnotenie_num	Dlzka	Rozpocet
IMDb_hodnotenie_num	1.00	0.32	0.63
Dlzka	0.32	1.00	0.09
Rozpocet	0.63	0.09	1.00

V našom prípade pracujeme iba s tromi premennými, ktoré typicky nebývajú extrémne silno korelované (hodnotenie, dĺžka a rozpočet), takže nie je nutné žiadnu z nich vylučovať.

Matica vzdialeností

Každému filmu zodpovedá jeden riadok pozorovaní. Vzdialenosť medzi filmami $i$ a $j$ je pri použití Euklidovskej vzdialenosti definovaná:

\[ d^{ij} = \sqrt{\sum_k (x^i_k - x^j_k)^2}, \]

kde $x^i_k$ je $k$-ta premenná (IMDb hodnotenie, dĺžka, rozpočet) pre film $i$. Tento typ vzdialenosti nazývame Euklidovská vzdialenosť. Vzdialenosti medzi jednotlivými filmami sa súhrnne vyjadrujú v matici vzdialenosti, ktorá je uvedená v Tab. 3.

Interpretácia je nasledovná: čím je hodnota v matici väčšia, tým sú si dva filmy z hľadiska zvolených premenných menej podobné (líšia sa napríklad v rozpočte, dĺžke alebo hodnotení). Naopak, malé vzdialenosti znamenajú podobnosť. V tabuľke používame ako identifikátor umiestnenie filmu v rebríčku (Umiestnenie), aby bola matica prehľadná.

# ============================
# 3) Distance matrix
# ============================

dist_mat <- round(dist(filmy_scaled, method = "euclidean"), 2)

as.matrix(dist_mat)[1:10, 1:10] %>%
  kable(caption = "Tab. 3: Časť matice Euklidovských vzdialeností medzi filmami (podľa umiestnenia v rebríčku)") %>%
  kable_styling(full_width = FALSE)

Tab. 3: Časť matice Euklidovských vzdialeností medzi filmami (podľa umiestnenia v rebríčku)
	159	207	179	14	195	170	50	118	43	229
159	0.00	2.63	3.76	5.88	4.90	3.57	2.62	2.80	2.48	3.57
207	2.63	0.00	1.51	5.83	2.30	0.98	2.15	1.07	2.20	0.96
179	3.76	1.51	0.00	5.17	1.71	0.93	2.75	1.07	2.90	1.06
14	5.88	5.83	5.17	0.00	6.37	5.81	4.66	4.79	4.78	5.90
195	4.90	2.30	1.71	6.37	0.00	1.35	3.48	2.44	3.62	1.34
170	3.57	0.98	0.93	5.81	1.35	0.00	2.54	1.25	2.65	0.12
50	2.62	2.15	2.75	4.66	3.48	2.54	0.00	1.85	0.21	2.54
118	2.80	1.07	1.07	4.79	2.44	1.25	1.85	0.00	1.97	1.32
43	2.48	2.20	2.90	4.78	3.62	2.65	0.21	1.97	0.00	2.65
229	3.57	0.96	1.06	5.90	1.34	0.12	2.54	1.32	2.65	0.00

(Pre prehľadnosť zobrazujeme len časť matice – prvých 10 filmov v riadkoch aj stĺpcoch.)

Princíp hierarchického zhlukovania (Wardova metóda)

Zhlukovanie pri Wardovej metóde prebieha zdola smerom nahor – začíname s jednočlennými klastrami a postupne zlučujeme dvojice klastrov. Táto metóda patrí medzi aglomeratívne hierarchické metódy. Minimalizuje nárast vnútornej variability pri spojení dvoch klastrov, pričom využíva nasledovné výpočty:

Wardova metóda minimalizuje sumu štvorcov chýb (Error Sum of Squares – ESS):

\[ ESS(C) = \sum_{i \in C} \lVert x_i - \bar{x}_C \rVert^2, \]

kde $C$ je zvažovaný klaster. V každom kroku zlučovania dvoch klastrov Wardova metóda hľadá minimálny prírastok sumy štvorcov chýb:

\[ \Delta ESS = ESS(A \cup B) - ESS(A) - ESS(B). \]

Obr. 2. Hierarchické zhlukovanie – dendrogram. Červená čiara určuje rez definujúci tri klastre.

# ============================
# 4) Hierarchical clustering
# ============================

hc <- hclust(dist_mat, method = "ward.D2")

plot(hc,
     labels = rownames(filmy_scaled),   # zobrazujeme len umiestnenie filmu
     main = "Hierarchické zhlukovanie filmov (Ward.D2)",
     xlab = "Umiestnenie filmu v rebríčku IMDb",
     sub = "")

k <- 3  # počet klastrov
h_cut <- hc$height[length(hc$height) - (k - 1)]
abline(h = h_cut, col = "red", lwd = 2, lty = 2)

klaster_membership <- cutree(hc, k = k)

filmy_klasters <- data.frame(
  Umiestnenie = rownames(filmy_num_complete),
  filmy_num_complete,
  klaster = factor(klaster_membership)
)

Tab. 4. Príslušnosť filmov do klastrov (podľa umiestnenia v rebríčku)

kable(filmy_klasters[, c("Umiestnenie", "klaster")],
      caption = "Tab. 4: Príslušnosť filmov do klastrov (identifikátorom je umiestnenie v rebríčku)") %>%
  kable_styling(full_width = FALSE)

Tab. 4: Príslušnosť filmov do klastrov (identifikátorom je umiestnenie v rebríčku)
	Umiestnenie	klaster
159	159	1
207	207	2
179	179	2
14	14	3
195	195	2
170	170	2
50	50	1
118	118	2
43	43	1
229	229	2
247	247	2
243	243	2
153	153	2
90	90	2
91	91	2
197	197	2
236	236	2
185	185	2
92	92	2
137	137	2

Vzhľadom na to, že filmy označujeme iba ich umiestnením v rebríčku, sú tabuľky a grafy prehľadné, pričom podrobné názvy filmov je možné dohľadať v Tab. 1.

Deskriptívne štatistiky – rozklad variability

Zaujíma nás, aká je variabilita jednotlivých premenných vo vnútri a medzi klastrami. Použijeme rozklad variability na:

TSS – celkovú sumu štvorcov odchýlok,
WSS – vnútroklastrovú variabilitu,
BSS – medzi-klastrovú variabilitu.

Na základe Tab. 5 môžeme posúdiť, ako dobre jednotlivé premenné prispievajú k odlíšeniu klastrov. Čím je podiel medzi-klastrovej variability na celkovej variabilite vyšší, tým lepšie daná premenná pomáha oddeľovať zhluky.

# ============================
# 5) Variability measures
# ============================

ssq <- function(x, m) sum((x - m)^2)

var_names <- colnames(filmy_scaled)

TSS <- sapply(var_names, function(v) ssq(filmy_scaled[, v], mean(filmy_scaled[, v])))

WSS <- sapply(var_names, function(v) {
  x <- filmy_scaled[, v]
  tapply(x, klaster_membership, function(z) ssq(z, mean(z))) |> sum()
})

BSS <- TSS - WSS

ss_table <- data.frame(
  Variable     = var_names,
  TSS          = TSS,
  WSS          = WSS,
  BSS          = BSS,
  Prop_Between = BSS / TSS
)

# Zaokrúhľujeme iba numerické stĺpce
ss_table_round <- ss_table
ss_table_round[, sapply(ss_table_round, is.numeric)] <-
  round(ss_table_round[, sapply(ss_table_round, is.numeric)], 3)

kable(ss_table_round,
      caption = "Tab. 5: Rozklad variability (celková, vnútri klastrov a medzi klastrami)") %>%
  kable_styling(full_width = FALSE)

Tab. 5: Rozklad variability (celková, vnútri klastrov a medzi klastrami)
	Variable	TSS	WSS	BSS	Prop_Between
IMDb_hodnotenie_num	IMDb_hodnotenie_num	19	4.942	14.058	0.740
Dlzka	Dlzka	19	10.216	8.784	0.462
Rozpocet	Rozpocet	19	2.005	16.995	0.894

Na základe tabuľky vidíme, že jednotlivé premenné sa líšia v tom, akú časť variability dokážeme vysvetliť rozdielmi medzi klastrami. Premenné s vyšším podielom medzi-klastrovej variability (vyššia hodnota Prop_Between) lepšie separujú jednotlivé skupiny filmov.

Centroidy – priemerné hodnoty premenných v klastroch

Záverečnú interpretáciu klastrov urobíme na základe tzv. centroidov, teda priemerných hodnôt sledovaných premenných v jednotlivých klastroch.

Tab. 6. Centroidy – priemerné hodnoty IMDb hodnotenia, dĺžky a rozpočtu podľa klastrov

descriptives <- filmy_klasters %>%
  group_by(klaster) %>%
  summarise(
    IMDb_hodnotenie_mean = mean(IMDb_hodnotenie_num, na.rm = TRUE),
    Dlzka_mean           = mean(Dlzka, na.rm = TRUE),
    Rozpocet_mean        = mean(Rozpocet, na.rm = TRUE)
  )

# zaokrúhľujeme iba numerické stĺpce, nie stĺpec "klaster"
descriptives_round <- descriptives
descriptives_round[, sapply(descriptives_round, is.numeric)] <-
  round(descriptives_round[, sapply(descriptives_round, is.numeric)], 2)

kable(descriptives_round,
      caption = "Tab. 6: Priemerné hodnoty sledovaných premenných v jednotlivých klastroch") %>%
  kable_styling(full_width = FALSE)

Tab. 6: Priemerné hodnoty sledovaných premenných v jednotlivých klastroch
klaster	IMDb_hodnotenie_mean	Dlzka_mean	Rozpocet_mean
1	8.40	185.00	3283333
2	8.15	112.56	12774875
3	8.80	148.00	160000000

Porovnaním centroidov môžeme charakterizovať jednotlivé klastre. Klaster, ktorý má najvyšší priemerný rozpočet a dĺžku, môžeme interpretovať ako skupinu veľkých, produkčne náročných filmov (často ide o tzv. „blockbustery“). Hodnotenie IMDb nám zároveň umožňuje rozlíšiť, ktorý klaster má v priemere vyššiu divácku odozvu.

Záver

Predložená analýza sa zaoberá zhlukovaním filmov na základe troch ukazovateľov: IMDb hodnotenia, dĺžky filmu a výšky rozpočtu. Na základe hierarchickej zhlukovej analýzy (Wardova metóda) sme identifikovali tri klastre, ktoré združujú filmy s podobnými charakteristikami.

Z praktických dôvodov sme pracovali iba s náhodne vybranou vzorkou 20 filmov z rebríčka IMDb a filmy sme v grafoch a tabuľkách označovali iba ich umiestnením v rebríčku. Tým sme dosiahli dobrú prehľadnosť výstupov, pričom plné názvy filmov je možné spätne dohľadať v úvodnej tabuľke.

Analýza ukázala, že:

zvolená kombinácia premenných umožňuje vytvoriť zhluky filmov, ktoré sa líšia najmä z hľadiska rozpočtu a dĺžky,
priemerné hodnoty (centroidy) klastrov umožňujú interpretovať jednotlivé skupiny ako viac či menej nákladné a dlhé filmy s rôznym priemerným hodnotením,
rozklad variability naznačuje, ktoré premenné najviac prispievajú k oddeleniu klastrov.

Autokorelácia rezíduí

1. Úvod a údaje

Databáza obsahuje najlepšie hodnotené filmy podľa IMDb s informáciou o roku uvedenia, žánri, dĺžke trvania, režisérovi, filmovom štúdiu a rozpočte. Nie všetky údaje budú použité, preto sa zameriame len na filmy jednej filmovej spoločnosti, aby sme mali časové usporiadanie – konkrétne na filmy štúdia Warner Bros. Pictures.

1.1 Úvod do problému, stanovenie hypotéz

Pracovná hypotéza je, že:

vyšší rozpočet by mohol viesť k lepšiemu filmu a vyššiemu IMDb hodnoteniu (očakávame kladné znamienko pri koeficiente rozpočtu),
primerane dlhší film (väčší priestor na rozvinutie deja) môže byť vnímaný lepšie, takže aj pri dĺžke očakávame skôr kladný vplyv,
pri roku uvedenia nepredpokladáme jednoznačný smer – ide skôr o kontrolnú premennú, ktorá zachytáva možné zmeny preferencií divákov v čase.

1.2 Príprava databázy, úprava údajov

Najprv načítame celý súbor, upravíme formát čísel a následne vyberieme len filmy štúdia Warner Bros. Pictures:

# kontrola, či súbor existuje v aktuálnom pracovnom priečinku
if (!file.exists("data_ekonometria.csv")) {
  stop("Súbor 'data_ekonometria.csv' sa nenašiel v aktuálnom priečinku. 
Skontroluj, prosím, že Rmd a CSV sú v tom istom priečinku.")
}

# načítanie údajov – bodkočiarka, desatinná čiarka
udaje <- read.csv2("data_ekonometria.csv",
                   header = TRUE,
                   stringsAsFactors = FALSE)

# premenovanie stĺpcov do "R-friendly" tvaru
colnames(udaje) <- c("Umiestnenie", "Nazov", "Rok", "Zaner",
                     "IMDb_raw", "Dlzka_min", "Rezia",
                     "Spolocnost", "Rozpocet_raw")

# úprava IMDb hodnotenia (z formátu "9,3" na numerické 9.3)
udaje$IMDb <- as.numeric(gsub(",", ".", udaje$IMDb_raw))

# úprava rozpočtu (odstránenie medzier a prevod na číslo)
udaje$Rozpocet <- as.numeric(gsub(" ", "", udaje$Rozpocet_raw))

# výber filmov Warner Bros. Pictures a relevantných premenných
udajeWB <- udaje[udaje$Spolocnost == "Warner Bros. Pictures",
                 c("Rok", "IMDb", "Dlzka_min", "Rozpocet")]

# zoradenie podľa roku
udajeWB <- udajeWB[order(udajeWB$Rok), ]

udajeWB

2. Lineárna regresia v základnom tvare

Budeme odhadovať rovnicu:

\[ \text{IMDb}_t = \beta_0 + \beta_1 \text{Dlzka}_t + \beta_2 \text{Rozpocet}_t + \beta_3 \text{Rok}_t + e_t, \]

kde:

$\text{IMDb}_t$ – IMDb hodnotenie vybraného filmu Warner Bros. v roku $t$,
$\text{Dlzka}_t$ – dĺžka filmu v minútach,
$\text{Rozpocet}_t$ – rozpočet filmu v dolároch,
$\text{Rok}_t$ – rok uvedenia filmu,
$e_t$ – náhodná zložka.

# základný regresný model
model <- lm(IMDb ~ Dlzka_min + Rozpocet + Rok, data = udajeWB)
summary(model)

## 
## Call:
## lm(formula = IMDb ~ Dlzka_min + Rozpocet + Rok, data = udajeWB)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33337 -0.13538 -0.05900  0.09991  0.53793 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 6.269e+00  5.168e+00   1.213   0.2356  
## Dlzka_min   3.413e-03  1.583e-03   2.156   0.0402 *
## Rozpocet    6.779e-10  7.355e-10   0.922   0.3648  
## Rok         7.713e-04  2.586e-03   0.298   0.7678  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2257 on 27 degrees of freedom
## Multiple R-squared:  0.2249, Adjusted R-squared:  0.1388 
## F-statistic: 2.612 on 3 and 27 DF,  p-value: 0.07183

3. Autokorelácia rezíduí

V tejto časti sa pozrieme na ďalší dôležitý predpoklad klasického lineárneho regresného modelu – nezávislosť rezíduí. Pri údajoch usporiadaných v čase (napr. filmy jedného štúdia podľa roku uvedenia) sa môže stať, že chyba v čase $t$ je systematicky spätá s chybou v čase $t-1$. Vtedy hovoríme o autokorelácii rezíduí.

3.1 Čo je autokorelácia rezíduí?

Autokorelácia rezíduí je situácia, keď platí:

\[ \operatorname{Corr}(e_t, e_{t-k}) \neq 0 \quad \text{pre niektoré } k \neq 0. \]

Najčastejšie sa skúma autokorelácia prvého rádu:

\[ e_t = \rho e_{t-1} + \nu_t,\quad |\rho| < 1. \]

3.2 Dôsledky autokorelácie rezíduí

odhady koeficientov $\hat{\beta}$ sú síce nestranné,
ale neefektívne,
štandardné chyby môžu byť podhodnotené → p-hodnoty sa javia menšie, než by mali byť,
môže dochádzať k falošnému zdaniu štatistickej významnosti,
teda t-testy a F-testy sú skreslené.

3.3 Detekcia autokorelácie

Grafická informácia

Najprv sa pozrieme na to, ako model vysvetlí vývoj IMDb hodnotenia filmov v čase. Pridáme k údajom vyrovnané (fitted) hodnoty z regresie a zobrazíme ich spolu s empirickými bodmi.

# pridáme vyrovnané hodnoty do dát
udajeWB$fitted <- fitted(model)

# scatterplot + vyrovnaná línia
ggplot(udajeWB, aes(x = Rok, y = IMDb)) +
  geom_point(size = 2) +
  geom_line(aes(y = fitted), size = 1) +
  labs(
    title = "IMDb hodnotenie filmov Warner Bros. Pictures: empirické vs. vyrovnané hodnoty",
    x = "Rok",
    y = "IMDb hodnotenie"
  ) +
  theme_minimal()

Z grafu vidíme, že hodnotenia filmov Warner Bros. sú rozmiestnené okolo vyrovnanej línie modelu relatívne pravidelne a nevidíme dlhé súvislé úseky, kde by boli rezíduá iba kladne alebo iba záporne. Už vizuálne teda nemáme silné podozrenie na výraznú autokoreláciu.

# uložíme si reziduá z pôvodného modelu
res <- residuals(model)

ACF graf (Autocorrelation Function)

Táto funkcia priraďuje odhad korelácie medzi rezíduami v aktuálnom období a rezíduami posunutými o $k$ období späť.

acf(res, lag.max = 4, main = "Autokorelačná funkcia rezíduí")

Na grafe je prerušovanou čiarou vyjadrený aj 95 % interval spoľahlivosti pre hodnotu autokorelačného koeficientu pri danom posune. V našich údajoch žiadny z odhadnutých autokorelačných koeficientov výrazne nevybočuje z tohto intervalu, takže ACF graf nenaznačuje výraznú sériovú koreláciu rezíduí.

Durbin–Watsonov test

Durbin–Watsonov koeficient (DW) je vypočítaný z rezíduí podľa vzorca

\[ DW = \frac{\sum_{t=2}^{n} (e_t - e_{t-1})^{2}}{\sum_{t=1}^{n} e_t^{2}}, \]

kde $n$ je počet pozorovaní. Medzi koeficientom autokorelácie dvoch susedných rezíduí a DW platí približný vzťah

\[ \hat{\rho} \approx 1 - \frac{DW}{2}. \]

Hodnoty DW sú v intervale od 0 po 4.

hodnoty blízke 0 → silná pozitívna autokorelácia,
hodnoty blízke 4 → silná negatívna autokorelácia,
hodnoty okolo 2 → bez problémov s autokoreláciou prvého rádu.

dwtest(model)

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 2.3437, p-value = 0.7856
## alternative hypothesis: true autocorrelation is greater than 0

V našom prípade sa Durbin–Watsonov koeficient pohybuje blízko hodnoty 2 a p-hodnota testu nenaznačuje štatisticky významnú autokoreláciu prvého rádu rezíduí.

Breusch–Godfreyov test (BG test)

BG test je formálnym testom autokorelácie (sériovej korelácie) rezíduí regresného modelu až do zvoleného rádu $p$.

Pre test autokorelácie až do rádu $p$ platí:

$H_0$: $\rho_1 = \rho_2 = \cdots = \rho_p = 0$ – žiadna sériová korelácia,
$H_1$: aspoň jedna z $\rho_j \neq 0$ – sériová korelácia prítomná.

bgtest(model, order = 1)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model
## LM test = 1.381, df = 1, p-value = 0.2399

V našom prípade BG test rovnako ako Durbin–Watsonov test nepotvrdzuje prítomnosť autokorelácie prvého rádu – p-hodnota je relatívne vysoká, takže nulovú hypotézu o „žiadnej autokorelácii“ nezamietame. Napriek tomu si demonštratívne ukážeme, ako by bolo možné model dynamizovať pomocou Koyckovho prístupu.

4. Dynamizácia modelu – Koyckov model

Teoreticky vychádzame z modelu s distribuovaným oneskorením, kde efekty vysvetľujúcich premenných v čase klesajú geometricky. Koyckova transformácia vedie k modelu:

\[ \text{IMDb}_t = \alpha(1-\lambda) + \beta_0 \text{Rozpocet}_t + \gamma_0 \text{Dlzka}_t + \lambda \text{IMDb}_{t-1} + v_t, \]

kde $\lambda$ reprezentuje rýchlosť prispôsobenia a $\text{IMDb}_{t-1}$ je oneskorená hodnota vysvetľovanej premennej.

V našom prípade:

$y_t$ = IMDb hodnotenie filmu Warner Bros. v roku $t$,
vysvetľujúce premenné: rozpočet a dĺžka filmu.

4.1 Odhad Koyckovho modelu v R

udajeWB <- udajeWB %>%
  arrange(Rok) %>%
  mutate(
    IMDb_lag1 = dplyr::lag(IMDb)
  )

model_koyck <- lm(IMDb ~ Rozpocet + Dlzka_min + IMDb_lag1,
                  data = udajeWB)

summary(model_koyck)

## 
## Call:
## lm(formula = IMDb ~ Rozpocet + Dlzka_min + IMDb_lag1, data = udajeWB)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32467 -0.13251 -0.00605  0.07865  0.52386 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.671e+00  1.693e+00   5.713 5.18e-06 ***
## Rozpocet     8.675e-10  5.726e-10   1.515    0.142    
## Dlzka_min    2.258e-03  1.751e-03   1.289    0.209    
## IMDb_lag1   -2.071e-01  1.887e-01  -1.097    0.283    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2231 on 26 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2568, Adjusted R-squared:  0.171 
## F-statistic: 2.994 on 3 and 26 DF,  p-value: 0.04905

dwtest(model_koyck)

## 
##  Durbin-Watson test
## 
## data:  model_koyck
## DW = 1.986, p-value = 0.4658
## alternative hypothesis: true autocorrelation is greater than 0

V dynamizovanom modeli je koeficient pri oneskorenej IMDb ($\text{IMDb}_{t-1}$) v našich údajoch relatívne malý a štatisticky nevýznamný, čo je v súlade s predchádzajúcimi testami – výrazný zotrvačný efekt v hodnoteniach filmov sa tu nepotvrdil. Porovnanie ukazovateľa Adjusted R-squared medzi pôvodným a dynamickým modelom tiež ukazuje, že pridaním oneskorenej hodnoty IMDb sa kvalita modelu výrazne nezlepšila.

4.2 Newey–West robustné (alebo aspoň robustné) štandardné chyby

Aj keď sme v našich údajoch nenašli silný dôkaz o autokorelácii, v praxi je často vhodné pracovať s robustnými štandardnými chybami. Pri krátkych časových radoch však môže výpočet Newey–Westovho odhadu zlyhať. Preto použijeme bezpečnú verziu:

pokúsime sa vypočítať Newey–Westov odhad,
ak to zlyhá, použijeme aspoň heteroskedasticky robustný odhad typu HC0.

vcov_nw <- tryCatch(
  {
    # konzervatívny malý lag, aby to nepadalo na krátkych radoch
    NeweyWest(model, prewhite = FALSE, adjust = TRUE, lag = 1)
  },
  error = function(e) {
    message("Newey-West zlyhal, používam vcovHC (HC0). Chyba bola: ", e$message)
    sandwich::vcovHC(model, type = "HC0")
  }
)

coeftest(model, vcov = vcov_nw)

## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 6.2692e+00 3.4929e+00  1.7948  0.08388 .
## Dlzka_min   3.4126e-03 1.4549e-03  2.3455  0.02660 *
## Rozpocet    6.7791e-10 9.1655e-10  0.7396  0.46591  
## Rok         7.7129e-04 1.7633e-03  0.4374  0.66529  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5. Záver

Na dátach filmov štúdia Warner Bros. Pictures sme odhadli lineárny regresný model, v ktorom IMDb hodnotenie závisí od dĺžky filmu, rozpočtu a roku uvedenia. Následne sme sa zamerali na diagnostiku autokorelácie rezíduí pomocou grafického posúdenia, ACF grafu, Durbin–Watsonovho testu a Breusch–Godfreyho testu.

V našej databáze sa výrazná autokorelácia rezíduí nepotvrdila, čo je prirodzené aj vzhľadom na to, že ide o filmy zoradené v čase, nie o klasický ekonomický časový rad s ročnou frekvenciou makroekonomických ukazovateľov. Napriek tomu sme demonštratívne ukázali, ako možno model dynamizovať pomocou Koyckovho modelu a ako použiť Newey–Westove (alebo aspoň heteroskedasticky robustné) štandardné chyby.

Multikolinearita v regresných modeloch

1. Úvod

Po autokorelácii a heteroskedasticite rezíduí je multikolinearita tretím závažným porušením predpokladov použitia metódy najmenších štvorcov.

Pri klasickej lineárnej regresii predpokladáme okrem iného, že matica vysvetľujúcich premenných X má lineárne nezávislé stĺpce. To znamená, že žiadny stĺpec (žiadna vysvetľujúca premenná) sa nedá presne vyjadriť ako lineárna kombinácia ostatných. Tento predpoklad zabezpečuje, že matica X’X je regulárna (má inverziu), a teda vieme vypočítať odhady regresných koeficientov beta = (X’X)^(-1) X’y.

Multikolinearita teda neskreslí priemer odhadov, ale zhorší ich presnosť a stabilitu.

2. Dôsledky multikolinearity

Hlavné dôsledky multikolinearity:

Nezavádza bias – odhady koeficientov sú v priemere správne.
Zvyšuje štandardné chyby koeficientov – intervaly sú širokejšie, p-hodnoty väčšie.
Odhady sú nestabilné – malé zmeny v dátach vedú k veľkým zmenám koeficientov.
Zhoršuje interpretáciu modelu – ťažko sa vysvetľuje samostatný vplyv jednej premennej.

3. Východiskový model a údaje

Budeme pracovať s databázou 250 najlepšie hodnotených filmov z IMDb (súbor data_ekonometria.csv). Zaujíma nás, ako sa hodnotenie IMDb (0–10) dá vysvetliť pomocou týchto premenných:

rok uvedenia filmu,
dĺžka filmu v minútach,
rozpočet filmu v dolároch.

Model:

IMDb = beta0 + beta1 * rok + beta2 * dlzka + beta3 * rozpocet + chyba

Najprv načítame a pripravíme údaje:

udaje <- read.csv("data_ekonometria.csv",
                  dec = ",", sep = ";", header = TRUE)

# premenovanie stĺpcov na jednoduchšie názvy
names(udaje) <- c("umiestnenie", "nazov", "rok", "zaner",
                  "imdb", "dlzka", "rezia", "spolocnost", "rozpocet_text")

# prevod rozpočtu z textu s medzerami na číselnú premennú v dolároch
udaje$rozpocet <- as.numeric(gsub(" ", "", udaje$rozpocet_text))

# vyber premenných, s ktorými budeme pracovať v regresii
udaje <- udaje[, c("imdb", "rok", "dlzka", "rozpocet")]

# jednoduchá imputácia – prípadné chýbajúce hodnoty nahradíme mediánom
column_medians <- sapply(udaje, median, na.rm = TRUE)
for (col in names(udaje)) {
  udaje[[col]][is.na(udaje[[col]])] <- column_medians[col]
}

nrow(udaje)

## [1] 250

head(udaje)

Komentár:

v databáze je 250 filmov – ide o prierezové údaje,
imdb je hodnotenie, typicky medzi 8 a 9.5,
rok sa pohybuje približne od polovice 20. storočia po súčasnosť,
dlzka je v minútach, väčšinou 90–200,
rozpocet je v dolároch, rádovo milióny až stovky miliónov.

4. Odhad základného regresného modelu

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

## 
## Call:
## lm(formula = imdb ~ rok + dlzka + 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 ***
## rok         -3.391e-04  6.267e-04  -0.541    0.589    
## dlzka        2.126e-03  4.689e-04   4.535 9.00e-06 ***
## 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 základného modelu:

R-squared je približne 0.087 – model vysvetľuje asi 8.7 % rozdielov v hodnotení IMDb. Nie je to veľa, ale vzhľadom na to, že máme len tri jednoduché vysvetľujúce premenné, je to realistické.
F-test pre celý model má veľmi malú p-hodnotu – aspoň jedna z premenných (rok, dlzka, rozpocet) má štatisticky významný vplyv na IMDb.
Koeficient pri dlzka je kladný, štatisticky významný: dlhšie filmy majú v priemere o niečo vyššie hodnotenie. Napríklad nárast dĺžky o 60 minút zvýši očakávané hodnotenie približne o 0.12 bodu.
Koeficient pri rok je blízko nule a štatisticky nevýznamny – po zohľadnení dĺžky a rozpočtu rok uvedenia filmu už nehrá veľkú rolu.
Koeficient pri rozpocet je extrémne malý (lebo pracujeme v dolároch) a tiež nie je štatisticky významny – vo vzorke top 250 filmov nevidíme jasný vzťah, že vyšší rozpočet znamená lepšie hodnotenie.
Intercept (konštanta) má hodnotu okolo 8.7 – berieme ho ako technický parameter, nie ako reálne „hodnotenie filmu z roku 0, s dĺžkou 0 a rozpočtom 0“.

Zhrnutie: najdôležitejšou premennou v tomto modeli je dĺžka filmu, rok a rozpočet už veľa navyše nevysvetľujú.

5. Korelačná matica a párové vzťahy

xvars <- udaje[, c("rok", "dlzka", "rozpocet")]
round(cor(xvars), 3)

##            rok dlzka rozpocet
## rok      1.000 0.111    0.452
## dlzka    0.111 1.000    0.110
## rozpocet 0.452 0.110    1.000

Interpretácia korelačnej matice:

rok a rozpocet: mierne kladná korelácia – novšie filmy majú často vyšší rozpočet,
rok a dlzka: korelácia je nízka – novšie filmy nie sú zásadne dlhšie či kratšie,
dlzka a rozpocet: korelácia je skôr slabá – drahšie filmy nie sú automaticky dlhšie.

Nevidíme extrémne korelácie typu 0.9, takže multikolinearita nie je na prvý pohľad zrejmá len z korelačnej matice.

pairs(xvars,
      main = "Scatterplotová matica – rok, dlzka, rozpocet")

Zo scatterplotov vidíme:

dĺžka filmu má mierne rastúci vzťah k IMDb (dlhšie filmy majú mierne vyššie hodnotenie),
vzťah rozpočtu k IMDb je veľmi rozptýlený – existujú lacnejšie, ale veľmi dobre hodnotené filmy.

6. VIF – Variance Inflation Factor

library(car)
vif(model)

##      rok    dlzka rozpocet 
## 1.262283 1.017010 1.262082

Interpretácia VIF:

VIF približne 1–2 znamená, že multikolinearita nie je silná z pohľadu danej premennej.
Pri našom modeli VIF vychádza okolo 1–1.3 – to je veľmi nízke.
To naznačuje, že žiadna premenná nie je extrémne lineárne závislá od ostatných.

Záver: VIF nám neindikuje vážny problém multikolinearity. To však ešte neznamená, že nemôže existovať problém s číselnou stabilitou kvôli rozdielnym rádom premenných – to zachytíme Condition Number.

7. Condition Number

X <- model.matrix(model)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)

condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number

## [1] 2345673

Pravidlo:

menej ako 10 – nízka multikolinearita,
10–30 – mierna,
30–100 – silná,
viac ako 100 – veľmi vážna.

V našom prípade je Condition Number veľmi vysoký (rádovo státisíce až milióny). To znamená:

matica X’X je zle podmienená,
malé zmeny údajov môžu spôsobiť výrazné zmeny v koeficientoch,
rozdielne rády premenných (rozpočet v desiatkach až stovkách miliónov, rok okolo 2000, dlzka okolo 100) spôsobujú numerickú nestabilitu.

VIF teda nie je problém, ale Condition Number ukazuje na číselnú nestabilitu modelu.

8. Riešenia multikolinearity a nestability

8.1 Vynechanie niektorej premennej

model_noRok <- lm(imdb ~ dlzka + rozpocet, data = udaje)
summary(model_noRok)

## 
## Call:
## lm(formula = imdb ~ dlzka + rozpocet, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41837 -0.16151 -0.06026  0.12843  0.98413 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.009e+00  6.216e-02 128.859  < 2e-16 ***
## dlzka       2.109e-03  4.671e-04   4.515 9.81e-06 ***
## rozpocet    2.793e-10  2.347e-10   1.190    0.235    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2214 on 247 degrees of freedom
## Multiple R-squared:  0.08607,    Adjusted R-squared:  0.07867 
## F-statistic: 11.63 on 2 and 247 DF,  p-value: 1.488e-05

model_noDlzka <- lm(imdb ~ rok + rozpocet, data = udaje)
summary(model_noDlzka)

## 
## Call:
## lm(formula = imdb ~ rok + rozpocet, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29853 -0.18305 -0.08277  0.11290  1.01132 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.565e+00  1.286e+00   6.661 1.75e-10 ***
## rok         -1.438e-04  6.495e-04  -0.221    0.825    
## rozpocet     4.228e-10  2.720e-10   1.554    0.121    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2304 on 247 degrees of freedom
## Multiple R-squared:  0.01084,    Adjusted R-squared:  0.00283 
## F-statistic: 1.353 on 2 and 247 DF,  p-value: 0.2603

model_noRozpocet <- lm(imdb ~ rok + dlzka, data = udaje)
summary(model_noRozpocet)

## 
## Call:
## lm(formula = imdb ~ rok + dlzka, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42572 -0.16175 -0.05905  0.12769  0.98021 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.963e+00  1.111e+00   7.165 8.94e-12 ***
## rok         2.476e-05  5.620e-04   0.044    0.965    
## dlzka       2.168e-03  4.685e-04   4.627 5.99e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2221 on 247 degrees of freedom
## Multiple R-squared:  0.08084,    Adjusted R-squared:  0.0734 
## F-statistic: 10.86 on 2 and 247 DF,  p-value: 3.011e-05

Typický obraz:

model bez dĺžky má výrazne nižší R-squared – dĺžka je kľúčová premenná,
model bez roku má podobný R-squared ako pôvodný model – rok veľa nepridáva,
model bez rozpočtu má tiež podobný R-squared – rozpočet tiež nepridáva veľa informácie navyše k dĺžke a roku.

Ak chceme jednoduchší model, dáva zmysel ponechať dĺžku (najdôležitejšia premenná) a pokojne vynechať rok alebo rozpočet, pretože ich prínos je malý.

8.2 Škálovanie premenných (standardizácia)

Teraz odstránime problém rozdielnych rádov premenných tým, že ich štandardizujeme (odčítame priemer a delíme smerodajnou odchýlkou).

udaje$rok_c       <- scale(udaje$rok,      center = TRUE, scale = TRUE)
udaje$dlzka_c     <- scale(udaje$dlzka,    center = TRUE, scale = TRUE)
udaje$rozpocet_c  <- scale(udaje$rozpocet, center = TRUE, scale = TRUE)

model_centered <- lm(imdb ~ rok_c + dlzka_c + rozpocet_c,
                     data = udaje)
summary(model_centered)

## 
## Call:
## lm(formula = imdb ~ rok_c + dlzka_c + rozpocet_c, 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.294400   0.014025 591.400   <2e-16 ***
## rok_c       -0.008543   0.015789  -0.541    0.589    
## dlzka_c      0.064272   0.014172   4.535    9e-06 ***
## rozpocet_c   0.020600   0.015788   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

vif(model_centered)

##      rok_c    dlzka_c rozpocet_c 
##   1.262283   1.017010   1.262082

Interpretácia:

R-squared ostáva prakticky rovnaké – vysvetlená variabilita sa nezmenila,
koeficienty teraz hovoria o zmene IMDb pri zmene premennej o 1 smerodajnú odchýlku,
VIF sú stále nízke (okolo 1–1.3).

Skontrolujeme Condition Number po škálovaní:

X <- model.matrix(model_centered)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)

condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number

## [1] 1.654121

Teraz je Condition Number blízko hodnoty 1–2 → matica je dobre podmienená, číselná stabilita je výrazne lepšia.

Nevýhoda: koeficienty sú menej intuitívne (pracujeme v jednotkách „1 smerodajná odchýlka“).

8.3 Transformácia rozpočtu na milióny dolárov

Ak chceme zachovať interpretáciu v reálnych jednotkách, môžeme upraviť len rozpočet: namiesto dolárov použijeme milióny dolárov.

udaje$rozpocet_mil <- udaje$rozpocet / 1e6
head(udaje)

model_rozp_mil <- lm(imdb ~ rok + dlzka + rozpocet_mil,
                     data = udaje)
summary(model_rozp_mil)

## 
## Call:
## lm(formula = imdb ~ rok + dlzka + rozpocet_mil, 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.6783808  1.2379574   7.010 2.28e-11 ***
## rok          -0.0003391  0.0006267  -0.541    0.589    
## dlzka         0.0021265  0.0004689   4.535 9.00e-06 ***
## rozpocet_mil  0.0003424  0.0002624   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

vif(model_rozp_mil)

##          rok        dlzka rozpocet_mil 
##     1.262283     1.017010     1.262082

Interpretácia:

koeficient pri rozpocet_mil hovorí, o koľko bodov sa zmení IMDb pri zvýšení rozpočtu o 1 milión dolárov,
VIF sú podobné ako predtým (multikolinearita medzi rokom a dĺžkou môže mierne pretrvávať),
Condition Number je nižší ako pri pôvodnom rozpočte v dolároch, ale stále vyšší ako po úplnom škálovaní všetkých premenných.

Je to kompromis: stále lepšie čísla ako pri pôvodnom rozpočte a koeficienty sú dobre interpretovateľné.

library(dplyr)
udaje <- udaje %>%
  dplyr::select(imdb, rok, dlzka, rozpocet, rozpocet_mil,
                rok_c, dlzka_c, rozpocet_c)

8.4 PCA (voliteľné)

PCA (Principal Component Analysis) vytvorí nové premenné – hlavné komponenty – ako lineárne kombinácie pôvodných premenných. Môžeme napríklad spojiť rok a rozpocet do jednej komponenty (PC1).

X_pca <- scale(udaje[, c("rok", "rozpocet")],
               center = TRUE, scale = TRUE)
pca_res <- prcomp(X_pca)

summary(pca_res)

## Importance of components:
##                           PC1    PC2
## Standard deviation     1.2049 0.7405
## Proportion of Variance 0.7258 0.2742
## Cumulative Proportion  0.7258 1.0000

udaje$PC1 <- pca_res$x[,1]

model_pca <- lm(imdb ~ dlzka + PC1, data = udaje)
summary(model_pca)

## 
## Call:
## lm(formula = imdb ~ dlzka + PC1, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42087 -0.16565 -0.05762  0.12685  0.98014 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.0172014  0.0627437 127.777  < 2e-16 ***
## dlzka        0.0021261  0.0004691   4.533 9.08e-06 ***
## PC1         -0.0085266  0.0117667  -0.725    0.469    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2218 on 247 degrees of freedom
## Multiple R-squared:  0.08278,    Adjusted R-squared:  0.07536 
## F-statistic: 11.15 on 2 and 247 DF,  p-value: 2.319e-05

vif(model_pca)

##    dlzka      PC1 
## 1.017009 1.017009

Interpretácia:

PC1 zachytáva väčšinu variability dvojice (rok, rozpocet) – je to akýsi index „modernosti a finančnej náročnosti“ filmu,
dĺžka ostáva významným vysvetľujúcim faktorom,
PC1 nemusí byť významný – kombinácia roku a rozpočtu stále nehrá veľkú rolu pri vysvetľovaní IMDb,
VIF sú nízke, pretože PC1 a dlzka sú nekorelované (ortogonálne).

Nevýhoda: PC1 je ťažšie interpretovateľná ako rok alebo rozpocet.

8.5 Hrebeňová regresia (ridge regression – voliteľné)

Hrebeňová regresia je modifikácia OLS, pri ktorej pridávame do matice X’X tzv. penalizáciu lambda * I. Odhady koeficientov sú mierne skreslené, ale stabilnejšie pri multikolinearite.

library(MASS)

X <- as.matrix(udaje[, c("rok", "dlzka", "rozpocet")])
y <- udaje$imdb

ridge_fit <- lm.ridge(y ~ X, lambda = seq(0, 10, 1))
ridge_fit

##                      Xrok      Xdlzka    Xrozpocet
##  0 8.678381 -0.0003390826 0.002126497 3.424405e-10
##  1 8.670991 -0.0003347735 0.002118001 3.407312e-10
##  2 8.663724 -0.0003305316 0.002109575 3.390442e-10
##  3 8.656578 -0.0003263555 0.002101216 3.373792e-10
##  4 8.649551 -0.0003222438 0.002092924 3.357357e-10
##  5 8.642640 -0.0003181951 0.002084698 3.341132e-10
##  6 8.635843 -0.0003142083 0.002076538 3.325112e-10
##  7 8.629156 -0.0003102819 0.002068443 3.309294e-10
##  8 8.622578 -0.0003064147 0.002060411 3.293672e-10
##  9 8.616106 -0.0003026056 0.002052443 3.278244e-10
## 10 8.609738 -0.0002988534 0.002044537 3.263005e-10

Interpretácia:

pri lambda = 0 dostaneme klasickú OLS regresiu,
ako lambda rastie, koeficienty sa zmenšujú smerom k nule a model je stabilnejší,
dôležitejšie premenné (ako dlzka) si držia nenulové koeficienty,
menej dôležité (rok, rozpocet) sa pri väčších lambda blížia k nule.

Ridge regresia je už pokročilejší nástroj, ale pekne ilustruje, ako sa dá technicky bojovať s multikolinearitou a nestabilitou odhadov.

9. Zhrnutie

Na databáze 250 najlepšie hodnotených filmov z IMDb sme skúmali vplyv roku, dĺžky a rozpočtu na hodnotenie IMDb.
Základný model ukázal, že najdôležitejšou premennou je dĺžka filmu – dlhšie filmy majú mierne vyššie hodnotenie.
Rok a rozpočet nie sú po zohľadnení dĺžky štatisticky významné.
VIF neukazuje vážnu multikolinearitu, ale Condition Number poukazuje na číselnú nestabilitu modelu kvôli rozdielnym rádom premenných.
Riešenia:
- vynechanie málo dôležitých premenných (rok, rozpocet),
- škálovanie premenných (výrazne zlepší Condition Number),
- transformácia rozpočtu na milióny dolárov,
- PCA a hrebeňová regresia (pokročilejšie metódy).

Celkovo vidíme, že multikolinearita a číselná nestabilita nie sú len teoretický problém, ale reálne ovplyvňujú interpretáciu regresného modelu aj pri praktických dátach, ako sú filmy z IMDb.

	Umiestnenie	klaster
159	159	1
207	207	2
179	179	2
14	14	3
195	195	2
170	170	2
50	50	1
118	118	2
43	43	1
229	229	2
247	247	2
243	243	2
153	153	2
90	90	2
91	91	2
197	197	2
236	236	2
185	185	2
92	92	2
137	137	2

	Umiestnenie	klaster
159	159	1
207	207	2
179	179	2
14	14	3
195	195	2
170	170	2
50	50	1
118	118	2
43	43	1
229	229	2
247	247	2
243	243	2
153	153	2
90	90	2
91	91	2
197	197	2
236	236	2
185	185	2
92	92	2
137	137	2

Finálne zadanie

Diana Ďuriančiková

December 2025

Obsah

Práca s databázou - import údajov, grafy, štatistiky

Úvod

Nastavenie

Načítanie a čistenie dát

Výber dvoch skupín

Grafy

Stĺpcový graf: priemerné IMDb podľa dekády a skupiny

Čiarový graf: medián rozpočtu podľa roku a skupiny

Základné štatistiky

Testovanie hypotéz

t-test: rozdiel IMDb medzi 2 skupinami v najnovšej dekáde

ANOVA: IMDb podľa žánru (top 5 žánrov)

Lineárna regresia

Lineárna regresia

Úvod do problému a hypotézy

Lineárna regresia

Čo objekt lm() obsahuje

Diagnostické grafy

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

Q–Q plot

Scale–Location

Residuals vs. Leverage

Testy normality a odľahlostí

Zhrnutie (Conclusion)

Heteroskedasticita

Testovanie prítomnosti heteroskedasticity

Špecifikácia modelu – hodnotenie filmov

1. Základná regresia

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

3. Grafická analýza

3.1 Graf Residuals vs. Fitted

3.2 C+R grafy (Component + Residual plots)

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

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

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

6.1 Zlom v autonómnom člene

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

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

Zhluková analýza filmov

Úvod

Škálovanie premenných a boxploty

Korelačná matica premenných

Matica vzdialeností

Princíp hierarchického zhlukovania (Wardova metóda)

Deskriptívne štatistiky – rozklad variability

Centroidy – priemerné hodnoty premenných v klastroch

Záver

Autokorelácia rezíduí

1. Úvod a údaje

1.1 Úvod do problému, stanovenie hypotéz

1.2 Príprava databázy, úprava údajov

2. Lineárna regresia v základnom tvare

3. Autokorelácia rezíduí

3.1 Čo je autokorelácia rezíduí?

3.2 Dôsledky autokorelácie rezíduí

3.3 Detekcia autokorelácie

Grafická informácia

ACF graf (Autocorrelation Function)

Durbin–Watsonov test

Breusch–Godfreyov test (BG test)

4. Dynamizácia modelu – Koyckov model

4.1 Odhad Koyckovho modelu v R

4.2 Newey–West robustné (alebo aspoň robustné) štandardné chyby

5. Záver

Multikolinearita v regresných modeloch

1. Úvod

2. Dôsledky multikolinearity

3. Východiskový model a údaje

4. Odhad základného regresného modelu

5. Korelačná matica a párové vzťahy

6. VIF – Variance Inflation Factor

7. Condition Number

8. Riešenia multikolinearity a nestability

8.1 Vynechanie niektorej premennej

8.2 Škálovanie premenných (standardizácia)

8.3 Transformácia rozpočtu na milióny dolárov

Čo objekt `lm()` obsahuje

	Umiestnenie	klaster
159	159	1
207	207	2
179	179	2
14	14	3
195	195	2
170	170	2
50	50	1
118	118	2
43	43	1
229	229	2
247	247	2
243	243	2
153	153	2
90	90	2
91	91	2
197	197	2
236	236	2
185	185	2
92	92	2
137	137	2