# Spustite raz v konzole, iba ak tieto balíky chýbajú.
install.packages(c("knitr", "rmarkdown", "broom", "kableExtra", "ggplot2"))cvičenie 4
Linearna regresia
1 Predmet výskumu
Pracujeme s údajmi o predaji umeleckých diel na aukciách (súbor transactions.csv).
Postup bude takýto:
- Zisťujeme, či väčšie dielo (väčšia rozloha) znamená vyššiu predajnú cenu.
- Graf závislosti ceny od rozlohy.
- Odhad regresnej priamky a jej interpretácia.
- Rozšírenie modelu o ďalšie vysvetľujúce premenné.
- Predikcia.
Modely odhadujeme funkciou lm(), výsledky pretvárame na tabuľky pomocou balíka broom a zobrazujeme ich cez kable() a kable_styling().
library(broom)
library(knitr)
library(kableExtra)
library(ggplot2)
# Pomocná funkcia: číslo so slovenskou desatinnou čiarkou pre texty v dokumente.
fmt <- function(x, d = 2) format(round(x, d), nsmall = d, decimal.mark = ",")1.1 Načítanie a upratanie CSV
Súbor transactions.csv nie je bežné CSV. Čiarka v ňom oddeľuje stĺpce, ale zároveň sa používa pri tisícoch (190,000), v dátumoch (March 31, 2006) a nachádza sa aj v názvoch diel či aukčných domov. Riadky majú preto rôzny počet stĺpcov (39 až 50 namiesto 38) a read.csv() by ich rozdelilo nesprávne.
Údaje preto čistíme ručne. Každý riadok rozdelíme podľa čiarok a potrebné hodnoty hľadáme podľa pevných bodov, ktoré ostávajú rovnaké aj pri posune stĺpcov:
- políčko
Lot ...(číslo lotu) – hneď za ním nasledujúprice_usd,cpi,real_price_usdaln_real_price, - políčko s obrázkom (
<div class=image ...) – o dve políčka za ním je rozloha diela v štvorcových palcoch (size (inch by inch)), - rok aukcie – štvorciferné číslo, ktoré sa opakuje hneď po texte dátumu,
- rok narodenia umelca a rok vzniku diela – dve políčka hneď za národnosťou umelca.
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.")
}
# Načítame surové riadky a odstránime nečitateľné znaky.
riadky <- readLines(cesta, warn = FALSE)
riadky <- iconv(riadky, from = "latin1", to = "ASCII", sub = "")
riadky <- riadky[-1] # vynecháme hlavičku
riadky <- riadky[nzchar(trimws(riadky))] # vynecháme prázdne riadky
# Text -> číslo; ak políčko nie je číslo (napr. "." alebo prázdne), vrátime NA.
ako_cislo <- function(x) {
if (length(x) == 1 && !is.na(x) && grepl("^-?[0-9]+(\\.[0-9]+)?$", x)) {
as.numeric(x)
} else {
NA_real_
}
}
# Z jedného riadku vytiahne päť potrebných hodnôt.
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] # kotva: číslo lotu
K <- which(startsWith(p, "<div class=image"))[1] # kotva: obrázok
if (is.na(L) || is.na(K)) return(vysledok)
# Reálna cena v USD a rozloha diela
vysledok["real_cena"] <- ako_cislo(p[L + 3])
vysledok["plocha"] <- ako_cislo(p[K + 2])
# Rok aukcie: štvorciferné číslo rovnaké ako predchádzajúce políčko (" 2014")
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
}
}
# Rok narodenia umelca a rok vzniku diela: za národnosťou (text)
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")Po upratovaní máme 114 283 záznamov o predaji. Jeden riadok je jedno dielo (lot) predané na aukcii.
1.2 Výber roka
Rovnako ako v pôvodnom zadaní pracujeme iba s jedným rokom. Zvolili sme rok aukcie 2014, pretože je to posledný rok s dostatkom pozorovaní, ktoré majú vyplnenú reálnu cenu (v roku 2015 reálna cena chýba úplne). Ide o prierezové údaje (cross-sectional data).
udaje.2014 <- transakcie[!is.na(transakcie$rok_aukcie) &
transakcie$rok_aukcie == 2014, ]
# Nezmyselné hodnoty nahradíme NA:
# cena a rozloha musia byť kladné; dielo nemohlo vzniknúť pred rokom 1800
# ani po roku aukcie.
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
# Nové premenné pre model
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")]
head(udaje.2014) |>
kable(format = "html", digits = 2, row.names = FALSE,
caption = "Prvých šesť diel vo vybraných údajoch (rok aukcie 2014)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| real_cena | ln_cena | plocha | ln_plocha | vek_diela | rok_narodenia |
|---|---|---|---|---|---|
| 16449 | 9.71 | 320.42 | 5.77 | 2 | 1957 |
| 60905 | 11.02 | 985.66 | 6.89 | 7 | 1957 |
| 83439 | 11.33 | 169.29 | 5.13 | 5 | 1957 |
| 293747 | 12.59 | 7430.69 | 8.91 | 6 | 1957 |
| 39368 | 10.58 | 182.27 | 5.21 | 4 | 1957 |
| 847598 | 13.65 | 2996.00 | 8.01 | NA | 1957 |
| Premenná | Význam |
|---|---|
| ln_cena | Prirodzený logaritmus reálnej ceny diela; vysvetľovaná premenná. |
| real_cena | Reálna cena diela v USD (očistená o infláciu pomocou CPI). |
| ln_plocha | Prirodzený logaritmus rozlohy diela. |
| plocha | Rozloha diela (výška × šírka) v štvorcových palcoch. |
| vek_diela | Vek diela v čase aukcie v rokoch (rok aukcie − rok vzniku diela). |
| rok_narodenia | Rok narodenia umelca. |
Ceny sú veľmi šikmo rozdelené (väčšina diel sa predá lacno, niekoľko za státisíce až milióny dolárov), preto logaritmujeme cenu aj rozlohu. Koeficient pri ln_plocha potom interpretujeme v percentách.
1.3 Chýbajúce hodnoty
Najprv spočítame, koľko hodnôt chýba (NA):
chybajuce <- data.frame(
Premenna = names(udaje.2014),
Pocet_NA = colSums(is.na(udaje.2014))
)
chybajuce |>
kable(format = "html", col.names = c("Premenná", "Počet chýbajúcich hodnôt"),
row.names = FALSE, caption = "Chýbajúce hodnoty pred vyradením") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Premenná | Počet chýbajúcich hodnôt |
|---|---|
| real_cena | 4085 |
| ln_cena | 4085 |
| plocha | 351 |
| ln_plocha | 351 |
| vek_diela | 1099 |
| rok_narodenia | 114 |
Chýba pomerne veľa hodnôt (hlavne cena alebo rok vzniku diela). Dopĺňanie mediánom by pri takom podiele výsledky skreslilo, preto diela s chýbajúcou hodnotou v niektorej z použitých premenných vyradíme.
premenne <- c("ln_cena", "ln_plocha", "vek_diela", "rok_narodenia")
udaje.2014.povodne <- udaje.2014 # ponecháme pôvodné údaje pre prípadnú neskoršiu prácu
udaje.2014 <- udaje.2014[complete.cases(udaje.2014[premenne]), ]
opis <- data.frame(
Premenna = premenne,
Minimum = sapply(udaje.2014[premenne], min),
Median = sapply(udaje.2014[premenne], median),
Priemer = sapply(udaje.2014[premenne], mean),
Maximum = sapply(udaje.2014[premenne], max)
)
opis |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Premenná", "Minimum", "Medián", "Priemer", "Maximum"),
caption = "Opisná štatistika po vyradení chýbajúcich hodnôt") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Premenná | Minimum | Medián | Priemer | Maximum |
|---|---|---|---|---|
| ln_cena | 2.20 | 8.91 | 9.13 | 17.08 |
| ln_plocha | -1.51 | 6.40 | 6.38 | 12.21 |
| vek_diela | 0.00 | 21.00 | 23.14 | 184.00 |
| rok_narodenia | 1919.00 | 1944.00 | 1946.39 | 2004.00 |
Z 14044 diel predaných v roku 2014 nám pre odhad ostalo 9015 diel.
2 Graf
Predpokladáme, že väčšie diela sa predávajú za vyššiu cenu. Zatiaľ je to len hypotéza.
graf_plocha <- ggplot(udaje.2014,
aes(x = ln_plocha, y = ln_cena)) +
geom_point(color = "steelblue", size = 1.5, alpha = 0.4) +
labs(x = "ln(rozloha diela v štvorcových palcoch)",
y = "ln(reálna cena v USD)") +
theme_minimal(base_size = 13)
graf_plocha
Z grafu vidno rastúcu závislosť, body sú však dosť rozptýlené: diela s podobnou rozlohou nemajú rovnakú cenu.
3 Jednoduchá lineárna regresia
3.1 Zápis modelu
Začneme modelom s jednou vysvetľujúcou premennou:
\[ \ln Cena_i = \beta_0 + \beta_1 \ln Plocha_i + u_i. \]
- \(i\) označuje dielo.
- \(\beta_0\) je konštanta, teda priesečník priamky s osou Y.
- \(\beta_1\) je smernica: o koľko sa zmení predikovaná hodnota \(\ln Cena\), keď \(\ln Plocha\) vzrastie o jednotku.
- \(u_i\) je náhodná zložka, teda to, čo jednoduchý model nevysvetľuje.
Priamka nemusí prechádzať všetkými bodmi. Metóda najmenších štvorcov ju určí tak, aby bol súčet štvorcov vertikálnych odchýlok bodov od priamky čo najmenší.
3.2 Odhad pomocou lm()
model1 <- lm(ln_cena ~ ln_plocha, data = udaje.2014)
vysledky1 <- tidy(model1)
vysledky1[c("term", "estimate")] |>
kable(format = "html", digits = 3, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Jednoduchá regresia: cena a rozloha diela") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 5.466 |
| ln_plocha | 0.575 |
Funkcia tidy(model1) z balíka broom vráti výsledky ako tabuľku, v ktorej jeden riadok zodpovedá jednému koeficientu (term je názov člena modelu, estimate jeho odhad, (Intercept) je konštanta). Zatiaľ zobrazujeme iba tieto dva stĺpce.
Funkcie kable() a kable_styling() upravujú len vzhľad tabuľky, odhad sa tým nemení.
Zápis ~ čítame ako „vysvetľované pomocou“: vľavo je vysvetľovaná premenná, vpravo vysvetľujúce. Konštantu R pridá automaticky.
Odhadnutá rovnica je:
\[ \widehat{\ln Cena}_i = 5,466 + 0,575\,\ln Plocha_i. \]
Smernica: Cena aj rozloha sú v logaritmoch, takže smernicu čítame v percentách: dielo s o 1 % väčšou rozlohou má v priemere o približne 0,57 % vyššiu predikovanú cenu.
Konštanta: Pri \(\ln Plocha = 0\) (rozloha 1 štvorcový palec) model predikuje \(\ln Cena\) približne 5,47. Takéto malé dielo je mimo typického rozsahu údajov, preto konštanta slúži hlavne na určenie polohy priamky.
Poznámka
Kladná závislosť neznamená, že zväčšením konkrétneho diela by sa jeho cena zvýšila. Diela sa líšia aj inými vlastnosťami (umelec, technika, kvalita, aukčný dom).
3.3 Regresná priamka v grafe
K bodovému grafu pridáme priamku pomocou geom_abline(), pričom konštantu aj smernicu berieme priamo z modelu model1.
graf_plocha +
geom_abline(intercept = coef(model1)["(Intercept)"],
slope = coef(model1)["ln_plocha"],
color = "firebrick", linewidth = 1)
4 Viacnásobná lineárna regresia
Cena diela môže závisieť od viacerých vecí naraz, preto model rozšírime o vek diela a rok narodenia umelca:
\[ \ln Cena_i = \beta_0 + \beta_1 \ln Plocha_i + \beta_2 VekDiela_i + \beta_3 RokNarodenia_i + u_i. \]
model <- lm(ln_cena ~ ln_plocha + vek_diela + rok_narodenia,
data = udaje.2014)
vysledky <- tidy(model)
vysledky[c("term", "estimate")] |>
kable(format = "html", digits = 4, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Viacnásobná regresia: rozloha, vek diela a rok narodenia umelca") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | -28.3380 |
| ln_plocha | 0.5620 |
| vek_diela | 0.0071 |
| rok_narodenia | 0.0173 |
Odhadnutá rovnica je:
\[ \widehat{\ln Cena}_i = -28,338 + 0,562\,\ln Plocha_i + 0,0071\,VekDiela_i + 0,0173\,RokNarodenia_i. \]
Koeficient pri ln_plocha teraz platí pri rovnakom veku diela a rovnakom roku narodenia umelca. Dielo s o 1 % väčšou rozlohou má pri nezmenených ostatných premenných o približne 0,56 % vyššiu predikovanú cenu.
Koeficient pri vek_diela vyzerá malý, lebo ide o zmenu za jeden rok. Pri logaritmovanej cene ho čítame približne v percentách: každý ďalší rok veku diela sa pri nezmenených ostatných premenných spája s rozdielom v cene približne 0,71 %.
Podobne pri rok_narodenia: umelec narodený o rok neskôr sa pri nezmenených ostatných premenných spája s rozdielom v cene približne 1,73 %. Je to len popis rozdielov medzi dielami v našej vzorke, nie dôkaz, že rok narodenia cenu spôsobuje.
Koeficient pri rozlohe sa po pridaní ďalších premenných zmenil (z 0,575 v jednoduchom na 0,562 vo viacnásobnom modeli), pretože viacnásobný model zohľadňuje aj vek diela a rok narodenia umelca.
5 Stručný súhrn modelu pomocou glance()
Zatiaľ čo tidy() dáva riadok pre každý koeficient, glance() dáva jeden riadok so súhrnom celého modelu. Vyberáme počet pozorovaní, \(R^2\) a odhad smerodajnej odchýlky chyby.
suhrn <- glance(model)
suhrn[c("nobs", "r.squared", "sigma")] |>
kable(format = "html", digits = c(0, 3, 3), row.names = FALSE,
col.names = c("Počet diel", "R²", "Odhad SD chyby (v log. jednotkách)"),
caption = "Základné charakteristiky viacnásobného modelu") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Počet diel | R² | Odhad SD chyby (v log. jednotkách) |
|---|---|---|
| 9015 | 0.176 | 1.826 |
nobsje počet diel použitých pri odhade.r.squaredje podiel výberovej variability \(\ln Cena\) zachytený modelom.sigmaje odhad štandardnej odchýlky modelovej chyby; keďže vysvetľovaná premenná je v logaritmoch, nie je v dolároch, ale v jednotkách \(\ln Cena\).
Hodnota \(R^2\) vo viacnásobnom modeli je približne 0,176, takže model vysvetľuje približne 17,6 % výberovej variability logaritmu ceny. Zvyšok pripadá na faktory mimo modelu (napríklad umelec, technika alebo aukčný dom). Neznamená to, že toľko percent predikcií je správnych, ani to nedokazuje kauzalitu.
Úplný výpis dá summary(model). tidy(model) obsahuje aj std.error, statistic a p.value, tie však zatiaľ nezobrazujeme ani neinterpretujeme, lebo patria k štatistickému testovaniu.
6 Skutočná hodnota, predikcia a rezíduum
Pre každé dielo porovnáme skutočnú hodnotu \(\ln Cena\) s hodnotou na regresnej ploche:
\[ e_i = y_i - \widehat y_i. \]
Rozdiel nazývame rezíduum. Kladné rezíduum znamená, že skutočná hodnota je vyššia ako predikcia, záporné opak. Rezíduum vieme spočítať, skutočnú chybu \(u_i\) však nepoznáme.
porovnanie <- data.frame(
Dielo = seq_len(nrow(udaje.2014)),
Skutocnost = udaje.2014$ln_cena,
Predikcia = fitted(model),
Reziduum = residuals(model)
)
head(porovnanie, 8) |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Dielo (poradie)", "Skutočnosť", "Predikcia", "Rezíduum"),
caption = "Skutočné a predikované hodnoty ln(ceny)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Dielo (poradie) | Skutočnosť | Predikcia | Rezíduum |
|---|---|---|---|
| 1 | 9.71 | 8.82 | 0.89 |
| 2 | 11.02 | 9.49 | 1.53 |
| 3 | 11.33 | 8.49 | 2.85 |
| 4 | 12.59 | 10.62 | 1.97 |
| 5 | 10.58 | 8.52 | 2.06 |
| 6 | 12.29 | 10.28 | 2.01 |
| 7 | 9.52 | 8.07 | 1.45 |
| 8 | 6.12 | 7.80 | -1.68 |
7 Predikcia pre zadanú hodnotu
Predikciu ukážeme na jednoduchom modeli s rozlohou. Pre diela s rozlohou 100, 600 a 3 000 štvorcových palcov (približne 25 × 4 palce, 24 × 25 palcov a 55 × 55 palcov) vyjde:
nove_udaje <- data.frame(plocha = c(100, 600, 3000))
nove_udaje$ln_plocha <- log(nove_udaje$plocha)
nove_udaje$Predikcia_ln <- predict(model1, newdata = nove_udaje)
nove_udaje$Predikcia_USD <- exp(nove_udaje$Predikcia_ln)
nove_udaje[c("plocha", "Predikcia_ln", "Predikcia_USD")] |>
kable(format = "html", digits = c(0, 2, 0), row.names = FALSE,
col.names = c("Rozloha (štvorcové palce)", "Predikcia ln(ceny)",
"Predikcia ceny (USD)"),
caption = "Predikcie jednoduchého regresného modelu") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Rozloha (štvorcové palce) | Predikcia ln(ceny) | Predikcia ceny (USD) |
|---|---|---|
| 100 | 8.11 | 3334 |
| 600 | 9.14 | 9333 |
| 3000 | 10.07 | 23532 |
Pomocou exp() prevedieme predikovaný logaritmus ceny späť na dolárovú hodnotu.