cvičenie 5

Základná diagnostika modelu

\[ \ln Cena_i = \beta_0 + \beta_1 \ln Plocha_i + \beta_2 VekDiela_i + \beta_3 RokNarodenia_i + u_i. \]

V Ccvičení 4 sme model len odhadli a popísali. Teraz sa pozrieme na to, či sa mu dá veriť:

• uloženie a obnovenie pracovnej pamäte (workspace),
• t-testy jednotlivých koeficientov, F-test celého modelu a kvalita vyrovnania (\(R^2\)),
• detailná grafická analýza rezíduí (plot(model)),
• test normality rezíduí a test odľahlých pozorovaní.

library(broom)
library(knitr)
library(kableExtra)
library(dplyr)
library(car)

# Číslo so slovenskou desatinnou čiarkou pre texty v dokumente.
fmt <- function(x, d = 2) format(round(x, d), nsmall = d, decimal.mark = ",")

1 Príprava databázy

Príprava údajov je rovnaká ako v cvičení 4

cesta <- "udaje/transactions.csv"
if (!file.exists(cesta)) cesta <- "transactions.csv"
if (!file.exists(cesta)) {
  stop("Vložte transactions.csv do podpriečinka udaje alebo vedľa QMD dokumentu.")
}

riadky <- readLines(cesta, warn = FALSE)
riadky <- iconv(riadky, from = "latin1", to = "ASCII", sub = "")
riadky <- riadky[-1]
riadky <- riadky[nzchar(trimws(riadky))]

ako_cislo <- function(x) {
  if (length(x) == 1 && !is.na(x) && grepl("^-?[0-9]+(\\.[0-9]+)?$", x)) {
    as.numeric(x)
  } else {
    NA_real_
  }
}

rozparsuj_riadok <- function(riadok) {
  vysledok <- c(real_cena = NA_real_, rok_aukcie = NA_real_,
                rok_narodenia = NA_real_, rok_diela = NA_real_,
                plocha = NA_real_)
  p <- trimws(strsplit(riadok, ",", fixed = TRUE)[[1]])
  n <- length(p)

  L <- which(startsWith(p, "Lot") & seq_len(n) >= 12)[1]
  K <- which(startsWith(p, "<div class=image"))[1]
  if (is.na(L) || is.na(K)) return(vysledok)

  vysledok["real_cena"] <- ako_cislo(p[L + 3])
  vysledok["plocha"]    <- ako_cislo(p[K + 2])

  for (i in (L + 5):min(n, L + 40)) {
    if (grepl("^[0-9]{4}$", p[i]) && p[i] == p[i - 1]) {
      vysledok["rok_aukcie"] <- as.numeric(p[i])
      break
    }
  }

  for (i in 2:(L - 3)) {
    if (grepl("^[0-9]{4}$", p[i]) && grepl("^([0-9]{4})?$", p[i + 1]) &&
        grepl("[A-Za-z]", p[i - 1])) {
      vysledok["rok_narodenia"] <- as.numeric(p[i])
      vysledok["rok_diela"]     <- ako_cislo(p[i + 1])
      break
    }
  }
  vysledok
}

matica <- vapply(riadky, rozparsuj_riadok, numeric(5), USE.NAMES = FALSE)
transakcie <- as.data.frame(t(matica))
names(transakcie) <- c("real_cena", "rok_aukcie", "rok_narodenia",
                       "rok_diela", "plocha")

udaje.2014 <- transakcie[!is.na(transakcie$rok_aukcie) &
                           transakcie$rok_aukcie == 2014, ]

udaje.2014$real_cena[!is.na(udaje.2014$real_cena) &
                       udaje.2014$real_cena <= 0] <- NA
udaje.2014$plocha[!is.na(udaje.2014$plocha) & udaje.2014$plocha <= 0] <- NA
zle.rok <- !is.na(udaje.2014$rok_diela) &
  (udaje.2014$rok_diela < 1800 | udaje.2014$rok_diela > udaje.2014$rok_aukcie)
udaje.2014$rok_diela[zle.rok] <- NA

udaje.2014$ln_cena   <- log(udaje.2014$real_cena)
udaje.2014$ln_plocha <- log(udaje.2014$plocha)
udaje.2014$vek_diela <- udaje.2014$rok_aukcie - udaje.2014$rok_diela

udaje.2014 <- udaje.2014[, c("real_cena", "ln_cena", "plocha", "ln_plocha",
                             "vek_diela", "rok_narodenia")]

# Diela s chýbajúcou hodnotou v použitých premenných vyradíme.
premenne <- c("ln_cena", "ln_plocha", "vek_diela", "rok_narodenia")
udaje.2014 <- udaje.2014[complete.cases(udaje.2014[premenne]), ]

Po príprave máme 9 015 diel predaných v roku 2014.

1.1 Uloženie a čistenie pamäte (workspace)

Aby sme v každom cvičení neopakovali načítanie a čistenie údajov, uložíme si všetky objekty z pamäte do súboru a v ďalšom cvičení ich na začiatku dokumentu načítame.

save.image("pamat.Rdata")

Vybrané objekty z pamäte vymažeme príkazom rm():

aaa <- 6
bbb <- "text"
rm(list = c("aaa", "bbb"))
"aaa" %in% ls()

[1] FALSE

Celú pamäť vyčistíme rm(list = ls()) a uložený obsah vrátime späť funkciou load():

rm(list = ls())
length(ls())

[1] 0

load("pamat.Rdata")
c("udaje.2014", "fmt") %in% ls()

[1] TRUE TRUE

load() budeme používať na začiatku dokumentu a save.image() na konci, aby sme si prenášali obsah pamäte z cvičenia do cvičenia. Treba však priebežne mazať nepotrebné objekty.

Poznámka

Balíky načítané cez library() sa uložením pamäte neukladajú, tie treba načítať v každom dokumente znova.

2 Odhad modelu a testy významnosti koeficientov

model <- lm(ln_cena ~ ln_plocha + vek_diela + rok_narodenia,
            data = udaje.2014)

tidy(model) |>
  mutate(
    Significance = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      p.value < 0.10  ~ ".",
      TRUE ~ ""
    )
  ) |>
  kable(
    format = "html", digits = 4,
    col.names = c("Premenná", "Odhad", "Štandardná chyba", "t-hodnota",
                  "p-hodnota", "Významnosť"),
    caption = "Odhady koeficientov viacnásobného modelu"
  ) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Odhady koeficientov viacnásobného modelu

Premenná	Odhad	Štandardná chyba	t-hodnota	p-hodnota	Významnosť
(Intercept)	-28.3380	3.2280	-8.7789	0	***
ln_plocha	0.5620	0.0138	40.8325	0	***
vek_diela	0.0071	0.0016	4.3531	0	***
rok_narodenia	0.0173	0.0016	10.5156	0	***

2.1 Ako čítať tabuľku

Stĺpec	Vysvetlenie
Odhad	Bodový odhad \(\hat\beta_j\): zmena očakávanej \(\ln Cena\) pri jednotkovej zmene premennej, ostatné premenné nezmenené (ceteris paribus).
Štandardná chyba	\(SE(\hat\beta_j)\), miera neistoty odhadu. Menšia chyba znamená presnejší odhad.
t-hodnota	\(t_j = \hat\beta_j / SE(\hat\beta_j)\), testovacia štatistika pre \(H_0:\beta_j=0\).
p-hodnota	Pravdepodobnosť získania aspoň takej extrémnej t-hodnoty, ak platí \(H_0\). Malá hodnota je dôkaz proti \(H_0\).
Významnosť	Hviezdičky: *** \(p<0{,}001\), ** \(p<0{,}01\), * \(p<0{,}05\), . \(p<0{,}10\).

Pre každý koeficient testujeme

\[ H_0:\beta_j=0 \qquad H_1:\beta_j\neq 0. \]

Za klasických predpokladov vrátane normality chýb má \(t_j\) pri platnosti \(H_0\) Studentovo \(t\)-rozdelenie s \(n-k\) stupňami voľnosti (\(k\) je počet koeficientov vrátane konštanty). Ak \(p<\alpha\) (zvyčajne \(\alpha=0{,}05\)), nulovú hypotézu zamietame.

Čo vidíme v našom modeli. Všetky tri smernice majú p-hodnotu pod 0,001, takže pri 5 % hladine zamietame \(H_0\) pre každú z nich. Pri \(n=\) 9015 pozorovaniach však aj malé efekty vychádzajú štatisticky významné; význam si preto treba overiť aj vecne:

• ln_plocha: o 1 % väčšia rozloha sa spája s približne 0,56 % vyššou cenou (ceteris paribus).
• vek_diela: jeden rok navyše sa spája s približne 0,71 % vyššou cenou.
• rok_narodenia: umelec narodený o rok neskôr sa spája s približne 1,73 % vyššou cenou.

Významnosť koeficientu nie je dôkaz kauzality. Pozor aj na to, že t-testy sú spoľahlivé len vtedy, keď model spĺňa predpoklady, čo si overíme v ďalších častiach.

3 Kvalita vyrovnania modelu

glance() dáva súhrn celého modelu. Stĺpec deviance (reziduálny súčet štvorcov) vynecháme.

glance(model) |>
  select(-deviance) |>
  kable(
    format = "html", digits = 4,
    col.names = c("R²", "Upravené R²", "Sigma", "F štatistika", "p-hodnota F",
                  "df", "logLik", "AIC", "BIC", "df.residual", "N"),
    caption = "Súhrnné charakteristiky modelu"
  ) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Súhrnné charakteristiky modelu

R²	Upravené R²	Sigma	F štatistika	p-hodnota F	df	logLik	AIC	BIC	df.residual	N
0.1763	0.176	1.8256	642.9807	0	3	-18216.14	36442.27	36477.81	9011	9015

Charakteristika	Význam	Hodnota v našom modeli
\(R^2 = 1 - RSS/TSS\)	Podiel variability \(\ln Cena\) zachytený modelom.	0,176
Upravené \(\bar R^2\)	\(R^2\) penalizované za počet premenných; vhodné pri porovnávaní modelov s rovnakou vysvetľovanou premennou a rovnakými pozorovaniami.	0,176
\(\hat\sigma=\sqrt{RSS/(n-k)}\)	Odhad smerodajnej odchýlky chyby, v jednotkách \(\ln Cena\).	1,826
F štatistika	Test spoločnej významnosti: \(H_0:\beta_1=\beta_2=\beta_3=0\).	643,0
df, df.residual	Počet vysvetľujúcich premenných a reziduálne stupne voľnosti \(n-k\).	3 a 9011
AIC, BIC	\(-2\ln L + 2k\) a \(-2\ln L + k\ln n\). Menšia hodnota je lepšia; BIC penalizuje zložitosť prísnejšie. Zmysel majú len pri porovnaní modelov na rovnakých údajoch s rovnakou vysvetľovanou premennou.	36442,3 a 36477,8

Model vysvetľuje len približne 17,6 % variability logaritmu ceny. F-test zamieta \(H_0\) o nulovosti všetkých smerníc, takže aspoň jedna premenná cenu štatisticky významne vysvetľuje. Malé \(R^2\) však ukazuje, že väčšina variability ceny ostáva nevysvetlená (umelec, technika, aukčný dom). Typická chyba modelu je približne 1,83 jednotky \(\ln Cena\), čo zodpovedá rozdielu približne 521 % v cene.

Vzťah medzi deviance, sigma a stupňami voľnosti si môžeme overiť:

c(sigma = sigma(model),
  z_deviance = sqrt(deviance(model) / df.residual(model)))

     sigma z_deviance 
   1.82564    1.82564 

3.1 F-test celkovej významnosti modelu

Zatiaľ čo t-test sa týka jedného koeficientu, F-test posudzuje model ako celok. Testujeme, či aspoň jedna vysvetľujúca premenná cenu vysvetľuje:

\[ H_0:\beta_1=\beta_2=\beta_3=0 \qquad H_1:\text{aspoň jeden z koeficientov }\beta_1,\beta_2,\beta_3\neq 0. \]

Testovacia štatistika je

\[ F=\frac{R^2/(k-1)}{(1-R^2)/(n-k)}, \]

kde \(k\) je počet koeficientov vrátane konštanty. Pri platnosti \(H_0\) a klasických predpokladoch má F-rozdelenie s \(k-1\) a \(n-k\) stupňami voľnosti.

f <- summary(model)$fstatistic
data.frame(F = f[["value"]], df1 = f[["numdf"]], df2 = f[["dendf"]],
           p = pf(f[["value"]], f[["numdf"]], f[["dendf"]], lower.tail = FALSE)) |>
  kable(format = "html", digits = 4,
        col.names = c("F štatistika", "df (čitateľ)", "df (menovateľ)", "p-hodnota"),
        caption = "F-test celkovej významnosti modelu") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

F-test celkovej významnosti modelu

F štatistika	df (čitateľ)	df (menovateľ)	p-hodnota
642.9807	3	9011	0

F-štatistika je 643,0 pri 3 a 9011 stupňoch voľnosti a p-hodnota je prakticky nulová (menej ako 0,001). Na 5 % hladine preto \(H_0\) zamietame: aspoň jedna z troch premenných má nenulový koeficient a model ako celok vysvetľuje cenu lepšie než model len s konštantou. Zamietnutie \(H_0\) ešte neznamená, že model je dobrý. Hodnotu \(R^2\) pozri vyššie.

4 Grafická diagnostika modelu

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

Diagnostické grafy modelu

par(mfrow = c(1, 1))

Červená krivka v grafoch je LOWESS (lokálne vážená regresia), ktorá vyhladzuje približne 2/3 najbližších pozorovaní, pričom najbližšie body majú najväčšiu váhu.

4.1 Residuals vs. Fitted

Graf ukazuje rezíduá proti vyrovnaným hodnotám. Pri správnej špecifikácii majú kolísať náhodne okolo nuly bez zjavného tvaru.

• Centrovanie: rezíduá kolíšu okolo nuly a červená čiara sa od nuly vzďaľuje len mierne, takže model nemá výrazné systematické skreslenie.
• Tvar: čiara je mierne zakrivená (vľavo vyššie, uprostred mierne pod nulou, vpravo opäť stúpa). Ide o náznak miernej nelinearity. Najvýraznejšie je vľavo, kde je málo pozorovaní.
• Rozptyl: rozptyl rezíduí je približne rovnaký v celom rozsahu vyrovnaných hodnôt, nevidno lievik.
• Odľahlé body: pozorovanie 58732 má rezíduum približne 8,7, čo je výrazne viac ako ostatné (rezíduá väčšiny pozorovaní sú medzi −5 a +6).

4.2 Q-Q plot

Štandardizované rezíduá sa porovnávajú s kvantilmi normálneho rozdelenia. Pri normalite ležia body na prerušovanej priamke.

• Väčšina bodov, najmä v strede (približne medzi −2 a +2), leží blízko priamky.
• Na pravom konci body ležia nad priamkou: pravý chvost je ťažší, než pri normálnom rozdelení. Na ľavom konci sú body mierne nad priamkou (ľahší ľavý chvost). Spolu to zodpovedá kladnej šikmosti rezíduí.
• Pozorovanie 58732 je najvýraznejšia odchýlka vpravo hore.

Zhoda s normalitou je teda približná, nie dokonalá. Formálny test je nižšie.

4.3 Scale-Location

Na osi y je druhá odmocnina absolútnych štandardizovaných rezíduí. Vodorovná červená čiara znamená konštantný rozptyl chýb.

• Červená čiara je takmer vodorovná, s miernym poklesom vľavo. Neukazuje systematický rast ani pokles rozptylu.
• Body netvoria lievik. Z grafu teda nevidíme závažné porušenie homoskedasticity. Je to ale len vizuálny dojem. Formálne sme rozptyl netestovali.
• Najvyššie body (58732, 47100, 47099) sú len jednotlivé odľahlé pozorovania.

4.4 Residuals vs. Leverage

Páka (leverage) meria, ako ďaleko je pozorovanie od priemeru vysvetľujúcich premenných. Pozorovanie je vplyvné, ak má veľkú páku aj veľké rezíduum. Čiarkované kontúry označujú Cookovu vzdialenosť 0,5 a 1.

• Takmer všetky pozorovania majú veľmi nízku páku (pod 0,005). Pri 9015 pozorovaniach je to očakávané.
• Niekoľko pozorovaní má vyššiu páku, najvyššiu (približne 0,022) má pozorovanie 76043. Jeho štandardizované rezíduum je ale malé (okolo −2), takže nie je vplyvné.
• Pozorovanie 58732 má veľké rezíduum, ale nízku páku.
• Žiadny bod sa nepribližuje ku kontúram Cookovej vzdialenosti 0,5 alebo 1. Maximálna Cookova vzdialenosť je 0,023. Žiadne pozorovanie teda neovplyvňuje odhady koeficientov neprimerane.

Grafy sú len orientačné. Pri veľkom počte pozorovaní sa body zlievajú, preto ich dopĺňame testami.

5 Testy rezíduí

5.1 Normalita rezíduí: Jarqueov-Berov test

Rezíduá sú lineárnou kombináciou chýb, \(\mathbf e=(\mathbf I-\mathbf X(\mathbf X^T\mathbf X)^{-1}\mathbf X^T)\mathbf u\), takže normalitu nepozorovateľných chýb \(\mathbf u\) posudzujeme nepriamo cez normalitu rezíduí \(\mathbf e\).

JB test využíva šikmosť \(S\) a špicatosť \(K\) (pri normálnom rozdelení \(S=0\), \(K=3\)):

\[ JB = \frac{n}{6}\left(S^2+\frac{(K-3)^2}{4}\right) \overset{a}{\sim}\chi^2(2). \]

\(H_0\): rezíduá sú normálne rozdelené. Ak \(p<0{,}05\), \(H_0\) zamietame. Štatistiku počítame priamo, aby sme nepotrebovali ďalší balík.

jb_test <- function(x) {
  x <- x - mean(x)
  n <- length(x)
  S <- mean(x^3) / mean(x^2)^1.5
  K <- mean(x^4) / mean(x^2)^2
  JB <- n / 6 * (S^2 + (K - 3)^2 / 4)
  data.frame(Sikmost = S, Spicatost = K, JB = JB,
             p = pchisq(JB, df = 2, lower.tail = FALSE))
}

jb1 <- jb_test(residuals(model))
jb1 |>
  kable(format = "html", digits = 4,
        col.names = c("Šikmosť", "Špicatosť", "JB štatistika", "p-hodnota"),
        caption = "Jarqueov-Berov test normality rezíduí") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Jarqueov-Berov test normality rezíduí

Šikmosť	Špicatosť	JB štatistika	p-hodnota
0.4743	3.2451	360.6361	0

JB test zamieta normalitu (JB = 360,6, p-hodnota <1e-16). Šikmosť 0,47 a špicatosť 3,25 sú však od hodnôt 0 a 3 odchýlené len mierne. Pri n = 9015 má test veľkú silu a zachytí aj malé odchýlky.

Poznámka

Existujú aj iné testy normality (Shapirov-Wilkov, Andersonov-Darlingov, Kolmogorovov-Smirnovov). JB test je v ekonometrii obvyklý a má jednoduchú interpretáciu cez šikmosť a špicatosť. Je asymptotický, preto ho dopĺňame Q-Q grafom.

6 Odľahlé pozorovania

6.1 Test odľahlých pozorovaní outlierTest()

Funkcia car::outlierTest() používa externe študentizované rezíduá

\[ t_i = \frac{e_i}{\hat\sigma_{(i)}\sqrt{1-h_{ii}}}, \]

kde \(h_{ii}\) je páka (leverage) a \(\hat\sigma_{(i)}\) je odhad \(\sigma\) po vynechaní \(i\)-teho pozorovania. Pri normalite chýb má \(t_i\) rozdelenie \(t\) s \(n-k-1\) stupňami voľnosti. Keďže testujeme všetky pozorovania, p-hodnoty sa upravujú Bonferroniho korekciou \(p_{Bonf}=\min(1, n\cdot p_i)\).

\(H_0\): pozorovanie nie je odľahlé, \(H_1\): pozorovanie je odľahlé. Rozhoduje Bonferroniho p-hodnota.

tabulka_outlier <- function(m) {
  ot <- outlierTest(m)
  data.frame(
    Pozorovanie = names(ot$rstudent),
    rstudent = unname(ot$rstudent),
    p = unname(ot$p),
    bonf = unname(ot$bonf.p)
  ) |>
    kable(format = "html", digits = 4,
          col.names = c("Pozorovanie (riadok)", "Študentizované rezíduum",
                        "Nekorigovaná p-hodnota", "Bonferroniho p-hodnota"),
          caption = "Test odľahlých pozorovaní s Bonferroniho korekciou") |>
    kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
}
tabulka_outlier(model)

Test odľahlých pozorovaní s Bonferroniho korekciou

Pozorovanie (riadok)	Študentizované rezíduum	Nekorigovaná p-hodnota	Bonferroniho p-hodnota
58732	4.8146	0	0.0135

Funkcia zobrazuje iba najextrémnejšie pozorovania. Ak je Bonferroniho p-hodnota pod 0,05, pozorovanie je štatisticky významne odľahlé. Odľahlé pozorovanie nemusí byť chybné a nemáme ho mazať automaticky: môže ísť o skutočne mimoriadny predaj (napríklad dielo slávneho autora).

V našom modeli identifikuje test pozorovanie 58732 ako štatisticky významne odľahlé (Bonferroniho p-hodnota je približne 0,014, teda menej ako 0,05). Je to jediné také pozorovanie z 9 015. Nemusí ísť o chybu, mohlo ísť o skutočne mimoriadny predaj.

6.2 Čo to znamená pre naše testy

Zistili sme dve veci: rezíduá nie sú presne normálne a jedno pozorovanie je významne odľahlé. Ovplyvňuje to závery t-testov a F-testu?

• Normalita. Presné t- a F-rozdelenie testových štatistík vyžaduje normálne chyby. Pri 9 015 pozorovaniach platia testy približne aj bez normality (asymptoticky) a odchýlka od normality je len mierna (pozri šikmosť a špicatosť v JB teste).
• Odľahlé pozorovanie. Overíme, či zmení odhady. Model odhadneme znova bez neho a porovnáme koeficienty.

odlahle <- names(outlierTest(model)$rstudent)[1]
model_bez <- update(model, data = udaje.2014[rownames(udaje.2014) != odlahle, ])

data.frame(
  Premenna = names(coef(model)),
  S = coef(model),
  Bez = coef(model_bez),
  p_s = summary(model)$coefficients[, 4],
  p_bez = summary(model_bez)$coefficients[, 4]
) |>
  kable(format = "html", digits = c(0, 4, 4, 4, 4), row.names = FALSE,
        col.names = c("Premenná", "Odhad (všetky diela)", "Odhad (bez odľahlého)",
                      "p-hodnota (všetky)", "p-hodnota (bez odľahlého)"),
        caption = "Model s odľahlým pozorovaním a bez neho") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Model s odľahlým pozorovaním a bez neho

Premenná	Odhad (všetky diela)	Odhad (bez odľahlého)	p-hodnota (všetky)	p-hodnota (bez odľahlého)
(Intercept)	-28.3380	-28.4992	0	0
ln_plocha	0.5620	0.5625	0	0
vek_diela	0.0071	0.0071	0	0
rok_narodenia	0.0173	0.0174	0	0

Po vynechaní pozorovania 58732 sa koeficienty zmenia len minimálne a všetky tri zostávajú štatisticky významné. \(R^2\) sa zmení z 0,176 na 0,177. Jediné odľahlé dielo teda závery neovplyvňuje, a preto ho v údajoch ponechávame.

Závery t-testov aj F-testu môžeme pri našej veľkej vzorke považovať za spoľahlivé z hľadiska normality a odľahlých hodnôt. Neoverovali sme však ostatné predpoklady (napríklad konštantný rozptyl chýb). Z grafov sme videli len vizuálny dojem.

save.image("pamat.Rdata")
rm(list = ls())