library(broom)
library(knitr)Ekonometria v R – cvičenie 5
Základná diagnostika modelu
1 Cieľ cvičenia
Pokračujeme v analýze hier vydaných v roku 2023, ktorú sme začali vo 4. cvičení. Používame rovnaké údaje a podobný postup analýzy ako minule, no model zjednodušíme. Namiesto údajov WHO z Gazdovho podkladu skúmame súvislosť medzi cenou hry a pozitívnymi a negatívnymi hodnoteniami. Žánre do modelu nezaraďujeme.
Na konci cvičenia budeme vedieť:
- Uložiť a obnoviť pracovné prostredie v R (workspace).
- Interpretovať t-testy koeficientov, celkový F-test a koeficient determinácie \(R^2\).
- Otestovať normalitu rezíduí a nájsť potenciálne odľahlé pozorovania.
- Prečítať štyri diagnostické grafy pomocou
plot(model). - Upraviť model, porovnať ho s pôvodným a zdôvodniť výsledok.
PoznámkaMetodická poznámka
Počet pozitívnych hodnotení nie je počet predajov ani priamo meraná kvalita hry. Ide o opisné asociácie, nie dôkaz, že zvýšenie ceny spôsobí viac hodnotení. Naším cieľom je najmä naučiť sa regresné testy a diagnostiku.
2 Príprava databázy
2.1 Načítanie údajov a balíkov
Súbor games_upravene.csv musí byť v rovnakom priečinku ako tento .qmd dokument. Nie je potrebné používať setwd() ani podpriečinok udaje.
cesta <- "games_upravene.csv"
if (!file.exists(cesta)) {
stop("Súbor games_upravene.csv vložte vedľa dokumentu QMD.")
}
udaje <- read.csv(cesta, sep = ",", dec = ".",
na.strings = c("", "NA"), stringsAsFactors = FALSE)
names(udaje) <- make.names(trimws(names(udaje)), unique = TRUE)
# Rovnako ako v predchádzajúcom cvičení vyberieme rok 2023.
udaje.2023 <- udaje[!is.na(udaje$Rok) & udaje$Rok == 2023,
c("ID", "Nazov", "Cena", "Pozitivne", "Negativne", "Rok")]
# Na odhad použijeme všetky hry s platnými nezápornými údajmi.
potrebne <- c("Cena", "Pozitivne", "Negativne")
udaje.2023 <- udaje.2023[complete.cases(udaje.2023[potrebne]) &
udaje.2023$Cena >= 0 &
udaje.2023$Pozitivne >= 0 &
udaje.2023$Negativne >= 0, ]
rownames(udaje.2023) <- NULL
# Logaritmus z (1 + počet): funguje aj pri nulovom počte hodnotení.
udaje.2023$LogPozitivne <- log1p(udaje.2023$Pozitivne)
udaje.2023$LogNegativne <- log1p(udaje.2023$Negativne)
# Iba na ukážku vyberieme rozmanité, skutočné hry s nenulovými hodnotami.
# Toto NIE JE filter pre regresiu: model zostáva odhadnutý na celej vzorke.
ukazkove_nazvy <- c("Summon", "Sea of Survivors", "Storyteller",
"Evil Dead: The Game", "System Shock",
"The Great War: Western Front™")
ukazka <- udaje.2023[match(ukazkove_nazvy, udaje.2023$Nazov),
c("Nazov", "Cena", "Pozitivne", "Negativne")]
ukazka <- ukazka[complete.cases(ukazka) & ukazka$Cena > 0 &
ukazka$Pozitivne > 0 & ukazka$Negativne > 0, ]
kable(ukazka, digits = 2, row.names = FALSE,
col.names = c("Hra", "Cena", "Pozitívne", "Negatívne"),
caption = "Ilustračný výber hier z roku 2023 (nenulové hodnoty)")| Hra | Cena | Pozitívne | Negatívne |
|---|---|---|---|
| Summon | 0.99 | 1140 | 42 |
| Sea of Survivors | 4.99 | 148 | 30 |
| Storyteller | 13.49 | 524 | 201 |
| Evil Dead: The Game | 19.79 | 221 | 99 |
| System Shock | 39.99 | 2278 | 174 |
| The Great War: Western Front™ | 34.99 | 455 | 123 |
Vo vzorke pracujeme s 15542 hrami. Premenná LogPozitivne je vysvetľovaná; Cena a LogNegativne sú vysvetľujúce premenné. Ukážková tabuľka obsahuje vybrané hry, regresia však využíva celú pripravenú vzorku, vrátane platných nulových hodnôt.
2.2 Uloženie všetkých objektov, čistenie a obnovenie workspace
Workspace je súbor objektov, ktoré máme vytvorené v R. Príkaz save.image() ho uloží; príkaz load() ho znova načíta.
# Uložíme aktuálny workspace do súboru.
save.image("pamat_cvicenie5.RData")
# Vypíšeme niektoré objekty, ktoré v prostredí existujú.
head(ls(), 12)[1] "cesta" "potrebne" "udaje" "udaje.2023"
[5] "ukazka" "ukazkove_nazvy"Mazanie objektov si ukážeme bezpečne na dvoch pomocných premenných. Neodstraňujeme dátovú tabuľku, aby mohli pokračovať ďalšie výpočty.
aaa <- 6
bbb <- "text"
rm(aaa, bbb)
ls()[1:min(12, length(ls()))][1] "cesta" "potrebne" "udaje" "udaje.2023"
[5] "ukazka" "ukazkove_nazvy"Úplné vyčistenie prostredia a obnovenie pamäte sa robí takto. Tento blok je ukážka na samostatné spustenie v konzole, pri renderovaní ho nevykonávame, aby neprerušil nadväzujúci dokument.
rm(list = ls())
load("pamat_cvicenie5.RData")
ls()
TipPraktický postup
Na konci práce môžeš použiť save.image("pamat_cvicenie5.RData"). Na ďalšom cvičení obnovíš uložené objekty príkazom load("pamat_cvicenie5.RData"). V Quarte je však spoľahlivejšie ponechať aj kód na načítanie CSV, aby sa dokument dal zostaviť úplne od začiatku.
3 Lineárna regresia
Použijeme zjednodušený viacnásobný model, ktorý nadväzuje na postup zo 4. týždňa, ale nezahŕňa žáner:
\[ \operatorname{LogPozitivne}_i = \beta_0 + \beta_1 \operatorname{Cena}_i + \beta_2 \operatorname{LogNegativne}_i + u_i. \]
model <- lm(LogPozitivne ~ Cena + LogNegativne,
data = udaje.2023)
summary(model)
Call:
lm(formula = LogPozitivne ~ Cena + LogNegativne, data = udaje.2023)
Residuals:
Min 1Q Median 3Q Max
-4.2574 -0.4534 -0.4408 0.2499 5.0054
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4407961 0.0083493 52.795 < 2e-16 ***
Cena 0.0025203 0.0004442 5.674 1.42e-08 ***
LogNegativne 1.2217800 0.0072157 169.322 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9224 on 15539 degrees of freedom
Multiple R-squared: 0.6666, Adjusted R-squared: 0.6665
F-statistic: 1.553e+04 on 2 and 15539 DF, p-value: < 2.2e-163.1 Odhady koeficientov a t-testy
Funkcia tidy() prehľadne zobrazí odhad koeficientu, jeho štandardnú chybu, t-štatistiku a p-hodnotu.
tab_koef <- tidy(model)
# Veľmi malé p-hodnoty zobrazíme ako nerovnosť, nie ako nulu.
tab_koef_vypis <- tab_koef
tab_koef_vypis$p.value <- format.pval(tab_koef$p.value, digits = 3, eps = 0.001)
kable(tab_koef_vypis, digits = 5,
col.names = c("Premenná", "Odhad", "Štandardná chyba", "t-štatistika", "p-hodnota"),
caption = "Výsledky t-testov jednotlivých koeficientov")| Premenná | Odhad | Štandardná chyba | t-štatistika | p-hodnota |
|---|---|---|---|---|
| (Intercept) | 0.44080 | 0.00835 | 52.79457 | <0.001 |
| Cena | 0.00252 | 0.00044 | 5.67429 | <0.001 |
| LogNegativne | 1.22178 | 0.00722 | 169.32197 | <0.001 |
Pre jednotlivý koeficient testujeme \(H_0: \beta_j=0\) proti \(H_1: \beta_j\neq0\). Ak \(p<0{,}05\), zamietame \(H_0\) na 5 % hladine významnosti. t-štatistika je odhad koeficientu / štandardná chyba.
koef_cena <- tab_koef[tab_koef$term == "Cena", ]
cat("Koeficient ceny:", round(koef_cena$estimate, 6), "\n")Koeficient ceny: 0.00252 cat("t-štatistika:", round(koef_cena$statistic, 3), "\n")t-štatistika: 5.674 cat("p-hodnota:", format.pval(koef_cena$p.value, digits = 3, eps = 0.001), "\n")p-hodnota: <0.001 Interpretácia nášho modelu: Pri nezmenenom počte negatívnych hodnotení sa zvýšenie ceny o 1 cenovú jednotku spája so zmenou LogPozitivne približne o 0.0025. Ide o zmenu logaritmickej transformácie, nie priamo o 0.0025 ďalších hodnotení. Súvislosť je v tomto modeli štatisticky významná na 5 % hladine.
4 Kvalita vyrovnania modelu
4.1 Koeficient determinácie, F-test a ďalšie charakteristiky
suhrn <- glance(model)
vyber <- data.frame(
Ukazovatel = c("Počet pozorovaní", "R²", "Upravené R²", "Reziduálna chyba (sigma)",
"F-štatistika", "p-hodnota F-testu", "AIC", "BIC"),
Hodnota = c(format(nobs(model), big.mark = " "),
round(suhrn$r.squared, 5), round(suhrn$adj.r.squared, 5),
round(suhrn$sigma, 5), round(suhrn$statistic, 3),
format.pval(suhrn$p.value, digits = 3, eps = 0.001),
round(suhrn$AIC, 2), round(suhrn$BIC, 2))
)
kable(vyber, row.names = FALSE, caption = "Základné charakteristiky pôvodného modelu")| Ukazovatel | Hodnota |
|---|---|
| Počet pozorovaní | 15 542 |
| R² | 0.66658 |
| Upravené R² | 0.66653 |
| Reziduálna chyba (sigma) | 0.92241 |
| F-štatistika | 15532.775 |
| p-hodnota F-testu | <0.001 |
| AIC | 41600.64 |
| BIC | 41631.25 |
Koeficient determinácie \(R^2\) ukazuje, akú časť variability vysvetľovanej premennej dokáže regresia zachytiť. Náš model má R² = 0.6666, teda vysvetľuje približne 66.7 % variability logaritmovaného počtu pozitívnych hodnotení (nie rovnaký podiel variability pôvodných počtov).
Upravené \(R^2\) zohľadňuje počet vysvetľujúcich premenných, a preto je vhodné pri porovnávaní modelov s rôznym počtom prediktorov. Sigma opisuje typickú veľkosť rezíduí na logaritmickej škále.
Celkový F-test skúma hypotézy:
\[H_0:\beta_1=\beta_2=0, \qquad H_1:\text{aspoň jeden koeficient nie je nulový}. \]
Pre náš model je F = 15 532.78 a p-hodnota < 0,001. Na hladine 5 % preto zamietame nulovú hypotézu a model je ako celok štatisticky významný.
PoznámkaČo neznamená významný F-test
Významný F-test nezaručuje, že model správne zachytáva všetky vzťahy, že sú rezíduá normálne alebo že existuje kauzálny účinok. Rovnako ani vyššie \(R^2\) neznamená, že sú splnené predpoklady regresie. To preverujeme diagnostikou.
5 Grafická analýza rezíduí
Rezíduum je rozdiel medzi pozorovanou a predikovanou hodnotou: \(e_i=y_i-\hat y_i\).
par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(model, which = 1:4, id.n = 3)
par(mfrow = c(1, 1))5.1 1. Residuals vs Fitted
Graf ukazuje rezíduá proti predikovaným hodnotám. Ideálne tvoria približne náhodný pás okolo nuly. Systematická krivka upozorňuje na možnú nelinearitu a meniaci sa rozptyl na heteroskedasticitu.
5.2 2. Normal Q-Q
Ak sú rezíduá približne normálne rozdelené, body ležia blízko referenčnej čiary. Odchýlky najmä na koncoch znamenajú, že rozdelenie rezíduí môže mať ťažšie chvosty alebo odľahlé hodnoty.
5.3 3. Scale-Location
Zobrazuje odmocninu absolútnej hodnoty štandardizovaných rezíduí proti predikciám. Približne vodorovný pás podporuje predpoklad podobného rozptylu chýb. Stúpajúci alebo klesajúci trend môže naznačovať heteroskedasticitu.
5.4 4. Cookova vzdialenosť
Tento graf zvýrazňuje pozorovania, ktoré môžu výraznejšie ovplyvňovať odhad regresie. V tomto dokumente používame štvrtý štandardný graf plot.lm, teda Cookovu vzdialenosť. Graf Residuals vs Leverage (spomínaný v Gazdovom podklade) si vykreslíme zvlášť:
plot(model, which = 5, id.n = 3)
Vysoká leverage znamená nezvyčajné kombinácie vysvetľujúcich premenných. Veľká Cookova vzdialenosť môže indikovať, že konkrétne pozorovanie výrazne mení odhady koeficientov.
6 Test normality rezíduí
6.1 Jarqueov–Berov test
Test vychádza zo šikmosti (\(S\)) a špicatosti (\(K\)):
\[JB=\frac{n}{6}\left(S^2+\frac{(K-3)^2}{4}\right).\]
Testujeme \(H_0\): rezíduá majú normálne rozdelenie; \(H_1\): rezíduá nemajú normálne rozdelenie. Na 5 % hladine zamietame \(H_0\), ak \(p<0{,}05\).
Výpočet je priamo v základnom R (bez potreby doplnkového balíka tseries):
jarque_bera <- function(x) {
x <- as.numeric(x)
n <- length(x)
s <- mean((x - mean(x))^3) / mean((x - mean(x))^2)^(3 / 2)
k <- mean((x - mean(x))^4) / mean((x - mean(x))^2)^2
JB <- n / 6 * (s^2 + (k - 3)^2 / 4)
data.frame(JB = JB, p_hodnota = pchisq(JB, df = 2, lower.tail = FALSE),
Sikmost = s, Spicatost = k)
}
test_normality <- jarque_bera(residuals(model))
test_normality_vypis <- test_normality
test_normality_vypis$p_hodnota <- format.pval(test_normality$p_hodnota,
digits = 3, eps = 0.001)
kable(test_normality_vypis, digits = 5,
col.names = c("JB štatistika", "p-hodnota", "Šikmosť", "Špicatosť"),
caption = "Jarqueov–Berov test rezíduí pôvodného modelu")| JB štatistika | p-hodnota | Šikmosť | Špicatosť |
|---|---|---|---|
| 7827.58 | <0.001 | 1.3439 | 5.20524 |
Výsledok: JB = 7827.58, p-hodnota < 0,001. Nulovú hypotézu o normalite rezíduí zamietame. Pri veľkej vzorke môžu byť aj menšie odchýlky od normality štatisticky významné, preto test interpretujeme spolu s Q-Q grafom. Odchýlky od predpokladov modelu treba zohľadniť pri interpretácii t-testov a F-testu.
7 Test odľahlých pozorovaní
Na podobnom princípe ako car::outlierTest() použijeme externe študentizované rezíduá a Bonferroniho korekciu p-hodnôt. Výpočet je opäť v základnom R.
\[t_i = \frac{e_i}{\hat\sigma_{(i)}\sqrt{1-h_{ii}}}, \qquad p_{i,Bonf}=\min(1,n p_i).\]
stud_res <- rstudent(model)
# rstudent() používa reziduálne stupne voľnosti mínus 1.
p_raw <- 2 * pt(abs(stud_res), df = df.residual(model) - 1,
lower.tail = FALSE)
p_bonf <- p.adjust(p_raw, method = "bonferroni")
odlahle <- data.frame(
Hra = udaje.2023$Nazov,
Studentizovane_reziduum = stud_res,
p_hodnota = p_raw,
p_Bonferroni = p_bonf,
Cookova_vzdialenost = cooks.distance(model)
)
# Ukážeme päť hier s najväčším absolútnym študentizovaným rezíduom.
poradie <- order(abs(odlahle$Studentizovane_reziduum), decreasing = TRUE)
odlahle_vypis <- head(odlahle[poradie, ], 5)
odlahle_vypis$p_hodnota <- format.pval(odlahle_vypis$p_hodnota,
digits = 3, eps = 0.001)
odlahle_vypis$p_Bonferroni <- format.pval(odlahle_vypis$p_Bonferroni,
digits = 3, eps = 0.001)
names(odlahle_vypis) <- c("Hra", "Rezíduum", "p-hodnota", "p-Bonf.", "Cookova vzd.")
kable(odlahle_vypis, digits = 5, row.names = FALSE,
caption = "Päť najvýraznejších odľahlých pozorovaní")| Hra | Rezíduum | p-hodnota | p-Bonf. | Cookova vzd. |
|---|---|---|---|---|
| Bzzzt | 5.43171 | <0.001 | < 0.001 | 0.00081 |
| Cattails: Wildwood Story | 5.28826 | <0.001 | 0.00195 | 0.00106 |
| Poosh XL | 4.94275 | <0.001 | 0.01210 | 0.00062 |
| MareQuest: An Interactive Tail | 4.76367 | <0.001 | 0.02981 | 0.00056 |
| Larcin Lazer | 4.72675 | <0.001 | 0.03577 | 0.00057 |
pocet_odlahlych <- sum(odlahle$p_Bonferroni < 0.05)Pri \(p_{Bonf}<0{,}05\) považujeme rezíduum za štatisticky významne odľahlé. V pôvodnom modeli sme tak našli 5 pozorovaní. Hry však automaticky neodstraňujeme: môžu to byť skutočné úspešné či neobvyklé tituly, nielen chyby databázy. Navyše presnosť testu závisí od predpokladov modelu vrátane rozdelenia chýb.
8 Úprava modelu
Gazdov podklad skúša logaritmickú transformáciu pri hľadaní vhodnejšieho modelu. Keďže počet hodnotení sme už logaritmovali v 4. týždni, teraz namiesto pôvodnej ceny použijeme log1p(Cena). Toto je vhodné aj pre bezplatné hry s cenou nula.
\[\operatorname{LogPozitivne}_i=\beta_0 + \beta_1\ln(1+\operatorname{Cena}_i) +\beta_2\operatorname{LogNegativne}_i+u_i.\]
model2 <- lm(LogPozitivne ~ log1p(Cena) + LogNegativne,
data = udaje.2023)
tab_model2 <- tidy(model2)
tab_model2$p.value <- format.pval(tab_model2$p.value, digits = 3, eps = 0.001)
kable(tab_model2, digits = 5,
col.names = c("Premenná", "Odhad", "Štandardná chyba", "t-štatistika", "p-hodnota"),
caption = "Koeficienty upraveného modelu")| Premenná | Odhad | Štandardná chyba | t-štatistika | p-hodnota |
|---|---|---|---|---|
| (Intercept) | 0.14396 | 0.01215 | 11.84409 | <0.001 |
| log1p(Cena) | 0.23158 | 0.00700 | 33.06888 | <0.001 |
| LogNegativne | 1.13955 | 0.00732 | 155.68740 | <0.001 |
8.1 Porovnanie modelov
porovnanie <- data.frame(
Model = c("Pôvodný (Cena)", "Upravený (log1p(Cena))"),
R2 = c(glance(model)$r.squared, glance(model2)$r.squared),
R2_upravene = c(glance(model)$adj.r.squared, glance(model2)$adj.r.squared),
AIC = c(AIC(model), AIC(model2)),
BIC = c(BIC(model), BIC(model2)),
JB = c(jarque_bera(residuals(model))$JB,
jarque_bera(residuals(model2))$JB)
)
kable(porovnanie, digits = 4, row.names = FALSE,
caption = "Porovnanie pôvodného a upraveného modelu")| Model | R2 | R2_upravene | AIC | BIC | JB |
|---|---|---|---|---|---|
| Pôvodný (Cena) | 0.6666 | 0.6665 | 41600.64 | 41631.25 | 7827.580 |
| Upravený (log1p(Cena)) | 0.6879 | 0.6878 | 40575.82 | 40606.43 | 6335.384 |
Vyhodnotenie: Pôvodný model má R² = 0.6666 a AIC 41600.64. Upravený model má R² = 0.6879 a AIC 40575.82. Nižšie AIC pri rovnakých údajoch uprednostňuje upravený model. Takéto porovnanie však nedokazuje, že všetky predpoklady modelu platia.
8.2 Diagnostické grafy upraveného modelu
par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(model2, which = 1:4, id.n = 3)
par(mfrow = c(1, 1))8.3 Normalita a odľahlé hodnoty upraveného modelu
test2 <- jarque_bera(residuals(model2))
stud2 <- rstudent(model2)
p2 <- 2 * pt(abs(stud2), df = df.residual(model2) - 1,
lower.tail = FALSE)
p2_bonf <- p.adjust(p2, method = "bonferroni")
porovnanie_testov <- data.frame(
Model = c("Pôvodný", "Upravený"),
JB = c(test_normality$JB, test2$JB),
p_JB = c(test_normality$p_hodnota, test2$p_hodnota),
Pocet_odlahlych = c(pocet_odlahlych, sum(p2_bonf < 0.05))
)
porovnanie_testov$p_JB <- format.pval(porovnanie_testov$p_JB,
digits = 3, eps = 0.001)
kable(porovnanie_testov, digits = 5, row.names = FALSE,
caption = "Normalita a počet Bonferroniho významných odľahlých pozorovaní")| Model | JB | p_JB | Pocet_odlahlych |
|---|---|---|---|
| Pôvodný | 7827.580 | <0.001 | 5 |
| Upravený | 6335.384 | <0.001 | 6 |
Interpretácia: Aj pri upravenom modeli test zamieta normalitu rezíduí. Logaritmická transformácia ceny znížila AIC, no treba osobitne posúdiť diagnostické grafy a odľahlé prípady. Odľahlé hry nemožno bez vysvetlenia vymazať.
9 Malé úlohy – vypracované
9.1 Úloha 1: Interpretuj t-test ceny
Riešenie: V tabuľke tidy(model) vyhľadáme premennú Cena. Jej odhad je 0.00252 a p-hodnota <0.001. Na 5 % hladine je koeficient štatisticky významný. Výsledok nepreukazuje príčinnosť.
9.2 Úloha 2: Interpretuj F-test a R²
Riešenie: F-test skúma spoločnú významnosť dvoch vysvetľujúcich premenných. Nulovú hypotézu zamietame. Model vysvetľuje približne 66.7 % variability LogPozitivne.
9.3 Úloha 3: Over normalitu a nájdi odľahlé prípady
Riešenie: Jarqueov–Berov test dáva p-hodnotu < 0,001. Po Bonferroniho korekcii má 5 hier p-hodnotu menšiu než 0,05. Odľahlosť nie je automatickým dôvodom na odstránenie hry.
9.4 Úloha 4: Ktorý model je lepší podľa AIC?
Riešenie: Nižšie AIC dosahuje upravený model s log1p(Cena). Porovnávame modely odhadnuté na tej istej vzorke a s rovnakou vysvetľovanou premennou. Treba však pozrieť aj grafy rezíduí.
9.5 Úloha 5: Ulož a znova otvor workspace
Riešenie: save.image("pamat_cvicenie5.RData") uloží aktuálne objekty. Následne ich možno po vyčistení konzoly načítať príkazom load("pamat_cvicenie5.RData"). Pri automatickom renderovaní ponecháme bezpečný ukážkový blok s eval: false.
10 Záver
V 5. cvičení sme na rovnakých dátach ako minule odhadli viacnásobnú regresiu, vyhodnotili t-testy, F-test a \(R^2\), skontrolovali diagnostické grafy, vykonali Jarqueov–Berov test a identifikovali odľahlé pozorovania. Skúsili sme tiež transformáciu ceny a modely porovnali pomocou AIC.
Z výsledkov nemôžeme vyvodzovať príčinné vzťahy; navyše pozitívne a negatívne hodnotenia sú úzko prepojené. Pri ďalšom výskume by bolo vhodné skúmať aj ďalšie premenné a predpoklady regresie.
# Uložíme všetky objekty vrátane modelov a diagnostických výsledkov.
save.image("pamat_cvicenie5.RData")