UJEP materiály

Pokročilé statistické metody

Užitečné odkazy

Statistika

Jednorozměrná statistika

Mnohorozměrná statistika

Metody

Metoda hlavních komponent

X <- wine[, setdiff(names(wine), "Cultivar")]
pca <- prcomp(X, scale. = TRUE)

Složení hlavní komponenty

pca$rotation

interpret_pca_loadings <- function(pca, components = 3, top_n = 5) {
  loadings <- pca$rotation

  for (pc in seq_len(components)) {
    cat(sprintf("PC%d:\n", pc))

    pc_loadings <- loadings[, pc]
    pc_loadings <- pc_loadings[order(abs(pc_loadings), decreasing = TRUE)]
    pc_loadings <- pc_loadings[seq_len(top_n)]

    for (var in names(pc_loadings)) {
      cat(sprintf("  %s: %.3f\n", var, pc_loadings[var]))
    }

    cat("\n")
  }

  invisible(loadings)
}

Každá komponenta má vlastní váhy a ty váhy říkají, které proměnné jsou pro danou komponentu více důležité a které méně. Komponenta funguje jako funkce, která z původních vlastností vína vypočte jedno číslo – skóre – určující jeho polohu vzhledem k dané komponentě.

PC_k = a1·x1 + a2·x2 + … + ap·xp

kde:

Konkrétní příklad:

PC1 =
-0.144·Alcohol
+0.245·Malic acid
+0.002·Ash
+0.239·Alcalinity of ash
-0.142·Magnesium
-0.395·Total phenols
-0.423·Flavanoids
+0.299·Nonflavanoid phenols
-0.313·Proanthocyanin
+0.089·Color intensity
-0.297·Hue
-0.376·OD280/OD315 of diluted wines
-0.287·Proline

Výstup a interpretace

interpret_pca <- function(pca, target_variance) {
  summary_pca <- summary(pca)
  explained <- summary_pca$importance["Proportion of Variance", ]
  cumulative <- summary_pca$importance["Cumulative Proportion", ]
  needed <- which(cumulative >= target_variance)[1]

  cat(sprintf("Požadavek: zachovat alespoň %.0f %% variability v datech.\n", target_variance * 100))
  cat(sprintf("Nejmenší počet komponent, který tento požadavek splňuje: %d.\n\n", needed))

  for (i in seq_len(needed)) {
    cat(sprintf("PC%d: %.1f %% (kumulativně %.1f %%)\n", i, explained[i] * 100, cumulative[i] * 100))
  }

  invisible(list(needed_components = needed, explained = explained, cumulative = cumulative))
}

Faktorová analýza

Diskriminační analýza

Shluková analýza

Kanonické korelace

Testy

Volba testu

Počet proměnných Počet skupin Test
1 2 t-test
1 >2 ANOVA
>1 2 Hotellingovo T²
>1 >2 MANOVA

t-test

ANOVA

Hotellingovo T²

MANOVA

Regrese

Jednoduchá lineární regrese

\[y = \beta_0 + \beta_1 x + e\] 
# x = Rozpočet na reklamu (Kč)
# y = Počet prodaných alb
x <- c(1000, 2000, 3000, 4000, 5000)
y <- c(130, 150, 170, 210, 240)
data <- data.frame(x, y)

model <- lm(y ~ x, data = data)
summary(model)
coef(model)

plot(
  data$x,
  data$y,
  pch = 19,
  xlab = "Rozpočet na reklamu",
  ylab = "Počet prodaných alb",
  main = "Lineární regrese: reklama vs. prodej"
)
abline(model, col = "red", lwd = 2)

Mnohonásobná lineární regrese

\[y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p + e\] 
x1 <- c(1000, 2000, 3000, 4000, 5000)
x2 <- c(2, 2, 4, 5, 5)
y  <- c(130, 145, 180, 215, 235)
data <- data.frame(y, x1, x2)

model <- lm(y ~ x1 + x2, data = data)
summary(model)
coef(model)

# Test homoskedasticity
library(lmtest)
bptest(model)

# Test multikolinearity
library(car)
vif(model)

# Test normálního rozdělení reziduí
shapiro.test(residuals(model))

# Cookova vzdálenost
cooks <- cooks.distance(model)
cooks[which.max(cooks.distance(model))]
plot(cooks, type = "h", lwd = 3, col = "darkred",
     main = "Cookova vzdálenost (Vlivné body)",
     xlab = "Index pozorování (řádek v datech)",
     ylab = "Vzdálenost (D)")
n <- nrow(data)
abline(h = 4/n, col = "blue", lty = 2)
which.max(cooks)

# install.packages("scatterplot3d")
library(scatterplot3d)
s3d <- scatterplot3d(
  data$x1, data$x2, data$y,
  pch = 16,
  color = "steelblue",
  grid = TRUE,
  box = TRUE,
  angle = 55,
  main = "Regresní rovina prodeje alb",
  xlab = "Reklama (x1)",
  ylab = "Hity (x2)",
  zlab = "Alba (y)"
)
s3d$plane3d(model, col = "red", lty.box = "dashed")

predict(model, newdata = data.frame(x1 = 6000, x2 = 3))

Koeficient determinace

Adjusted R-squared

model_summary <- summary(model)
model_summary$adj.r.squared

Shapiro-Wilkův test

Breusch-Paganův test

Cookova vzdálenost

Multikolinearita

\[VIF_{i} = \frac{1}{1 - R_{i}^2}\]
Hodnota VIF Význam Interpretace
1 Žádná korelace Proměnné spolu nesouvisí.
(1, 5> Střední korelace Proměnné jsou trochu podobné, ale model je stále stabilní.
>5 Vysoká korelace Proměnné jsou příliš podobné. Model je nespolehlivý.

Kroková regrese

data(mtcars)

model_full <- lm(mpg ~ ., data = mtcars)
model_null <- lm(mpg ~ 1, data = mtcars)

step(model_null, scope = formula(model_full), direction = "forward")
step(model_full, direction = "backward")
step(model_full, direction = "both")

Metoda maximální věrohodnosti

Regrese s kategorickou proměnnou

x <- c(1000, 2000, 3000, 4000, 5000)
y <- c(130, 150, 170, 210, 240)

# Dvě kapely jsou domácí (0), tři jsou zahraniční (1)
typ <- factor(c("domaci", "domaci", "zahranicni", "zahranicni", "zahranicni"))
data <- data.frame(x, y, typ)

# y = beta0 + beta1*x + beta2*typ
model_dummy <- lm(y ~ x + typ, data = data)
summary(model_dummy)

# Zkoumáme, jestli reklama funguje u zahraničních kapel jinak než u domácích
# y = beta0 + beta1*x + beta2*typ + beta3*(x * typ)
model_interakce <- lm(y ~ x * typ, data = data)
summary(model_interakce)

plot(
  data$x, data$y,
  pch = 19, col = as.numeric(data$typ),
  xlab = "Rozpočet na reklamu",
  ylab = "Počet prodaných alb",
  main = "Regrese s kategorickou proměnnou"
)

# Přímka pro domácí (černá)
abline(a = coef(model_interakce)[1], b = coef(model_interakce)[2], col = "black", lwd = 2)
# Přímka pro zahraniční (červená)
abline(a = coef(model_interakce)[1] + coef(model_interakce)[3], b = coef(model_interakce)[2] + coef(model_interakce)[4], col = "red", lwd = 2)

legend("topleft", legend = levels(data$typ), col = 1:2, pch = 19)
\[y = \beta_0 + \beta_1 x + \beta_2 D + e\]

Dummy variable

Pořadí (Index) (Intercept) sezonaLeto sezonaPodzim sezonaZima
1 (Jaro) 1 0 0 0
2 (Léto) 1 1 0 0
3 (Podzim) 1 0 1 0
4 (Zima) 1 0 0 1 
sezona <- factor(c("Jaro", "Jaro", "Leto", "Leto", "Podzim", "Podzim", "Zima", "Zima"))
prodeje <- c(100, 110, 150, 140, 120, 130, 80, 90)
data <- data.frame(prodeje, sezona)

model <- lm(prodeje ~ sezona, data = data)
summary(model)
model.matrix(model) # Tento příkaz vykreslí tabulku výše

Zápis modelu: \(y = \beta_0 + \beta_1 D_{Leto} + \beta_2 D_{Podzim} + \beta_3 D_{Zima} + \epsilon\)

Pokud chceme spočítat třeba prodeje pro Zimu, dosadíme do rovnice takto: \(D_{Leto} = 0\) \(D_{Podzim} = 0\) \(D_{Zima} = 1\) Rovnice se ti zjednoduší na $y = \beta_0 + \beta_3$

Logistická regrese

# Hodiny učení a výsledek (0 = neuspěl, 1 = uspěl)
data_test <- data.frame(
  hodiny = c(2, 5, 8, 10, 12, 15, 18, 20, 25, 30),
  uspech = c(0, 0, 0, 0, 1, 0, 1, 1, 1, 1)
)
model <- glm(uspech ~ hodiny, data = data_test, family = "binomial")
bod_zlomu <- -coef(model)[1] / coef(model)[2]

plot(uspech ~ hodiny, data = data_test, pch = 16, xlab = "Hodin učení", ylab = "Pravděpodobnost úspěchu", main = "Logistická regrese")
grid(col = "lightgray", lty = "dotted")
abline(h = 0.5, col = "red", lty = 2)
abline(v = bod_zlomu, col = "darkgreen", lty = 2)
# Přidání sigmoidy (S-křivky)
curve(predict(model, data.frame(hodiny = x), type = "response"), add = TRUE, col = "blue", lwd = 2)

# # Kolik % šance mám, když se učím 13 hodin?
predict(model, newdata = data.frame(hodiny = 13), type = "response")

vysledky <- function(m, d) {
  d$Sance <- round(predict(m, type = "response"), 3)
  d$Predpoved <- ifelse(d$Sance > 0.5, 1, 0)
  return(d)
}

vysledky(model, data_test)

Sigmoida

\[\sigma(z) = \frac{1}{1 + e^{-z}}\]

Poměr šancí

\[\log\!\left(\frac{p}{1 - p}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p\]

Ordinální regrese

\[\log\!\left(\frac{P(y \le k)}{P(y > k)}\right) = \theta_k - \beta_1 x_1 - \beta_2 x_2 - \cdots - \beta_p x_p\] 
library(MASS)
set.seed(1)

# Počet hodin, kolik student studoval a známku, kterou dostal
# C < B < A
data_skola <- data.frame(
  hodiny = c(2, 5, 8, 10, 11, 12, 15, 18, 20, 22, 25, 30),
  znamka = factor(c("C", "C", "B", "C", "B", "B", "B", "A", "B", "A", "A", "A"), levels = c("C", "B", "A"), ordered = TRUE)
)

# Ordinální regrese
m <- polr(znamka ~ hodiny, data = data_skola, Hess = TRUE)

ctable <- coef(summary(m))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = round(p, 4)))
#             Value Std. Error  t value p value
# hodiny  0.6018547  0.2966476 2.028854  0.0425
# C|B     5.0492378  2.9428805 1.715747  0.0862
# B|A    11.3597420  5.7131212 1.988360  0.0468

# Koeficient (0.6019): Je kladný, což znamená, že s každou další hodinou učení se zvyšuje šance na lepší známku.
# P-hodnota (0.0425): Je menší než 0.05, takže vliv učení je statisticky významný.

# Tady také model říká, jak těžké je "přeskočit" z jedné kategorie do druhé:
# - Práh C|B (5.049): Bod, kde se láme šance mezi C a lepšími známkami.
# - Práh B|A (11.360): Bod, kde se láme šance na čisté A.
# Závěr je tedy takový, že rozdíl mezi těmito čísly je velký (skoro 6 bodů na logitové stupnici).
# To znamená, že dostat se z B na A vyžaduje mnohem intenzivnější studium než se jen "odlepit" od Cčka.

# 4. Predikce pro studenta, který studoval 16 hodinam
predict(m, newdata = data.frame(hodiny = 16), type = "probs")
#          C          B          A
# 0.01014641 0.83927457 0.15057902
# Šance na C: 1,0 % (skoro nemožné, aby to takhle zkazil).
# Šance na B: 83,9 % (téměř jistota, že dostane béčko).
# Šance na A: 15,1 % (existuje naděje, ale model ho vidí spíše jako typického B studenta).

Poissonova regrese

\[\log(y) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p\] 
# Počet absencí studenta podle jeho průměru známek
data_absence <- data.frame(
  absence = c(0, 1, 2, 2, 3, 5, 8, 10),
  znamka = c(1.0, 1.2, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0)
)

# Poissonova regrese
mod_pois <- glm(absence ~ znamka, data = data_absence, family = "poisson")
exp(coef(mod_pois))

plot(absence ~ znamka, data = data_absence, pch = 16, col = "darkblue", main = "Poissonova regrese: Absence vs. Známka", xlab = "Průměr známek", ylab = "Počet absencí")
# Exponenciální křivka
curve(predict(mod_pois, data.frame(znamka = x), type = "response"), add = TRUE, col = "red", lwd = 2)

Polynomická regrese

\[y = \beta_0 + \beta_1 x + \beta_2 x^2 + \cdots + \beta_k x^k + e\] 
# Stres (0-100) a Výkon (0-100)
# Příliš málo stresu vede k nudě a nízkému výkonu, příliš mnoho stresu k vyhoření. Nejlepší výkon je někde uprostřed.
# Data tedy tvoří obrácené U (kopec).
data_vykon <- data.frame(
  stres = c(10, 20, 30, 40, 50, 60, 70, 80, 90, 100),
  vykon = c(20, 45, 65, 85, 95, 90, 75, 50, 30, 10)
)

# Polynomická regrese (2. řádu, tedy kvadratická)
mod_poly <- lm(vykon ~ poly(stres, degree = 2, raw = TRUE), data = data_vykon)

plot(vykon ~ stres, data = data_vykon, pch = 16, main = "Polynomická regrese (Kvadratická)")
curve(predict(mod_poly, data.frame(stres = x)), add = TRUE, col = "red", lwd = 2)