Robustní a neparametrické statistické metody – projekt
Tento projekt se zabývá robustními a neparametrickými statistickými metodami, které slouží jako alternativa ke klasickým statistickým metodám. Tyto metody jsou velmi často založené na silných předpokladech, například o rozdělení pravděpodobnosti, ze kterého pochází náhodný výběr, proto je nutné v praxi při analýze tyto předpoklady ověřit a pokud nejsou splněny, použít jiné (robustní či neparametrické) metody, které jsou vůči nesplnění předpokladů odolnější.
Cílem projektu je porovnat klasické a robustní, resp. neparametrické metody na simulovaných datech. Nejprve se budeme věnovat metodám pro odhad regresních koeficientů v lineárním modelu, přičemž bude porušen předpoklad normality náhodných chyb a také se v datech budou vyskytovat odlehlá pozorování.
Ve druhé části projektu se poté zaměříme na pořadové testy a jejich porovnání s klasickými statistickými metodami. Použijeme stejná data jako v situaci lineárního modelu, budeme však analyzovat pouze jednu vybranou dimenzi.
Dokument je rozčleněn do následujících sekcí:
- Vytvoření datového souboru – sekce se zabývá simulováním dat, která následně použijeme pro porovnání klasických a robustních metod pro odhad regresních koeficientů v lineárním regresním modelu,
- Klasické statistické metody – v této sekci aplikujeme klasický lineární regresní model a ověříme splnění všech předpokladů modelu,
- Robustní statistické metody – zde použijeme robustní odhady v lineárním regresním modelu a porovnáme s klasickým odhadem,
- Simulační studie – sekce se věnuje všem třem předchozím bodům dohromady, přičemž nás zajímají závěry a porovnání metod na základě více než jednoho nasimulovaného náhodného výběru,
- Pořadové testy – následně se budeme věnovat neparametrickým testům jakožto alternativě ke klasickým, přičemž se podíváme na jedno i dvouvýběrové testy a následně na porovnání výsledků s klasickými testy.
Začněme tedy nejprve vytvořením datového souboru.
Vytvoření datového souboru
# nacteni knihoven
library(MASS)
library(HSAUR)
library(ggplot2)
library(rrcov)
library(ddalpha)
library(ICSNP)
library(aplpack)
library(dplyr)
library(tidyr)
library(tidyverse)
library(gtsummary)
library(quantreg)
library(mblm)
library(nortest)
library(BSDA)Pro nasimulování dat pro účely lineární regrese si nejprve musíme
definovat funkci, resp. funkce, ze které budou naše data pocházet.
Uvažujme lineární regresní model s jednou spojitou kovariátou a jednou
kategoriální se dvěma kategoriemi. Ty označme například Muž
a Žena, aby se nám snáze interpretovali výsledky a data
byla lépe uchopitelná. Uvažujme pro každé pohlaví jinou regresní
parabolu, která může představovat například vývoj výšky jedince v
závislosti na věku.
V tomto projektu budeme pracovat s následujícími funkcemi. Pro muže se bude jednat o funkci ve tvaru
\[ f_{male}(x) = 50 + 9x - 0,\!1 x^2 \]
a pro ženy o funkci ve tvaru
\[ f_{female}(x) = 78 + 6x - 0,\!05 x^2. \]
# regresni funkce pro muze
male_lm <- function(x) {return(50 + 9 * x - 0.1 * x^2)}
female_lm <- function(y) {return(78 + 6 * y - 0.05 * y^2)}Tyto křivky si pro lepší představu vykresleme na intervalu \(I = [5; 18]\), přičemž barevně odližíme křivky pro muže (modrá barva) a ženy (červená barva). Toto barevné značení budeme využívat v celém dokumentu.
# vykresleni generovanych primek
xx <- seq(5, 18, length = 501)
data.frame(x = c(xx, xx),
y = c(male_lm(xx), female_lm(xx)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi)) |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line() +
theme_bw() +
scale_color_manual(values=c("deepskyblue2", "tomato3"))
Simulace dat
Nyní již můžeme přistoupit k simulaci dat. Aby byla analýza
reproducibilní, nastavíme generátor pseudonáhodných čísel
(set.seed(21)). Dále zvolíme hodnoty mu a
n, které představují střed (střední hodnotu rozdělení, ze
kterého budeme generovat hodnoty kovariáty), resp. počet pozorování na
ose \(x\) pro muže a ženy. Hodnota
sigma_x je směrodatnou odchylkou pro generování hodnot
kovariáty a df_y jsou stupně volnosti Studentova rozdělení,
ze kterého budeme generovat náhodné chyby. Toto rozdělení jsme zvolili
proto, aby byl porušen jeden z předpokladů klasického lineárního modelu.
Volíme relativně malou hodnotu stupňů volnosti, aby se rozdělení
výrazněji lišilo od normálního.
Nejprve tedy generujeme hodnoty kovariáty v závoslosti na hodnotě kategoriální proměnné (muž, resp. žena) z normálního rozdělení \[ X|(Gender = g) \sim \text N(\mu_g, \sigma_x^2). \] Následně generujeme hodnoty závislé proměnné pomocí výše vytvořených funkcí. Pokud bychom chtěli matematicky zapsat, jak vypadá model našich dat, mohli bychom uvažovat například následující zápis: \[ Y_i = \beta_0 + \beta_1x_i + \beta_2x_i^2 + \beta_3 \times \mathbb I\{Gender_i = F\} +\qquad\qquad\qquad\ \ \ \ \\ \beta_4x_i \times \mathbb I\{Gender_i = F\} + \beta_5x_i^2 \times \mathbb I\{Gender_i = F\} + \varepsilon_i,\\ \varepsilon_i \stackrel{iid}{\sim} 2\cdot t_{df_y}, \quad i = 1, 2, \dots, n_M + n_F. \]
# simulace dat
set.seed(21)
mu_male <- 12.5
mu_female <- 11.5
n_male <- 85
n_female <- 94
sigma_x <- 3.5
df_y <- 2
x_male <- rnorm(n_male, mean = mu_male, sd = sigma_x)
x_female <- rnorm(n_female, mean = mu_female, sd = sigma_x)
y_male <- male_lm(x_male) + rt(n_male, df_y) * 2
y_female <- female_lm(x_female) + rt(n_female, df_y) * 2
data <- data.frame(
x = c(x_male, x_female),
y = c(y_male, y_female),
Pohlavi = rep(c('Muz', 'Zena'), c(n_male, n_female))
) |>
mutate(Pohlavi = factor(Pohlavi))Vygenerovaná data si vykresleme a přidejme i křivky, ze kterých se data generovala.
# vykresleni dat
xx_male <- seq(min(x_male), max(x_male), length = 501)
xx_female <- seq(min(x_female), max(x_female), length = 501)
data.frame(x = c(xx_male, xx_female),
y = c(male_lm(xx_male), female_lm(xx_female)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx_male))) |>
mutate(Pohlavi = factor(Pohlavi)) |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.3) +
theme_bw() +
geom_point(data = data, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3"))
Přidání odlehlých hodnot
Abychom lépe viděli výhody použití robustních odhadů a zároveň porušili ještě více předpoklady klasického lineárního modelu, přidáme do datového souboru odlehlá pozorování. Nebudeme je však přidávat pevně, ale opět si je nasimulujeme, nyní z dvourozměrného normálního rozdělení s danou kovarianční maticí a vektory středních hodnot. Pro muže vygenerujeme dva a pro ženy tři odlehlá pozorování.
# pridame odlehle hodnoty
center_male <- c(21, 140)
center_female <- c(20, 145)
Sigma <- matrix(c(0.5, -0.1, -0.1, 10), nrow = 2, byrow = F)
df <- rbind(mvrnorm(n = 2, mu = center_male, Sigma = Sigma),
mvrnorm(n = 3, mu = center_female, Sigma = Sigma)) |>
data.frame()
colnames(df) <- c('x', 'y')
df <- df |> mutate(Pohlavi = rep(levels(data$Pohlavi), c(2, 3)) |> factor())
df <- rbind(data, df)Opět si výsledná data znázorněme graficky. Tato data jsou již připravena pro použití na klasické i robustní metody.
# vykresleni dat
xx_male <- seq(min(x_male), max(x_male), length = 501)
xx_female <- seq(min(x_female), max(x_female), length = 501)
data.frame(x = c(xx_male, xx_female),
y = c(male_lm(xx_male), female_lm(xx_female)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx_male))) |>
mutate(Pohlavi = factor(Pohlavi)) |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.3) +
theme_bw() +
geom_point(data = df, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3"))
Klasické statistické metody
Podívejme se nejprve na to, jak si s našimi nasimulovanými daty poradí klasická statistika, tedy klasický lineární regresní model, který předpokládá normalitu náhodných chyb, nezávislost a homoskedasticitu.
Lineární regresní model
Jelikož v praxi nevíme, jak daná funkční závislost vypadala, postupujme nyní podobně (i když my původní funkce známe). Definujme několik modelů a porovnejme je mezi sebou, čímž vybereme ten pro naše data nejlepší. Jelikož máme jednu spojitou a jednu kategoriální kovariátu, uvažujme následující modely:
- lineární aditivní model,
- lineární model s interakcí,
- kvadratický aditivní model,
- kvadratický model s interakcí u lineárního členu,
- kvadratický model s interakcí u kvadratického členu a konečně
- kvadratický model s interakcí u obou členů (nejobecnější kvadratický model, tento model je “správný”).
mod0 <- lm(y ~ x + Pohlavi, data = df)
mod1 <- lm(y ~ x * Pohlavi, data = df)
mod2 <- lm(y ~ x + I(x^2) + Pohlavi, data = df)
mod3 <- lm(y ~ x * Pohlavi + I(x^2), data = df)
mod4 <- lm(y ~ x + I(x^2) * Pohlavi, data = df)
mod5 <- lm(y ~ (x + I(x^2)) * Pohlavi, data = df)
cr <- data.frame(Model = paste0('model ', 0:5),
AIC = c(AIC(mod0), AIC(mod1), AIC(mod2),
AIC(mod3), AIC(mod4), AIC(mod5)),
BIC = c(BIC(mod0), BIC(mod1), BIC(mod2),
BIC(mod3), BIC(mod4), BIC(mod5)),
R_sq_adj = c(summary(mod0)$adj.r.squared, summary(mod1)$adj.r.squared,
summary(mod2)$adj.r.squared, summary(mod3)$adj.r.squared,
summary(mod4)$adj.r.squared, summary(mod5)$adj.r.squared))Pro porovnání modelů použijeme kritéria \(AIC\), \(BIC\) a \(R^2_{adj}\). V tabulce níže jsou uvedeny hodnoty jednotlivých kritérií pro všech šest výše zmíněných modelů.
| Model | AIC | BIC | \(R^2_{adj}\) |
|---|---|---|---|
| model 0 | 1302.614 | 1315.474 | 0.8436 |
| model 1 | 1292.374 | 1308.448 | 0.8529 |
| model 2 | 1270.193 | 1286.268 | 0.8696 |
| model 3 | 1229.218 | 1248.508 | 0.8962 |
| model 4 | 1237.024 | 1256.313 | 0.8917 |
| model 5 | 1226.360 | 1248.865 | 0.8983 |
Vidíme, že nejlépe vychází model 5, tedy ten nejobecnější a ten,
který odpovídá realitě. Modely také můžeme porovnat pomocí funkce
anova()
print(paste0('Analysis of Variance between model 3 and model 5: p-value = ', round(anova(mod3, mod5)$'Pr(>F)'[2], 4)))## [1] "Analysis of Variance between model 3 and model 5: p-value = 0.0304"print(paste0('Analysis of Variance between model 4 and model 5: p-value = ', round(anova(mod4, mod5)$'Pr(>F)'[2], 8)))## [1] "Analysis of Variance between model 4 and model 5: p-value = 0.00047381"Opět vychází nejlépe poslední model. Vypišme si odhady regresních koeficientů společně s intervaly spolehlivosti a příslušnými p-hodnotami do tabulky.
mod <- mod5
tbl_regression(mod,
pvalue_fun = function(x) style_pvalue(x, digits = 3),
estimate_fun = function(x) style_number(x, digits = 3)) %>%
bold_labels() | Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| x | 13.555 | 11.451, 15.658 | <0.001 |
| I(x^2) | -0.306 | -0.386, -0.227 | <0.001 |
| Pohlavi | |||
| Muz | — | — | |
| Zena | 39.129 | 22.771, 55.486 | <0.001 |
| x * Pohlavi | |||
| x * Zena | -4.888 | -7.596, -2.179 | <0.001 |
| I(x^2) * Pohlavi | |||
| I(x^2) * Zena | 0.120 | 0.011, 0.229 | 0.030 |
| 1 CI = Confidence Interval | |||
Můžeme si také vykreslit obě regresní paraboly společně s bodovými intervaly spolehlivosti. Všimněme si, jak si oněch pět odlehlých pozorování (úplně vpravo s hodnotami závislé proměnné pod 150) k sobě přitáhlo regresní křivky.
xx <- seq(min(df$x), max(df$x), length = 501)
df_plot <- data.frame(x = c(xx, xx),
y = c(predict(mod, newdata = data.frame(x = xx,
Pohlavi = 'Muz')),
predict(mod, newdata = data.frame(x = xx,
Pohlavi = 'Zena'))),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi))
CI <- rbind(predict(mod,
newdata = data.frame(x = xx, Pohlavi = 'Muz'),
interval = 'confidence'),
predict(mod,
newdata = data.frame(x = xx, Pohlavi = 'Zena'),
interval = 'confidence')) |> as.data.frame()
df_plot$lwr <- CI$lwr
df_plot$upr <- CI$upr
df_plot |> ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.9, linewidth = 0.8) +
theme_bw() +
geom_point(data = df, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3")) +
geom_ribbon(aes(ymin = lwr, ymax = upr, fill = Pohlavi), alpha = 0.3,
colour = NA) +
scale_fill_manual(values=c("deepskyblue2", "tomato3"))
Diagnostika modelu
Podívejme se nyní na diagnostiku výsledného modelu. Zajímá nás zejména normalita, odlehlá pozorování a homoskedasticita. Z grafů níže můžeme vidět, že ani jeden diagnostický graf nevypadá dobře. Cookova vzdálenost dvou pozorování je nad hranicí 0,5 a jednoho dokonce nad hodnotou 1. Q-Q graf reziduí vykazuje silné porušení normality (těžké chvosty) a graf standardizovaných reziduí proti fitovaným hodnotám ukazuje na heteroskedasticitu.
Pro formálnost ještě proveďme Shapirův-Wilkův test normality. P-hodnota je velmi malá, tedy i test zamítá normalitu reziduí.
## [1] "Shapiro-Wilk normality test: p-value = 5.59938297005324e-14"Jednou z možností, jak taková data uchopit, je odstranit odlehlá pozorování. Viděli jsme (ať už z Cookovy vzdálenosti, ostatních diagnostických grafů nebo naší znalosti vzniku dat), že se v datech nacházejí odlehlé hodnoty. Odstraňme je proto z další analýzy a proveďme celý postup znovu.
df_reduced <- df[-c(180:184), ]
mod0 <- lm(y ~ x + Pohlavi, data = df_reduced)
mod1 <- lm(y ~ x * Pohlavi, data = df_reduced)
mod2 <- lm(y ~ x + I(x^2) + Pohlavi, data = df_reduced)
mod3 <- lm(y ~ x * Pohlavi + I(x^2), data = df_reduced)
mod4 <- lm(y ~ x + I(x^2) * Pohlavi, data = df_reduced)
mod5 <- lm(y ~ (x + I(x^2)) * Pohlavi, data = df_reduced)
cr <- data.frame(Model = paste0('model ', 0:5),
AIC = c(AIC(mod0), AIC(mod1), AIC(mod2),
AIC(mod3), AIC(mod4), AIC(mod5)),
BIC = c(BIC(mod0), BIC(mod1), BIC(mod2),
BIC(mod3), BIC(mod4), BIC(mod5)),
R_sq_adj = c(summary(mod0)$adj.r.squared, summary(mod1)$adj.r.squared,
summary(mod2)$adj.r.squared, summary(mod3)$adj.r.squared,
summary(mod4)$adj.r.squared, summary(mod5)$adj.r.squared))| Model | AIC | BIC | \(R^2_{adj}\) |
|---|---|---|---|
| model 0 | 1064.243 | 1076.993 | 0.9510 |
| model 1 | 1023.830 | 1039.766 | 0.9611 |
| model 2 | 1065.848 | 1081.785 | 0.9509 |
| model 3 | 1009.358 | 1028.483 | 0.9643 |
| model 4 | 1020.193 | 1039.318 | 0.9621 |
| model 5 | 1004.692 | 1027.003 | 0.9655 |
print(paste0('Analysis of Variance between model 3 and model 5: p-value = ', round(anova(mod3, mod5)$'Pr(>F)'[2], 4)))## [1] "Analysis of Variance between model 3 and model 5: p-value = 0.0113"print(paste0('Analysis of Variance between model 4 and model 5: p-value = ', round(anova(mod4, mod5)$'Pr(>F)'[2], 8)))## [1] "Analysis of Variance between model 4 and model 5: p-value = 4.011e-05"Opět se jako nejvíce vhodný pro modelování našich dat jeví model 5, který je nejvíce obecný. Vypišme si znovu odhady regresních koeficientů.
mod <- mod5
tbl_regression(mod,
pvalue_fun = function(x) style_pvalue(x, digits = 3),
estimate_fun = function(x) style_number(x, digits = 3)) %>%
bold_labels() | Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| x | 9.524 | 8.179, 10.868 | <0.001 |
| I(x^2) | -0.127 | -0.179, -0.074 | <0.001 |
| Pohlavi | |||
| Muz | — | — | |
| Zena | 30.794 | 20.462, 41.126 | <0.001 |
| x * Pohlavi | |||
| x * Zena | -3.820 | -5.608, -2.031 | <0.001 |
| I(x^2) * Pohlavi | |||
| I(x^2) * Zena | 0.098 | 0.023, 0.173 | 0.011 |
| 1 CI = Confidence Interval | |||
Nyní již odhadnuté regresní paraboly vypadají mnohem lépe.
xx <- seq(min(df$x), max(df$x), length = 501)
df_plot <- data.frame(x = c(xx, xx),
y = c(predict(mod, newdata = data.frame(x = xx,
Pohlavi = 'Muz')),
predict(mod, newdata = data.frame(x = xx,
Pohlavi = 'Zena'))),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi))
CI <- rbind(predict(mod,
newdata = data.frame(x = xx, Pohlavi = 'Muz'),
interval = 'confidence'),
predict(mod,
newdata = data.frame(x = xx, Pohlavi = 'Zena'),
interval = 'confidence')) |> as.data.frame()
df_plot$lwr <- CI$lwr
df_plot$upr <- CI$upr
df_plot |> ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.9, linewidth = 0.8) +
theme_bw() +
geom_point(data = df_reduced, aes(x, y, colour = Pohlavi)) +
geom_point(data = df, aes(x, y, colour = Pohlavi), alpha = 0.25) +
scale_color_manual(values=c("deepskyblue2", "tomato3")) +
geom_ribbon(aes(ymin = lwr, ymax = upr, fill = Pohlavi), alpha = 0.3,
colour = NA) +
scale_fill_manual(values=c("deepskyblue2", "tomato3"))
Podívejme se na diagnostiku modelu s odtraněnými odlehlými pozorováními. Všimněme si, že většina grafů – zejména první graf reziduí proti fitovaným hodnotám a graf standardizovaných reziduí – se výrazně zlepšily. Naopak Q-Q graf reziduí vypadá velmi podobně, opět rezidua vykazují výrazné těžké chvosty.
## [1] "Shapiro-Wilk normality test: p-value = 1.21527977505141e-09"P-hodnota Shapirova-Wilkova testu normality je výrazně pod hladinou významnosti, tedy není splněn předpoklad normality. Bez tohoto předpokladu se nemůžeme spolehnout na výsledky provedených testů o regresních koeficientech a na intervaly spolehlivosti.
Diskuze
Podívejme se na srovnání odhadů regresních koeficientů s jejich skutečnými hodnotami. Z tabulky níže vidíme, že závěry, které bychom udělali na základě lineárního modelu s odlehlými pozorováními, by byly velmi nepřesné a zavádějící. Naopak odhady po odstranění odlehlých pozorování jsou výrazně lepší, avšak musíme mít na paměti, že jsme jednak využili naši znalost o tom, která jsou odlehlá pozorování a také to, že ani jeden z modelů nesplňuje předpoklad normality.
mod_lm_out <- lm(y ~ (x + I(x^2)) * Pohlavi, data = df)
mod_lm_red <- lm(y ~ (x + I(x^2)) * Pohlavi, data = df_reduced)
# skutecne odhady zname ze simulace
estimates <- data.frame(beta_0_male = 50, beta_1_male = 9,
beta_2_male = -0.1, beta_0_female = 78,
beta_1_female = 6, beta_2_female = -0.05)
# odhady v modelu s outliery
beta_out <- mod_lm_out |> coef()
estimates <- rbind(estimates,
c(beta_out[1], beta_out[2], beta_out[3],
beta_out[1] + beta_out[4], beta_out[2] + beta_out[5],
beta_out[3] + beta_out[6]))
# odhady v modelu s odstranenymi outliery
beta_red <- mod_lm_red |> coef()
estimates <- rbind(estimates,
c(beta_red[1], beta_red[2], beta_red[3],
beta_red[1] + beta_red[4], beta_red[2] + beta_red[5],
beta_red[3] + beta_red[6]))
rownames(estimates) <- c('Skutečné hodnoty', 'LM s outliery', 'LM bez outlierů')| \(\beta_0^{(male)}\) | \(\beta_1^{(male)}\) | \(\beta_2^{(male)}\) | \(\beta_0^{(female)}\) | \(\beta_1^{(female)}\) | \(\beta_2^{(female)}\) | |
|---|---|---|---|---|---|---|
| Skutečné hodnoty | 50.0000 | 9.0000 | -0.1000 | 78.0000 | 6.0000 | -0.0500 |
| LM s outliery | 26.7907 | 13.5549 | -0.3065 | 65.9193 | 8.6671 | -0.1863 |
| LM bez outlierů | 47.3872 | 9.5238 | -0.1266 | 78.1811 | 5.7040 | -0.0286 |
Robustní statistické metody
Vzhledem k porušeným předpokladům klasického lineárního regresního modelu přistupme k alternativním metodám, tedy k robustním metodám pro odhady v lineárním regresním modelu. Těch existuje celá řada, nás bude zajímat i jejich porovnání mezi sebou, proto uvažujeme následující metody:
- kvantilová regrese,
- Huberův M-odhad,
- Huberův GM-odhad,
- nejmenší useknuté čtverce a
- nejmenší medián čtverců.
Nejprve pomocí všech výše zmíněných metod odhadneme vektor regresních koeficientů a vykreslíme pro každou metodu zvlášť odhad regresních parabol pro muže a ženy zvlášť. Následně metody mezi sebou porovnáme a také je porovnáme s klasickými metodami.
Kvantilová regrese
Nejprve se podívejme na kvantilovou regresi, přičemž budeme uvažovat \(\alpha = 0,\!5\), tedy regresní medián (least absolute deviation). Ten je definován následovně: \[ \widehat{\boldsymbol\beta}_{0,5} = \text{argmin}_{{\boldsymbol\beta}} \sum_{i = 1}^n |Y_i - \boldsymbol x_i^\top {\boldsymbol\beta}|. \] Tento odhad je robustní v odchylkách \(y\), ale ne v odchylkách \(x\). Proto je jeho bod selhání nulový.
V R pro výpočet regresního mediánu využijeme funkci
rq() s parametrem tau = 0.5.
Můžeme se podívat na odhady jednotlivých regresních koeficientů a také si vykreslit odhadnuté regresní paraboly. Na obrázku níže můžeme vidět, že ačkoli byl model natrénovaný na celkovém počtu pozorování včetně odlehlých pozorování, odhady jsou výrazně lepší v porovnání s klasickou lineární regresí. Na obrázku jsou přikresleny průsvitnější čarou skutečné paraboly, ze kterých byla data generována.
| coefficients | lower bd | upper bd | |
|---|---|---|---|
| (Intercept) | 48.8397276 | 35.4224421 | 52.3056607 |
| x | 9.3858089 | 8.6109549 | 11.9115663 |
| I(x^2) | -0.1244017 | -0.2247797 | -0.0802067 |
| PohlaviZena | 26.2739897 | 22.5877810 | 39.6723116 |
| x:PohlaviZena | -2.6659646 | -4.4978943 | -2.0744122 |
| I(x^2):PohlaviZena | 0.0350744 | -0.0425463 | 0.0794125 |
xx <- seq(min(df$x), max(df$x), length = 501)
df_plot <- data.frame(x = c(xx, xx),
y = c(predict(mod_LAD, newdata = data.frame(x = xx,
Pohlavi = 'Muz')),
predict(mod_LAD, newdata = data.frame(x = xx,
Pohlavi = 'Zena'))),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi))
data.frame(x = c(xx, xx),
y = c(male_lm(xx), female_lm(xx)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi),
Metoda = 'Skutecnost') |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.3) +
geom_line(data = df_plot, alpha = 0.9, linewidth = 0.8) +
theme_bw() +
geom_point(data = df, aes(x, y, colour = Pohlavi), alpha = 0.25) +
geom_point(data = df_reduced, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3"))
M-odhad s Huberovou funkcí \(\psi\)
Další možností je využít M-odhad parametru \(\boldsymbol \beta\), konkrétně s Huberovou
funkcí \(\psi\). Tento odhad je
definován jako řešení soustavy \[
\sum_{i = 1}^n \psi\left(Y_i - \boldsymbol x_i^\top
{\boldsymbol\beta}\right) \boldsymbol x_i=0.
\] Bod selhání tohoto odhadu je nula. V R pro
výpočet M-odhadu využijeme funkci rlm() s parametrem
method = 'M' a scale.est = 'Huber'.
Opět se můžeme podívat na odhady jednotlivých regresních koeficientů a také si vykreslit odhadnuté regresní paraboly. Na obrázku níže můžeme vidět odhady, které jsou výrazně lepší v porovnání s klasickou lineární regresí.
| Value | Std. Error | t value | |
|---|---|---|---|
| (Intercept) | 45.5474737 | 3.1097937 | 14.646461 |
| x | 9.8269348 | 0.4913386 | 20.000330 |
| I(x^2) | -0.1392501 | 0.0186087 | -7.483056 |
| PohlaviZena | 28.6273128 | 3.8206253 | 7.492834 |
| x:PohlaviZena | -2.9806622 | 0.6326192 | -4.711621 |
| I(x^2):PohlaviZena | 0.0459044 | 0.0253744 | 1.809085 |
xx <- seq(min(df$x), max(df$x), length = 501)
df_plot <- data.frame(x = c(xx, xx),
y = c(predict(mod_Hub, newdata = data.frame(x = xx,
Pohlavi = 'Muz')),
predict(mod_Hub, newdata = data.frame(x = xx,
Pohlavi = 'Zena'))),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi))
data.frame(x = c(xx, xx),
y = c(male_lm(xx), female_lm(xx)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi),
Metoda = 'Skutecnost') |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.3) +
geom_line(data = df_plot, alpha = 0.9, linewidth = 0.8) +
theme_bw() +
geom_point(data = df, aes(x, y, colour = Pohlavi), alpha = 0.25) +
geom_point(data = df_reduced, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3"))
GM-odhad s Huberovou funkcí \(\psi\)
Zobecněním předchozí metody jsou GM-odhady parametru \(\boldsymbol \beta\), konkrétně opět využijeme Huberovu funkci \(\psi\). Tento odhad je definován jako řešení soustavy \[ \widehat{\boldsymbol\beta} = \text{argmin}_{{\boldsymbol\beta \in \mathbb R^p}} \sum_{i = 1}^n \rho\left( Y_i - \boldsymbol x_i^\top {\boldsymbol\beta}\right) \cdot w(\boldsymbol x_i), \] kde \(w:\mathbb R^p\rightarrow[0,\infty)\) je váhová funkce. Bod selhání tohoto odhadu je \(\varepsilon^* \leq \frac{1}{p+1}\).
V R pro výpočet GM-odhadu využijeme funkci
rlm() s parametrem method = 'M' a
scale.est = 'Huber'. Narozdíl od předchozího M-odhadu ale
nastavíme i parametr weights =, který nese informaci o
váhách. Volba vah není jednoznačná, po vyzkoušení několika možností se
zdají být nejlepší váhy ve tvaru \[
w(\boldsymbol x_i) = \frac{1}{\|\boldsymbol x_i\|}.
\]
w1 <- 1 / hat(df$x)
w2 <- sqrt(1 - hat(df$x))
w3 <- c()
for(i in 1:length(df$x)) {
w3 <- c(w3,
1 / sqrt(sum(model.matrix(mod_lm_out)[i, ]^2)))
}
w4 <- 1 / abs(df$x - mean(df$x)) # nedava dobre vysledky
w5 <- 1 / (df$x - mean(df$x))^2 # nedava dobre vysledky
mod_GM <- rlm(y ~ (x + I(x^2)) * Pohlavi, data = df, weights = w3,
method = 'M', scale.est = 'Huber')| Value | Std. Error | t value | |
|---|---|---|---|
| (Intercept) | 48.9232562 | 1.9823274 | 24.679706 |
| x | 9.3060394 | 0.3800257 | 24.487924 |
| I(x^2) | -0.1207209 | 0.0167996 | -7.185954 |
| PohlaviZena | 26.1972475 | 2.0555593 | 12.744584 |
| x:PohlaviZena | -2.5953575 | 0.4162222 | -6.235510 |
| I(x^2):PohlaviZena | 0.0312939 | 0.0200962 | 1.557210 |
xx <- seq(min(df$x), max(df$x), length = 501)
df_plot <- data.frame(x = c(xx, xx),
y = c(predict(mod_GM, newdata = data.frame(x = xx,
Pohlavi = 'Muz')),
predict(mod_GM, newdata = data.frame(x = xx,
Pohlavi = 'Zena'))),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi))
data.frame(x = c(xx, xx),
y = c(male_lm(xx), female_lm(xx)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi),
Metoda = 'Skutecnost') |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.3) +
geom_line(data = df_plot, alpha = 0.9, linewidth = 0.8) +
theme_bw() +
geom_point(data = df, aes(x, y, colour = Pohlavi), alpha = 0.25) +
geom_point(data = df_reduced, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3"))
Z obrázku i hodnot regresních koeficientů vidíme, že tato metoda je (alespoň opticky) zatím nejblíže skutečným regresním parabolám.
Nejmenší useknuté čtverce (LTS)
Předposlední metodou, kterou použijeme k odhadu parametru regresních koeficientů, jsou nejmenší useknuté čtverce (least trimmed squares). Odhad touto metodou je definován jako řešení optimalizační úlohy \[ \widehat{\boldsymbol\beta} = \text{argmin}_{{\boldsymbol\beta \in \mathbb R^p}} \sum_{i = 1}^n r_{(i)}^2\left({\boldsymbol\beta}\right), \] kde \(r_{i}\left({\boldsymbol\beta}\right) = Y_i - \boldsymbol x_i^\top\boldsymbol\beta\) a \(r_{(1)}^2\left({\boldsymbol\beta}\right) \leq r_{(2)}^2\left({\boldsymbol\beta}\right) \leq \dots \leq r_{(n)}^2\left({\boldsymbol\beta}\right)\). Tato metoda dosahuje bodu selhání 0,5.
V R pro výpočet LTS využijeme funkci lqs()
s parametrem method = 'lts'.
xx <- seq(min(df$x), max(df$x), length = 501)
df_plot <- data.frame(x = c(xx, xx),
y = c(predict(mod_LTS, newdata = data.frame(x = xx,
Pohlavi = 'Muz')),
predict(mod_LTS, newdata = data.frame(x = xx,
Pohlavi = 'Zena'))),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi))
data.frame(x = c(xx, xx),
y = c(male_lm(xx), female_lm(xx)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi),
Metoda = 'Skutecnost') |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.3) +
geom_line(data = df_plot, alpha = 0.9, linewidth = 0.8) +
theme_bw() +
geom_point(data = df, aes(x, y, colour = Pohlavi), alpha = 0.25) +
geom_point(data = df_reduced, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3"))
Z grafu výše vidíme, že odhady regresních parabol jsou velmi blízko skutečným křivkám.
Nejmenší medián čtverců (LMS)
Poslední metodou je nejmenší medián čtverců (least median of
squares). Odhad touto metodou je definován jako řešení
minimalizační úlohy \[
\widehat{\boldsymbol\beta} = \text{argmin}_{{\boldsymbol\beta \in
\mathbb R^p}} \ \text{median}\big\{(Y_1 - \boldsymbol
x_1^\top\boldsymbol\beta)^2, (Y_2 - \boldsymbol
x_2^\top\boldsymbol\beta)^2, \dots, (Y_n - \boldsymbol
x_n^\top\boldsymbol\beta)^2\big\}.
\] Tato metoda dosahuje také bodu selhání 0,5. V R
pro výpočet LMS využijeme funkci lqs() s parametrem
method = 'lms'.
xx <- seq(min(df$x), max(df$x), length = 501)
df_plot <- data.frame(x = c(xx, xx),
y = c(predict(mod_LMS, newdata = data.frame(x = xx,
Pohlavi = 'Muz')),
predict(mod_LMS, newdata = data.frame(x = xx,
Pohlavi = 'Zena'))),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi))
data.frame(x = c(xx, xx),
y = c(male_lm(xx), female_lm(xx)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi),
Metoda = 'Skutecnost') |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.3) +
geom_line(data = df_plot, alpha = 0.9, linewidth = 0.8) +
theme_bw() +
geom_point(data = df, aes(x, y, colour = Pohlavi), alpha = 0.25) +
geom_point(data = df_reduced, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3"))
Porovnání robustních odhadů a klasického LRM
Konečně se podívejme na porovnání jednotlivých rubustních metod mezi sebou a také s klasickým lineárním regresním modelem s a bez odstranění outlierů. Nejprve si spočtěme odhady regresních koeficientů, resp. jejich přepočet na hodnoty, které byly použity v definici generujících funkcí.
# odhady v modelu kvantilove regrese
beta <- mod_LAD |> coef()
estimates <- rbind(estimates,
c(beta[1], beta[2], beta[3],
beta[1] + beta[4], beta[2] + beta[5],
beta[3] + beta[6]))
# odhady v modelu Huberova M-odhadu
beta <- mod_Hub |> coef()
estimates <- rbind(estimates,
c(beta[1], beta[2], beta[3],
beta[1] + beta[4], beta[2] + beta[5],
beta[3] + beta[6]))
# odhady v modelu Huberova GM-odhadu
beta <- mod_GM |> coef()
estimates <- rbind(estimates,
c(beta[1], beta[2], beta[3],
beta[1] + beta[4], beta[2] + beta[5],
beta[3] + beta[6]))
# odhady v modelu LTS
beta <- mod_LTS |> coef()
estimates <- rbind(estimates,
c(beta[1], beta[2], beta[3],
beta[1] + beta[4], beta[2] + beta[5],
beta[3] + beta[6]))
# odhady v modelu LMS
beta <- mod_LMS |> coef()
estimates <- rbind(estimates,
c(beta[1], beta[2], beta[3],
beta[1] + beta[4], beta[2] + beta[5],
beta[3] + beta[6]))
# odhady v modelu s odstranenymi outliery
rownames(estimates)[4:8] <- c('Kvantilová regrese',
"Huberův M-odhad",
"Huberův GM-odhad",
"Nejmenší useknuté čtverce",
"Nejmenší medián čtverců")| \(\beta_0^{(male)}\) | \(\beta_1^{(male)}\) | \(\beta_2^{(male)}\) | \(\beta_0^{(female)}\) | \(\beta_1^{(female)}\) | \(\beta_2^{(female)}\) | |
|---|---|---|---|---|---|---|
| Skutečné hodnoty | 50.0000 | 9.0000 | -0.1000 | 78.0000 | 6.0000 | -0.0500 |
| LM s outliery | 26.7907 | 13.5549 | -0.3065 | 65.9193 | 8.6671 | -0.1863 |
| LM bez outlierů | 47.3872 | 9.5238 | -0.1266 | 78.1811 | 5.7040 | -0.0286 |
| Kvantilová regrese | 48.8397 | 9.3858 | -0.1244 | 75.1137 | 6.7198 | -0.0893 |
| Huberův M-odhad | 45.5475 | 9.8269 | -0.1393 | 74.1748 | 6.8463 | -0.0933 |
| Huberův GM-odhad | 48.9233 | 9.3060 | -0.1207 | 75.1205 | 6.7107 | -0.0894 |
| Nejmenší useknuté čtverce | 52.0929 | 8.6637 | -0.0858 | 73.7700 | 6.8101 | -0.0862 |
| Nejmenší medián čtverců | 52.1976 | 8.6637 | -0.0858 | 73.8747 | 6.8101 | -0.0862 |
Z tabulky výše můžeme vidět, že v porovnání s klasickým lineárním modelem s outliery jsou všechny ostatní odhady výrazně blíže skutečnosti. Nejlépe vychází nejspíše GM-odhad s Huberovou funkcí. Můžeme si také vykreslit všechny odhady do jednoho grafu.
xx <- seq(min(df$x), max(df$x), length = 501)
df_plot <- data.frame(
x = c(xx, xx),
y = c(predict(mod_lm_out, newdata = data.frame(x = xx, Pohlavi = 'Muz')),
predict(mod_lm_out, newdata = data.frame(x = xx, Pohlavi = 'Zena'))),
Pohlavi = factor(rep(c('Muz', 'Zena'), each = length(xx))),
Model = 'LM s outliery') %>%
bind_rows(data.frame(
x = c(xx, xx),
y = c(predict(mod_lm_red, newdata = data.frame(x = xx, Pohlavi = 'Muz')),
predict(mod_lm_red, newdata = data.frame(x = xx, Pohlavi = 'Zena'))),
Pohlavi = factor(rep(c('Muz', 'Zena'), each = length(xx))),
Model = 'LM bez outlierů')) %>%
bind_rows(data.frame(
x = c(xx, xx),
y = c(predict(mod_LAD, newdata = data.frame(x = xx, Pohlavi = 'Muz')),
predict(mod_LAD, newdata = data.frame(x = xx, Pohlavi = 'Zena'))),
Pohlavi = factor(rep(c('Muz', 'Zena'), each = length(xx))),
Model = 'Kvantilová regrese')) %>%
bind_rows(data.frame(
x = c(xx, xx),
y = c(predict(mod_Hub, newdata = data.frame(x = xx, Pohlavi = 'Muz')),
predict(mod_Hub, newdata = data.frame(x = xx, Pohlavi = 'Zena'))),
Pohlavi = factor(rep(c('Muz', 'Zena'), each = length(xx))),
Model = 'Huberův M-odhad')) %>%
bind_rows(data.frame(
x = c(xx, xx),
y = c(predict(mod_GM, newdata = data.frame(x = xx, Pohlavi = 'Muz')),
predict(mod_GM, newdata = data.frame(x = xx, Pohlavi = 'Zena'))),
Pohlavi = factor(rep(c('Muz', 'Zena'), each = length(xx))),
Model = 'Huberův GM-odhad')) %>%
bind_rows(data.frame(
x = c(xx, xx),
y = c(predict(mod_LTS, newdata = data.frame(x = xx, Pohlavi = 'Muz')),
predict(mod_LTS, newdata = data.frame(x = xx, Pohlavi = 'Zena'))),
Pohlavi = factor(rep(c('Muz', 'Zena'), each = length(xx))),
Model = 'Nejmenší useknuté čtverce')) %>%
bind_rows(data.frame(
x = c(xx, xx),
y = c(predict(mod_LMS, newdata = data.frame(x = xx, Pohlavi = 'Muz')),
predict(mod_LMS, newdata = data.frame(x = xx, Pohlavi = 'Zena'))),
Pohlavi = factor(rep(c('Muz', 'Zena'), each = length(xx))),
Model = 'Nejmenší medián čtverců')) %>%
mutate(Model = factor(Model,
levels = c('LM s outliery', 'LM bez outlierů',
'Kvantilová regrese', 'Huberův M-odhad',
'Huberův GM-odhad',
'Nejmenší useknuté čtverce',
'Nejmenší medián čtverců'),
ordered = T))
data.frame(x = c(xx, xx),
y = c(male_lm(xx), female_lm(xx)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi)) |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.3) +
geom_line(data = df_plot, aes(group = interaction(Model, Pohlavi),
linetype = Model),
alpha = 0.9, linewidth = 0.5) +
theme_bw() +
geom_point(data = df, aes(x, y, colour = Pohlavi), alpha = 0.25) +
geom_point(data = df_reduced, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3"))
Takový graf je poměrně nepřehledný, proto si přibližme pouze část nahoře vpravo, kde se odhady liší nejvíce.
data.frame(x = c(xx, xx),
y = c(male_lm(xx), female_lm(xx)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi)) |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.3) +
geom_line(data = df_plot, aes(group = interaction(Model, Pohlavi),
linetype = Model),
alpha = 0.9, linewidth = 0.5) +
theme_bw() +
geom_point(data = df, aes(x, y, colour = Pohlavi), alpha = 0.25) +
geom_point(data = df_reduced, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3")) +
ylim(c(145, 190)) +
xlim(c(13, 21))
Vidíme, že nejmenší useknuté čtverce téměř přesně kopírují skutečnou parabolu pro muže (modrá průsvitná plná čára). Pro ženy je většina metod o něco méně úspěšná v odhadu regresní paraboly, avšak všechny si vedou podstatně lépe než klasická lineární regrese (plná čára).
Protože je však graf stále málo přehledný, vykresleme si i každou metodu do grafu zvlášť.
data.frame(x = c(xx, xx),
y = c(male_lm(xx), female_lm(xx)),
Pohlavi = rep(c('Muz', 'Zena'), each = length(xx))) |>
mutate(Pohlavi = factor(Pohlavi)) |>
ggplot(aes(x = x, y = y, colour = Pohlavi)) +
geom_line(alpha = 0.3) +
geom_line(data = df_plot |> filter(Model != 'LM s outliery'),
aes(group = interaction(Model, Pohlavi), linetype = Model),
alpha = 0.9, linewidth = 0.5) +
theme_bw() +
geom_point(data = df, aes(x, y, colour = Pohlavi), alpha = 0.25) +
geom_point(data = df_reduced, aes(x, y, colour = Pohlavi)) +
scale_color_manual(values=c("deepskyblue2", "tomato3")) +
facet_wrap(~Model) +
theme(legend.position = 'none')
Simulační studie
Všechny předchozí závěry jsme dělali na základě pouze jednoho vygenerovaného datového souboru. Proto tyto závěry nemají velkou vypovídající hodnotu a chceme generování dat a celou analýzu zopakovat pro více opakování, přičemž závěry nyní uděláme na základě průměrů odhadů, resp. na základě střední čtvercové chyby (mean squared error, MSE).
Nejprve si definujeme funkce simulace(),
get_estimates() a norma2(), které nám
zpřehlední kód. Střední čtvercová chyba pro mnohorozměrný parametr \(\boldsymbol \beta\) a jeho odhad \(\widehat{\boldsymbol \beta}\) se definuje
jako \[
MSE(\widehat{\boldsymbol \beta},\boldsymbol \beta) = \frac{1}{N}\sum_{i
= 1}^N \Big\|\widehat{\boldsymbol \beta} - \boldsymbol \beta\Big\|^2_2.
\]
V naší situaci však není tato definice optimální, neboť každá ze složek \(\boldsymbol \beta\) nabývá hodnot na jiné škále (od 0,01 po 80), proto při použití definice uvedené výše budou odchylky v jednotlivých složkách neporovnatelné. Proto budeme jednotlivé složky rozdílu \(\widehat{\boldsymbol \beta} - \boldsymbol \beta\) škálovat příslušnými složkami \(\boldsymbol \beta\). Celkem tedy \[ \widetilde{MSE}(\widehat{\boldsymbol \beta},\boldsymbol \beta) = \frac{1}{N}\sum_{i = 1}^N \left\|\frac{\widehat{\boldsymbol \beta} - \boldsymbol \beta}{\boldsymbol \beta}\right\|^2_2. \]
# simulace dat
simulace <- function(mu_male = 12.5, mu_female = 11.5,
n_male = 85, n_female = 94, sigma_x = 3.5, df_y = 2,
center_male = c(21, 140), center_female = c(20, 145),
n_out_male = 2, n_out_female = 3) {
x_male <- rnorm(n_male, mean = mu_male, sd = sigma_x)
x_female <- rnorm(n_female, mean = mu_female, sd = sigma_x)
y_male <- male_lm(x_male) + rt(n_male, df_y) * 2
y_female <- female_lm(x_female) + rt(n_female, df_y) * 2
data <- data.frame(
x = c(x_male, x_female),
y = c(y_male, y_female),
Pohlavi = rep(c('Muz', 'Zena'), c(n_male, n_female))
) |>
mutate(Pohlavi = factor(Pohlavi))
Sigma <- matrix(c(0.5, -0.1, -0.1, 10), nrow = 2, byrow = F)
df <- rbind(mvrnorm(n = n_out_male, mu = center_male, Sigma = Sigma),
mvrnorm(n = n_out_female, mu = center_female, Sigma = Sigma)) |> data.frame()
colnames(df) <- c('x', 'y')
df <- df |> mutate(Pohlavi = rep(levels(data$Pohlavi), c(n_out_male, n_out_female)) |> factor())
df <- rbind(data, df)
return(df)
}get_estimates <- function(df) {
# nafitujeme modely
mod_lm_out <- lm(y ~ (x + I(x^2)) * Pohlavi, data = df)
df_reduced <- df[-c(180:184), ]
mod_lm_red <- lm(y ~ (x + I(x^2)) * Pohlavi, data = df_reduced)
mod_LAD <- rq(y ~ (x + I(x^2)) * Pohlavi, data = df, tau = 0.5)
mod_Hub <- rlm(y ~ (x + I(x^2)) * Pohlavi, data = df,
method = 'M', scale.est = 'Huber')
w3 <- c()
for(i in 1:length(df$x)) {
w3 <- c(w3,
1 / sqrt(sum(model.matrix(mod_lm_out)[i, ]^2)))
}
mod_GM <- rlm(y ~ (x + I(x^2)) * Pohlavi, data = df, weights = w3,
method = 'M', scale.est = 'Huber')
mod_LTS <- lqs(y ~ (x + I(x^2)) * Pohlavi, data = df, method = 'lts')
mod_LMS <- lqs(y ~ (x + I(x^2)) * Pohlavi, data = df, method = 'lms')
# urcime odhady
estimates <- data.frame(beta_0_male = 50, beta_1_male = 9,
beta_2_male = -0.1, beta_0_female = 78,
beta_1_female = 6, beta_2_female = -0.05)
beta_out <- mod_lm_out |> coef()
estimates <- rbind(estimates,
c(beta_out[1], beta_out[2], beta_out[3],
beta_out[1] + beta_out[4], beta_out[2] + beta_out[5],
beta_out[3] + beta_out[6]))
beta_red <- mod_lm_red |> coef()
estimates <- rbind(estimates,
c(beta_red[1], beta_red[2], beta_red[3],
beta_red[1] + beta_red[4], beta_red[2] + beta_red[5],
beta_red[3] + beta_red[6]))
rownames(estimates) <- c('Skutečné hodnoty', 'LM s outliery', 'LM bez outlierů')
beta <- mod_LAD |> coef()
estimates <- rbind(estimates,
c(beta[1], beta[2], beta[3],
beta[1] + beta[4], beta[2] + beta[5],
beta[3] + beta[6]))
beta <- mod_Hub |> coef()
estimates <- rbind(estimates,
c(beta[1], beta[2], beta[3],
beta[1] + beta[4], beta[2] + beta[5],
beta[3] + beta[6]))
beta <- mod_GM |> coef()
estimates <- rbind(estimates,
c(beta[1], beta[2], beta[3],
beta[1] + beta[4], beta[2] + beta[5],
beta[3] + beta[6]))
beta <- mod_LTS |> coef()
estimates <- rbind(estimates,
c(beta[1], beta[2], beta[3],
beta[1] + beta[4], beta[2] + beta[5],
beta[3] + beta[6]))
beta <- mod_LMS |> coef()
estimates <- rbind(estimates,
c(beta[1], beta[2], beta[3],
beta[1] + beta[4], beta[2] + beta[5],
beta[3] + beta[6]))
rownames(estimates)[4:8] <- c('Kvantilová regrese',
"Huberův M-odhad",
"Huberův GM-odhad",
"Nejmenší useknuté čtverce",
"Nejmenší medián čtverců")
return(estimates)
}norma2 <- function(est, real = c(50, 9, -0.1, 78, 6, -0.05)) {
norma <- norm(as.matrix((est - real)/real), type = 'f')^2
return(norma)
}Zvolme N = 1000 a set.seed(21). Nejprve
vždy nasimulejeme data pomocí funkce simulace() a následně
na její výstupní argument zavoláme funkci get_estimates(),
ketrá nám vrátí tabulku s hodnotami odhadů regresních koeficientů pro
jednotlivé metody.
N <- 1000
set.seed(21)
results <- array(NA, dim = c(7, 6, N))
for (i in 1:N) {
df <- simulace()
estimates <- get_estimates(df)
results[, , i] <- estimates[-1, ] |> as.matrix()
}
Mean <- apply(results, c(1, 2), mean)
Mean <- rbind(c(50, 9, -0.1, 78, 6, -0.05), Mean)
rownames(Mean) <- rownames(estimates)V tabulce níže jsou uvedeny průměry odhadů pro jednotlivé metody a jednotlivé složky vektoru regresních parametrů. V prvním řádku jsou pak připomenuty skutečné hodnoty parametru. Vidíme, že až na klasický lineární model s odlehlými pozorováními všechny metody dávají slušné odhady. Nejlépe vychází tyto metody:
- lineární model s odstraněnými outliery,
- nejmenší useknuté čtverce a
- nejmenší medián čtverců.
| \(\beta_0^{(male)}\) | \(\beta_1^{(male)}\) | \(\beta_2^{(male)}\) | \(\beta_0^{(female)}\) | \(\beta_1^{(female)}\) | \(\beta_2^{(female)}\) | |
|---|---|---|---|---|---|---|
| Skutečné hodnoty | 50.0000 | 9.0000 | -0.1000 | 78.0000 | 6.0000 | -0.0500 |
| LM s outliery | 27.1132 | 13.5192 | -0.3053 | 65.1709 | 8.8657 | -0.1942 |
| LM bez outlierů | 49.8413 | 9.0259 | -0.1010 | 78.2312 | 5.9587 | -0.0483 |
| Kvantilová regrese | 47.9368 | 9.4104 | -0.1188 | 75.4674 | 6.5507 | -0.0774 |
| Huberův M-odhad | 46.6060 | 9.6664 | -0.1302 | 74.1939 | 6.8394 | -0.0920 |
| Huberův GM-odhad | 48.1698 | 9.4069 | -0.1205 | 76.2962 | 6.4536 | -0.0763 |
| Nejmenší useknuté čtverce | 49.6165 | 9.0408 | -0.1008 | 78.0142 | 5.9932 | -0.0496 |
| Nejmenší medián čtverců | 50.0488 | 9.0010 | -0.1004 | 77.7704 | 6.0380 | -0.0514 |
Nakonec si ještě spočítejme hodnoty MSE pro jednotlivé metody.
MSE_df <- data.frame(MSE = rep(NA, 7),
row.names = rownames(estimates)[-1])
for(method in 1:7) {
MSE_df[method, 1] <- apply(results[method, , ], 2, norma2) |> mean()
}| \(MSE\) | |
|---|---|
| LM s outliery | 15.1704 |
| LM bez outlierů | 1.0477 |
| Kvantilová regrese | 0.7484 |
| Huberův M-odhad | 1.3316 |
| Huberův GM-odhad | 0.9342 |
| Nejmenší useknuté čtverce | 1.8946 |
| Nejmenší medián čtverců | 1.7187 |
Z tabulky výše vidíme, že kvantilová regrese je nejstabilnější co se týče odhadování regresních koeficientů. Dobře je na tom také Huberův GM-odhad a lineární model bez odlehlých pozorování.
Pořadové testy
Nyní se přesuňme ke druhé části projektu, kterou jsou pořadové
testy. Budeme pracovat se stejným typem dat jako v předchozím
případě, jen upravíme odlehlá pozorování tak, aby byla opravdu zřejmě
odlehlá (pozorování, která jsou odlehlá ve dvou rozměrech, se nemusí
jevit jako odlehlá v jedné dimenzi). Proto změníme střed, ze kterého
generujeme odlehlá pozorování ve funkci simulace(). Také
mírně zvýšíme jejich počet, aby byly závěry o testování hypotéz
výraznější.
set.seed(21)
data <- simulace(center_male = c(33, 140), center_female = c(33, 150),
n_out_male = 5, n_out_female = 7)Jednovýběrové pořadové testy
Začněme nejprve jednovýběrovou situací. Bude nás zajímat, zda
vygenerované hodnoty (na ose \(x\)) pro
mužskou populaci pocházejí z rozdělení s mediánem \(\Delta\) (střední hodnotou) rovným \(12,\!5\) (to je skutečná střední hodnota
normálního rozdělení, ze kterého jsme generovali hodnoty ve funkci
simulace()). Poznamenejme, že jsou splněny technické
předpoklady pro použití pořadových testů – spojité a symetrické
rozdělení kolem mediánu \(\Delta\)
(jelikož uvažujeme data z normálního rozdělení s přidáním odlehlých
pozorování). Budeme tedy testovat nulovou hypotézu \[
H_0: \Delta = \Delta_0\ \ vs\ \ H_1: \Delta \neq \Delta_0,
\] kde \(\Delta_0 =
12,\!5\).
Nejprve se na data podívejme graficky. Vykresleme si histogram \(x\)-ových hodnot pro muže superponovaný odhadnutou křivkou hustoty normálního rozdělení s parametry určenými z náhodného výběru (plná modrá čára). Přidáme také jádrový odhad hustoty (přerušovaná modrá čára).
data |> filter(Pohlavi == 'Muz') |>
ggplot(aes(x = x)) +
geom_histogram(aes(y = after_stat(density)),
fill = 'deepskyblue2', alpha = 0.5, colour = 'black',
bins = ceiling(1 + log2(length(data$x)))) +
theme_bw() +
geom_function(fun = dnorm, args = list(mean = mean(data$x), sd = sd(data$x)),
color = "darkblue", linewidth = 1.0) +
geom_density(linewidth = 0.7, color = 'darkblue', linetype = 'dashed',
adjust = 1.5)
Vidíme, že odlehlá pozorování výrazně vychylují oba odhady hustoty. Pro posouzení normality dat můžeme využít například Q-Q graf.
data |> filter(Pohlavi == 'Muz') |> ggplot(aes(sample = x)) +
stat_qq() +
theme_bw() +
geom_qq_line()
Odlehlá pozorování se opět výrazně odlišují od zbytku dat a tudíž od normality. Bez těchto pozorování bychom správně normalitu nezamítali. Formálně si můžeme předpoklad normality otestovat například pomocí Shapirova-Wilkova testu normality.
print(paste0('Shapiro-Wilk normality test: p-value = ',
shapiro.test(data$x[data$Pohlavi == 'Muz'])$p.value |> round(10)))## [1] "Shapiro-Wilk normality test: p-value = 2.75e-08"Jeho p-hodnota je výrazně menší než jakákoli rozumně zvolená hladina významnosti, proto zamítáme hypotézu o normalitě náhodného výběru. Nemůžeme tedy testovat hypotézu \(H_0\) pomocí t-testu, můžeme však využít pořadových testů, které normalitu nepředpokládají.
Jednou z možností je využít jednovýběrový Wilcoxonův
test. V R využijeme funkci
wilcox.test() s parametry exact = TRUE a
correct = TRUE. Výsledná p-hodnota je rovna 0,07, na
hladině významnosti 0,05 tedy nezamítáme nulovou hypotézu \(H_0\) o mediánu hodnot pro mužskou populaci
rovném hodnotě 12,5. Tento závěr je v souladu se skutečností.
p_val <- wilcox.test(data$x[data$Pohlavi == 'Muz'],
mu = 12.5, exact = TRUE, correct = TRUE)$p.value
print(paste0('Wilcoxon signed rank exact test: p-value = ',
p_val |> round(5)))## [1] "Wilcoxon signed rank exact test: p-value = 0.07006"Podobně můžeme testovat hypotézu pomocí znaménkového testu. Jeho p-hodnota je rovna 0,07255, opět tedy správně nezamítáme nulovou hypotézu \(H_0\).
U <- sum(data$x[data$Pohlavi == 'Muz'] > 12.5)
n <- length(data$x[data$Pohlavi == 'Muz'])
p_val <- binom.test(U, n)$p.value
print(paste0('Exact binomial test: p-value = ',
p_val |> round(5)))## [1] "Exact binomial test: p-value = 0.07255"p_val <- SIGN.test(data$x[data$Pohlavi == 'Muz'], md = 12.5)$p.value
print(paste0('One-sample Sign-Test: p-value = ',
p_val |> round(5)))## [1] "One-sample Sign-Test: p-value = 0.07255"Dvouvýběrové pořadové testy
Podívejme se nyní na dvouvýběrový problém, tedy budeme uvažovat kromě mužské populace také ženskou a bude nás zajímat, zda se liší jejich distribuce posunutím. Víme, že jsme generovali \(x\)-ové hodnoty z normálních rozdělení se stejnými rozptyly, ale s hodnotami \(\mu_{male} = 12,\!5\) a \(\mu_{female} = 11,\!5\). Budeme tedy testovat hypotézu \[ H_0: F(x) = G(x) \ \ vs \ \ H_1: G(x) = F(x - \Delta),\ \Delta \neq 0, \] kde \(F\), resp. \(G\) jsou distribuční funkce \(x\)-ových hodnot pro mužskou, resp. ženskou populaci. V našem případě je \(\Delta = -1\), tedy platí alternativa \(H_1\).
Nejprve si data vykresleme graficky. Histogramy zvlášť pro muže a ženy jsou na obrázku níže s přikresleným jádrovým odhadem hustoty. Vidíme, že odlehlá pozorování opět výrazně vychylují data od normality.
data |>
ggplot(aes(x = x)) +
geom_histogram(aes(y = after_stat(density)),
fill = 'deepskyblue2', alpha = 0.5, colour = 'black',
bins = ceiling(1 + log2(length(data$x[data$Pohlavi == 'Muz'])))) +
theme_bw() +
#geom_function(fun = dnorm, args = list(mean = mean(data$x), sd = sd(data$x)),
# color = "darkblue", linewidth = 1.0) +
geom_density(linewidth = 0.7, color = 'darkblue', linetype = 'solid',
adjust = 1.5) +
facet_wrap(~Pohlavi)
Stejné závěry o normalitě usoudíme i z pohledu na Q-Q grafy pro jednotlivé populace.
Formálně můžeme normalitu otestovat pomocí Shapirova-Wilkova testu normality. Pro obě populace jsou p-hodnoty téměř nulové, proto zamítáme předpoklad o normalitě obou náhodných výběrů.
## [1] "Shapiro-Wilk normality test: p-value = 2.75e-08"## [1] "Shapiro-Wilk normality test: p-value = 2.41e-11"Pro otestování hypotézy \(H_0\) využijeme opět Wilcoxonův test, nyní dvouvýběrový. Jeho p-hodnota je rovna 0,00516, zamítáme tedy nulovou hypotézu \(H_0\) ve prospěch alternativy \(H_1\), která tvrdí, že se distribuce liší posunutím. Víme, že tento závěr odpovídá skutečnosti. Nulovou hypotézu tedy zamítáme správně.
p_val <- wilcox.test(X, Y, exact = TRUE, correct = TRUE)$p.value
print(paste0('Wilcoxon rank sum exact test: p-value = ',
p_val |> round(5)))## [1] "Wilcoxon rank sum exact test: p-value = 0.00516"Porovnání s klasickými testy
Podívejme se nakonec na porovnání pořadových testů s klasickými testy. Zajímá nás, k jakým závěrům bychom dospěli, pokud bychom (nesprávně) použili klasické statistické metody, které předpokládají normalitu. Testy nulových hypotéz pro jednovýběrový i dvouvýběrový problém provedeme pomocí t-testu.
V případě jednoho výběru testujeme, zda je medián \(x\)-ových hodnot pro muže roven 12,5. V případě pořadových testů jsme nulovou hypotézu nezamítli. Klasický t-test nám dá p-hodnotu rovnu 0,01965, na hladině významnosti 0,05 bychom tedy zamítali nulovou hypotézu ve prospěch alternativy. Dopustili bychom se tedy chyby prvního typu, kdy zamítáme pravdivou nulovou hypotézu.
p_val <- t.test(X, mu = 12.5)$p.value
print(paste0('One Sample t-test: p-value = ',
p_val |> round(5)))## [1] "One Sample t-test: p-value = 0.01965"V případě dvouvýběrového problému jsme při použití pořadových testů zamítali nulovou hypotézu, která tvrdila, že distribuce \(x\)-ových hodnot jsou v mužské a ženské populaci stejné, proti alternativě, že se distribuce liší posunutím.
Při testování této hypotézy klasickým t-testem dostaneme p-hodnotu rovnu 0,1238, nezamítali bychom tedy nulovou hypotézu \(H_0\) a dopustili bychom se proto chyby druhého typu (neboť v tomto případě platí alternativa).
p_val <- t.test(X, Y, var.equal = T, paired = F)$p.value
print(paste0('Two Sample t-test: p-value = ',
p_val |> round(5)))## [1] "Two Sample t-test: p-value = 0.12381"Závěr
Cílem projektu bylo porovnat klasické statistické metody s těmi robustními, resp. neparametrickými. V první části projektu jsme se věnovali lineárními regresnímu modelu, ve druhé pak pořadovým testům a jejich porovnání s klasickým jednovýběrovým, resp. dvouvýběrovým t-testem.
Viděli jsme, že při přítomnosti odlehlých pozorování klasický lineární model může vést ke špatným závěrům o správném vztahu mezi kovariáty a závislou proměnnou. Jednou z možností, která však vyžaduje hlubší znalost problému a více práce, je odstranění odlehlých hodnot. Pomocí simulační studie jsme došli k závěru, že na naše data nejlépe fungují dvě robustní metody – nejmenší useknuté čtverce a nejmenší medián čtverců. Měli bychom mít na paměti, že ačkoli odstranění odlehlých pozorování výrazně vylepší odhady regresních koeficientů spočítané pomocí klasického lineárního modelu, stále není splněn předpoklad normality náhodných chyb a proto se nemůžeme spoléhat na závěry vyplývající z testů provedených pomocí tohoto lineárního modelu.
Ve druhé části projektu jsme ukázali, jak nesprávné použití klasických statistických metod (t-test) může vést ke zcela opačným závěrům, než jaké nám poskytují pořadové testy. Navíc závěry udělané na základě klasických metod nejsou správné a můžeme se dopustit chyby prvního, resp. druhého typu.
Viděli jsme tedy, že nesprávné použití klasických statistických metod může vést ke špatným závěrům o datech. Měli bychom tedy vždy přistupovat k předpokladům klasickým metod obezřetně a pokud data vykazují jejich porušení, použít pro analýzu, resp. testování hypotéz odpovídající robustní či neparametrickou alternativu.