library(gginference)
library(ggpval)
library(ggpubr)
library(ggsignif)
library(ggstatsplot)
library(ggdist)
library(infer)
library(ggplot2)
library(tidyverse)
library(data.table)
library(ggmosaic)4 Wizualizacja wyników wnioskowania statystycznego - wybrane biblioteki programu R
W tym rozdziale przedstawiono wybrane biblioteki programu R, które wspierają graficzną prezentację wyników wnioskowania statystycznego. Omówione pakiety rozszerzają możliwości ggplot2, umożliwiając wizualizację rezultatów testów (wartości statystyk, obszary krytyczne, \(p\)‑wartości, przedziały ufności), a także łączenie klasycznych procedur testowych z nowoczesnymi metodami symulacyjnymi. Szczególny nacisk położono na to, jak wyniki testowania hipotez można w sposób czytelny i zwięzły zaprezentować na wykresach.
W kolejnych podrozdziałach zaprezentowano zastosowania wybranych pakietów programu R do wizualizacji wyników testów: chi‑kwadrat (zgodności i niezależności), testów dla proporcji, współczynnika korelacji, testów t (dla prób niezależnych i zależnych), testów dla wariancji oraz analizy wariancji (ANOVA). Zwieńczeniem rozdziału jest przykład „wizualizacji wnioskowania” w pełni zdefiniowanej przez użytkownika – implementacja testu pustych cel jako testu wielowymiarowej normalności wraz z graficzną prezentacją rozkładu statystyki, obszaru krytycznego i decyzji testowej.
4.1 Wnioskowanie statystyczne - wybrane biblioteki programu R
W Kończak (2024) wskazano na różnorodne biblioteki pozwalające rozszerzyć możliwości pakietu graficznego ggplot2. W tej części skupiono się na pakietach, które mogą być zastosowane w graficznych prezentacjach rezultatów wnioskowania statystycznego. Wykaz pakietów, które zostały wykorzystane w tym rozdziale przedstawiono w tabeli tabela 4.1. Przedstawione w tej tabeli pakiety rozszerzeń w różny sposób wspomagają przeprowadzenie wnioskowania statystycznego. Pakiety gginference, ggpval, ggpubr, ggsignif, ggstatsplot, ggdist i infer umożliwiają przeprowadzenie wnioskowania statystycznego wzbogacając prezentacje o wyświetlenie wyników testowania hipotez na wykresie. Mogą to być na przykład wartości statystyk testowych, zaznaczone obszary krytyczne i wskazane \(p\)-wartości.
| Biblioteka | Opis |
|---|---|
| gginference | Wizualizacja wyników wybranych testów statystycznych. |
| ggpval | Przeprowadzanie testów statystycznych i wprowadzanie \(p\)-wartości na wykres. |
| ggpubr | Funkcje do tworzenia i dostosowywania gotowych do publikacji wykresów opartych na ggplot2. |
| ggsignif | Wskazanie na wykresie, czy pomiędzy grupami istnieją istotne statystycznie różnice. |
| ggstatsplot | Konstrukcja wykresów ze szczegółami wyników testów statystycznych. |
| ggdist | Zestaw funkcji zaprojektowanych specjalnie do wizualizacji rozkładów i niepewności. |
| infer | Pakiet do przeprowadzania testów statystycznych zgodnie ze składnią tidyverse . |
Na wstępie konieczne jest załadowanie bibliotek, które będą wykorzystywane w tym rozdziale. W rozdziale zostaną wykorzystane następujące biblioteki:
4.2 Wizualizacja wyników wnioskowania statystycznego
Wnioskowanie statystyczne to dział statystyki zajmujący się metodami umożliwiającymi uogólnianie wyników uzyskanych na podstawie próby na całą populację. Obejmuje ono zarówno szacowanie nieznanych parametrów populacji (np. średniej, wariancji, współczynnika korelacji), jak i testowanie hipotez dotyczących tej populacji. Kluczową cechą wnioskowania statystycznego jest uwzględnianie niepewności i błędów wynikających z tego, że badania przeprowadza się na próbie, a nie na całej populacji. W tym punkcie przedstawione będą możliwości wizualizacji wyników klasycznych testów statystycznych w oparciu o wybrane biblioteki programu R.
4.2.1 Wizualizacja w ggplot2
Bardzo często przed wykonaniem testu statystycznego konieczne jest sprawdzenie wymaganych założeń, a najczęściej jest to wymaganie dotyczące normalności rozkładu badanej zmiennej losowej. Wymóg ten pojawia się przy stosowaniu wielu testów parametrycznych. W związku z powyższym w pierwszym etapie weryfikowana jest stosowna hipoteza. Przykład graficznego wglądu w proces weryfikacji hipotezy głoszącej, że zmienna mpg (zbiór mtcars) ma rozkład normalny realizuje następujący kod:
ggplot(mtcars, aes(mpg)) +
geom_histogram(aes(y =after_stat(density)),
breaks = seq(0, 40, by = 2),
colour = "red",
fill = "white") +
stat_function(fun = dnorm, args = list(mean = mean(mtcars$mpg), sd = sd(mtcars$mpg)))+theme_minimal()
Wynik realizacji powyższego kodu został przedstawiony na rysunku rys. 4.1. Na tym rysunku poza histogramem, który jest oszacowaniem teoretycznego rozkładu gęstości badanej zmiennej umieszczona została krzywa gęstości rozkładu normalnego z parametrami, które zostały oszacowane na podstawie danych.
W badaniach statystycznych bardzo często zachodzi konieczność porównania dwóch populacji na podstawie pobranych prób losowych. Zwykle intersujące jest porównanie gęstości rozkładów, których estymatorem są odpowiednie histogramy. Wstępną analizę graficzną rozkładów dla dwóch grup umożliwia następujący kod.
df_mtcars <- mtcars
df_mtcars$group <- factor(df_mtcars$am, labels = c("Automatyczna", "Ręczna"))
bw=3
grid_data <- df_mtcars %>%
group_by(group) %>%
summarise(
mean = mean(mpg),
sd = sd(mpg),
n = n(),
.groups = "drop") %>%
mutate(
mpg = map(group, ~ seq(min(df_mtcars$mpg), max(df_mtcars$mpg), length.out = 100)),
y = pmap(list(mpg, mean, sd, n), function(x, m, s, n) dnorm(x, m, s) * bw * n)) %>%
unnest(c(mpg, y))
ggplot() +
geom_histogram(data = df_mtcars, aes(x = mpg), binwidth = bw, colour = "blue", fill = "lightblue", alpha = 0.6) +
geom_line(data = grid_data, aes(x = mpg, y = y), color = "red", linewidth = 1.2) +
facet_wrap(~ group, ncol = 1) +
theme_minimal()
Wynik realizacji powyższego kodu został przedstawiony na rysunku rys. 4.2. Zauważalne jest że wartość oczekiwana mpg jest mniejsza dla samochodów z automatyczną skrzynią biegów, a jednocześnie odchylenie standardowe mpg w tej grupie jest większe.
4.2.2 Pakiet gginference
Przedstawiona wcześniej biblioteka gginference pozwala nie tylko na wykreślenie teoretycznych postaci gęstości lub dystrybuant zmiennych losowych, ale również na przeprowadzenie testów statystycznych i zaprezentowanie odpowiednich wyników w formie graficznej. W tej części zostanie wskazanych kilka przykładów takich realizacji z wykorzystaniem tej biblioteki.
We wnioskowaniu statystycznym stosuje się różne testy odwołujące się do statystyki chi-kwadrat. Do najprostszych testów opartych na tej statystyce zaliczany jest test zgodności dla proporcji.
W badaniach pytano o preferencje koloru samochodu. Wśród 320 badanych otrzymano wyniki, które przedstawia tabela 4.2
| Kolor samochodu | Srebrny | Czarny | Biały | Niebieski | Czerwony |
|---|---|---|---|---|---|
| Liczba wskazań | 82 | 70 | 65 | 55 | 48 |
Test zgodności dla proporcji pozwala na weryfikację hipotezy głoszącej, że wartości z wyróżnionych kategorii są przyjmowane w zadanych proporcjach. Na podstawie powyższych danych weryfikowana jest hipoteza \(H_0\) głosząca, że wszystkie kolory samochodów są tak samo preferowane, wobec hipotezy alternatywnej \(H_1\), głoszącej że kolory nie są jednakowo preferowane.
x <- c(Srebrny = 82, Czarny = 70, Biały = 65, Niebieski = 55, Czerwony = 48)
ggchisqtest(chisq.test(x))
Rezultat realizacji powyższych komend przedstawiono na rys. 4.3. Na wykresie pionowa linia wskazuje wartość obserwowanej statystyki testowej \((\chi^2 = 10,906)\). Przedstawiono także wartość krytyczną \((\chi^2_{0,05} = 9,488)\) i obszar krytyczny dla tego testu (szare pole). Graficzna prezentacja wyników przy przyjętym poziomie istotności \(\alpha=0,05\) prowadzi do odrzucenia hipotezy \(H_0\), a więc można twierdzić, że kolory samochodów nie cieszą się jednakowym zainteresowaniem badanych.
Kolejnym bardzo często wykorzystywanym testem w badaniach naukowych jest test niezależności chi-kwadrat. Test ten pozwala na weryfikację hipotezy o niezależności dwóch zmiennych jakościowych, dla których zostały wyróżnione pewne kategorie. W poniżej przedstawionym przykładzie weryfikacji została poddana hipoteza \(H_0\), głosząca, że przeżycie katastrofy statku Titanic (zbiór danych Titanic) nie zależało od klasy, którą podróżowała dana osoba. Weryfikacji tej hipotezy służy następujący kod
tab_Titanic <- margin.table(Titanic, margin = c(1, 4))
ggchisqtest(chisq.test(tab_Titanic))
Na rys. 4.4 przedstawiono wizualizację wyników weryfikacji hipotezy głoszącej, że klasa, w której przebywał pasażer statku Titanic nie miała wpływu na przeżycie katastrofy. Zaznaczony na tym rysunku obszar krytyczny i wartość krytyczna \(\chi^2_{0,05} = 7,815\) oraz wartość statystyki testowej \(\chi^2 = 190,401\) prowadzą do odrzucenia hipotezy o niezależności dwóch rozważanych zmiennych przy poziomie istotności \(\alpha=0,05\). Bardzo duża wartość statystyki testowej świadczy przeciw tej hipotezie. Można twierdzić, że klasa miała wpływ na przeżycie katastrofy.
Do wizualizacji wyników testu na podstawie danych z tablic wielodzielczych można wykorzystać także funkcję mosaicplot dostępną w podstawowej wersji programu R.
plot(Titanic,shade=TRUE, main='')
Wykres mozaikowy rys. 6.38 przedstawia zależności pomiędzy klasą podróży, płcią oraz wiekiem pasażera a przeżyciem katastrofy. Każdy prostokąt na wykresie odpowiada określonej grupie pasażerów, a jego wielkość odzwierciedla liczebność tej grupy. Kolorowanie (shade=TRUE) wskazuje na istotność statystyczną odchyleń od oczekiwanej liczebności – intensywniejsze kolory oznaczają większe różnice.
Najwyższy odsetek przeżywalności dotyczył kobiet i dzieci, zwłaszcza podróżujących pierwszą i drugą klasą. Mężczyźni, szczególnie z niższych klas, mieli znacznie mniejsze szanse na przeżycie. Klasa podróży miała istotny wpływ na przeżywalność – pasażerowie pierwszej klasy przeżywali częściej niż pasażerowie trzeciej klasy. Wiek również odgrywał rolę – dzieci miały większe szanse na przeżycie niż dorośli, zwłaszcza w wyższych klasach. Wykres mozaikowy wyraźnie pokazuje, że przeżywalność na Titanicu była silnie związana z płcią, wiekiem i klasą podróży pasażerów.
Jednym z najczęściej stawianych pytań w badaniach statystycznych jest pytanie o występowanie zależności pomiędzy zmiennymi. Jeżeli dotyczy to zmiennych mierzonych na skalach mocnych, to zazwyczaj wykorzystuje się test dla współczynnika korelacji liniowej Pearsona.
Weryfikowana hipoteza \(H_0\) ma postać:
\[H_0: \rho = 0\] wobec hipotezy alternatywnej
\[H_1: \rho \ne 0.\]
Przykład zastosowania takiego testu dla danych o szerokości i długości kielicha kwiatu iris na podstawie danych ze zbioru iris przedstawia poniższy kod.
corr_test <- cor.test(iris$Sepal.Length, iris$Sepal.Width)
ggcortest(corr_test)
Na rys. 4.6 przedstawiono graficznie wyniki testu. Wizualizacja wskazuje na brak podstaw do odrzucenia hipotezy \(H_0\) głoszącej, że nie ma zależności pomiędzy długością i szerokością kielicha kwiatu iris (dane ze zbioru iris).
Kolejny przykład przedstawia zastosowanie testu dla równości proporcji w trzech populacjach. Weryfikowana jest hipoteza:
\[H_0: p_1 = p_2 = p_3\]
wobec hipotezy alternatywnej
\[H_1: \text{nie wszystkie proporcje są jednakowe }\]
m <- c(52, 25, 50)
n <- c(80, 36, 64)
pr_test <- prop.test(m, n)
ggproptest(pr_test)
Wynik realizacji testu dla równości proporcji został przedstawiony na rysunku rys. 4.7. Widoczne jest, że wartość statystyki testowej \((\chi^2=2,975)\) jest znacznie mniejsza od wartości krytycznej \((\chi^2_\alpha=5,991)\), co prowadzi do stwierdzenia braku podstaw do odrzucenia hipotezy zerowej, głoszącej że proporcje w trzech populacjach są jednakowe.
Do najczęściej stosowanych testów statystycznych należy zaliczyć test \(t\)-Studenta dla równości dwóch wartości oczekiwanych. Przykład zastosowania takiego testu, pozwalającego na weryfikację hipotezy, głoszącej, że przeciętna liczba mil przejechanych na galonie paliwa jest taka sama w samochodach z automatyczną i ręczną skrzynią biegów (zbiór danych mtcars) przedstawia kolejny kod.
t_test <- t.test(mtcars$mpg ~ mtcars$am)
ggttest(t_test)
Wynik tego testu został przedstawiony na rys. 4.8. Wartość statystyki testowej \((t = -3,761)\) znajduje się w obszarze krytycznym, a więc można twierdzić, że liczba mil przejechanych na galonie paliwa w samochodach z ręczną i automatyczną skrzynią biegów nie jest jednakowa.
Inną formą testu \(t\)-Studenta dla dwóch prób jest test dla danych sparowanych (skojarzonych). Przykład zastosowania takiego testu przedstawia poniższy kod.
t_test2 <- t.test(x = rnorm(10), y = rnorm(10,1,1), paired=TRUE)
ggttest(t_test2)
Wynik realizacji testu \(t\)-Studenta dla danych sparowanych został przedstawiony na rysunku rys. 4.9. Wartość statystyki testowej \((t = -1,605)\) nie znajduje się w obszarze krytycznym, a więc brak podstaw do odrzucenia hipotezy \(H_0\).
Niekiedy podczas prowadzenia wnioskowania statystycznego zachodzi potrzeba sprawdzenia równości wariancji w dwóch populacjach. Przykład realizacji tego testu dla umownych danych przedstawia następujący kod:
x <- rnorm(50, mean = 0, sd = 2)
y <- rnorm(30, mean = 1, sd = 1)
var_test <- var.test(x, y)
ggvartest(var_test)
Wynik realizacji kodu dla weryfikacji hipotezy o równości dwóch wariancji został przedstawiony na rysunku rys. 4.10. Wartość statystyki testowej \((F = 2,07)\) znajduje się w obszarze krytycznym. W związku z tym można twierdzić, że wariancje w obu populacjach nie są jednakowe.
Dotychczas przedstawione testy pozwalały na porównanie dwóch populacji. W badaniach statystycznych takie porównania często odnoszą się do większej liczby populacji. Porównanie wartości oczekiwanych dla trzech lub większej liczby populacji można przeprowadzić z wykorzystaniem testu analizy wariancji (ANOVA). Przykład kodu realizującego taki test pozwalającego weryfikować hipotezę, głoszącą, że przeciętna liczba przejechanych mil na jednym galonie paliwa (mpg) jest jednakowa w trzech grupach samochodów wyróżnionych ze względu na liczbę biegów (gear), przedstawiono poniżej
mt_aov <- aov(mpg~factor(gear), data = mtcars)
ggaov(mt_aov, colaccept = "grey89", colreject = "black")
Wynik realizacji powyższego kodu został przedstawiony na rysunku rys. 4.11. Wartość statystyki testowej \((F=10,9007)\) znajduje się w obszarze krytycznym i na tej podstawie można stwierdzić, że przeciętna liczba mil przejechanych na jednym galonie paliwa (mpg) nie jest jednakowa dla samochodów o różnej liczbie biegów (gear).
4.2.3 Pakiet ggpval
Kolejną biblioteką umożliwiającą nie tylko przeprowadzenie testów statystycznych ale również graficzną prezentację wyników takich testów jest biblioteka ggpval. Biblioteka ta pozwala na czytelną wizualizację występujących różnic pomiędzy porównywalnymi populacjami.
Do najczęściej stosowanych testów statystycznych należy test \(t\)-Studenta dla dwóch prób. Pakiet ggpval umożliwia wizualizację wyników tego testu. Możliwości te zostaną przedstawione na przykładzie danych umownych, które są utworzone w następujacy sposób:
A <- rnorm(40, 0, 3)
B <- rnorm(40, 2, 4)
G <- rep(c("G1", "G2"), each = 20)
dt <- data.table(A, B, G)
dt <- melt(dt, id.vars = "G")Prezentację graficzną dwóch zbiorów wraz z zaznaczeniem różnic przedstawia poniższy kod:
p <- ggplot(dt, aes(variable, value)) +
geom_boxplot() +
geom_jitter()
add_pval(p, pairs = list(c(1, 2)), fold_change=TRUE)+theme_minimal()
Na rys. 4.12 przedstawiono wykres pudełkowy pozwalający na porównanie rozkładu dwóch prób. Na wykresie dodatkowo wskazano \(p\)-wartość, która potwierdza, że na poziomie istotności \(\alpha = 0,05\) należy odrzucić hipotezę \(H_0\), a więc można twierdzić, że wartości oczekiwane w populacjach nie są jednakowe.
W poprzednim przykładzie porównano dwie próby. Poniższy kod prowadzi do wizualizacji porównania parami dla 6 prób. Wynik został przedstawiony na rys. 4.13. Na wykresie zamieszczone zostały \(p\)-wartości dla porównań parami zmiennych 1 i 2 oraz 3 i 5. Przyjmując poziom istotności \(\alpha=0,05\) w pierwszym przypadku różnice są istotne, natomiast w drugim nie są istotne, co zostało zilustrowane na rys. 4.13.
p <- ggplot(chickwts, aes(feed, weight)) +
geom_boxplot() +
geom_jitter()
add_pval(p, pairs = list(c(1, 2),c(3,5)), fold_change=TRUE)+theme_minimal()
Na rys. 4.14 przedstawiono porównanie trzech grup w innej formie - z wykorzystaniem paneli. Dla porównywanych zmiennych na wykresach zostały zamieszczone \(p\)-wartości.
p <- ggplot(mtcars, aes(factor(am), mpg)) +
geom_boxplot() +
geom_jitter() +
facet_wrap(~cyl)
add_pval(p, pairs = list(c(1, 2)))+theme_minimal()
Na rys. 4.15 przedstawiono porównania, nie jak poprzednio w formie wykresów pudełkowych, ale w formie wykresów słupkowych. Podobnie jak w poprzednich przykładach na wykresie wskazano \(p\)-wartości, dla par zmiennych statystycznie istotnie różniących się.
mt=as.data.table(mtcars)
mt[, mu := mean(mpg), by = c("am", "cyl")]
mt[, se := sd(mpg) / sqrt(.N), by = c("am", "cyl")]
p_bar <- ggplot(mt, aes(x=factor(cyl), y=mu, fill = factor(cyl))) +
geom_bar(stat = "identity", position = 'dodge') +
geom_errorbar(aes(ymin=mu-3*se, ymax=mu+3*se),width = .2) +
facet_wrap(~am)
add_pval(p_bar, pairs = list(c(1, 2)), response = 'mpg')+theme_minimal()
4.2.4 Pakiet ggpubr
Kolejną biblioteką, która może być wykorzystana do graficznej prezentacji wyników testowania hipotez jest biblioteka ggpubr. Możliwości tej biblioteki są zbliżone do wcześniej opisanej biblioteki ggpval. Funkcje tej biblioteki pozwalają tradycyjne wykresy w prosty sposób wzbogacić o wyniki testów statystycznych.
W poniższym przykładzie została przedstawiona graficzna prezentacja mocy silnika (hp) dla trzech grup wyróżnionych ze względu na liczbę cylindrów (cyl) na podstawie danych pochodzących ze zbioru mtcars.
p <- ggboxplot(mtcars, x = "cyl", y = "hp",
color = "cyl", palette =c("red", "green", "blue"),
add = "jitter", shape = "cyl")
porownanie <- list( c("4", "6"), c("6", "8"), c("4", "8") )
p + stat_compare_means(comparisons = porownanie)+
stat_compare_means(label.y = 320)
rys. 4.16 przedstawia wykres pudełkowy dla zmiennej hp ze zbioru mtcars. Na wykresie zostały wskazane \(p\)-wartości dla porównań wartości oczekiwanych hp dla trzech wyróżnionych grup ze względu na liczbę cylindrów.
Inną formę prezentacji, nie w formie wykresów pudełkowych a wiolinowych z pudełkami, wyników wnioskowania dla tych samych zmiennych przedstawia rys. 4.17.
ggviolin(mtcars, x = "cyl", y = "hp", fill = "cyl",
palette = c("red", "green", "blue"),
add = "boxplot", add.params = list(fill = "white"))+
stat_compare_means(comparisons = porownanie, label = "p.signif")+
stat_compare_means(label.y = 50) 
4.2.5 Pakiet ggsignif
Kolejnym pakietem umożliwiającym wizualizację wyników wnioskowania statystycznego jest pakiet ggsignif. Umożliwia on dodawanie oznaczeń istotności statystycznej (np. \(p\)-wartości, gwiazdek, linii łączących grupy) bezpośrednio na wykresach. Jest szczególnie przydatny przy prezentacji wyników testów statystycznych na wykresach typu boxplot, barplot, violin plot itp.
Na rys. 4.18 w sposób czytelny zaprezentowano różnice w zużyciu paliwa pomiędzy samochodami 4-, 6- i 8-cylindrowymi oraz informację, które z tych różnic, i na jakim poziomie, są statystycznie istotne.
ggplot(mtcars, aes(factor(cyl), mpg)) +
geom_boxplot() +
geom_signif(
comparisons = list(c("4", "6"), c("6", "8"),c("4", "8")),
map_signif_level = TRUE, textsize = 6) +
ylim(10, 40)+theme_minimal()
rys. 4.19 łączy wizualną prezentację rozkładów (boxplot) z oceną statystyczną jednej, konkretnej pary gatunków, uwzględniając dodatkowo informację o „szerokości płatka” za pomocą koloru.
test <- t.test(
iris[iris$Petal.Width > 1 & iris$Species == "versicolor", "Sepal.Width"],
iris[iris$Species == "virginica", "Sepal.Width"])$p.value
ggplot(iris, aes(x = Species, y = Sepal.Width, fill = Petal.Width > 1)) +
geom_boxplot(position = "dodge") +
geom_signif(
annotation = formatC(test, digits = 1),
y_position = 4.05, xmin = 2.2, xmax = 3,
tip_length = c(0.2, 0.04))+theme_minimal()
Na rys. 4.19 przedstawiono porównanie szerokości działki kielicha dla trzech gatunków z wyróżnieniem szerokości płatka.
Kolejny kod i ?fig-g30porównanie masy piskląt w zależności od rodzaju paszy. Na wykresie wskazano istotne różnice dla wybranych par paszy.
ggplot(data = chickwts,mapping = aes(x = weight,y = feed,fill=feed)) +
geom_boxplot(orientation = "y") +
geom_signif(
comparisons = list(c("horsebean", "sunflower"), c("horsebean", "linseed"), c("casein", "meatmeal")),
map_signif_level = TRUE,
textsize = 6,
margin_top = 0.08,
step_increase = 0.05,
tip_length = 0.01,
orientation = "y")+theme_minimal()
Wykres rys. 4.20 potwierdza, że rodzaj paszy wpływa znacząco na masę kurcząt: casein i sunflower prowadzą do najwyższej wagi, a horsebean do wyraźnie najniższej. Dla wybranych par paszy wskazano istotne różnice.
4.2.6 Pakiet ggstatsplot
Pakiet ggstatsplot to narzędzie do wizualizacji danych w programie R, które łączy prezentację graficzną z przedstawieniem wyników testów statystycznych. Pakiet automatycznie dobiera i wykonuje odpowiednie testy statystyczne (np. t-test, ANOVA, testy nieparametryczne, testy dla danych jakościowych) w zależności od typu danych i liczby grup.
Poniższy kod poddaje analizie przeżywalność katastrofy statku Titanic.
ggpiestats(
data = titanic,
x = Class,
y = Survived,
legend.title = "Klasa")
rys. 4.21 przedstawia udział pasażerów, którzy przeżyli i nie przeżyli katastrofy Titanica, w podziale na klasy podróży (1st, 2nd, 3rd, Crew). Oba wykresy kołowe pokazują proporcje pasażerów według klas. Pasażerowie pierwszej klasy mieli najwyższy odsetek przeżycia – większość osób z tej grupy przeżyła katastrofę. Pasażerowie trzeciej klasy mieli najniższy odsetek przeżycia – większość osób z tej grupy nie przeżyła. Wyniki testu statystycznego (np. test chi-kwadrat), prezentowane na wykresie, wskazują, że różnice w przeżywalności pomiędzy klasami są statystycznie istotne (\(p\)-wartość < 0.05).
ggbarstats(
data = titanic,
x = Class,
y = Survived,
xlab = "Przeżycie katastrofy",
legend.title = "Klasa",
ggplot.component = list(ggplot2::scale_x_discrete(guide = ggplot2::guide_axis(n.dodge = 2))),
palette = "Set2")
Wykres słupkowy (rys. 4.22) przedstawia liczbę pasażerów, którzy przeżyli i nie przeżyli katastrofy Titanica, w podziale na klasy podróży (1st, 2nd, 3rd, Crew). Każdy słupek reprezentuje liczbę osób w danej klasie, a kolory odpowiadają poszczególnym klasom według legendy. Pasażerowie pierwszej klasy mieli najwyższy odsetek przeżycia – w tej grupie liczba osób, które przeżyły, jest zbliżona do liczby tych, które nie przeżyły, a nawet może być wyższa. Pasażerowie trzeciej klasy mieli najniższy odsetek przeżycia – w tej grupie zdecydowanie przeważają osoby, które nie przeżyły katastrofy. Na wykresie widoczne są także wyniki testu chi-kwadrat niezależności, które potwierdzają, że różnice w przeżywalności pomiędzy klasami są istotne statystycznie (\(p\)-wartość < 0.05).
Kolejne dwa wykresy dotyczą wizualizacji wyników testowania dla danych, z wcześniej wspomnianego zbioru mtcars.
model <- stats::lm(formula = mpg ~ am * cyl, data = mtcars)
ggcoefstats(model)
Wykres rys. 4.23 pokazuje punkty reprezentujące oszacowane wartości współczynników regresji. Te współczynniki określają ilościowo wpływ każdej zmiennej predykcyjnej na zmienną zależną (mpg). Paski błędu wokół punktów reprezentują przedziały ufności (zwykle 95%) dla współczynników. Te przedziały wskazują zakres, w którym z dużą pewnością znajduje się rzeczywista wartość współczynnika korelacji w populacji. Wykres wykorzystuje kolory do wskazania istotności statystycznej każdego współczynnika. Współczynniki z \(p-wartościami\) poniżej wybranego poziomu istotności (np. 0,05) są uważane za statystycznie istotne, co oznacza, że jest mało prawdopodobne, aby obserwowany efekt wystąpił przypadkowo.
ggbetweenstats(data = mtcars,x = cyl, y = mpg )
rys. 4.24 przedstawia porównanie średniego zużycia paliwa (mpg) w samochodach z różną liczbą cylindrów (4, 6 i 8) w zbiorze danych mtcars. Każda grupa na osi OX odpowiada innej liczbie cylindrów, a na osi OY przedstawione są wartości mpg dla poszczególnych samochodów. Samochody z 4 cylindrami mają najwyższe średnie wartości mpg (czyli są najbardziej ekonomiczne). Samochody z 8 cylindrami mają najniższe wartości mpg, co oznacza, że zużywają najwięcej paliwa. Na wykresie prezentowane sa wyniki testu ANOVA i porównania parami. Różnice pomiędzy grupami są istotne statystycznie (\(p\)-wartości < 0.05).
4.2.7 Pakiet ggdist
Pakiet ggdist umożliwia konstrukcję elastycznych i bogatych w informacje wizualizacji, które pokazują rozkłady danych oraz niepewność estymacji. Poza prezentacją gęstości umożliwia pokazanie np. średnich i przedziałów ufności.
mtcars$gear <- as.factor(mtcars$gear)
ggplot(mtcars, aes(x = mpg, y = gear, fill = gear)) +
stat_halfeye(
adjust = 0.4,
.width = c(0.90),
point_interval = mean_qi,
alpha = 0.7) +
labs(
x = "mpg",
y = "Liczba biegów") +
theme_minimal() +
theme(legend.position = "none")
Na rys. 4.25 pokazano jednocześnie oszacowanie gęstości, średnią oraz przedziały ufności zmiennej mpg dla samochodów z różną liczbą biegów (gear).
Poniższy kod pozwala porównać rozkład mpg dla samochodów z automatyczną i ręczną skrzynią biegów. Używa gradientu kolorów do pokazania gęstości prawdopodobieństwa, co jest alternatywą dla klasycznego wykresu skrzypcowego.
ggplot(mtcars, aes(x = mpg, y = factor(am), fill = after_stat(x))) +
stat_gradientinterval(
show.legend = FALSE ,
.width = c(0.90)) +
scale_fill_viridis_c() +
labs(
x = "Liczba mil na galonie paliwa",
y = "Typ skrzyni biegów") +
theme_minimal()
Na podstawie rys. 4.26 można zauważyć, że samochody z ręczną skrzynią biegów mają generalnie przejechaną większą liczbę mil na galonie paliwa. Dla obu wariantów skrzyni biegów na wykresie zostały zaznaczone średnia wartość mpg i przedział ufności (poziom ufności = 0,90).
4.2.8 Pakiet infer
Pakiet infer w R został skonstruowany z myślą o prostym i intuicyjnym przeprowadzaniu wnioskowania statystycznego, szczególnie za pomocą metod symulacyjnych, takich jak testy permutacyjne czy bootstrap. Dzięki zastosowaniu składni zgodnej z tidyverse, infer pozwala użytkownikom na łatwe definiowanie hipotez, generowanie rozkładów statystyk testowych oraz wyznaczenie \(p\)-wartości i przedziałów ufności. Pakiet ten ułatwia stosowanie nowoczesnych nieklasycznych metod statystycznych, oferując czytelne i elastyczne podejście do analizy danych. Pakiet ten umożliwia także przeprowadzenie klasycznego testowania i wizualizację wyników. Przykład realizacji takiego testu dotyczącego porównania mpg dla samochodów z 4 i 6 cylindrami przedstawia poniższy kod, a wizualizacja wyników jest zamieszczona na rys. 4.27.
df <- mtcars %>%
filter(cyl %in% c(4, 6)) %>%
mutate(cyl = factor(cyl))
obs_stat <- df %>%
specify(mpg ~ cyl) %>%
calculate(stat = "diff in means", order = c("6", "4"))
t_test <- t.test(mpg ~ cyl, data = df)
n1 <- sum(df$cyl == 6)
n2 <- sum(df$cyl == 4)
s2 <- var(df$mpg) # lwspólna wariancja z obu grup
se <- sqrt(s2 * (1/n1 + 1/n2))
null_dist <- tibble(
stat = seq(obs_stat$stat - 4*se, obs_stat$stat + 8*se, length.out = 1000)) %>%
mutate(density = dnorm(stat, mean = 0, sd = se))
ggplot(null_dist, aes(x = stat, y = density)) +
geom_line(color = "blue") +
xlim(-8,8)+
geom_vline(xintercept = obs_stat$stat, color = "red", linetype = "dashed") +
labs(subtitle = paste0("p-value = ", signif(t_test$p.value, 3)),x = "Różnica średnich", y = "Gęstość") +
theme_minimal()
Rezultat realizacji powyższego kodu przedstawiono na rys. 4.27. Na wykresie zaprezentowano gęstość rozkładu różnicy średnich przy założeniu słuszności \(H_0\). Zaznaczona jest także wartość statystyki testowej i odpowiadająca jej \(p\)-wartość. W omawianym przypadku należy odrzucić hipotezę \(H_0\), a więc można twierdzić, że przeciętna liczba mil przejechanych na galonie paliwa dla samochodów z 4 i 6 cylindrami nie jest jednakowa.
Z wykorzystaniem pakietu infer można uzyskać symulacyjne przybliżenie rozkładu teoretycznego przy założeniu hipotezy \(H_0\) (rys. 4.28).
null_dist <- df %>%
specify(mpg ~ cyl) %>%
hypothesize(null = "independence") %>%
generate(reps = 1000, type = "permute") %>%
calculate(stat = "diff in means", order = c("6", "4"))
visualize(null_dist) + labs(title = NULL)+theme_minimal()
Możliwości pakietu infer zostaną wykorzystane w kolejnym rozdziale do weryfikacji i wizualizacji testów permutacyjnych.
4.3 Wizualizacja wyników wnioskowania - implementacja użytkownika
W poprzednim punkcie przedstawiono zastosowanie różnych bibliotek programu R do wizualizacji wyników wnioskowania statystycznego. Praktyka badań naukowych prowadzi często do konieczności samodzielnego opracowania procedur wnioskowania i niekiedy też odpowiednich ich wizualizacji. Przykład takiego rozwiązania w zastosowaniu do testu pustych cel zostanie przedstawiony w tym punkcie.
4.3.1 Test pustych cel - idea, założenia i hipotezy
Test pustych cel (ang. empty cells test) jest jednym z testów zgodności. Za pomocą tego testu można weryfikować hipotezę dotyczącą postaci rozkładu zmiennej losowej. Test pustych cel to statystyczna metoda weryfikacji losowości, wywodząca się z klasycznego problemu rozmieszczeń (ang. occupancy problem). W problemie tym badane jest czy \(n\) elementów zostało rozmieszczonych w \(m\) celach w sposób całkowicie losowy. Ideę testu oraz własności przedstawiają m.in. David (1950); Hellwig (1965), Domański (2011) oraz Domański (2012). Test pustych cel pozwala na bardzo różnorodne zastosowania. Modyfikacje testu dla przypadków jedno- i dwuwymiarowego przedstawił Kończak (2008) i Kończak (2009). Z kolei zastosowanie testu pustych cel dla weryfikacji wielowymiarowej normalności przedstawił Kończak (2025).
Cały obszar zmienności zmiennej losowej dzieli się na \(m\) cel i zlicza się liczbę elementów z losowej próby w każdej celi. Statystyką testową jest liczba cel, które po zakończeniu procesu pozostają puste \((K_n)\). Liczba ta jest porównywana z wartością krytyczną. Domański & Pruska (2000) przedstawili interpolowane wartości krytyczne dla statystyki testu pustych komórek. Załóżmy, że z populacji pobrano próbę losową składającą się z \(n\) elementów. Weryfikacji poddana zostanie hipoteza dotycząca postaci rozkładu.
Dane: liczba elementów \(n\) oraz stała liczba cel \(m\).
Niezależność: każdy element trafia do celu niezależnie od pozostałych elementów.
Jednakowe prawdopodobieństwo: Prawdopodobieństwo trafienia elementu do dowolnej z \(m\) cel jest jednakowe i wynosi: \(p = \frac{1}{m}\).
Brak ograniczeń pojemności: Każda cela może pomieścić dowolną liczbę elementów.
Niech będzie pobrana próba \(n\)-elementowa z populacji o dystrybuancie \(F(x)\). Weryfikowana jest hipoteza o zgodności rozkładu z pewnym ustalonym rozkładem teoretycznym \(F_0(x)\).
\[H_0: F(x) = F_0(x)\]
\[H_1: F(x) \ne F_0(x)\] gdzie \(F_0\) jest daną dystrybuantą.
Obszar zmienności zmiennej losowej \(X\) dzielony jest na \(m\) rozłącznych cel \(M_i (i = 1, 2, \ldots, m)\), takich że
\[\mathbb{X} = \bigcup_{i=1}^{m} M_i\] gdzie \(\mathbb{X}\) jest obszarem zmiennosci zmiennej losowej \(X\) oraz
\(M_i \cap M_j = \emptyset\) dla \(i, j \in \{1, 2, \ldots, m\}, i \ne j\),
\(P(x \in M_i) = \frac{1}{m}\) dla \(i \in \{1, 2, \ldots, m\}\), jeśli hipoteza \(H_0\) jest prawdziwa.
Statystyka testowa (liczba pustych cel) ma następująca postać
\[K= card \{j: m_j=0\}.\]
Dla przyjętego poziomu istotności \(\alpha\) obszar krytyczny jest następujący:
\[\mathbf{K}=\{k: k \ge K_{n, \alpha} \}\] gdzie \(K_{n, \alpha}\) są wartościami krytycznymi uzyskanymi z tablic (Domański & Pruska, 2000) lub uzyskanymi symulacyjnie.
4.3.2 Test pustych cel - implementacja testu normalności
Test pustych cel, to ogólna idea. Na tej podstawie mogą być konstruowane specyficzne testy statystyczne. Kończak (2005) przedstawił modyfikację testu na przypadek wyznaczania nie tylko pustych cel, ale także cel jednoelementowych. Z kolei Kończak (2025) przedstawia konstrukcję wielowymiarowego testu normalności z zadanymi parametrami.
Poniżej przedstawiono konstrukcję testu pustych cel dla weryfikacji hipotezy o normalności rozkładu. Przed zastosowaniem testu zostaną określone dwie funkcje. Funkcja ect_Kn służy do obliczania wartości statystyki testowej (\(K_n\)), czyli liczby pustych cel.
ect_Kn <- function(data, mean_val, sd_val, num_bins) {
p_vals <- pnorm(data, mean = mean_val, sd = sd_val)
p_vals <- pmax(1e-10, pmin(1 - 1e-10, p_vals))
bins <- ceiling(p_vals * num_bins)
return(num_bins - length(unique(bins)))
}Funkcja ect_norm realizuje procedurę testu zgodności z rozkładem normalnym. Dla zadanego poziomu istotności \(\alpha\) w oparciu o symulacje komputerowe wyznacza wartość krytyczną (\(K_{n, \alpha}\)) \(p\)-wartość. Na podstawie \(R\) losowych prób z rozkładu normalnego konstruowany jest rozkład odniesienia statystyki przy nieznanych parametrach. Wynikiem jest decyzja statystyczna oparta na porównaniu z symulacjami oraz wizualizacja graficzna wyniku na tle histogramu.
ect_norm <- function(x, alpha=0.05, R = 2000, plot = TRUE) {
m <- n <- length(x)
mu_obs <- mean(x, na.rm = TRUE)
sd_obs <- sd(x, na.rm = TRUE)
K0 <- ect_Kn(x, mu_obs, sd_obs, m)
K_sim <- numeric(R)
for (i in 1:R) {
x_sim <- rnorm(n, mean = mu_obs, sd = sd_obs)
sim_mu <- mean(x_sim)
sim_sd <- sd(x_sim)
K_sim[i] <- ect_Kn(x_sim, sim_mu, sim_sd, m)
}
p_value <- (sum(K_sim >= K0) + 1) / (R + 1)
if (plot) {
K_crit <- quantile(K_sim, 1 - alpha)
hist_res <- hist(K_sim, breaks = seq(min(K_sim)-0.5, max(K_sim)+0.5, by=1), plot = FALSE)
col_hist <- if (p_value < 0.05) "mistyrose" else "lightgreen"
plot(hist_res, col = col_hist, freq = FALSE, xlim = c(0,12),
main = "", xlab = "Liczba pustych cel", ylab="Gęstość")
abline(v = K_crit, col = "red", lwd = 2, lty = 2)
abline(v = K0, col = "blue", lwd = 3)
axis(side = 1, at = K_crit, labels = expression(K[n ~ "," ~ alpha]), font = 2, padj = 1.5)
axis(side = 1, at = K0, labels = expression(K[0]), font = 2, padj = 1.5)
}
return(list(n_empty = K0, p.value = p_value, estimated_params = list(mean=mu_obs, sd=sd_obs)))
}Poniżej przedstawiono zastosowanie wprowadzonego powyżej testu pustych cel dla weryfikacji hipotezy o normalności rozkładu analizowanej zmiennej. Przeprowadzono weryfikację hipotezy dla 2 prób o liczebnościach \(n = 15\). Pierwsza z tych prób została wygenerowana z rozkładu normalnego, a druga z rozkładu wykładniczego. Tylko w pierwszym z tych przypadków należy oczekiwać, że hipoteza \(H_0\) nie zostanie odrzucona.
Na rys. 4.29 i rys. 4.30 przedstawiono graficznie wyniki testowania hipotezy o normalności rozkładu dla tych dwóch rozkładów.
Na wykresach zaznaczono liczbę pustych cel \((K_0)\) (kolor niebieski), wyznaczoną symulacyjnie wartość krytyczną \((K_{n, \alpha})\) (czerwona linia przerywana). Dodatkowo kolorem słupków rozróżniono przypadki braku podstaw do odrzucenia hipotezy \(H_0\) (kolor jasnozielony) i odrzucenia tej hipotezy (kolor jasnoczerwony). Na wykresach rys. 4.29 i rys. 4.30 przedstawiono najważniejsze wyniki przeprowadzonego testu, a dodatkowo decyzja została zakodowana w kolorze słupków wykresu, co może być pomocne zwłaszcza w praktycznych zastosowaniach biznesowych.
dane1 <- rnorm(15, mean = 10, sd = 2)
dane2 <- rexp(15,rate=2)
ect_norm(dane1)$p.value
ect_norm(dane2)$p.value# [1] 0.8850575
# [1] 0.01549225
4.3.3 Dwuwymiarowy test normalności - wizualizacja wyników
Kończak (2025) przedstawił proponozycję wielomiarowego testu normalności opartego na idei pustych cel. Weryfikowana jest hipoteza o rozkładzie normalnym z zadanymi parametrami. Poniżej zostanie przytoczona idea testu dla rozkładu dwuwymiarowego.
Konstrukcja testu opiera się na partycjonowaniu przestrzeni dwuwymiarowej na rozłączne obszary (cele) o jednakowym prawdopodobieństwie teoretycznym. Podział ten przebiega dwuetapowo:
Przestrzeń dzielona jest na \(k\) koncentrycznych pierścieni eliptycznych, dobranych tak, aby każdy z nich gromadził równą masę prawdopodobieństwa \((\frac{1}{k})\) zgodnie z teoretycznym rozkładem normalnym.
Każdy pierścień jest następnie dzielony na cztery sektory przez główne osie elipsoidy, wyznaczone przez wektory własne macierzy kowariancji. W rezultacie otrzymuje się układ cel, do których – przy założeniu hipotezy zerowej – obserwacje trafiają z jednakowym prawdopodobieństwem. Idea konstrukcji cel została przedstawiona na rys. 4.31.
Źródło: opracowanie własne na podstawie Kończak (2025)
Wyniki testu normalności zawierają w szczególności informacje o wartości statystyki testowej \((K)\), wartości krytycznej \((K_{n, \alpha})\), \(p\)-wartość, decyzja. Wszystkie te wyniki warto przedstawić w formie wizualizacji jak np. na rys. 4.32.
Źródło: opracowanie własne
rys. 4.32 przedstawia graficzną ilustrację wyniku testu wielowymiarowej normalności opartego na liczbie pustych cel. Propozycja i własności testu zostały przedstawione w Kończak (2025). Wykres ten pokazuje relację między wartością obliczoną z danych \(K\) a rozkładem teoretycznym (symulacyjnym) przy założeniu prawdziwości hipotezy zerowej. Całość składa się z czterech wykresów tworzących kompleksową wizualizację wyników testu dwuwymiarowej normalności. W kolejnych panelach przedstawiono następujące informacje.
Panel 1:
Wykres przedstawia obserwacje po transformacji do przestrzeni standaryzowanej, w której elipsoidy stają się okręgami. Wykres umożliwia ocenę, czy obserwacje rozkładają się równomiernie w celach. Występowanie skupień jest argumentem przeciw hipotezie \(H_0\).
Panel 2:
Wykres typu mapa ciepła pokazuje liczbę obserwacji w każdej celi w formie macierzy. Wykres przedstawia macierz, gdzie wiersze odpowiadają kolejnym elipsoidom, a kolumny ćwiartkom \((Q1–Q4)\). Każdy prostokąt reprezentuje jedną celę. Wewnątrz każdego pola podana jest liczba obserwacji, które wpadły do danej celi. Cele puste zostały wyróżnione kolorem czerwonym.
Panel 3:
Wykres słupkowy przedstawia rozkład empiryczny statystyki testowej (liczby pustych cel \(K\)) przy założeniu, że hipoteza zerowa jest prawdziwa. Został on wyznaczony symulacyjnie. Kolor histogramu jest wizualnym wskaźnikiem decyzji testowej. Jasnozielony oznacza brak podstaw do odrzucenia hipotezy zerowej. Różowy oznacza odrzucenie hipotezy zerowej. Czerwona linia przerywana wskazuje interpolowaną wartość krytyczną \((K_{n, \alpha})\) dla ustalonego poziomu istotności \(\alpha\). Jest to granica obszaru krytycznego. W tym teście obszar krytyczny znajduje się na prawo od tej linii. Jeśli liczba pustych cel przekroczy tę linię, oznacza to, że należy odrzucić hipotezę \(H_0\). Niebieska linia ciągła Wskazuje wartość statystyki testowej \((K)\) obliczoną dla badanej próby.
Panel 4:
Jest to panel tekstowy, który przedstawia kluczowe wyniki testu w jednym miejscu.
Statystyka \(K\): Liczba pustych cel znaleziona w danych.
Wartość krytyczna \(K_{n, \alpha}\): Próg liczby pustych cel, powyżej którego przy ustalonym poziomie istotności α należy odrzucić hipotezę \(H_0\).
\(E[K]\): Oczekiwana liczba pustych cel.
\(P-value\): Prawdopodobieństwo uzyskania takiego lub gorszego wyniku przy założeniu \(H0\).
Decyzja: Komunikat (“ODRZUĆ H0” czerwoną czcionką lub “BRAK PODSTAW DO ODRZUCENIA H0” zieloną czcionką).
David, F. N. (1950). Two combinatorial tests whether a sample has come from a given population. W Biometrika (T. 37, s. 97–110).
Domański, C. (2011). Własności testu Davida-Hellwiga. Prace Naukowe Uniwersytetu Ekonomicznego we Wrocławiu, 165, 40–49.
Domański, C. (2012). Statistical Tests based on Empty Cells. Acta Universitatis Lodziensis. Folia Oeconomica, 269, 39–47.
Domański, C., & Pruska, K. (2000). Nieklasyczne metody statystyczne (s. 368). Polskie Wydawnictwo Ekonomiczne.
Hellwig, Z. (1965). Test zgodności dla małej próbyt. Przegląd Statystyczny, 12, 99–112.
Kończak, G. (2005). On the Modification of David-Hellwig Test. Innovation in Classification, Data Science and Information Systems. Proceedings of the 27th Annual Conference of the Gesellschaft für Klassifikation e. V., Brandenburg University of Technology, Cottbus, March 12–14, 2003, 138–145.
Kończak, G. (2008). On the Multivariate Goodness-of-Fit test. Compstat 2008: Proceedings in Computational Statistics, 751–758.
Kończak, G. (2009). On the Modification of the Empty Cells Test. Acta Universitatis Lodziensis. Folia Oeconomica, 228, 228.
Kończak, G. (2024). Wizualizacja wyników badań naukowych. Zasady, metody i narzędzia. Uniwersytet Ekonomiczny w Katowicach.
Kończak, G. (2025). On a new goodness-of-fit test for multivariate normality with fixed parameters based on the David-Hellwig test idea. Statistics in Transition New Series, 26(4), 151–166. https://doi.org/10.59139/stattrans-2025-044
Wells, H. G. (1903). Mankind in the Making. Macmillan; Co.