Histogram w R: hist() z przykładami

Co pokazuje histogram

Histogram odpowiada na jedno pytanie o jedną zmienną liczbową: jak rozkładają się wartości? Dzieli zakres danych na kolejne przedziały (ang. bins), liczy, ile obserwacji trafia do każdego, i rysuje dla każdego przedziału słupek o wysokości równej tej liczbie. Tam, gdzie danych jest dużo, słupki są wysokie, a tam, gdzie mało, niskie.

Ten jeden obrazek mówi rzeczy, których nie powie średnia ani odchylenie standardowe: czy dane mają jeden szczyt, czy dwa, czy są symetryczne, czy przesunięte w stronę jednego ogona, i czy są pojedyncze wartości daleko od reszty. Zwykle to pierwszy wykres, jaki robisz przy nowym zbiorze danych, zaraz obok liczb ze statystyki opisowej.

W R służy do tego funkcja hist() i wystarczy jedna linia.

Pierwszy histogram z hist()

Zasymulujmy wzrost 200 osób o rozkładzie normalnym wokół 170 cm z odchyleniem standardowym 10 (co robią set.seed() i rnorm(), wyjaśnia artykuł o liczbach losowych):

set.seed(42)
heights <- rnorm(200, mean = 170, sd = 10)

hist(heights)

R sam wybiera przedziały i rysuje dzwonowaty stos słupków: niskie słupki w okolicach 145 cm, rosnące przez sto sześćdziesiąt kilka, szczyt około 170 i symetryczny spadek za 190. Oś pionowa ma etykietę "Frequency" (surowe liczebności), a tytuł i etykieta osi pochodzą od nazwy zmiennej.

Ustawienia domyślne są celowo skromne. Działają tu wszystkie zwykłe argumenty upiększające:

hist(heights,
     main   = "Distribution of heights",
     xlab   = "Height (cm)",
     col    = "steelblue",
     border = "white")

col wypełnia słupki, border koloruje ich obrys, a border = "white" daje czysty wygląd rozdzielonych słupków, który widać w większości publikacji. Wszystko, co przewodnik po plot() mówi o kolorach i tytułach, obowiązuje bez zmian.

Szerokość przedziału: argument breaks

Najważniejszym argumentem hist() jest breaks, bo szerokość przedziału decyduje o tym, jaką historię opowiada wykres. Te same dane, dwa ustawienia:

hist(heights, breaks = 5)    # five wide bins: a crude, blocky bell
hist(heights, breaks = 30)   # thirty narrow bins: detail, plus noise

Przy breaks = 5 histogram jest tak zgrubny, że wszystko wygląda jak jeden gładki garb: drugie skupisko albo luka w danych byłyby niewidoczne. Przy breaks = 30 widać drobną strukturę, ale losowe drgania zaczynają udawać cechy danych. Żadne z ustawień nie jest "poprawne"; uczciwie jest wypróbować kilka i sprawdzić, które cechy przetrwają.

breaks przyjmuje trzy formy:

• Liczba, np. breaks = 30: sugestia liczby przedziałów. R dopasowuje ją tak, żeby trafić w równe granice, więc możesz dostać 28 albo 33 przedziały. Każdego to kiedyś zaskakuje.
• Wektor punktów podziału, np. breaks = seq(140, 200, by = 5): dokładne granice przedziałów, bez negocjacji. Używaj tego, gdy przedziały muszą się pokrywać w kilku porównywanych histogramach.
• Nazwa reguły, np. breaks = "FD" dla reguły Freedmana-Diaconisa, która dobiera szerokość na podstawie rozrzutu danych i liczebności próby i dobrze radzi sobie z danymi skośnymi.

Histogram gęstości i krzywa normalna

Domyślnie wysokości słupków to liczebności (freq = TRUE). Ustawienie freq = FALSE przeskalowuje słupki tak, że ich łączne pole wynosi 1: to skala gęstości. Kształt się nie zmienia, zmienia się oś pionowa. Chodzi o to, że histogram gęstości jest w tej samej skali co funkcja gęstości prawdopodobieństwa, więc możesz nałożyć teoretyczną krzywą bezpośrednio na dane:

hist(heights,
     freq = FALSE,
     col = "gray90",
     main = "Heights vs. a normal curve",
     xlab = "Height (cm)")

curve(dnorm(x, mean = mean(heights), sd = sd(heights)),
      add = TRUE, col = "tomato", lwd = 2)

curve() z add = TRUE rysuje krzywą dzwonową na istniejącym histogramie (x w dnorm(x, ...) to symbol zastępczy, który wypełnia curve(), a nie twoja zmienna). Jeśli słupki przylegają do krzywej, model normalny jest rozsądny; tam, gdzie od niej odstają (ciężki ogon, drugi szczyt), modelowi czegoś brakuje. Pomiń freq = FALSE, a krzywa będzie bezużytecznie pełzać po dole wykresu, bo liczebności i gęstości są w różnych skalach.

Liczby stojące za słupkami

hist() nie tylko rysuje: zwraca też podział na przedziały jako listę. Z plot = FALSE całkiem pomija rysowanie i tylko liczy, więc ten przykład możesz uruchomić od razu tutaj:

set.seed(42)
heights <- rnorm(200, mean = 170, sd = 10)

h <- hist(heights, breaks = 5, plot = FALSE)
print(h$breaks)   # the bin edges
print(h$counts)   # observations per bin
print(h$mids)     # the midpoint of each bin

# The same idea with table() + cut():
print(table(cut(heights, breaks = 5)))

h$counts to histogram w postaci danych: każda liczba to wysokość jednego słupka. Linia z table(cut(...)) to tekstowy odpowiednik: cut() dzieli wartości na przedziały, table() liczy każdy przedział. To dobry test poprawności, który warto wypisać przed wykresem albo zamiast niego. Jeśli przypiszesz wynik bez plot = FALSE (h <- hist(heights)), R narysuje wykres i zwróci listę; przypisanie nie wyłącza rysowania.

Histogramy w ggplot2

Wersja w ggplot2 zamienia breaks na binwidth, co często jest bardziej naturalnym pokrętłem: podajesz, jak szeroki jest przedział w jednostkach danych, zamiast ile przedziałów chcesz:

library(ggplot2)

ggplot(data.frame(heights), aes(x = heights)) +
    geom_histogram(binwidth = 5, fill = "steelblue", color = "white") +
    labs(title = "Distribution of heights", x = "Height (cm)", y = "Count")

binwidth = 5 oznacza, że każdy słupek obejmuje 5 cm. ggplot2 ostrzega, jeśli nie ustawisz binwidth, i po cichu przyjmuje 30 przedziałów. Posłuchaj ostrzeżenia, bo domyślna wartość rzadko jest właściwą szerokością dla twoich danych. Do szybkiego podglądu jednej zmiennej hist() wymaga mniej pisania; ggplot2 opłaca się, gdy histogram potrzebuje paneli, grup albo spójnego motywu z resztą twoich wykresów.

Co warto zapamiętać

• Histogram pokazuje rozkład jednej zmiennej liczbowej: hist(x) i gotowe.
• breaks to argument, który ma znaczenie: liczba jest tylko sugestią, wektor ustala dokładne granice, a "FD" wybiera rozsądną szerokość. Zawsze wypróbuj więcej niż jedno ustawienie.
• freq = FALSE przełącza na skalę gęstości i właśnie dzięki temu curve(dnorm(...), add = TRUE) sensownie nakłada krzywą normalną.
• hist(x, plot = FALSE) zwraca przedziały jako dane ($breaks, $counts, $mids); table(cut(x, breaks)) to odpowiednik czysto tekstowy.
• W ggplot2 użyj geom_histogram(binwidth = ...) i ustaw binwidth samodzielnie.

Dalej: wykres pudełkowy, czyli wykres rozkładu, który najlepiej sprawdza się przy porównywaniu grup obok siebie.