1 Statystyki opisowe
Język R głównie wykorzystuj się do przetwarzania, analizy i
wizualizacji danych.
Po pierwsze musimy ustalić z jakim typem danych mamy do czynienia.
Wyróżniamy:
dane ilościowe - są to liczby odpowiadające mierzalnym
wielkościom, np. temperatura, liczba lat, waga w kilogramach,
dane jakościowe - są to cechy jakościowe obiektów, np. płeć,
marka samochodu, poziom wykształcenia.
Aby określić jakie dane mamy w naszym zbiorze danych najlepiej jest
wywołać funkcje str() oraz summary(), np:
library("Przewodnik") #biblioteka, z której wykorzystamy przykładowe dane
dane <- daneSoc #daneSoc - jedne z danych dostępne w z tej biblioteki
str(dane)
## 'data.frame': 204 obs. of 7 variables:
## $ wiek : int 70 66 71 57 45 48 53 61 63 47 ...
## $ wyksztalcenie : Factor w/ 4 levels "podstawowe","srednie",..: 4 4 4 2 2 2 2 3 2 3 ...
## $ st.cywilny : Factor w/ 2 levels "singiel","w zwiazku": 2 2 1 2 2 2 2 2 2 2 ...
## $ plec : Factor w/ 2 levels "kobieta","mezczyzna": 2 1 1 2 1 2 2 2 1 2 ...
## $ praca : Factor w/ 2 levels "nie pracuje",..: 2 2 2 2 2 1 2 2 2 2 ...
## $ cisnienie.skurczowe : int 143 123 167 150 130 138 149 103 125 158 ...
## $ cisnienie.rozkurczowe: int 83 80 80 87 83 75 80 63 78 89 ...
## wiek wyksztalcenie st.cywilny plec
## Min. :22.00 podstawowe:93 singiel :120 kobieta : 55
## 1st Qu.:30.00 srednie :55 w zwiazku: 84 mezczyzna:149
## Median :45.00 wyzsze :34
## Mean :43.16 zawodowe :22
## 3rd Qu.:53.00
## Max. :75.00
## praca cisnienie.skurczowe cisnienie.rozkurczowe
## nie pracuje : 52 Min. : 93.0 Min. : 57.00
## uczen lub pracuje:152 1st Qu.:126.0 1st Qu.: 77.00
## Median :137.5 Median : 80.00
## Mean :137.0 Mean : 80.43
## 3rd Qu.:148.0 3rd Qu.: 86.25
## Max. :178.0 Max. :107.00
Widzimy zatem, że wykształcenie, stan cywilny, płeć i praca są
zmiennymi jakościowymi, natomiast wiek, ciśnienie skurczowe i ciśnienie
rozkurczowe są zmiennymi ilościowymi.
2 Dane jakościowe
2.1 Statystyki dla
zmiennych jakościowych
Dane jakościowe są podzielone na kategorie i zwykle są analizowane
przy użyciu liczności.
2.1.1 Tabela
liczności
Zliczenie liczby obserwacji w każdej kategorii dokonuje się za pomocą
funkcji table().
##
## kobieta mezczyzna
## 55 149
Możesz także analizować dwie zmienne jakościowe jednocześnie tzn.
macierz liczności:
##
## singiel w zwiazku
## kobieta 38 17
## mezczyzna 82 67
Można również tworzyć table wielowymiarowe:
## , , = singiel
##
##
## podstawowe srednie wyzsze zawodowe
## kobieta 18 9 8 3
## mezczyzna 46 17 12 7
##
## , , = w zwiazku
##
##
## podstawowe srednie wyzsze zawodowe
## kobieta 4 7 2 4
## mezczyzna 25 22 12 8
2.2 Wykresy dla danych
jakościowych
Wykresy dla danych jakościowych pozwalają na szybką interpretację i
porównanie różnic między kategoriami. Oto kilka podstawowych typów
wykresów używanych do wizualizacji danych jakościowych:
2.2.1 Wykresy
słupkowe
Na jednej osi wykresów słupkowych mamy częstotliwość występowania
danej kategorii, na drugiej osi natomiast kategorię. Możemy mieć wykresy
słupkowe pionowe i poziome.
Służy nam do tego funkcja barplot() do której przekazuje
się table(dane_jakościowe), wektor częstotliwości danych
jakościowych. Na przykład:
Dodatkowe argumenty, takie jak horiz=TRUE sprawia, że
słupki są poziome, las=1 wskazuje w jakim kierunku powinny
być etykiety przy osiach, a cex.names kontroluje rozmiar
etykiet przy osi pionowej.
Za pomocą biblioteki ggplot:
# Przetestuj poniższy kod
library(ggplot2)
ggplot(dane, aes(x = wyksztalcenie)) +
geom_bar() +
coord_flip() + # Odwrócenie osi
theme(axis.text.y = element_text(size = 10),
axis.title.y = element_text(size = 20)) # Dostosowanie rozmiaru tekstu etykiet
Możemy również wykonać bardziej złożone wykresy dla macierzy
liczności, tj. table(dane$plec, dane$wyksztalcenie).
t <- table(dane$plec, dane$wyksztalcenie)
barplot(t,las=1, beside=TRUE, col=c("cyan", "darkblue")) # beside=TRUE sprawia że słupki są obok siebie, dodajemy kolory i legendę
legend("topright", c("kobieta", "mezczyzna" ), fill=c("cyan", "darkblue"))
ggplot(dane, aes(x = wyksztalcenie, fill = plec)) +
geom_bar(position = position_dodge()) +
scale_fill_manual(values = c("cyan", "darkblue")) +
labs(x = "Wykształcenie", y = "Liczba", fill = "Płeć")
2.2.2 Wykresy kołowe
Zamiast udziału liczbowego możemy również określić udział procentowy
zmiennych jakościowych. Do tego świetnie nadaje się wykres kołowy.
Z użyciem ggplot:
ggplot(dane, aes(x="", y="wyksztalcenie", fill=wyksztalcenie))+
geom_bar(stat="identity",width = 1) +
coord_polar("y") + # pytanie: jakie działanie ma ta linijka?
labs(title = "Proporcje dotyczące wykształacenia")
Możemy również przekazać więcej zmiennych do tego typu wykresu.
Należy jednak pamiętać, że zbyt duża ilość zaciemni wykres.
pie(table(dane$wyksztalcenie, dane$plec), labels = c("kobieta, podstawowe","kobieta, średnie", "kobieta, wyższe", "kobieta, zadowe", "meżczyzna, podstawowe", "meżczyzna, srednie", "mężczyzna, wyższe", "meżczyzna, wyższe" ))
Domyślnie funkcja pie() przy macierzach nieliczność nie
podpisuje wycinków koła. Należy nazwy przekazać do funkcji za pomocą
argumentu labels.
Z użyciem ggplot:
dane_summary <- as.data.frame(table(dane$plec, dane$wyksztalcenie))
colnames(dane_summary) <- c("plec", "wyksztalcenie", "count")
dane_summary$label <- paste(dane_summary$plec, dane_summary$wyksztalcenie, sep = ", ")
ggplot(dane_summary, aes(x = "", y = count, fill = label)) +
geom_bar(width = 1, stat = "identity", color = "black") +
coord_polar(theta = "y")
2.2.3 Wykres
mozaikowy
Wykres mozaikowy jest wykorzystywany, aby zobrazować zależność
pomiędzy dwoma lub większą ilością zmiennych jakościowych.
Używając funkcji mozaicplot()można przedstawić
liczebność wartości w danym zbiorze za pomocą pola prostokątów. Na
podstawie tego wykresy można stwierdzić, czy dwie zmienne od siebie
zależą oraz jak często występują dane wartości/pary wartości.
Spójrzmy na poniższy wykres, który przedstawia częstotliwość
występowania różnych grup wykształcenia:
Uwaga Jeżeli dwie zmienne są niezależne, to mozaika
przypomina regularną kratę. Gdy pojawia się zależność pomiędzy
zmiennymi, zazwyczaj łatwo ją zauważyć, np:
Zauważmy, że najwięcej osób pracujących występuje wśród osób z
wykształceniem wyższym. Zwróćmy również uwagę, że zachowała się
szerokość prostokątów z wcześniejszego polecenia, zatem widzimy, że w
naszym zbiorze danych najwięcej jest osób z wykształceniem
podstawowym.
Można za pomocą wykresu mosaicplot() zaprezentować
wyniki dla większej ilości zmiennych jakościowych. Należy jednak uważać,
ponieważ umieszczenie bardzo dużej ilości zmiennych zmniejszy czytelność
wykresu.
mosaicplot(~wyksztalcenie+praca+st.cywilny, data = daneSoc, col=c("darkorange", "grey"))
legend("bottomleft", c("single", "w związku" ), fill=c("darkorange", "grey"))
Można spróbować dodać kolory i legendę w celu zachowania czytelności.
Na podstawie wykresu jesteśmy w takie powiedzieć, że najmniejsze
bezrobocie jest wśród singli.
2.3 Zadania - zmienne
jakościowe
- Załaduj poniższy zbiór danych:
set.seed(111111) # W miejsce "jedynek" wpisz swój numer indeksu
dane_j <- data.frame(
Płeć = sample(c("Mężczyzna", "Kobieta"), 100, replace = TRUE),
Rok = sample(c("1", "2", "3", "4"), 100, replace = TRUE),
Matematyka = sample(c("Lubi matematykę", "Nie lubi matematyki"), 100, replace = TRUE),
Programowanie = sample(c("Lubi programowanie", "Nie lubi programowania"), 100, replace = TRUE),
Sztuka = sample(c("Lubi sztukę", "Nie lubi sztuki"), 100, replace = TRUE)
)
# Podgląd danych
head(dane_j)
Zacznij od wyświetlenia podsumowania dla danych przy użyciu funkcji
summary(). Jeżeli wyświetlona dane są typu
character należy wymusić na programie, aby interpretować
odpowiednie kolumny jako zmienne ilościowe lub jakościowe, używając do
tego np. funkcji as.number() lub as.numeric().
Przykładowo,
Następnie dla powyższego zestawu danych:
Stwórz tabelę liczności i proporcji dla zmiennej Programowanie i
dla zmiennej Sztuka.
Stwórz wykresy słupkowe i kołowe dla obu zmiennych.
Oblicz tabele liczności i proporcji dla zmiennych Płeć i
Programowanie oraz Płeć i Sztuka.
Stwórz wykresy mozaikowe dla obu zestawów zmiennych.
Na podstawie uzyskanych tabel liczności, proporcji i wykresów,
zinterpretuj wyniki. Jakie są preferencje studentów dotyczące różnych
zajęć w zależności od płci? Jakie wnioski można wyciągnąć na podstawie
analizy?
Uwaga: Wykresy powinny mieć tytuły, podpisane osie i
zmienione kolory.
- Wykonaj wykresy z powyższego zadania w wykorzystaniem biblioteki
ggplot2.
3 Dane ilościowe
3.1 Wykresy dla danych
ilościowych
Przedstawimy również kilka podstawowych wykresów wykonywanych dla
zmiennych ilościowych, które wykonuje się przy wstępnej analizie
danych.
3.1.1 Histogram
Jest jednym z najpopularniejszych statystyk graficznych. Przedstawia
rozkład wartości zmiennych ilościowych.
Jak działa histogram?
Histogram składa się z słupków, które ilustrują częstość występowania
danych w określonych przedziałach.
Konstrukcja histogramu:
Wybór ilości przedziałów
Dla każdego przedziału liczona jest liczba obserwacji, które
wpadają w dany zakres.
Każdy przedział jest reprezentowany przez słupek na wykresie,
gdzie wysokość słupka oznacza liczbę obserwacji w danym
przedziale.
Idea histogramu: \[f(x)=\frac{1}{n\cdot
\text{długość przedziału}}\#\{ \text{obserwacje które wpadają w dany
przedział zawierający x}\},\] gdzie \(n\) to ilość obserwacji w naszym zbiorze
danych.
Histogram wyznacza rozkład empiryczny
Oznacza to, że histogram pokazuje, jak często poszczególne wartości
(lub ich przedziały) występują w badanej próbce. Histogram pozwala
zaobserwować kluczowe cechy rozkładu, takie jak modę, skośność, kurtozę,
obecność wartości odstających oraz potencjalną multimodalność (więcej o
tych pojęciach w trakcie przykładów).
Uwaga W praktyce, gdy liczba danych i liczba
przedziałów są odpowiednio duże, histogram może być używany jako
estymator gęstości prawdopodobieństwa.
Histogram generujemy za pomocą funkcji hist().
Opcjonalnymi argumentami jest np. breaks, dzięki któremu
możemy określić liczbę lub sposób podziału na przedziały (może być to
liczba lub nazwa algorytmu). Jeżeli nie określimy liczby przedziałów, to
zostanie ona dobrana w zależności od liczby obserwacji. Za pomocą
argumentów freq i probability możemy określić,
czy na histogramie mają bym zaznaczone frakcje, czy liczebności
wystąpień.
Przykład:
hist(dane$cisnienie.skurczowe, main="Histogram jednomodalny", xlab="Ciśnienie skurczowe", ylab="Frakcja wystepowania", probability = TRUE)
W ggplot
ggplot(dane, aes(x=cisnienie.skurczowe,after_stat(density)))+
geom_histogram()+
labs(title="Histogram jednomodalny")
Zauważmy, że kształt tego histogramu jest w przybliżeniu symetryczny,
ma również jedno maksimum zwane również modą. Jest to
przykład histogramu jednomodalnego.
Na poniższym rysunku są zobrazowane inne kształty histogramów.
Przykłady kształtów histogramów
Średnia, mediana i moda w symetrycznych,
prawoskośnych i lewoskośnych danych
W przypadku symetrycznych histogramów moda, średnia i mediana
pokrywają się (tak jak na panelu A). Jeżeli histogram ma ogon, który
wydaje się przeciągnięty w prawo to nazywamy go
prawoskośnym (panel B). Dodanie większej ilości
wartości w pobliżu maksimum silniej wpływa na średnią niż na modę, czy
medianę. Analogiczny przypadek jest, gdy ogon pojawia się po lewej
stronie, to znaczy mamy więcej wartości w pobliżu minimum. Taki
histogram nazywamy lewoskośnym.
Histogram do porównania zbiorów danych dla różnych
grup
Wykres dwóch zestawów danych w postaci dwóch histogramów może być
bardzo łatwym pierwszym krokiem do porównania środka, rozrzutu i
kształtu każdego z zestawu danych.
library(RVAideMemoire)# Biblioteka do obliczania mody za pomocą funkcji mod().
par(mfrow=c(2,1))
kobieta=dane$plec=="kobieta"
ciskobiet=dane[kobieta,]$cisnienie.skurczowe
mezczyzna=dane$plec=="mezczyzna"
cismezczyzna=dane[mezczyzna,]$cisnienie.skurczowe
hist(ciskobiet, breaks=20, xlim = c(93,178), main="Ciśnienie skurczowe dla kobiet")
abline(v=mean(ciskobiet), col="red")
abline(v=median(ciskobiet),col="blue")
abline(v=mod(ciskobiet), col="green")
hist(cismezczyzna, breaks = 25, xlim = c(93,178), main = "Ciśnienie skurczowe dla mężczyzn")
abline(v=mean(cismezczyzna), col="red")
abline(v=median(cismezczyzna),col="blue")
abline(v=mod(cismezczyzna), col="green")
UWAGA: aby porównywać dwa histogramy ważne jest by
oś x dla oby wykresów była tak samo wyskalowana, można to ustalić za
pomocą argumentu xlim.
Histogram identyfikuje wartości odstające
Na powyższym przykładzie większość
obserwacji wpada w przedział pomiędzy 0 a 200 mm, natomiast jedna
wartość mieści się w przedział pomiędzy 550, a 600. Taką wartość, która
odbiega mocno od pozostałych wartości znajdujących się w zbiorze danych,
nazywamy wartością odstającą.
3.1.2 Jądrowy estymator
gęstości
Można o nim myśleć jako o gładkiej wersji histogramu. Jest bardziej
zaawansowanym narzędziem do przedstawienia empirycznego rozkładu
gęstości dla posiadanych danych.
Zauważmy, że histogram mocno zależy od wyboru ilości przedziałów, ale
również od tego, w jakich miejscach dokonujemy dyskretyzacji, co może
wpływać na jego kształt i prowadzić do błędnych wniosków.
Jądrowy estymator gęstości rozwiązuje problemy jakie pojawiają się
przy histogramach.
Idea jądrowego estymatora gęstości: \[\hat{f}(x)=\frac{1}{n\cdot\text{długość
przedziału}}\#\{\text{obserwacje które wpadają w małe otoczenie punktu
x}\}.\] Definicja jądrowego estymatora
gęstości
Niech \(X_1, \dots, X_n\) będzie
próbą danych, wtedy \[\hat{f}(x)=\frac{1}{nh}\sum_{i=1}^n
K\left(\frac{x-X_i}{h}\right),\] gdzie:
- \(\hat{f}(x)\) to estymowana
gęstość w punkcie \(x\),
- \(h\) to szerokość okna
(bandwidth),
- \(K\) to funkcja jądra (kernel
function),
nazywamy jądrowym estymatorem gęstości.
Jako funkcję jądra bardzo często przyjmuje się tak zwane jądro
Gaussa tj.
\[K(u)=\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}u^2}.\]
Szerokość okna h
Szerokość okna \(h\) kontroluje
stopień wygładzenia estymacji. Małe \(h\) prowadzi do estymacji o wysokiej
wariancji (szczegółowej, “szumowej”), natomiast duże \(h\) prowadzi do estymacji o wysokim
obciążeniu (gładkiej, “rozmytej”). Dla funkcji density()
(patrz przykłady niżej) można ją ustalić ręcznie lub wskazać regułę jej
wyboru, która określi najodpowiedniejszą szerokość, np:
- reguła kciuka -
bw="nrd0" - jest ona używana
domyślnie,
- reguła Scotta -
bw="nrd",
- cross-validation -
bw="ucv",
- estymator Sheathera-Jonesa -
bw="SJ".
Intuicyjnie:
Każdy punkt danych \(X_i\) jest
zastępowany przez funkcję jądra \(K\),
która jest przesunięta tak, aby była wycentrowana na \(X_i\). A następnie te funkcje są sumowane i
znormalizowane przez liczbę punktów i szerokość okna, aby uzyskać
końcową estymację gęstości w każdym punkcie \(x\).
Jądrowy estymator gęstości
W praktyce, aby użyć jądrowego estymatora gęstości wystarczy użyć
funkcji density().
par(mfrow=c(3,1))
plot(density(dane$wiek)) # domyślna wartość okna
plot(density(dane$wiek, bw=0.5)) #mała szerokość okna
plot(density(dane$wiek, bw=10)) #duża szerokość okna
Powyższy przykład ukazuje, jak wykres zmienia się w zależności od
szerokości okna. Jeżeli jest ona mała to na jądrowym estymatorze
gęstości widzimy bardzo dużo szczegółów. Często w praktyce nie jest to
korzystne, ponieważ chcemy wyciągnąć ogólne wnioski i zniwelować
obecność szumu wśród naszych obserwacji. Jeżeli okno jest zbyt duże to
tracimy informacje o kształcie rozkładu naszych danych. Dlatego tak
ważne jest znalezienie kompromisu pomiędzy tymi aspektami.
Jądrowy estymator gęstości często można przedstawiać na tym samym
wykresie co histogram. Pamiętajmy jednak, że w takim przypadku na
histogramie muszą być zaznaczone frakcje, a nie ilości wystąpień.
Użyjemy wcześniejszego przykładu:
par(mfrow=c(2,1))
hist(ciskobiet,breaks = 25, xlim = c(93,178), ylim = c(0,0.07), main="Ciśnienie skurczowe dla kobiet", probability = TRUE)
lines(density(ciskobiet))
hist(cismezczyzna, breaks = 25, xlim = c(93,178), ylim = c(0,0.07), main = "Ciśnienie skurczowe dla meżczyzn", probability = TRUE)
lines(density(cismezczyzna))
Przy użyciu ggplot:
ggplot(dane, aes(x = cisnienie.skurczowe, y = ..density..)) +
geom_histogram(bins = 25, fill = "white", color = "black") +
geom_density(color = "blue", size = 1) +
facet_wrap(~plec, ncol = 1) +
labs(title = "Ciśnienie skurczowe", x = "Ciśnienie skurczowe", y = "Density")
3.1.3 Dystrybuanta
empiryczna
Jest to funkcja schodkowa, która reprezentuje proporcję obserwacji w
próbce, które są mniejsze lub równe \(t\).
To znaczy dla obserwacji \(X_1, X_2, \dots,
X_n\) dystrybuanta empiryczna definiowana jest w następujący
sposób: \[\hat{F}_n(t)=\frac{1}{n}\sum_{i=1}^n\mathbb{I}\{X_i<t\}.\]
Można ją wyznaczyć używając komendy ecdf(). Wynikiem
funkcji ecdf() jest obiekt klasy ecdf, dla
której zaimplementowana jest przeciążona funkcja plot()
prezentująca graficzną reprezentację obliczonej dystrybuanty
empirycznej.
plot(ecdf(ciskobiet), main="Ciśnienie skurczowe dla kobiet")
plot(ecdf(cismezczyzna), main = "Ciśnienie skurczowe dla meżczyzn", add=TRUE,col="red")
Dystrybuantę można wykorzystać do porównania kilka grup obserwacji.
Wystarczy na jednym wykresie przedstawić dystrybuanty dla różnych grup,
wykorzystując argument add=TRUE dla funkcji
plot() (wersja plot przeciążona, w zwykłym
plot nie ma tego argumentu).
Możemy również użyć biblioteki ggplot:
ggplot(dane, aes(x = cisnienie.skurczowe, color = plec)) +
stat_ecdf(geom = "step") +
labs(title = "Ciśnienie skurczowe dla kobiet i meżczyzn", x = "Ciśnienie skurczowe", y = "ECDF") +
scale_color_manual(values = c("kobieta" = "blue", "mezczyzna" = "red"))
3.1.4 Wykres
punktowy
Wykres punktowy, często zwany również wykresem rozrzutu jest jednym
podstawowych wykresów pokazujących zależność pomiędzy dwoma zmiennymi
ilościowymi.
W R najprościej taki wykres jest uzyskać poprzez użycie funkcji
plot(), ale można również zastosować funkcję
scatterplot() z biblioteki car, która posiada
kilka interesujących dodatkowych funkcji:
- prezentacja danych z podziałem na grupy,
reg.line=TRUE - dodanie trendu liniowego,smooth=TRUE- dodanie trendu nieliniowego.
library(car)
scatterplot(dane$cisnienie.rozkurczowe ~ dane$wiek | dane$st.cywilny,
xlab = "Wiek (lata)",
ylab = "Ciśnienie rozkurczowe (mmHg)",
main = "Ciśnienie rozkurczowe względem wieku w zależności od stanu cywilnego",
legend=list(coords="topleft", title="Stan cywilny"))
W formule po znaku |, możemy określić zmienną
warunkującą służąca do podziału na podgrupy. W powyższym przykładzie
pleć jest wyróżniona za pomocą równych kształtów i kolorów.
Przy użyciu ggplot:
ggplot(dane, aes(x = wiek, y = cisnienie.rozkurczowe, color = st.cywilny, shape = st.cywilny)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Scatter Plot of Diastolic Blood Pressure vs Age by Marital Status",
x = "Wiek",
y = "Ciśnienie rozkurczowe",
color = "Stan cywilny",
shape = "Stan cywilny")
Wykresy rozrzutu są przydatne do wykrywania i opisywania zależności
liniowych lub monotonicznych. Na przykład, na podstawie powyższego
wykresu można powiedzieć, że ciśnienie rozkurczowe dla singli z wiekiem
nie zmienia się, a dla osób w związku maleje wraz z wiekiem.
3.1.5 Wykresy
pudełkowe
Dla porównania rozkładów wartości w grupach popularną metodą jest
wykres pudełkowy, zwany też boxplotem lub wykresem ramka-wąsy. Za jego
pomocą można zobrazować położenie, rozmieszczenie i kształt rozkładu
empirycznego badanej cechy.
Jak się tworzy wykres pudełkowy
Elementy wykresu ramka-wąsy.
Pudełko - środkowa część wykresu, rozciąga się od pierwszego
kwartylu (Q1) do trzeciego kwartylu (Q3).
Wewnątrz pudełka jest linia, która wskazuje na medianę.
Wąsy rozciągają się od najmniejszej do maksymalnej wartości, lecz
nie dalej niż \(1.5 \cdot IQR\), gdzie
IQR =Q3-Q1(rozstęp międzykwartylowy). Pozostałe wartości są uznawane za
odstające.
Do uzyskania wykresy pudełkowego służy funkcja
boxplot(). Jako pierwszy argument możemy podać wektor,
listę wektorów, ramkę danych lub formułę ze zmiennymi, które mają się
znaleźć na wykresie.
boxplot(dane$cisnienie.skurczowe,dane$cisnienie.rozkurczowe, horizontal=TRUE, names=c("Skurczowe", "Rozkurczowe"))
lub za pomocą biblioteki ggplot:
dane_long <- data.frame(
Value = c(dane$cisnienie.skurczowe, dane$cisnienie.rozkurczowe),
Type = rep(c("Skurczowe", "Rozkurczowe"), each = nrow(dane))
)
ggplot(dane_long, aes(x = Value, y = Type)) +
geom_boxplot() +
labs(title = "Boxplot of Blood Pressure Types",
x = "Blood Pressure (mmHg)",
y = "")
Na powyższym wykresie można między innymi zauważyć, że ciśnienie
rozkurczowe przyjmuje generalnie wartości z zakresu od 65 do 98 i jest
mniejsze niż ciśnienie skurczowe.
boxplot(dane$wiek~dane$wyksztalcenie, varwidth=TRUE, col="lightblue", ylab="wiek", xlab="Wykształcenie", las=1)
Za pomocą ~ możemy porównać wiek w różnych
grupach wykształcenia. Dodatkowo, argument
varwidth=TRUE powoduje, że szerokość pudełka odpowiadać
będzie liczebności porównywanej grupy.
Za pomocą biblioteki ggplot:
ggplot(dane, aes(x = wyksztalcenie, y = wiek, fill = wyksztalcenie)) +
geom_boxplot(varwidth = TRUE) +
labs(x = "Wykształcenie",
y = "Wiek") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
3.1.6 Zadania - wykresy
dla danych jakościowych
Załaduj zbiór danych diamonds z biblioteki
ggplot. Następnie wykonaj poniższe zadania:
Dla zmiennej price utwórz histogram o 30
przedziałach i zaznacz na nim pionowymi kreskami średnią i medianę.
Nadaj tytuł, podpisz osie i dodaj legendę. Określ cechy otrzymanego
histogramu? Dorysuj na histogramie estymator gęstości
Wykonaj powyższy histogram w ggplocie.
Za pomocą wykresu punktowego zilustruj związek między ceną a
ilością karatów.
Wykonaj wykresy pudełkowe i określ, jak zmienia się cena w
zależności od szlifu (zmienna cut).
Wykonaj wykres dystrybuanty empirycznej dla ceny.
3.2 Statystyki dla danych
ilościowych
Są to liczbowe miary, opisujące podstawowe własności cech rozkładów.
Często chcielibyśmy mieć liczbową odpowiedź, gdzie leży “centrum” lub
“środek” próby, jak duże jest rozproszenie cechy w próbie wokół owego
“centrum” oraz jaki jest kształt rozkładu?
Wskaźniki odnoszące się do “centrum” lub “środka” próby nazywamy
wskaźnikami położenia, wskaźniki określające rozproszenie cech
nazywamy wskaźnikami rozproszenia. Natomiast wartości
określające kształt nazywamy wskaźnikami kształtu.
3.2.1 Wskaźniki
położenia
3.2.1.1 Średnia
Niech \(x_1, \dots, x_n\) będzie
próbą o liczności \(n\). Wtedy \[\bar{x}=\frac{1}{n}\sum_{i=1}^n x_i\]
nazywamy średnią arytmetyczną wartości cechy w próbie.
Aby obliczyć średnią wystarczy użyć funkcji mean().
Przykład: Rozważmy rozkład miesięcznych zasadniczych
wynagrodzeń pracowników z wyższym wykształceniem zatrudnionych w pewnej
firmie:
salary <- c(rep(3500,6), rep(4000,8), rep(4100, 7), rep(4500, 4), rep(5000, 3), rep(6000,2), 13000 )
mean(salary) #średnie wynagrodzenie
## [1] 4506.452
Zauważmy, że histogram jest prawostronnie skośny i posiada wartość
odstającą na prawo. W rezultacie wartość średnia jest znacznie bardziej
przesunięta na prawo od mody. Wyobraźmy sobie sytuację, że młody
absolwent wyższej uczelni zgłasza się na rozmowę kwalifikacyjną do pracy
w tej firmie. Kandydat dowiaduje się że będzie zarabiać 3500 zł, ale że
średnie miesięczne wynagrodzenie w wynosi 4500 zł. Zatem, powiada
wiceprezes firmy, ma pan przed sobą wspaniałe możliwości awansu i
znaczne wyższego wynagrodzenia. Czego kandydat nie słyszy, to tego, że
około 2/3 pracowników tej firmy zarabia mniej niż 4100 zł. Wysoka
średnia jest konsekwencją bardzo wysokich zarobków kierownictwa.
Kandydat dobrze by uczynił, gdyby zapytał się o inny wskaźnik położenia
zwany medianą.
3.2.1.2 Mediana
Medianą jest środkowa wartość (drugi kwartyl) próby uporządkowanej
niemalejąco, od wartości najmniejszej do wartości największej. Niech
\(x_{(1)}, \dots, x_{(n)}\) będzie
uporządkowana niemalejąco próba naszych danych. Wtedy medianę w próbie
nazywamy następującą wielkość: \[x_{\text
{med }}= \begin{cases}x_{((n+1) / 2)}, & \text { gdy } n \text {
nieparzyste } \\ \frac{1}{2}\left(x_{(n / 2)}+x_{(n / 2+1)}\right),
& \text { gdy } n \text { parzyste. }\end{cases}\]
## [1] 4100
Ciąg dalszy przykładu:
Mediana znacznie lepiej oddaje zarobkowe perspektywy kandydata.
Wszystko dlatego, że mediana jest niewrażliwa na wartości
odstające.
Warto jednak zwrócić uwagę, że jeżeli próba jest mała i kolejne
uporządkowane elementy próby są od siebie dość odległe, to mediana jest
niestabilnym wskaźnikiem położenia. W celu uniknięcia niestabilności
można użyć jeszcze innego wskaźnika położenia jakim jest średnia
ucinana.
Ogólniej, można używać statystyki kwartyli, które dla
uporządkowanego zbioru danych wyznaczają miejsca gdzie dane są dzielone
na poszczególne segmenty. Na przykład:
pierwszy kwartyl (Q1), to wartość poniżej której znajduje się 25%
danych,
trzeci kwartyl (Q3), to wartość poniżej której znajduje się 75%
danych.
3.2.1.3 Średnia
ucinana
Powstaje ona przez obliczenie średniej próby powstałej przez
usunięcie z próby oryginalnej \(k\)
wartości najmniejszych i \(k\) wartości
największych:
\[\bar{x}_{tk}=\frac{1}{n-2k}\sum_{i=k+1}^{n-k}
x_{(i)}.\] Aby obliczyć średnia ucinaną w R, wystarczy użyć
argumentu trim do funkcji mean() w którym
podamy frakcje elementów, o które chcemy uciąć naszą średnią.
## [1] 4168
Podamy przykład jeszcze jednego wskaźnika położenia. Którym warto
posłużyć się na przykład wtedy, gdy skrajne obserwacje charakteryzują
się dużą niepewnością co do ich rzeczywistych wartości.
3.2.1.4 Średnia
Winsorowska
Podobnie jak w przypadku średniej ucinanej, obliczanie średniej
winsorowskiej opiera się na \(n-2k\)
środkowych elementach próby, ale uwzględniając fakt pojawienia się w
próbie \(k\) wartości większych niż
\(x_{(k+1)}\) oraz \(k\) wartości mniejszych niż \(x_{(n-k)}\) przy obliczeniu średniej tak
jakby \(x_{(k+1)}\) oraz \(x_{(n-k)}\) wystąpiły dodatkowo \(k\) razy.
Formalnie, można wyrazić ją za pomocą poniższego wzoru: \[\bar{x}_{w k}=\frac{1}{n}\left[(k+1) \cdot
x_{(k+1)}+\sum_{i=k+2}^{n-k-1} x_{(i)}+(k+1)\cdot x_{(n-k)}\right]
.\] W R należy użyć funkcji winsor.mean() z pakietu
psych.
## [1] 4183.871
3.2.2 Wskaźniki
rozproszenia
Zredukowanie informacji o rozkładzie cech w próbie do wartości
średniej jest zdecydowanie zbyt drastyczne. Potrzebujemy jeszcze
informacji o stopniu rozproszenia próby wokół średniej.
3.2.2.1 Rozstęp
Rozstęp to różnica pomiędzy największą, a najmniejsza
wartością w próbie.
## [1] 9500
3.2.2.2 Wariancja (drugi
moment centralny)
Gdy wariancja jest mała, to oznacza to, że obserwacje są bardziej
skupione wokół centrum. Natomiast, gdy wariancja jest duża to rozstęp
jest duży i obserwacje znajdują się szeroko wokół swojego centrum.
Wariancję (var()) w próbie nazwiemy następującą
wielkość: \[s^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})^2.\]
Pierwiastek z wariancji nazywamy odchyleniem standardowym
(funkcja w R to sd()).
Warto zauważyć, że wariancja jest wskaźnikiem, który nie jest odporny
na wartości odstające.
## [1] 60.32462
## [1] 246.904
par(mfrow=c(2,1))
hist(dane$cisnienie.rozkurczowe, xlim = c(57,178))
hist(dane$cisnienie.skurczowe, xlim = c(57,178))
Na podstawie histogramów również jesteśmy w stanie ocenić, czy dane
mają dużą czy małą wariancję.
3.2.3 Wskaźniki
kształtu
Informują nasz jaki jest kształt rozkładu naszych danych
3.2.3.1 Skośność (Trzeci
moment centralny)
Jest miarą asymetrii. Wyraża się poniższym wzorem \[\gamma=
\frac{n\sum_{i=1}^n(x_i-\bar{x})^3}{(n-1)(n-2)s^3}.\] Może
przyjmować wartości dodatnie, ujemne i zerowe mając wpływ na to jak
wygląda rozkład:
- gdy \(\gamma=0\), to rozkład jest
symetryczny,
- gdy \(\gamma<0\), to rozkład
jest lewo skośny,
- gdy \(\gamma>0\), to rozkład
jest prawo skośny.
Wpływ skośności na kształt rozkładu
Aby obliczyć skośność w R wystarczy użyć funkcji
skewness().
3.2.3.2 Kurtoza
Wyraża się poniższym wzorem:
\[\text{kurtoza}=
\frac{n(n+1)\sum_{i=1}^n(x_i-\bar{x})^4}{(n-1)(n-2)(n-3)s^4}-3\frac{(n-1)^2}{(n-2)(n-3)}.\]
Również może przyjmować wartości dodatnie, ujemne i zerowe. W zależności
od tego wpływa na kształt ogonów rozkładu.
Wpływ kurtozy na rozkład
Jeżeli kurtoza jest ujemna, to wtedy ogony rozkładu stają się
ciężkie, coraz więcej obserwacji przyjmuje wartości skrajne.
Dla rozkładu normalnego kurtoza=3. Dlatego też często stosuje się
skorygowaną kurtozę=kurtoza-3, która dla rozkładu normalnego
wynosi 0. Następnie względem niego porównujemy ciężar ogonów
rozkładu.
Kurtozę można obliczyć za pomocą funkcji kurtosis z
pakietu moments.
4 Wstępna analiza
danych
W większości przypadków zebrane dane wymagają wstępnego przetwarzania
zanim zostanie na nich wykonana analiza statystyczna.
4.1 Kwestia brakujących
obserwacji
Dość powszechnym zjawiskiem jest brak co najmniej jednej wartości
charakteryzującej analizowane próbki. Przyczyną może być np. jakiś błąd
w procesie gromadzenia informacji, niemożliwość wykorzystania pewnych
form pomiarów lub podczas wypełniania formularza pewne pola mogły
pozostać puste.
W R brakujące dane oznaczane są symbolem NA.
Zignorowanie tych wartości może prowadzić do trudności z obliczeniem
podstawowych statystyk lub do błędnych wnioskowań.
Wczytajmy zbiór danych airquality z brakującymi
obserwacjami.
## [1] 153 6
W celu sprawdzenia czy w zbiorze danych są i ile jest brakujących
obserwacji możemy użyć funkcji sum(is.na()):
## [1] 44
Zatem, posiadamy 44 brakujące obserwacje. Możemy również użyć funkcji
summary(), która również może służyć do wykrycia
brakujących wartości:
## Ozone Solar.R Wind Temp
## Min. : 1.00 Min. : 7.0 Min. : 1.700 Min. :56.00
## 1st Qu.: 18.00 1st Qu.:115.8 1st Qu.: 7.400 1st Qu.:72.00
## Median : 31.50 Median :205.0 Median : 9.700 Median :79.00
## Mean : 42.13 Mean :185.9 Mean : 9.958 Mean :77.88
## 3rd Qu.: 63.25 3rd Qu.:258.8 3rd Qu.:11.500 3rd Qu.:85.00
## Max. :168.00 Max. :334.0 Max. :20.700 Max. :97.00
## NA's :37 NA's :7
## Month Day
## Min. :5.000 Min. : 1.0
## 1st Qu.:6.000 1st Qu.: 8.0
## Median :7.000 Median :16.0
## Mean :6.993 Mean :15.8
## 3rd Qu.:8.000 3rd Qu.:23.0
## Max. :9.000 Max. :31.0
##
Jednym z najprostszych sposobów jest usunięcie całych wierszy (lub
kolumn) ze zbioru danych. Można to zrobić za pomocą funkcji
na.omit()
## [1] 111 6
Pomysł ten jest zdecydowanie nie najlepszym rozwiązaniem. Usunięcie
kilkudziesięciu procent zebranych wyników, ponieważ brakuje jednego
pomiaru, prowadzi do niepotrzebnego pogorszenia właściwości procedur
statystycznych, zwiększenia wariancji estymatorów, zmniejszenia mocy
testów, itp. Zatem usunięcie wierszy w których pojawiają się wartości
NA ma sens, gdy jest ich mało i pojawiają się one losowo
(Missing at Random). Zamiast usuwać całe wiersze, można w ich miejsce
wstawić sztucznie inny pomiar. Można to zrobić na dwa sposoby:
W miejsce brakującej wartości wstawia się wartość
charakterystyczną dla danej zmiennej. Najczęściej medianę, średnią,
ewentualnie modę dla danej zmiennej. Bardziej wyrafinowanym sposobem
jest dopasowanie modelu regresji i wstawienie wartości wyznaczonej z
modelu regresji (o tym w przyszłości).
Wygenerowanie wielu replikacji (powtórzeń) zbioru danych z
powtórzeniami. Dla każdego powtórzenia brakujące obserwacje zastępowane
są przez obserwacje wygenerowane z rozkładu określonego dla danej
zmiennej. Kolejne kroki analizy są wykonywane dla każdej z replikacji
zbioru danych, a następnie wyniki poszczególnych kroków są ze sobą
zestawiane (więcej na temat tej metody w kolejnych
laboratoriach).
Uwaga: W przypadku, gdy analizujemy zbiór danych
najlepiej dowiedzieć się o nim jak najwięcej, ponieważ w przypadku, gdy
brakujące wartości powstają w inny sposób na przykład dlatego, że
wychodzą poza skalę czujnika pomiaru, to wstawianie tam średniej/mediany
jest błędne.
Teraz skupimy się na pierwszej metodzie. W celu podstawienia za
brakujące wartości innych wygodnie jest posłużyć się funkcją
inpute() z pakiety Hmisc. Za pomocą argumentu
fun= możemy wybrać, czy podstawiamy średnią
(mean), czy medianę (jest ona domyślnym parametrem tej
funkcji). Jeżeli fun=random, to brakujące wartości będą
zastąpione przez pozostałe wartości z próby, wybrane metodą losowania ze
zwracaniem.
4.2 Zadania - brakujące
dane
Załaduj zbiór danych mieszkania z pakietu
Przewodnik. Następnie dla kolumny cena dla indeksów
c(1,2,3,50, 60, 80, 100,105, 110, 120, 180, 181,188, 190, 195,196, 197, 198, 199, 200)
wpisz wartość NA. Następnie, usuń i uzupełnij brakujące
dane poznanymi metodami. Porównaj otrzymany zestaw danych do
oryginalnego. Między innymi oblicz błąd tj. wartość bezwzględną różnicy
średniej i mediany przed usunięciem danych i po zastosowaniu metody
radzenia sobie z wartościami NA. Dla której metody błąd był
największy, a dla której najmniejszy?
5 Skalowanie
Skalowanie zmiennych jest ważnym krokiem w przygotowaniu danych do
analizy, ale nie zawsze jest konieczne.
Skalujemy zmienne gdy chcemy użyć ich:
do algorytmów na odległość (\(k\)-means, KNN, SVD),
regresja (liniowa, logistyczna, grzbietowa, LASSO),
analiza składowych głównych (PCA),
sieci neuronowe.
Nie skalujemy zmiennych: gdy chcemy użyć ich:
drzew losowych, lasów losowych, Gradient Boosting Machines
(GBM),
algorytmów, które działają na bazie rang lub
klasyfikacji.
Pojawiło się tutaj kilka nowych skrótów od nazw algorytmów, które w
przyszłość może większość z Was pozna (np. na II stopniu na Statistical
Learning, Machine Learning, Wstęp do analizy danych, Deep Learning). Na
razie jest to informacja tylko poglądowa. Niestety, nie damy rady omówić
ich na tym przedmiocie.
5.1 Normalizacja
Zakresy zmiennych na ogół się różnią (np. z powodu
jednostek/zakresów), co może wpływać na to, że niektóre zmienne w
większym zakresie będą miały wpływ na wyniki.
Normalizacja min-max skaluje wartości do przedziału \([0,1]\) według poniższego wzoru: \[x_{\text{norm}}^{(i)}=\frac{x^{(i)}-x_{\min}}{x_{\max}-x_{\min}},\]
gdzie \(x^{(i)}\) oznacza daną próbkę,
\(x_{\max}\), \(x_{\min}\) to odpowiednio maksimum i
minimum dla wartości próbki w danej kolumnie.
Aby znormalizować dane w poniższym przykładzie użyjemy funkcji
lapply():
df <- data.frame("Wiek"=c(40, 13, 15, 45),
"Wzrost[cm]"=c(168, 146, 170, 165),
"Waga[kg]"=c(85,27, 42, 67))
# Normalizacja
normalize <- function(x) {
return((x - min(x)) / (max(x) - min(x)))
}
df_normalized <- as.data.frame(lapply(df, normalize))
df_normalized
## Wiek Wzrost.cm. Waga.kg.
## 1 0.84375 0.9166667 1.0000000
## 2 0.00000 0.0000000 0.0000000
## 3 0.06250 1.0000000 0.2586207
## 4 1.00000 0.7916667 0.6896552
Normalizacja utrudnia interpretację danych, ponieważ są one oderwane
od jednostek i skali. W szczególności najmłodsza osoba teraz ma wiek 0,
a najstarsza 1.
Normalizacja jest bardzo wrażliwa na wartości odstające, jeżeli
chociaż jedna obserwacja znacznie odbiega od pozostałych, to wpłynie to
na proces normalizacji.
5.2 Standaryzacja
Standaryzacja skaluje wartości tak, żeby miały średnią \(0\) i odchylenie standardowe \(1\).
Procedurę standaryzacji można zapisać za pomocą poniższego wzoru:
\[x_{\text{std}}^{(i)}=\frac{x^{(i)}-\bar{x}^{(i)}}{s^{(i)}},\]
gdzie \(\bar{x}^{(i)}\) jest średnia
próbek z kolumny \(i\), a \(s^{(i)}\) jest związanym z nią odchyleniem
standardowym.
Aby zestandaryzować dane możemy użyć funkcji
scale():
## Wiek Wzrost.cm. Waga.kg.
## 1 0.7078014 0.5214718 1.1531499
## 2 -0.9186358 -1.4737246 -1.0950079
## 3 -0.7981590 0.7028533 -0.5135878
## 4 1.0089934 0.2493996 0.4554457
Standaryzacja zachowuje informacje o wartościach odstających, ale nie
jest na nie tak podatna. Pokazuje również jak daleko w jakim kierunku od
średniej (tj. \(0\)) są obserwacje.
5.3 Zadania
Zadanie 1 Za pomocą funkcji cov()
wykonaj macierz kowariancji dla zbioru danych daneSoc z
pakietu Przewodnik dla danych nie skalowanych i po
standaryzacji. Do macierz kowariancji weź wszystkie zmienne ilościowe.
Przypomnij co to jest macierz kowariancji i jakie ma własności. Porównaj
dwie otrzymacie macierze, czym się różnią i jakie wnioski na ich
podstawie możny wyciągnąć?
6 Zadanie - wstępna
analiza danych
W zależności od ostatniej cyfry indeksu wybierz odpowiedni zbiór
danych. Zbiory danych znajdują się na stronie przedmiotu. Pod linkami
niżej znajdują się opisy zawartości zbiorów danych:
0. Australian Health Service Utilization Data - dane0.csv - dane
0;
1. House Prices in the city of Windsor, Canada - dane1.csv - dane
1;
2. Doctoral Publications - dane2. csv - dane
2;
3. Travel Mode Choice Data - dane3.csv - dane
3;
4. Economic Journals Data Set - dane4.csv - dane
4;
5. Wife Working Hours - dane5.csv - dane
5;
6. Epilepsy Data - dane6.csv - dane
6;
7. Days not Spent at School - dane7.csv - dane
7;
8. Porsche and Jaguar Prices - dane8.csv - dane
8;
9. Overdrawn Checking Account - dane9.csv - dane
9.
Załaduj zbiór danych i wykonaj następujące kroki:
a. Opisz zbiór danych, jeśli jest taka informacja to opisz sposób, w
jaki dane zostały pozyskane.
b. Określ które zmienne są ilościowe, a które są jakościowe. Sprawdź,
czy są wartości brakujące, jeśli tak to użyj jednej z podanych metod
uzupełniania danych.
Uwaga: często zmienne jakościowe są przez język R
rozpoznane jako charakter. Aby R interpretował je jako
zmienne jakościowe musimy wykonać as.factor(zmiennej) i
przypisać ją do “starej” zmiennej.
c. Wykonaj przykładowe wykresy i tabele liczności dla wybranych
zmiennych jakościowych i opisz co na ich podstawie możemy powiedzieć o
posiadanym zbiorze danych. Nie trzeba wykonywać wykresów dla wszystkich
zmiennych jakościowych, możemy wybrać takie, aby móc wyciągnąć jakieś
logiczne wnioski.
d. Wykonaj wykresy poznane w tym laboratorium dla wybranej zmiennej
ilościowej. Oblicz wskaźniki położenia (średnia, mediana, moda),
rozrzutu (odchylenie standardowe, wariancja), skośności i kurtozy.
Porównaj obliczone wskaźniki z wyglądem histogramu, opisując, jak
wskaźniki te wpływają na kształt i interpretację wykresu.
1 Statystyki opisowe
Język R głównie wykorzystuj się do przetwarzania, analizy i wizualizacji danych.
Po pierwsze musimy ustalić z jakim typem danych mamy do czynienia. Wyróżniamy:
dane ilościowe - są to liczby odpowiadające mierzalnym wielkościom, np. temperatura, liczba lat, waga w kilogramach,
dane jakościowe - są to cechy jakościowe obiektów, np. płeć, marka samochodu, poziom wykształcenia.
Aby określić jakie dane mamy w naszym zbiorze danych najlepiej jest
wywołać funkcje str() oraz summary(), np:
library("Przewodnik") #biblioteka, z której wykorzystamy przykładowe dane
dane <- daneSoc #daneSoc - jedne z danych dostępne w z tej biblioteki
str(dane)## 'data.frame': 204 obs. of 7 variables:
## $ wiek : int 70 66 71 57 45 48 53 61 63 47 ...
## $ wyksztalcenie : Factor w/ 4 levels "podstawowe","srednie",..: 4 4 4 2 2 2 2 3 2 3 ...
## $ st.cywilny : Factor w/ 2 levels "singiel","w zwiazku": 2 2 1 2 2 2 2 2 2 2 ...
## $ plec : Factor w/ 2 levels "kobieta","mezczyzna": 2 1 1 2 1 2 2 2 1 2 ...
## $ praca : Factor w/ 2 levels "nie pracuje",..: 2 2 2 2 2 1 2 2 2 2 ...
## $ cisnienie.skurczowe : int 143 123 167 150 130 138 149 103 125 158 ...
## $ cisnienie.rozkurczowe: int 83 80 80 87 83 75 80 63 78 89 ...## wiek wyksztalcenie st.cywilny plec
## Min. :22.00 podstawowe:93 singiel :120 kobieta : 55
## 1st Qu.:30.00 srednie :55 w zwiazku: 84 mezczyzna:149
## Median :45.00 wyzsze :34
## Mean :43.16 zawodowe :22
## 3rd Qu.:53.00
## Max. :75.00
## praca cisnienie.skurczowe cisnienie.rozkurczowe
## nie pracuje : 52 Min. : 93.0 Min. : 57.00
## uczen lub pracuje:152 1st Qu.:126.0 1st Qu.: 77.00
## Median :137.5 Median : 80.00
## Mean :137.0 Mean : 80.43
## 3rd Qu.:148.0 3rd Qu.: 86.25
## Max. :178.0 Max. :107.00Widzimy zatem, że wykształcenie, stan cywilny, płeć i praca są zmiennymi jakościowymi, natomiast wiek, ciśnienie skurczowe i ciśnienie rozkurczowe są zmiennymi ilościowymi.
2 Dane jakościowe
2.1 Statystyki dla
zmiennych jakościowych
Dane jakościowe są podzielone na kategorie i zwykle są analizowane
przy użyciu liczności.
2.1.1 Tabela
liczności
Zliczenie liczby obserwacji w każdej kategorii dokonuje się za pomocą
funkcji table().
##
## kobieta mezczyzna
## 55 149
Możesz także analizować dwie zmienne jakościowe jednocześnie tzn.
macierz liczności:
##
## singiel w zwiazku
## kobieta 38 17
## mezczyzna 82 67
Można również tworzyć table wielowymiarowe:
## , , = singiel
##
##
## podstawowe srednie wyzsze zawodowe
## kobieta 18 9 8 3
## mezczyzna 46 17 12 7
##
## , , = w zwiazku
##
##
## podstawowe srednie wyzsze zawodowe
## kobieta 4 7 2 4
## mezczyzna 25 22 12 8
2.2 Wykresy dla danych
jakościowych
Wykresy dla danych jakościowych pozwalają na szybką interpretację i
porównanie różnic między kategoriami. Oto kilka podstawowych typów
wykresów używanych do wizualizacji danych jakościowych:
2.2.1 Wykresy
słupkowe
Na jednej osi wykresów słupkowych mamy częstotliwość występowania
danej kategorii, na drugiej osi natomiast kategorię. Możemy mieć wykresy
słupkowe pionowe i poziome.
Służy nam do tego funkcja barplot() do której przekazuje
się table(dane_jakościowe), wektor częstotliwości danych
jakościowych. Na przykład:
Dodatkowe argumenty, takie jak horiz=TRUE sprawia, że
słupki są poziome, las=1 wskazuje w jakim kierunku powinny
być etykiety przy osiach, a cex.names kontroluje rozmiar
etykiet przy osi pionowej.
Za pomocą biblioteki ggplot:
# Przetestuj poniższy kod
library(ggplot2)
ggplot(dane, aes(x = wyksztalcenie)) +
geom_bar() +
coord_flip() + # Odwrócenie osi
theme(axis.text.y = element_text(size = 10),
axis.title.y = element_text(size = 20)) # Dostosowanie rozmiaru tekstu etykiet
Możemy również wykonać bardziej złożone wykresy dla macierzy
liczności, tj. table(dane$plec, dane$wyksztalcenie).
t <- table(dane$plec, dane$wyksztalcenie)
barplot(t,las=1, beside=TRUE, col=c("cyan", "darkblue")) # beside=TRUE sprawia że słupki są obok siebie, dodajemy kolory i legendę
legend("topright", c("kobieta", "mezczyzna" ), fill=c("cyan", "darkblue"))
ggplot(dane, aes(x = wyksztalcenie, fill = plec)) +
geom_bar(position = position_dodge()) +
scale_fill_manual(values = c("cyan", "darkblue")) +
labs(x = "Wykształcenie", y = "Liczba", fill = "Płeć")
2.2.2 Wykresy kołowe
Zamiast udziału liczbowego możemy również określić udział procentowy
zmiennych jakościowych. Do tego świetnie nadaje się wykres kołowy.
Z użyciem ggplot:
ggplot(dane, aes(x="", y="wyksztalcenie", fill=wyksztalcenie))+
geom_bar(stat="identity",width = 1) +
coord_polar("y") + # pytanie: jakie działanie ma ta linijka?
labs(title = "Proporcje dotyczące wykształacenia")
Możemy również przekazać więcej zmiennych do tego typu wykresu.
Należy jednak pamiętać, że zbyt duża ilość zaciemni wykres.
pie(table(dane$wyksztalcenie, dane$plec), labels = c("kobieta, podstawowe","kobieta, średnie", "kobieta, wyższe", "kobieta, zadowe", "meżczyzna, podstawowe", "meżczyzna, srednie", "mężczyzna, wyższe", "meżczyzna, wyższe" ))
Domyślnie funkcja pie() przy macierzach nieliczność nie
podpisuje wycinków koła. Należy nazwy przekazać do funkcji za pomocą
argumentu labels.
Z użyciem ggplot:
dane_summary <- as.data.frame(table(dane$plec, dane$wyksztalcenie))
colnames(dane_summary) <- c("plec", "wyksztalcenie", "count")
dane_summary$label <- paste(dane_summary$plec, dane_summary$wyksztalcenie, sep = ", ")
ggplot(dane_summary, aes(x = "", y = count, fill = label)) +
geom_bar(width = 1, stat = "identity", color = "black") +
coord_polar(theta = "y")
2.2.3 Wykres
mozaikowy
Wykres mozaikowy jest wykorzystywany, aby zobrazować zależność
pomiędzy dwoma lub większą ilością zmiennych jakościowych.
Używając funkcji mozaicplot()można przedstawić
liczebność wartości w danym zbiorze za pomocą pola prostokątów. Na
podstawie tego wykresy można stwierdzić, czy dwie zmienne od siebie
zależą oraz jak często występują dane wartości/pary wartości.
Spójrzmy na poniższy wykres, który przedstawia częstotliwość
występowania różnych grup wykształcenia:
Uwaga Jeżeli dwie zmienne są niezależne, to mozaika
przypomina regularną kratę. Gdy pojawia się zależność pomiędzy
zmiennymi, zazwyczaj łatwo ją zauważyć, np:
Zauważmy, że najwięcej osób pracujących występuje wśród osób z
wykształceniem wyższym. Zwróćmy również uwagę, że zachowała się
szerokość prostokątów z wcześniejszego polecenia, zatem widzimy, że w
naszym zbiorze danych najwięcej jest osób z wykształceniem
podstawowym.
Można za pomocą wykresu mosaicplot() zaprezentować
wyniki dla większej ilości zmiennych jakościowych. Należy jednak uważać,
ponieważ umieszczenie bardzo dużej ilości zmiennych zmniejszy czytelność
wykresu.
mosaicplot(~wyksztalcenie+praca+st.cywilny, data = daneSoc, col=c("darkorange", "grey"))
legend("bottomleft", c("single", "w związku" ), fill=c("darkorange", "grey"))
Można spróbować dodać kolory i legendę w celu zachowania czytelności.
Na podstawie wykresu jesteśmy w takie powiedzieć, że najmniejsze
bezrobocie jest wśród singli.
2.3 Zadania - zmienne
jakościowe
- Załaduj poniższy zbiór danych:
set.seed(111111) # W miejsce "jedynek" wpisz swój numer indeksu
dane_j <- data.frame(
Płeć = sample(c("Mężczyzna", "Kobieta"), 100, replace = TRUE),
Rok = sample(c("1", "2", "3", "4"), 100, replace = TRUE),
Matematyka = sample(c("Lubi matematykę", "Nie lubi matematyki"), 100, replace = TRUE),
Programowanie = sample(c("Lubi programowanie", "Nie lubi programowania"), 100, replace = TRUE),
Sztuka = sample(c("Lubi sztukę", "Nie lubi sztuki"), 100, replace = TRUE)
)
# Podgląd danych
head(dane_j)
Zacznij od wyświetlenia podsumowania dla danych przy użyciu funkcji
summary(). Jeżeli wyświetlona dane są typu
character należy wymusić na programie, aby interpretować
odpowiednie kolumny jako zmienne ilościowe lub jakościowe, używając do
tego np. funkcji as.number() lub as.numeric().
Przykładowo,
Następnie dla powyższego zestawu danych:
Stwórz tabelę liczności i proporcji dla zmiennej Programowanie i
dla zmiennej Sztuka.
Stwórz wykresy słupkowe i kołowe dla obu zmiennych.
Oblicz tabele liczności i proporcji dla zmiennych Płeć i
Programowanie oraz Płeć i Sztuka.
Stwórz wykresy mozaikowe dla obu zestawów zmiennych.
Na podstawie uzyskanych tabel liczności, proporcji i wykresów,
zinterpretuj wyniki. Jakie są preferencje studentów dotyczące różnych
zajęć w zależności od płci? Jakie wnioski można wyciągnąć na podstawie
analizy?
Uwaga: Wykresy powinny mieć tytuły, podpisane osie i
zmienione kolory.
- Wykonaj wykresy z powyższego zadania w wykorzystaniem biblioteki
ggplot2.
2.1 Statystyki dla zmiennych jakościowych
Dane jakościowe są podzielone na kategorie i zwykle są analizowane przy użyciu liczności.
2.1.1 Tabela
liczności
Zliczenie liczby obserwacji w każdej kategorii dokonuje się za pomocą
funkcji table().
##
## kobieta mezczyzna
## 55 149
Możesz także analizować dwie zmienne jakościowe jednocześnie tzn.
macierz liczności:
##
## singiel w zwiazku
## kobieta 38 17
## mezczyzna 82 67
Można również tworzyć table wielowymiarowe:
## , , = singiel
##
##
## podstawowe srednie wyzsze zawodowe
## kobieta 18 9 8 3
## mezczyzna 46 17 12 7
##
## , , = w zwiazku
##
##
## podstawowe srednie wyzsze zawodowe
## kobieta 4 7 2 4
## mezczyzna 25 22 12 8
2.1.1 Tabela liczności
Zliczenie liczby obserwacji w każdej kategorii dokonuje się za pomocą
funkcji table().
##
## kobieta mezczyzna
## 55 149Możesz także analizować dwie zmienne jakościowe jednocześnie tzn. macierz liczności:
##
## singiel w zwiazku
## kobieta 38 17
## mezczyzna 82 67Można również tworzyć table wielowymiarowe:
## , , = singiel
##
##
## podstawowe srednie wyzsze zawodowe
## kobieta 18 9 8 3
## mezczyzna 46 17 12 7
##
## , , = w zwiazku
##
##
## podstawowe srednie wyzsze zawodowe
## kobieta 4 7 2 4
## mezczyzna 25 22 12 82.2 Wykresy dla danych jakościowych
Wykresy dla danych jakościowych pozwalają na szybką interpretację i porównanie różnic między kategoriami. Oto kilka podstawowych typów wykresów używanych do wizualizacji danych jakościowych:
2.2.1 Wykresy
słupkowe
Na jednej osi wykresów słupkowych mamy częstotliwość występowania
danej kategorii, na drugiej osi natomiast kategorię. Możemy mieć wykresy
słupkowe pionowe i poziome.
Służy nam do tego funkcja barplot() do której przekazuje
się table(dane_jakościowe), wektor częstotliwości danych
jakościowych. Na przykład:
Dodatkowe argumenty, takie jak horiz=TRUE sprawia, że
słupki są poziome, las=1 wskazuje w jakim kierunku powinny
być etykiety przy osiach, a cex.names kontroluje rozmiar
etykiet przy osi pionowej.
Za pomocą biblioteki ggplot:
# Przetestuj poniższy kod
library(ggplot2)
ggplot(dane, aes(x = wyksztalcenie)) +
geom_bar() +
coord_flip() + # Odwrócenie osi
theme(axis.text.y = element_text(size = 10),
axis.title.y = element_text(size = 20)) # Dostosowanie rozmiaru tekstu etykiet
Możemy również wykonać bardziej złożone wykresy dla macierzy
liczności, tj. table(dane$plec, dane$wyksztalcenie).
t <- table(dane$plec, dane$wyksztalcenie)
barplot(t,las=1, beside=TRUE, col=c("cyan", "darkblue")) # beside=TRUE sprawia że słupki są obok siebie, dodajemy kolory i legendę
legend("topright", c("kobieta", "mezczyzna" ), fill=c("cyan", "darkblue"))
ggplot(dane, aes(x = wyksztalcenie, fill = plec)) +
geom_bar(position = position_dodge()) +
scale_fill_manual(values = c("cyan", "darkblue")) +
labs(x = "Wykształcenie", y = "Liczba", fill = "Płeć")
2.2.2 Wykresy kołowe
Zamiast udziału liczbowego możemy również określić udział procentowy
zmiennych jakościowych. Do tego świetnie nadaje się wykres kołowy.
Z użyciem ggplot:
ggplot(dane, aes(x="", y="wyksztalcenie", fill=wyksztalcenie))+
geom_bar(stat="identity",width = 1) +
coord_polar("y") + # pytanie: jakie działanie ma ta linijka?
labs(title = "Proporcje dotyczące wykształacenia")
Możemy również przekazać więcej zmiennych do tego typu wykresu.
Należy jednak pamiętać, że zbyt duża ilość zaciemni wykres.
pie(table(dane$wyksztalcenie, dane$plec), labels = c("kobieta, podstawowe","kobieta, średnie", "kobieta, wyższe", "kobieta, zadowe", "meżczyzna, podstawowe", "meżczyzna, srednie", "mężczyzna, wyższe", "meżczyzna, wyższe" ))
Domyślnie funkcja pie() przy macierzach nieliczność nie
podpisuje wycinków koła. Należy nazwy przekazać do funkcji za pomocą
argumentu labels.
Z użyciem ggplot:
dane_summary <- as.data.frame(table(dane$plec, dane$wyksztalcenie))
colnames(dane_summary) <- c("plec", "wyksztalcenie", "count")
dane_summary$label <- paste(dane_summary$plec, dane_summary$wyksztalcenie, sep = ", ")
ggplot(dane_summary, aes(x = "", y = count, fill = label)) +
geom_bar(width = 1, stat = "identity", color = "black") +
coord_polar(theta = "y")
2.2.3 Wykres
mozaikowy
Wykres mozaikowy jest wykorzystywany, aby zobrazować zależność
pomiędzy dwoma lub większą ilością zmiennych jakościowych.
Używając funkcji mozaicplot()można przedstawić
liczebność wartości w danym zbiorze za pomocą pola prostokątów. Na
podstawie tego wykresy można stwierdzić, czy dwie zmienne od siebie
zależą oraz jak często występują dane wartości/pary wartości.
Spójrzmy na poniższy wykres, który przedstawia częstotliwość
występowania różnych grup wykształcenia:
Uwaga Jeżeli dwie zmienne są niezależne, to mozaika
przypomina regularną kratę. Gdy pojawia się zależność pomiędzy
zmiennymi, zazwyczaj łatwo ją zauważyć, np:
Zauważmy, że najwięcej osób pracujących występuje wśród osób z
wykształceniem wyższym. Zwróćmy również uwagę, że zachowała się
szerokość prostokątów z wcześniejszego polecenia, zatem widzimy, że w
naszym zbiorze danych najwięcej jest osób z wykształceniem
podstawowym.
Można za pomocą wykresu mosaicplot() zaprezentować
wyniki dla większej ilości zmiennych jakościowych. Należy jednak uważać,
ponieważ umieszczenie bardzo dużej ilości zmiennych zmniejszy czytelność
wykresu.
mosaicplot(~wyksztalcenie+praca+st.cywilny, data = daneSoc, col=c("darkorange", "grey"))
legend("bottomleft", c("single", "w związku" ), fill=c("darkorange", "grey"))
Można spróbować dodać kolory i legendę w celu zachowania czytelności.
Na podstawie wykresu jesteśmy w takie powiedzieć, że najmniejsze
bezrobocie jest wśród singli.
2.2.1 Wykresy słupkowe
Na jednej osi wykresów słupkowych mamy częstotliwość występowania danej kategorii, na drugiej osi natomiast kategorię. Możemy mieć wykresy słupkowe pionowe i poziome.
Służy nam do tego funkcja barplot() do której przekazuje
się table(dane_jakościowe), wektor częstotliwości danych
jakościowych. Na przykład:
Dodatkowe argumenty, takie jak horiz=TRUE sprawia, że
słupki są poziome, las=1 wskazuje w jakim kierunku powinny
być etykiety przy osiach, a cex.names kontroluje rozmiar
etykiet przy osi pionowej.
Za pomocą biblioteki ggplot:
# Przetestuj poniższy kod
library(ggplot2)
ggplot(dane, aes(x = wyksztalcenie)) +
geom_bar() +
coord_flip() + # Odwrócenie osi
theme(axis.text.y = element_text(size = 10),
axis.title.y = element_text(size = 20)) # Dostosowanie rozmiaru tekstu etykietMożemy również wykonać bardziej złożone wykresy dla macierzy
liczności, tj. table(dane$plec, dane$wyksztalcenie).
t <- table(dane$plec, dane$wyksztalcenie)
barplot(t,las=1, beside=TRUE, col=c("cyan", "darkblue")) # beside=TRUE sprawia że słupki są obok siebie, dodajemy kolory i legendę
legend("topright", c("kobieta", "mezczyzna" ), fill=c("cyan", "darkblue"))
ggplot(dane, aes(x = wyksztalcenie, fill = plec)) +
geom_bar(position = position_dodge()) +
scale_fill_manual(values = c("cyan", "darkblue")) +
labs(x = "Wykształcenie", y = "Liczba", fill = "Płeć") 2.2.2 Wykresy kołowe
Zamiast udziału liczbowego możemy również określić udział procentowy zmiennych jakościowych. Do tego świetnie nadaje się wykres kołowy.
Z użyciem ggplot:
ggplot(dane, aes(x="", y="wyksztalcenie", fill=wyksztalcenie))+
geom_bar(stat="identity",width = 1) +
coord_polar("y") + # pytanie: jakie działanie ma ta linijka?
labs(title = "Proporcje dotyczące wykształacenia")Możemy również przekazać więcej zmiennych do tego typu wykresu. Należy jednak pamiętać, że zbyt duża ilość zaciemni wykres.
pie(table(dane$wyksztalcenie, dane$plec), labels = c("kobieta, podstawowe","kobieta, średnie", "kobieta, wyższe", "kobieta, zadowe", "meżczyzna, podstawowe", "meżczyzna, srednie", "mężczyzna, wyższe", "meżczyzna, wyższe" ))
Domyślnie funkcja pie() przy macierzach nieliczność nie
podpisuje wycinków koła. Należy nazwy przekazać do funkcji za pomocą
argumentu labels.
Z użyciem ggplot:
dane_summary <- as.data.frame(table(dane$plec, dane$wyksztalcenie))
colnames(dane_summary) <- c("plec", "wyksztalcenie", "count")
dane_summary$label <- paste(dane_summary$plec, dane_summary$wyksztalcenie, sep = ", ")
ggplot(dane_summary, aes(x = "", y = count, fill = label)) +
geom_bar(width = 1, stat = "identity", color = "black") +
coord_polar(theta = "y") 2.2.3 Wykres mozaikowy
Wykres mozaikowy jest wykorzystywany, aby zobrazować zależność pomiędzy dwoma lub większą ilością zmiennych jakościowych.
Używając funkcji mozaicplot()można przedstawić
liczebność wartości w danym zbiorze za pomocą pola prostokątów. Na
podstawie tego wykresy można stwierdzić, czy dwie zmienne od siebie
zależą oraz jak często występują dane wartości/pary wartości.
Spójrzmy na poniższy wykres, który przedstawia częstotliwość występowania różnych grup wykształcenia:
Uwaga Jeżeli dwie zmienne są niezależne, to mozaika przypomina regularną kratę. Gdy pojawia się zależność pomiędzy zmiennymi, zazwyczaj łatwo ją zauważyć, np:
Zauważmy, że najwięcej osób pracujących występuje wśród osób z wykształceniem wyższym. Zwróćmy również uwagę, że zachowała się szerokość prostokątów z wcześniejszego polecenia, zatem widzimy, że w naszym zbiorze danych najwięcej jest osób z wykształceniem podstawowym.
Można za pomocą wykresu mosaicplot() zaprezentować
wyniki dla większej ilości zmiennych jakościowych. Należy jednak uważać,
ponieważ umieszczenie bardzo dużej ilości zmiennych zmniejszy czytelność
wykresu.
mosaicplot(~wyksztalcenie+praca+st.cywilny, data = daneSoc, col=c("darkorange", "grey"))
legend("bottomleft", c("single", "w związku" ), fill=c("darkorange", "grey"))
Można spróbować dodać kolory i legendę w celu zachowania czytelności. Na podstawie wykresu jesteśmy w takie powiedzieć, że najmniejsze bezrobocie jest wśród singli.
2.3 Zadania - zmienne jakościowe
- Załaduj poniższy zbiór danych:
set.seed(111111) # W miejsce "jedynek" wpisz swój numer indeksu
dane_j <- data.frame(
Płeć = sample(c("Mężczyzna", "Kobieta"), 100, replace = TRUE),
Rok = sample(c("1", "2", "3", "4"), 100, replace = TRUE),
Matematyka = sample(c("Lubi matematykę", "Nie lubi matematyki"), 100, replace = TRUE),
Programowanie = sample(c("Lubi programowanie", "Nie lubi programowania"), 100, replace = TRUE),
Sztuka = sample(c("Lubi sztukę", "Nie lubi sztuki"), 100, replace = TRUE)
)
# Podgląd danych
head(dane_j)Zacznij od wyświetlenia podsumowania dla danych przy użyciu funkcji
summary(). Jeżeli wyświetlona dane są typu
character należy wymusić na programie, aby interpretować
odpowiednie kolumny jako zmienne ilościowe lub jakościowe, używając do
tego np. funkcji as.number() lub as.numeric().
Przykładowo,
Następnie dla powyższego zestawu danych:
Stwórz tabelę liczności i proporcji dla zmiennej Programowanie i dla zmiennej Sztuka.
Stwórz wykresy słupkowe i kołowe dla obu zmiennych.
Oblicz tabele liczności i proporcji dla zmiennych Płeć i Programowanie oraz Płeć i Sztuka.
Stwórz wykresy mozaikowe dla obu zestawów zmiennych.
Na podstawie uzyskanych tabel liczności, proporcji i wykresów, zinterpretuj wyniki. Jakie są preferencje studentów dotyczące różnych zajęć w zależności od płci? Jakie wnioski można wyciągnąć na podstawie analizy?
Uwaga: Wykresy powinny mieć tytuły, podpisane osie i zmienione kolory.
- Wykonaj wykresy z powyższego zadania w wykorzystaniem biblioteki
ggplot2.
3 Dane ilościowe
3.1 Wykresy dla danych
ilościowych
Przedstawimy również kilka podstawowych wykresów wykonywanych dla
zmiennych ilościowych, które wykonuje się przy wstępnej analizie
danych.
3.1.1 Histogram
Jest jednym z najpopularniejszych statystyk graficznych. Przedstawia
rozkład wartości zmiennych ilościowych.
Jak działa histogram?
Histogram składa się z słupków, które ilustrują częstość występowania
danych w określonych przedziałach.
Konstrukcja histogramu:
Wybór ilości przedziałów
Dla każdego przedziału liczona jest liczba obserwacji, które
wpadają w dany zakres.
Każdy przedział jest reprezentowany przez słupek na wykresie,
gdzie wysokość słupka oznacza liczbę obserwacji w danym
przedziale.
Idea histogramu: \[f(x)=\frac{1}{n\cdot
\text{długość przedziału}}\#\{ \text{obserwacje które wpadają w dany
przedział zawierający x}\},\] gdzie \(n\) to ilość obserwacji w naszym zbiorze
danych.
Histogram wyznacza rozkład empiryczny
Oznacza to, że histogram pokazuje, jak często poszczególne wartości
(lub ich przedziały) występują w badanej próbce. Histogram pozwala
zaobserwować kluczowe cechy rozkładu, takie jak modę, skośność, kurtozę,
obecność wartości odstających oraz potencjalną multimodalność (więcej o
tych pojęciach w trakcie przykładów).
Uwaga W praktyce, gdy liczba danych i liczba
przedziałów są odpowiednio duże, histogram może być używany jako
estymator gęstości prawdopodobieństwa.
Histogram generujemy za pomocą funkcji hist().
Opcjonalnymi argumentami jest np. breaks, dzięki któremu
możemy określić liczbę lub sposób podziału na przedziały (może być to
liczba lub nazwa algorytmu). Jeżeli nie określimy liczby przedziałów, to
zostanie ona dobrana w zależności od liczby obserwacji. Za pomocą
argumentów freq i probability możemy określić,
czy na histogramie mają bym zaznaczone frakcje, czy liczebności
wystąpień.
Przykład:
hist(dane$cisnienie.skurczowe, main="Histogram jednomodalny", xlab="Ciśnienie skurczowe", ylab="Frakcja wystepowania", probability = TRUE)
W ggplot
ggplot(dane, aes(x=cisnienie.skurczowe,after_stat(density)))+
geom_histogram()+
labs(title="Histogram jednomodalny")
Zauważmy, że kształt tego histogramu jest w przybliżeniu symetryczny,
ma również jedno maksimum zwane również modą. Jest to
przykład histogramu jednomodalnego.
Na poniższym rysunku są zobrazowane inne kształty histogramów.
Przykłady kształtów histogramów
Średnia, mediana i moda w symetrycznych,
prawoskośnych i lewoskośnych danych
W przypadku symetrycznych histogramów moda, średnia i mediana
pokrywają się (tak jak na panelu A). Jeżeli histogram ma ogon, który
wydaje się przeciągnięty w prawo to nazywamy go
prawoskośnym (panel B). Dodanie większej ilości
wartości w pobliżu maksimum silniej wpływa na średnią niż na modę, czy
medianę. Analogiczny przypadek jest, gdy ogon pojawia się po lewej
stronie, to znaczy mamy więcej wartości w pobliżu minimum. Taki
histogram nazywamy lewoskośnym.
Histogram do porównania zbiorów danych dla różnych
grup
Wykres dwóch zestawów danych w postaci dwóch histogramów może być
bardzo łatwym pierwszym krokiem do porównania środka, rozrzutu i
kształtu każdego z zestawu danych.
library(RVAideMemoire)# Biblioteka do obliczania mody za pomocą funkcji mod().
par(mfrow=c(2,1))
kobieta=dane$plec=="kobieta"
ciskobiet=dane[kobieta,]$cisnienie.skurczowe
mezczyzna=dane$plec=="mezczyzna"
cismezczyzna=dane[mezczyzna,]$cisnienie.skurczowe
hist(ciskobiet, breaks=20, xlim = c(93,178), main="Ciśnienie skurczowe dla kobiet")
abline(v=mean(ciskobiet), col="red")
abline(v=median(ciskobiet),col="blue")
abline(v=mod(ciskobiet), col="green")
hist(cismezczyzna, breaks = 25, xlim = c(93,178), main = "Ciśnienie skurczowe dla mężczyzn")
abline(v=mean(cismezczyzna), col="red")
abline(v=median(cismezczyzna),col="blue")
abline(v=mod(cismezczyzna), col="green")
UWAGA: aby porównywać dwa histogramy ważne jest by
oś x dla oby wykresów była tak samo wyskalowana, można to ustalić za
pomocą argumentu xlim.
Histogram identyfikuje wartości odstające
Na powyższym przykładzie większość
obserwacji wpada w przedział pomiędzy 0 a 200 mm, natomiast jedna
wartość mieści się w przedział pomiędzy 550, a 600. Taką wartość, która
odbiega mocno od pozostałych wartości znajdujących się w zbiorze danych,
nazywamy wartością odstającą.
3.1.2 Jądrowy estymator
gęstości
Można o nim myśleć jako o gładkiej wersji histogramu. Jest bardziej
zaawansowanym narzędziem do przedstawienia empirycznego rozkładu
gęstości dla posiadanych danych.
Zauważmy, że histogram mocno zależy od wyboru ilości przedziałów, ale
również od tego, w jakich miejscach dokonujemy dyskretyzacji, co może
wpływać na jego kształt i prowadzić do błędnych wniosków.
Jądrowy estymator gęstości rozwiązuje problemy jakie pojawiają się
przy histogramach.
Idea jądrowego estymatora gęstości: \[\hat{f}(x)=\frac{1}{n\cdot\text{długość
przedziału}}\#\{\text{obserwacje które wpadają w małe otoczenie punktu
x}\}.\] Definicja jądrowego estymatora
gęstości
Niech \(X_1, \dots, X_n\) będzie
próbą danych, wtedy \[\hat{f}(x)=\frac{1}{nh}\sum_{i=1}^n
K\left(\frac{x-X_i}{h}\right),\] gdzie:
- \(\hat{f}(x)\) to estymowana
gęstość w punkcie \(x\),
- \(h\) to szerokość okna
(bandwidth),
- \(K\) to funkcja jądra (kernel
function),
nazywamy jądrowym estymatorem gęstości.
Jako funkcję jądra bardzo często przyjmuje się tak zwane jądro
Gaussa tj.
\[K(u)=\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}u^2}.\]
Szerokość okna h
Szerokość okna \(h\) kontroluje
stopień wygładzenia estymacji. Małe \(h\) prowadzi do estymacji o wysokiej
wariancji (szczegółowej, “szumowej”), natomiast duże \(h\) prowadzi do estymacji o wysokim
obciążeniu (gładkiej, “rozmytej”). Dla funkcji density()
(patrz przykłady niżej) można ją ustalić ręcznie lub wskazać regułę jej
wyboru, która określi najodpowiedniejszą szerokość, np:
- reguła kciuka -
bw="nrd0" - jest ona używana
domyślnie,
- reguła Scotta -
bw="nrd",
- cross-validation -
bw="ucv",
- estymator Sheathera-Jonesa -
bw="SJ".
Intuicyjnie:
Każdy punkt danych \(X_i\) jest
zastępowany przez funkcję jądra \(K\),
która jest przesunięta tak, aby była wycentrowana na \(X_i\). A następnie te funkcje są sumowane i
znormalizowane przez liczbę punktów i szerokość okna, aby uzyskać
końcową estymację gęstości w każdym punkcie \(x\).
Jądrowy estymator gęstości
W praktyce, aby użyć jądrowego estymatora gęstości wystarczy użyć
funkcji density().
par(mfrow=c(3,1))
plot(density(dane$wiek)) # domyślna wartość okna
plot(density(dane$wiek, bw=0.5)) #mała szerokość okna
plot(density(dane$wiek, bw=10)) #duża szerokość okna
Powyższy przykład ukazuje, jak wykres zmienia się w zależności od
szerokości okna. Jeżeli jest ona mała to na jądrowym estymatorze
gęstości widzimy bardzo dużo szczegółów. Często w praktyce nie jest to
korzystne, ponieważ chcemy wyciągnąć ogólne wnioski i zniwelować
obecność szumu wśród naszych obserwacji. Jeżeli okno jest zbyt duże to
tracimy informacje o kształcie rozkładu naszych danych. Dlatego tak
ważne jest znalezienie kompromisu pomiędzy tymi aspektami.
Jądrowy estymator gęstości często można przedstawiać na tym samym
wykresie co histogram. Pamiętajmy jednak, że w takim przypadku na
histogramie muszą być zaznaczone frakcje, a nie ilości wystąpień.
Użyjemy wcześniejszego przykładu:
par(mfrow=c(2,1))
hist(ciskobiet,breaks = 25, xlim = c(93,178), ylim = c(0,0.07), main="Ciśnienie skurczowe dla kobiet", probability = TRUE)
lines(density(ciskobiet))
hist(cismezczyzna, breaks = 25, xlim = c(93,178), ylim = c(0,0.07), main = "Ciśnienie skurczowe dla meżczyzn", probability = TRUE)
lines(density(cismezczyzna))
Przy użyciu ggplot:
ggplot(dane, aes(x = cisnienie.skurczowe, y = ..density..)) +
geom_histogram(bins = 25, fill = "white", color = "black") +
geom_density(color = "blue", size = 1) +
facet_wrap(~plec, ncol = 1) +
labs(title = "Ciśnienie skurczowe", x = "Ciśnienie skurczowe", y = "Density")
3.1.3 Dystrybuanta
empiryczna
Jest to funkcja schodkowa, która reprezentuje proporcję obserwacji w
próbce, które są mniejsze lub równe \(t\).
To znaczy dla obserwacji \(X_1, X_2, \dots,
X_n\) dystrybuanta empiryczna definiowana jest w następujący
sposób: \[\hat{F}_n(t)=\frac{1}{n}\sum_{i=1}^n\mathbb{I}\{X_i<t\}.\]
Można ją wyznaczyć używając komendy ecdf(). Wynikiem
funkcji ecdf() jest obiekt klasy ecdf, dla
której zaimplementowana jest przeciążona funkcja plot()
prezentująca graficzną reprezentację obliczonej dystrybuanty
empirycznej.
plot(ecdf(ciskobiet), main="Ciśnienie skurczowe dla kobiet")
plot(ecdf(cismezczyzna), main = "Ciśnienie skurczowe dla meżczyzn", add=TRUE,col="red")
Dystrybuantę można wykorzystać do porównania kilka grup obserwacji.
Wystarczy na jednym wykresie przedstawić dystrybuanty dla różnych grup,
wykorzystując argument add=TRUE dla funkcji
plot() (wersja plot przeciążona, w zwykłym
plot nie ma tego argumentu).
Możemy również użyć biblioteki ggplot:
ggplot(dane, aes(x = cisnienie.skurczowe, color = plec)) +
stat_ecdf(geom = "step") +
labs(title = "Ciśnienie skurczowe dla kobiet i meżczyzn", x = "Ciśnienie skurczowe", y = "ECDF") +
scale_color_manual(values = c("kobieta" = "blue", "mezczyzna" = "red"))
3.1.4 Wykres
punktowy
Wykres punktowy, często zwany również wykresem rozrzutu jest jednym
podstawowych wykresów pokazujących zależność pomiędzy dwoma zmiennymi
ilościowymi.
W R najprościej taki wykres jest uzyskać poprzez użycie funkcji
plot(), ale można również zastosować funkcję
scatterplot() z biblioteki car, która posiada
kilka interesujących dodatkowych funkcji:
- prezentacja danych z podziałem na grupy,
reg.line=TRUE - dodanie trendu liniowego,smooth=TRUE- dodanie trendu nieliniowego.
library(car)
scatterplot(dane$cisnienie.rozkurczowe ~ dane$wiek | dane$st.cywilny,
xlab = "Wiek (lata)",
ylab = "Ciśnienie rozkurczowe (mmHg)",
main = "Ciśnienie rozkurczowe względem wieku w zależności od stanu cywilnego",
legend=list(coords="topleft", title="Stan cywilny"))
W formule po znaku |, możemy określić zmienną
warunkującą służąca do podziału na podgrupy. W powyższym przykładzie
pleć jest wyróżniona za pomocą równych kształtów i kolorów.
Przy użyciu ggplot:
ggplot(dane, aes(x = wiek, y = cisnienie.rozkurczowe, color = st.cywilny, shape = st.cywilny)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Scatter Plot of Diastolic Blood Pressure vs Age by Marital Status",
x = "Wiek",
y = "Ciśnienie rozkurczowe",
color = "Stan cywilny",
shape = "Stan cywilny")
Wykresy rozrzutu są przydatne do wykrywania i opisywania zależności
liniowych lub monotonicznych. Na przykład, na podstawie powyższego
wykresu można powiedzieć, że ciśnienie rozkurczowe dla singli z wiekiem
nie zmienia się, a dla osób w związku maleje wraz z wiekiem.
3.1.5 Wykresy
pudełkowe
Dla porównania rozkładów wartości w grupach popularną metodą jest
wykres pudełkowy, zwany też boxplotem lub wykresem ramka-wąsy. Za jego
pomocą można zobrazować położenie, rozmieszczenie i kształt rozkładu
empirycznego badanej cechy.
Jak się tworzy wykres pudełkowy
Elementy wykresu ramka-wąsy.
Pudełko - środkowa część wykresu, rozciąga się od pierwszego
kwartylu (Q1) do trzeciego kwartylu (Q3).
Wewnątrz pudełka jest linia, która wskazuje na medianę.
Wąsy rozciągają się od najmniejszej do maksymalnej wartości, lecz
nie dalej niż \(1.5 \cdot IQR\), gdzie
IQR =Q3-Q1(rozstęp międzykwartylowy). Pozostałe wartości są uznawane za
odstające.
Do uzyskania wykresy pudełkowego służy funkcja
boxplot(). Jako pierwszy argument możemy podać wektor,
listę wektorów, ramkę danych lub formułę ze zmiennymi, które mają się
znaleźć na wykresie.
boxplot(dane$cisnienie.skurczowe,dane$cisnienie.rozkurczowe, horizontal=TRUE, names=c("Skurczowe", "Rozkurczowe"))
lub za pomocą biblioteki ggplot:
dane_long <- data.frame(
Value = c(dane$cisnienie.skurczowe, dane$cisnienie.rozkurczowe),
Type = rep(c("Skurczowe", "Rozkurczowe"), each = nrow(dane))
)
ggplot(dane_long, aes(x = Value, y = Type)) +
geom_boxplot() +
labs(title = "Boxplot of Blood Pressure Types",
x = "Blood Pressure (mmHg)",
y = "")
Na powyższym wykresie można między innymi zauważyć, że ciśnienie
rozkurczowe przyjmuje generalnie wartości z zakresu od 65 do 98 i jest
mniejsze niż ciśnienie skurczowe.
boxplot(dane$wiek~dane$wyksztalcenie, varwidth=TRUE, col="lightblue", ylab="wiek", xlab="Wykształcenie", las=1)
Za pomocą ~ możemy porównać wiek w różnych
grupach wykształcenia. Dodatkowo, argument
varwidth=TRUE powoduje, że szerokość pudełka odpowiadać
będzie liczebności porównywanej grupy.
Za pomocą biblioteki ggplot:
ggplot(dane, aes(x = wyksztalcenie, y = wiek, fill = wyksztalcenie)) +
geom_boxplot(varwidth = TRUE) +
labs(x = "Wykształcenie",
y = "Wiek") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
3.1.6 Zadania - wykresy
dla danych jakościowych
Załaduj zbiór danych diamonds z biblioteki
ggplot. Następnie wykonaj poniższe zadania:
Dla zmiennej price utwórz histogram o 30
przedziałach i zaznacz na nim pionowymi kreskami średnią i medianę.
Nadaj tytuł, podpisz osie i dodaj legendę. Określ cechy otrzymanego
histogramu? Dorysuj na histogramie estymator gęstości
Wykonaj powyższy histogram w ggplocie.
Za pomocą wykresu punktowego zilustruj związek między ceną a
ilością karatów.
Wykonaj wykresy pudełkowe i określ, jak zmienia się cena w
zależności od szlifu (zmienna cut).
Wykonaj wykres dystrybuanty empirycznej dla ceny.
3.2 Statystyki dla danych
ilościowych
Są to liczbowe miary, opisujące podstawowe własności cech rozkładów.
Często chcielibyśmy mieć liczbową odpowiedź, gdzie leży “centrum” lub
“środek” próby, jak duże jest rozproszenie cechy w próbie wokół owego
“centrum” oraz jaki jest kształt rozkładu?
Wskaźniki odnoszące się do “centrum” lub “środka” próby nazywamy
wskaźnikami położenia, wskaźniki określające rozproszenie cech
nazywamy wskaźnikami rozproszenia. Natomiast wartości
określające kształt nazywamy wskaźnikami kształtu.
3.2.1 Wskaźniki
położenia
3.2.1.1 Średnia
Niech \(x_1, \dots, x_n\) będzie
próbą o liczności \(n\). Wtedy \[\bar{x}=\frac{1}{n}\sum_{i=1}^n x_i\]
nazywamy średnią arytmetyczną wartości cechy w próbie.
Aby obliczyć średnią wystarczy użyć funkcji mean().
Przykład: Rozważmy rozkład miesięcznych zasadniczych
wynagrodzeń pracowników z wyższym wykształceniem zatrudnionych w pewnej
firmie:
salary <- c(rep(3500,6), rep(4000,8), rep(4100, 7), rep(4500, 4), rep(5000, 3), rep(6000,2), 13000 )
mean(salary) #średnie wynagrodzenie
## [1] 4506.452
Zauważmy, że histogram jest prawostronnie skośny i posiada wartość
odstającą na prawo. W rezultacie wartość średnia jest znacznie bardziej
przesunięta na prawo od mody. Wyobraźmy sobie sytuację, że młody
absolwent wyższej uczelni zgłasza się na rozmowę kwalifikacyjną do pracy
w tej firmie. Kandydat dowiaduje się że będzie zarabiać 3500 zł, ale że
średnie miesięczne wynagrodzenie w wynosi 4500 zł. Zatem, powiada
wiceprezes firmy, ma pan przed sobą wspaniałe możliwości awansu i
znaczne wyższego wynagrodzenia. Czego kandydat nie słyszy, to tego, że
około 2/3 pracowników tej firmy zarabia mniej niż 4100 zł. Wysoka
średnia jest konsekwencją bardzo wysokich zarobków kierownictwa.
Kandydat dobrze by uczynił, gdyby zapytał się o inny wskaźnik położenia
zwany medianą.
3.2.1.2 Mediana
Medianą jest środkowa wartość (drugi kwartyl) próby uporządkowanej
niemalejąco, od wartości najmniejszej do wartości największej. Niech
\(x_{(1)}, \dots, x_{(n)}\) będzie
uporządkowana niemalejąco próba naszych danych. Wtedy medianę w próbie
nazywamy następującą wielkość: \[x_{\text
{med }}= \begin{cases}x_{((n+1) / 2)}, & \text { gdy } n \text {
nieparzyste } \\ \frac{1}{2}\left(x_{(n / 2)}+x_{(n / 2+1)}\right),
& \text { gdy } n \text { parzyste. }\end{cases}\]
## [1] 4100
Ciąg dalszy przykładu:
Mediana znacznie lepiej oddaje zarobkowe perspektywy kandydata.
Wszystko dlatego, że mediana jest niewrażliwa na wartości
odstające.
Warto jednak zwrócić uwagę, że jeżeli próba jest mała i kolejne
uporządkowane elementy próby są od siebie dość odległe, to mediana jest
niestabilnym wskaźnikiem położenia. W celu uniknięcia niestabilności
można użyć jeszcze innego wskaźnika położenia jakim jest średnia
ucinana.
Ogólniej, można używać statystyki kwartyli, które dla
uporządkowanego zbioru danych wyznaczają miejsca gdzie dane są dzielone
na poszczególne segmenty. Na przykład:
pierwszy kwartyl (Q1), to wartość poniżej której znajduje się 25%
danych,
trzeci kwartyl (Q3), to wartość poniżej której znajduje się 75%
danych.
3.2.1.3 Średnia
ucinana
Powstaje ona przez obliczenie średniej próby powstałej przez
usunięcie z próby oryginalnej \(k\)
wartości najmniejszych i \(k\) wartości
największych:
\[\bar{x}_{tk}=\frac{1}{n-2k}\sum_{i=k+1}^{n-k}
x_{(i)}.\] Aby obliczyć średnia ucinaną w R, wystarczy użyć
argumentu trim do funkcji mean() w którym
podamy frakcje elementów, o które chcemy uciąć naszą średnią.
## [1] 4168
Podamy przykład jeszcze jednego wskaźnika położenia. Którym warto
posłużyć się na przykład wtedy, gdy skrajne obserwacje charakteryzują
się dużą niepewnością co do ich rzeczywistych wartości.
3.2.1.4 Średnia
Winsorowska
Podobnie jak w przypadku średniej ucinanej, obliczanie średniej
winsorowskiej opiera się na \(n-2k\)
środkowych elementach próby, ale uwzględniając fakt pojawienia się w
próbie \(k\) wartości większych niż
\(x_{(k+1)}\) oraz \(k\) wartości mniejszych niż \(x_{(n-k)}\) przy obliczeniu średniej tak
jakby \(x_{(k+1)}\) oraz \(x_{(n-k)}\) wystąpiły dodatkowo \(k\) razy.
Formalnie, można wyrazić ją za pomocą poniższego wzoru: \[\bar{x}_{w k}=\frac{1}{n}\left[(k+1) \cdot
x_{(k+1)}+\sum_{i=k+2}^{n-k-1} x_{(i)}+(k+1)\cdot x_{(n-k)}\right]
.\] W R należy użyć funkcji winsor.mean() z pakietu
psych.
## [1] 4183.871
3.2.2 Wskaźniki
rozproszenia
Zredukowanie informacji o rozkładzie cech w próbie do wartości
średniej jest zdecydowanie zbyt drastyczne. Potrzebujemy jeszcze
informacji o stopniu rozproszenia próby wokół średniej.
3.2.2.1 Rozstęp
Rozstęp to różnica pomiędzy największą, a najmniejsza
wartością w próbie.
## [1] 9500
3.2.2.2 Wariancja (drugi
moment centralny)
Gdy wariancja jest mała, to oznacza to, że obserwacje są bardziej
skupione wokół centrum. Natomiast, gdy wariancja jest duża to rozstęp
jest duży i obserwacje znajdują się szeroko wokół swojego centrum.
Wariancję (var()) w próbie nazwiemy następującą
wielkość: \[s^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})^2.\]
Pierwiastek z wariancji nazywamy odchyleniem standardowym
(funkcja w R to sd()).
Warto zauważyć, że wariancja jest wskaźnikiem, który nie jest odporny
na wartości odstające.
## [1] 60.32462
## [1] 246.904
par(mfrow=c(2,1))
hist(dane$cisnienie.rozkurczowe, xlim = c(57,178))
hist(dane$cisnienie.skurczowe, xlim = c(57,178))
Na podstawie histogramów również jesteśmy w stanie ocenić, czy dane
mają dużą czy małą wariancję.
3.2.3 Wskaźniki
kształtu
Informują nasz jaki jest kształt rozkładu naszych danych
3.2.3.1 Skośność (Trzeci
moment centralny)
Jest miarą asymetrii. Wyraża się poniższym wzorem \[\gamma=
\frac{n\sum_{i=1}^n(x_i-\bar{x})^3}{(n-1)(n-2)s^3}.\] Może
przyjmować wartości dodatnie, ujemne i zerowe mając wpływ na to jak
wygląda rozkład:
- gdy \(\gamma=0\), to rozkład jest
symetryczny,
- gdy \(\gamma<0\), to rozkład
jest lewo skośny,
- gdy \(\gamma>0\), to rozkład
jest prawo skośny.
Wpływ skośności na kształt rozkładu
Aby obliczyć skośność w R wystarczy użyć funkcji
skewness().
3.2.3.2 Kurtoza
Wyraża się poniższym wzorem:
\[\text{kurtoza}=
\frac{n(n+1)\sum_{i=1}^n(x_i-\bar{x})^4}{(n-1)(n-2)(n-3)s^4}-3\frac{(n-1)^2}{(n-2)(n-3)}.\]
Również może przyjmować wartości dodatnie, ujemne i zerowe. W zależności
od tego wpływa na kształt ogonów rozkładu.
Wpływ kurtozy na rozkład
Jeżeli kurtoza jest ujemna, to wtedy ogony rozkładu stają się
ciężkie, coraz więcej obserwacji przyjmuje wartości skrajne.
Dla rozkładu normalnego kurtoza=3. Dlatego też często stosuje się
skorygowaną kurtozę=kurtoza-3, która dla rozkładu normalnego
wynosi 0. Następnie względem niego porównujemy ciężar ogonów
rozkładu.
Kurtozę można obliczyć za pomocą funkcji kurtosis z
pakietu moments.
3.1 Wykresy dla danych ilościowych
Przedstawimy również kilka podstawowych wykresów wykonywanych dla zmiennych ilościowych, które wykonuje się przy wstępnej analizie danych.
3.1.1 Histogram
Jest jednym z najpopularniejszych statystyk graficznych. Przedstawia
rozkład wartości zmiennych ilościowych.
Jak działa histogram?
Histogram składa się z słupków, które ilustrują częstość występowania
danych w określonych przedziałach.
Konstrukcja histogramu:
Wybór ilości przedziałów
Dla każdego przedziału liczona jest liczba obserwacji, które
wpadają w dany zakres.
Każdy przedział jest reprezentowany przez słupek na wykresie,
gdzie wysokość słupka oznacza liczbę obserwacji w danym
przedziale.
Idea histogramu: \[f(x)=\frac{1}{n\cdot
\text{długość przedziału}}\#\{ \text{obserwacje które wpadają w dany
przedział zawierający x}\},\] gdzie \(n\) to ilość obserwacji w naszym zbiorze
danych.
Histogram wyznacza rozkład empiryczny
Oznacza to, że histogram pokazuje, jak często poszczególne wartości
(lub ich przedziały) występują w badanej próbce. Histogram pozwala
zaobserwować kluczowe cechy rozkładu, takie jak modę, skośność, kurtozę,
obecność wartości odstających oraz potencjalną multimodalność (więcej o
tych pojęciach w trakcie przykładów).
Uwaga W praktyce, gdy liczba danych i liczba
przedziałów są odpowiednio duże, histogram może być używany jako
estymator gęstości prawdopodobieństwa.
Histogram generujemy za pomocą funkcji hist().
Opcjonalnymi argumentami jest np. breaks, dzięki któremu
możemy określić liczbę lub sposób podziału na przedziały (może być to
liczba lub nazwa algorytmu). Jeżeli nie określimy liczby przedziałów, to
zostanie ona dobrana w zależności od liczby obserwacji. Za pomocą
argumentów freq i probability możemy określić,
czy na histogramie mają bym zaznaczone frakcje, czy liczebności
wystąpień.
Przykład:
hist(dane$cisnienie.skurczowe, main="Histogram jednomodalny", xlab="Ciśnienie skurczowe", ylab="Frakcja wystepowania", probability = TRUE)
W ggplot
ggplot(dane, aes(x=cisnienie.skurczowe,after_stat(density)))+
geom_histogram()+
labs(title="Histogram jednomodalny")
Zauważmy, że kształt tego histogramu jest w przybliżeniu symetryczny,
ma również jedno maksimum zwane również modą. Jest to
przykład histogramu jednomodalnego.
Na poniższym rysunku są zobrazowane inne kształty histogramów.
Przykłady kształtów histogramów
Średnia, mediana i moda w symetrycznych,
prawoskośnych i lewoskośnych danych
W przypadku symetrycznych histogramów moda, średnia i mediana
pokrywają się (tak jak na panelu A). Jeżeli histogram ma ogon, który
wydaje się przeciągnięty w prawo to nazywamy go
prawoskośnym (panel B). Dodanie większej ilości
wartości w pobliżu maksimum silniej wpływa na średnią niż na modę, czy
medianę. Analogiczny przypadek jest, gdy ogon pojawia się po lewej
stronie, to znaczy mamy więcej wartości w pobliżu minimum. Taki
histogram nazywamy lewoskośnym.
Histogram do porównania zbiorów danych dla różnych
grup
Wykres dwóch zestawów danych w postaci dwóch histogramów może być
bardzo łatwym pierwszym krokiem do porównania środka, rozrzutu i
kształtu każdego z zestawu danych.
library(RVAideMemoire)# Biblioteka do obliczania mody za pomocą funkcji mod().
par(mfrow=c(2,1))
kobieta=dane$plec=="kobieta"
ciskobiet=dane[kobieta,]$cisnienie.skurczowe
mezczyzna=dane$plec=="mezczyzna"
cismezczyzna=dane[mezczyzna,]$cisnienie.skurczowe
hist(ciskobiet, breaks=20, xlim = c(93,178), main="Ciśnienie skurczowe dla kobiet")
abline(v=mean(ciskobiet), col="red")
abline(v=median(ciskobiet),col="blue")
abline(v=mod(ciskobiet), col="green")
hist(cismezczyzna, breaks = 25, xlim = c(93,178), main = "Ciśnienie skurczowe dla mężczyzn")
abline(v=mean(cismezczyzna), col="red")
abline(v=median(cismezczyzna),col="blue")
abline(v=mod(cismezczyzna), col="green")
UWAGA: aby porównywać dwa histogramy ważne jest by
oś x dla oby wykresów była tak samo wyskalowana, można to ustalić za
pomocą argumentu xlim.
Histogram identyfikuje wartości odstające
Na powyższym przykładzie większość
obserwacji wpada w przedział pomiędzy 0 a 200 mm, natomiast jedna
wartość mieści się w przedział pomiędzy 550, a 600. Taką wartość, która
odbiega mocno od pozostałych wartości znajdujących się w zbiorze danych,
nazywamy wartością odstającą.
3.1.2 Jądrowy estymator
gęstości
Można o nim myśleć jako o gładkiej wersji histogramu. Jest bardziej
zaawansowanym narzędziem do przedstawienia empirycznego rozkładu
gęstości dla posiadanych danych.
Zauważmy, że histogram mocno zależy od wyboru ilości przedziałów, ale
również od tego, w jakich miejscach dokonujemy dyskretyzacji, co może
wpływać na jego kształt i prowadzić do błędnych wniosków.
Jądrowy estymator gęstości rozwiązuje problemy jakie pojawiają się
przy histogramach.
Idea jądrowego estymatora gęstości: \[\hat{f}(x)=\frac{1}{n\cdot\text{długość
przedziału}}\#\{\text{obserwacje które wpadają w małe otoczenie punktu
x}\}.\] Definicja jądrowego estymatora
gęstości
Niech \(X_1, \dots, X_n\) będzie
próbą danych, wtedy \[\hat{f}(x)=\frac{1}{nh}\sum_{i=1}^n
K\left(\frac{x-X_i}{h}\right),\] gdzie:
- \(\hat{f}(x)\) to estymowana
gęstość w punkcie \(x\),
- \(h\) to szerokość okna
(bandwidth),
- \(K\) to funkcja jądra (kernel
function),
nazywamy jądrowym estymatorem gęstości.
Jako funkcję jądra bardzo często przyjmuje się tak zwane jądro
Gaussa tj.
\[K(u)=\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}u^2}.\]
Szerokość okna h
Szerokość okna \(h\) kontroluje
stopień wygładzenia estymacji. Małe \(h\) prowadzi do estymacji o wysokiej
wariancji (szczegółowej, “szumowej”), natomiast duże \(h\) prowadzi do estymacji o wysokim
obciążeniu (gładkiej, “rozmytej”). Dla funkcji density()
(patrz przykłady niżej) można ją ustalić ręcznie lub wskazać regułę jej
wyboru, która określi najodpowiedniejszą szerokość, np:
- reguła kciuka -
bw="nrd0" - jest ona używana
domyślnie,
- reguła Scotta -
bw="nrd",
- cross-validation -
bw="ucv",
- estymator Sheathera-Jonesa -
bw="SJ".
Intuicyjnie:
Każdy punkt danych \(X_i\) jest
zastępowany przez funkcję jądra \(K\),
która jest przesunięta tak, aby była wycentrowana na \(X_i\). A następnie te funkcje są sumowane i
znormalizowane przez liczbę punktów i szerokość okna, aby uzyskać
końcową estymację gęstości w każdym punkcie \(x\).
Jądrowy estymator gęstości
W praktyce, aby użyć jądrowego estymatora gęstości wystarczy użyć
funkcji density().
par(mfrow=c(3,1))
plot(density(dane$wiek)) # domyślna wartość okna
plot(density(dane$wiek, bw=0.5)) #mała szerokość okna
plot(density(dane$wiek, bw=10)) #duża szerokość okna
Powyższy przykład ukazuje, jak wykres zmienia się w zależności od
szerokości okna. Jeżeli jest ona mała to na jądrowym estymatorze
gęstości widzimy bardzo dużo szczegółów. Często w praktyce nie jest to
korzystne, ponieważ chcemy wyciągnąć ogólne wnioski i zniwelować
obecność szumu wśród naszych obserwacji. Jeżeli okno jest zbyt duże to
tracimy informacje o kształcie rozkładu naszych danych. Dlatego tak
ważne jest znalezienie kompromisu pomiędzy tymi aspektami.
Jądrowy estymator gęstości często można przedstawiać na tym samym
wykresie co histogram. Pamiętajmy jednak, że w takim przypadku na
histogramie muszą być zaznaczone frakcje, a nie ilości wystąpień.
Użyjemy wcześniejszego przykładu:
par(mfrow=c(2,1))
hist(ciskobiet,breaks = 25, xlim = c(93,178), ylim = c(0,0.07), main="Ciśnienie skurczowe dla kobiet", probability = TRUE)
lines(density(ciskobiet))
hist(cismezczyzna, breaks = 25, xlim = c(93,178), ylim = c(0,0.07), main = "Ciśnienie skurczowe dla meżczyzn", probability = TRUE)
lines(density(cismezczyzna))
Przy użyciu ggplot:
ggplot(dane, aes(x = cisnienie.skurczowe, y = ..density..)) +
geom_histogram(bins = 25, fill = "white", color = "black") +
geom_density(color = "blue", size = 1) +
facet_wrap(~plec, ncol = 1) +
labs(title = "Ciśnienie skurczowe", x = "Ciśnienie skurczowe", y = "Density")
3.1.3 Dystrybuanta
empiryczna
Jest to funkcja schodkowa, która reprezentuje proporcję obserwacji w
próbce, które są mniejsze lub równe \(t\).
To znaczy dla obserwacji \(X_1, X_2, \dots,
X_n\) dystrybuanta empiryczna definiowana jest w następujący
sposób: \[\hat{F}_n(t)=\frac{1}{n}\sum_{i=1}^n\mathbb{I}\{X_i<t\}.\]
Można ją wyznaczyć używając komendy ecdf(). Wynikiem
funkcji ecdf() jest obiekt klasy ecdf, dla
której zaimplementowana jest przeciążona funkcja plot()
prezentująca graficzną reprezentację obliczonej dystrybuanty
empirycznej.
plot(ecdf(ciskobiet), main="Ciśnienie skurczowe dla kobiet")
plot(ecdf(cismezczyzna), main = "Ciśnienie skurczowe dla meżczyzn", add=TRUE,col="red")
Dystrybuantę można wykorzystać do porównania kilka grup obserwacji.
Wystarczy na jednym wykresie przedstawić dystrybuanty dla różnych grup,
wykorzystując argument add=TRUE dla funkcji
plot() (wersja plot przeciążona, w zwykłym
plot nie ma tego argumentu).
Możemy również użyć biblioteki ggplot:
ggplot(dane, aes(x = cisnienie.skurczowe, color = plec)) +
stat_ecdf(geom = "step") +
labs(title = "Ciśnienie skurczowe dla kobiet i meżczyzn", x = "Ciśnienie skurczowe", y = "ECDF") +
scale_color_manual(values = c("kobieta" = "blue", "mezczyzna" = "red"))
3.1.4 Wykres
punktowy
Wykres punktowy, często zwany również wykresem rozrzutu jest jednym
podstawowych wykresów pokazujących zależność pomiędzy dwoma zmiennymi
ilościowymi.
W R najprościej taki wykres jest uzyskać poprzez użycie funkcji
plot(), ale można również zastosować funkcję
scatterplot() z biblioteki car, która posiada
kilka interesujących dodatkowych funkcji:
- prezentacja danych z podziałem na grupy,
reg.line=TRUE - dodanie trendu liniowego,smooth=TRUE- dodanie trendu nieliniowego.
library(car)
scatterplot(dane$cisnienie.rozkurczowe ~ dane$wiek | dane$st.cywilny,
xlab = "Wiek (lata)",
ylab = "Ciśnienie rozkurczowe (mmHg)",
main = "Ciśnienie rozkurczowe względem wieku w zależności od stanu cywilnego",
legend=list(coords="topleft", title="Stan cywilny"))
W formule po znaku |, możemy określić zmienną
warunkującą służąca do podziału na podgrupy. W powyższym przykładzie
pleć jest wyróżniona za pomocą równych kształtów i kolorów.
Przy użyciu ggplot:
ggplot(dane, aes(x = wiek, y = cisnienie.rozkurczowe, color = st.cywilny, shape = st.cywilny)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Scatter Plot of Diastolic Blood Pressure vs Age by Marital Status",
x = "Wiek",
y = "Ciśnienie rozkurczowe",
color = "Stan cywilny",
shape = "Stan cywilny")
Wykresy rozrzutu są przydatne do wykrywania i opisywania zależności
liniowych lub monotonicznych. Na przykład, na podstawie powyższego
wykresu można powiedzieć, że ciśnienie rozkurczowe dla singli z wiekiem
nie zmienia się, a dla osób w związku maleje wraz z wiekiem.
3.1.5 Wykresy
pudełkowe
Dla porównania rozkładów wartości w grupach popularną metodą jest
wykres pudełkowy, zwany też boxplotem lub wykresem ramka-wąsy. Za jego
pomocą można zobrazować położenie, rozmieszczenie i kształt rozkładu
empirycznego badanej cechy.
Jak się tworzy wykres pudełkowy
Elementy wykresu ramka-wąsy.
Pudełko - środkowa część wykresu, rozciąga się od pierwszego
kwartylu (Q1) do trzeciego kwartylu (Q3).
Wewnątrz pudełka jest linia, która wskazuje na medianę.
Wąsy rozciągają się od najmniejszej do maksymalnej wartości, lecz
nie dalej niż \(1.5 \cdot IQR\), gdzie
IQR =Q3-Q1(rozstęp międzykwartylowy). Pozostałe wartości są uznawane za
odstające.
Do uzyskania wykresy pudełkowego służy funkcja
boxplot(). Jako pierwszy argument możemy podać wektor,
listę wektorów, ramkę danych lub formułę ze zmiennymi, które mają się
znaleźć na wykresie.
boxplot(dane$cisnienie.skurczowe,dane$cisnienie.rozkurczowe, horizontal=TRUE, names=c("Skurczowe", "Rozkurczowe"))
lub za pomocą biblioteki ggplot:
dane_long <- data.frame(
Value = c(dane$cisnienie.skurczowe, dane$cisnienie.rozkurczowe),
Type = rep(c("Skurczowe", "Rozkurczowe"), each = nrow(dane))
)
ggplot(dane_long, aes(x = Value, y = Type)) +
geom_boxplot() +
labs(title = "Boxplot of Blood Pressure Types",
x = "Blood Pressure (mmHg)",
y = "")
Na powyższym wykresie można między innymi zauważyć, że ciśnienie
rozkurczowe przyjmuje generalnie wartości z zakresu od 65 do 98 i jest
mniejsze niż ciśnienie skurczowe.
boxplot(dane$wiek~dane$wyksztalcenie, varwidth=TRUE, col="lightblue", ylab="wiek", xlab="Wykształcenie", las=1)
Za pomocą ~ możemy porównać wiek w różnych
grupach wykształcenia. Dodatkowo, argument
varwidth=TRUE powoduje, że szerokość pudełka odpowiadać
będzie liczebności porównywanej grupy.
Za pomocą biblioteki ggplot:
ggplot(dane, aes(x = wyksztalcenie, y = wiek, fill = wyksztalcenie)) +
geom_boxplot(varwidth = TRUE) +
labs(x = "Wykształcenie",
y = "Wiek") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
3.1.6 Zadania - wykresy
dla danych jakościowych
Załaduj zbiór danych diamonds z biblioteki
ggplot. Następnie wykonaj poniższe zadania:
Dla zmiennej price utwórz histogram o 30
przedziałach i zaznacz na nim pionowymi kreskami średnią i medianę.
Nadaj tytuł, podpisz osie i dodaj legendę. Określ cechy otrzymanego
histogramu? Dorysuj na histogramie estymator gęstości
Wykonaj powyższy histogram w ggplocie.
Za pomocą wykresu punktowego zilustruj związek między ceną a
ilością karatów.
Wykonaj wykresy pudełkowe i określ, jak zmienia się cena w
zależności od szlifu (zmienna cut).
Wykonaj wykres dystrybuanty empirycznej dla ceny.
3.1.1 Histogram
Jest jednym z najpopularniejszych statystyk graficznych. Przedstawia rozkład wartości zmiennych ilościowych.
Jak działa histogram?
Histogram składa się z słupków, które ilustrują częstość występowania danych w określonych przedziałach.
Konstrukcja histogramu:
Wybór ilości przedziałów
Dla każdego przedziału liczona jest liczba obserwacji, które wpadają w dany zakres.
Każdy przedział jest reprezentowany przez słupek na wykresie, gdzie wysokość słupka oznacza liczbę obserwacji w danym przedziale.
Idea histogramu: \[f(x)=\frac{1}{n\cdot \text{długość przedziału}}\#\{ \text{obserwacje które wpadają w dany przedział zawierający x}\},\] gdzie \(n\) to ilość obserwacji w naszym zbiorze danych.
Histogram wyznacza rozkład empiryczny
Oznacza to, że histogram pokazuje, jak często poszczególne wartości (lub ich przedziały) występują w badanej próbce. Histogram pozwala zaobserwować kluczowe cechy rozkładu, takie jak modę, skośność, kurtozę, obecność wartości odstających oraz potencjalną multimodalność (więcej o tych pojęciach w trakcie przykładów).
Uwaga W praktyce, gdy liczba danych i liczba przedziałów są odpowiednio duże, histogram może być używany jako estymator gęstości prawdopodobieństwa.
Histogram generujemy za pomocą funkcji hist().
Opcjonalnymi argumentami jest np. breaks, dzięki któremu
możemy określić liczbę lub sposób podziału na przedziały (może być to
liczba lub nazwa algorytmu). Jeżeli nie określimy liczby przedziałów, to
zostanie ona dobrana w zależności od liczby obserwacji. Za pomocą
argumentów freq i probability możemy określić,
czy na histogramie mają bym zaznaczone frakcje, czy liczebności
wystąpień.
Przykład:
hist(dane$cisnienie.skurczowe, main="Histogram jednomodalny", xlab="Ciśnienie skurczowe", ylab="Frakcja wystepowania", probability = TRUE)
W ggplot
ggplot(dane, aes(x=cisnienie.skurczowe,after_stat(density)))+
geom_histogram()+
labs(title="Histogram jednomodalny")Zauważmy, że kształt tego histogramu jest w przybliżeniu symetryczny, ma również jedno maksimum zwane również modą. Jest to przykład histogramu jednomodalnego.
Na poniższym rysunku są zobrazowane inne kształty histogramów.
W przypadku symetrycznych histogramów moda, średnia i mediana pokrywają się (tak jak na panelu A). Jeżeli histogram ma ogon, który wydaje się przeciągnięty w prawo to nazywamy go prawoskośnym (panel B). Dodanie większej ilości wartości w pobliżu maksimum silniej wpływa na średnią niż na modę, czy medianę. Analogiczny przypadek jest, gdy ogon pojawia się po lewej stronie, to znaczy mamy więcej wartości w pobliżu minimum. Taki histogram nazywamy lewoskośnym.
Histogram do porównania zbiorów danych dla różnych grup
Wykres dwóch zestawów danych w postaci dwóch histogramów może być bardzo łatwym pierwszym krokiem do porównania środka, rozrzutu i kształtu każdego z zestawu danych.
library(RVAideMemoire)# Biblioteka do obliczania mody za pomocą funkcji mod().
par(mfrow=c(2,1))
kobieta=dane$plec=="kobieta"
ciskobiet=dane[kobieta,]$cisnienie.skurczowe
mezczyzna=dane$plec=="mezczyzna"
cismezczyzna=dane[mezczyzna,]$cisnienie.skurczowe
hist(ciskobiet, breaks=20, xlim = c(93,178), main="Ciśnienie skurczowe dla kobiet")
abline(v=mean(ciskobiet), col="red")
abline(v=median(ciskobiet),col="blue")
abline(v=mod(ciskobiet), col="green")
hist(cismezczyzna, breaks = 25, xlim = c(93,178), main = "Ciśnienie skurczowe dla mężczyzn")
abline(v=mean(cismezczyzna), col="red")
abline(v=median(cismezczyzna),col="blue")
abline(v=mod(cismezczyzna), col="green")
UWAGA: aby porównywać dwa histogramy ważne jest by
oś x dla oby wykresów była tak samo wyskalowana, można to ustalić za
pomocą argumentu xlim.
Histogram identyfikuje wartości odstające
Na powyższym przykładzie większość
obserwacji wpada w przedział pomiędzy 0 a 200 mm, natomiast jedna
wartość mieści się w przedział pomiędzy 550, a 600. Taką wartość, która
odbiega mocno od pozostałych wartości znajdujących się w zbiorze danych,
nazywamy wartością odstającą.
3.1.2 Jądrowy estymator gęstości
Można o nim myśleć jako o gładkiej wersji histogramu. Jest bardziej zaawansowanym narzędziem do przedstawienia empirycznego rozkładu gęstości dla posiadanych danych.
Zauważmy, że histogram mocno zależy od wyboru ilości przedziałów, ale również od tego, w jakich miejscach dokonujemy dyskretyzacji, co może wpływać na jego kształt i prowadzić do błędnych wniosków.
Jądrowy estymator gęstości rozwiązuje problemy jakie pojawiają się przy histogramach.
Idea jądrowego estymatora gęstości: \[\hat{f}(x)=\frac{1}{n\cdot\text{długość przedziału}}\#\{\text{obserwacje które wpadają w małe otoczenie punktu x}\}.\] Definicja jądrowego estymatora gęstości
Niech \(X_1, \dots, X_n\) będzie próbą danych, wtedy \[\hat{f}(x)=\frac{1}{nh}\sum_{i=1}^n K\left(\frac{x-X_i}{h}\right),\] gdzie:
- \(\hat{f}(x)\) to estymowana gęstość w punkcie \(x\),
- \(h\) to szerokość okna (bandwidth),
- \(K\) to funkcja jądra (kernel function),
nazywamy jądrowym estymatorem gęstości.
Jako funkcję jądra bardzo często przyjmuje się tak zwane jądro Gaussa tj.
\[K(u)=\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}u^2}.\]
Szerokość okna h
Szerokość okna \(h\) kontroluje
stopień wygładzenia estymacji. Małe \(h\) prowadzi do estymacji o wysokiej
wariancji (szczegółowej, “szumowej”), natomiast duże \(h\) prowadzi do estymacji o wysokim
obciążeniu (gładkiej, “rozmytej”). Dla funkcji density()
(patrz przykłady niżej) można ją ustalić ręcznie lub wskazać regułę jej
wyboru, która określi najodpowiedniejszą szerokość, np:
- reguła kciuka -
bw="nrd0"- jest ona używana domyślnie, - reguła Scotta -
bw="nrd", - cross-validation -
bw="ucv", - estymator Sheathera-Jonesa -
bw="SJ".
Intuicyjnie:
Każdy punkt danych \(X_i\) jest zastępowany przez funkcję jądra \(K\), która jest przesunięta tak, aby była wycentrowana na \(X_i\). A następnie te funkcje są sumowane i znormalizowane przez liczbę punktów i szerokość okna, aby uzyskać końcową estymację gęstości w każdym punkcie \(x\).
W praktyce, aby użyć jądrowego estymatora gęstości wystarczy użyć
funkcji density().
par(mfrow=c(3,1))
plot(density(dane$wiek)) # domyślna wartość okna
plot(density(dane$wiek, bw=0.5)) #mała szerokość okna
plot(density(dane$wiek, bw=10)) #duża szerokość okna
Powyższy przykład ukazuje, jak wykres zmienia się w zależności od szerokości okna. Jeżeli jest ona mała to na jądrowym estymatorze gęstości widzimy bardzo dużo szczegółów. Często w praktyce nie jest to korzystne, ponieważ chcemy wyciągnąć ogólne wnioski i zniwelować obecność szumu wśród naszych obserwacji. Jeżeli okno jest zbyt duże to tracimy informacje o kształcie rozkładu naszych danych. Dlatego tak ważne jest znalezienie kompromisu pomiędzy tymi aspektami.
Jądrowy estymator gęstości często można przedstawiać na tym samym wykresie co histogram. Pamiętajmy jednak, że w takim przypadku na histogramie muszą być zaznaczone frakcje, a nie ilości wystąpień.
Użyjemy wcześniejszego przykładu:
par(mfrow=c(2,1))
hist(ciskobiet,breaks = 25, xlim = c(93,178), ylim = c(0,0.07), main="Ciśnienie skurczowe dla kobiet", probability = TRUE)
lines(density(ciskobiet))
hist(cismezczyzna, breaks = 25, xlim = c(93,178), ylim = c(0,0.07), main = "Ciśnienie skurczowe dla meżczyzn", probability = TRUE)
lines(density(cismezczyzna))
Przy użyciu ggplot:
ggplot(dane, aes(x = cisnienie.skurczowe, y = ..density..)) +
geom_histogram(bins = 25, fill = "white", color = "black") +
geom_density(color = "blue", size = 1) +
facet_wrap(~plec, ncol = 1) +
labs(title = "Ciśnienie skurczowe", x = "Ciśnienie skurczowe", y = "Density") 3.1.3 Dystrybuanta empiryczna
Jest to funkcja schodkowa, która reprezentuje proporcję obserwacji w próbce, które są mniejsze lub równe \(t\).
To znaczy dla obserwacji \(X_1, X_2, \dots,
X_n\) dystrybuanta empiryczna definiowana jest w następujący
sposób: \[\hat{F}_n(t)=\frac{1}{n}\sum_{i=1}^n\mathbb{I}\{X_i<t\}.\]
Można ją wyznaczyć używając komendy ecdf(). Wynikiem
funkcji ecdf() jest obiekt klasy ecdf, dla
której zaimplementowana jest przeciążona funkcja plot()
prezentująca graficzną reprezentację obliczonej dystrybuanty
empirycznej.
plot(ecdf(ciskobiet), main="Ciśnienie skurczowe dla kobiet")
plot(ecdf(cismezczyzna), main = "Ciśnienie skurczowe dla meżczyzn", add=TRUE,col="red")
Dystrybuantę można wykorzystać do porównania kilka grup obserwacji.
Wystarczy na jednym wykresie przedstawić dystrybuanty dla różnych grup,
wykorzystując argument add=TRUE dla funkcji
plot() (wersja plot przeciążona, w zwykłym
plot nie ma tego argumentu).
Możemy również użyć biblioteki ggplot:
ggplot(dane, aes(x = cisnienie.skurczowe, color = plec)) +
stat_ecdf(geom = "step") +
labs(title = "Ciśnienie skurczowe dla kobiet i meżczyzn", x = "Ciśnienie skurczowe", y = "ECDF") +
scale_color_manual(values = c("kobieta" = "blue", "mezczyzna" = "red"))3.1.4 Wykres punktowy
Wykres punktowy, często zwany również wykresem rozrzutu jest jednym podstawowych wykresów pokazujących zależność pomiędzy dwoma zmiennymi ilościowymi.
W R najprościej taki wykres jest uzyskać poprzez użycie funkcji
plot(), ale można również zastosować funkcję
scatterplot() z biblioteki car, która posiada
kilka interesujących dodatkowych funkcji:
- prezentacja danych z podziałem na grupy,
reg.line=TRUE- dodanie trendu liniowego,smooth=TRUE- dodanie trendu nieliniowego.
library(car)
scatterplot(dane$cisnienie.rozkurczowe ~ dane$wiek | dane$st.cywilny,
xlab = "Wiek (lata)",
ylab = "Ciśnienie rozkurczowe (mmHg)",
main = "Ciśnienie rozkurczowe względem wieku w zależności od stanu cywilnego",
legend=list(coords="topleft", title="Stan cywilny"))
W formule po znaku |, możemy określić zmienną
warunkującą służąca do podziału na podgrupy. W powyższym przykładzie
pleć jest wyróżniona za pomocą równych kształtów i kolorów.
Przy użyciu ggplot:
ggplot(dane, aes(x = wiek, y = cisnienie.rozkurczowe, color = st.cywilny, shape = st.cywilny)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Scatter Plot of Diastolic Blood Pressure vs Age by Marital Status",
x = "Wiek",
y = "Ciśnienie rozkurczowe",
color = "Stan cywilny",
shape = "Stan cywilny")Wykresy rozrzutu są przydatne do wykrywania i opisywania zależności liniowych lub monotonicznych. Na przykład, na podstawie powyższego wykresu można powiedzieć, że ciśnienie rozkurczowe dla singli z wiekiem nie zmienia się, a dla osób w związku maleje wraz z wiekiem.
3.1.5 Wykresy pudełkowe
Dla porównania rozkładów wartości w grupach popularną metodą jest wykres pudełkowy, zwany też boxplotem lub wykresem ramka-wąsy. Za jego pomocą można zobrazować położenie, rozmieszczenie i kształt rozkładu empirycznego badanej cechy.
Jak się tworzy wykres pudełkowy
Pudełko - środkowa część wykresu, rozciąga się od pierwszego kwartylu (Q1) do trzeciego kwartylu (Q3).
Wewnątrz pudełka jest linia, która wskazuje na medianę.
Wąsy rozciągają się od najmniejszej do maksymalnej wartości, lecz nie dalej niż \(1.5 \cdot IQR\), gdzie IQR =Q3-Q1(rozstęp międzykwartylowy). Pozostałe wartości są uznawane za odstające.
Do uzyskania wykresy pudełkowego służy funkcja
boxplot(). Jako pierwszy argument możemy podać wektor,
listę wektorów, ramkę danych lub formułę ze zmiennymi, które mają się
znaleźć na wykresie.
boxplot(dane$cisnienie.skurczowe,dane$cisnienie.rozkurczowe, horizontal=TRUE, names=c("Skurczowe", "Rozkurczowe"))
lub za pomocą biblioteki ggplot:
dane_long <- data.frame(
Value = c(dane$cisnienie.skurczowe, dane$cisnienie.rozkurczowe),
Type = rep(c("Skurczowe", "Rozkurczowe"), each = nrow(dane))
)
ggplot(dane_long, aes(x = Value, y = Type)) +
geom_boxplot() +
labs(title = "Boxplot of Blood Pressure Types",
x = "Blood Pressure (mmHg)",
y = "")Na powyższym wykresie można między innymi zauważyć, że ciśnienie rozkurczowe przyjmuje generalnie wartości z zakresu od 65 do 98 i jest mniejsze niż ciśnienie skurczowe.
boxplot(dane$wiek~dane$wyksztalcenie, varwidth=TRUE, col="lightblue", ylab="wiek", xlab="Wykształcenie", las=1)
Za pomocą ~ możemy porównać wiek w różnych
grupach wykształcenia. Dodatkowo, argument
varwidth=TRUE powoduje, że szerokość pudełka odpowiadać
będzie liczebności porównywanej grupy.
Za pomocą biblioteki ggplot:
ggplot(dane, aes(x = wyksztalcenie, y = wiek, fill = wyksztalcenie)) +
geom_boxplot(varwidth = TRUE) +
labs(x = "Wykształcenie",
y = "Wiek") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) 3.1.6 Zadania - wykresy dla danych jakościowych
Załaduj zbiór danych
diamondsz bibliotekiggplot. Następnie wykonaj poniższe zadania:Dla zmiennej
priceutwórz histogram o 30 przedziałach i zaznacz na nim pionowymi kreskami średnią i medianę. Nadaj tytuł, podpisz osie i dodaj legendę. Określ cechy otrzymanego histogramu? Dorysuj na histogramie estymator gęstościWykonaj powyższy histogram w ggplocie.
Za pomocą wykresu punktowego zilustruj związek między ceną a ilością karatów.
Wykonaj wykresy pudełkowe i określ, jak zmienia się cena w zależności od szlifu (zmienna
cut).Wykonaj wykres dystrybuanty empirycznej dla ceny.
3.2 Statystyki dla danych ilościowych
Są to liczbowe miary, opisujące podstawowe własności cech rozkładów. Często chcielibyśmy mieć liczbową odpowiedź, gdzie leży “centrum” lub “środek” próby, jak duże jest rozproszenie cechy w próbie wokół owego “centrum” oraz jaki jest kształt rozkładu?
Wskaźniki odnoszące się do “centrum” lub “środka” próby nazywamy wskaźnikami położenia, wskaźniki określające rozproszenie cech nazywamy wskaźnikami rozproszenia. Natomiast wartości określające kształt nazywamy wskaźnikami kształtu.
3.2.1 Wskaźniki
położenia
3.2.1.1 Średnia
Niech \(x_1, \dots, x_n\) będzie
próbą o liczności \(n\). Wtedy \[\bar{x}=\frac{1}{n}\sum_{i=1}^n x_i\]
nazywamy średnią arytmetyczną wartości cechy w próbie.
Aby obliczyć średnią wystarczy użyć funkcji mean().
Przykład: Rozważmy rozkład miesięcznych zasadniczych
wynagrodzeń pracowników z wyższym wykształceniem zatrudnionych w pewnej
firmie:
salary <- c(rep(3500,6), rep(4000,8), rep(4100, 7), rep(4500, 4), rep(5000, 3), rep(6000,2), 13000 )
mean(salary) #średnie wynagrodzenie
## [1] 4506.452
Zauważmy, że histogram jest prawostronnie skośny i posiada wartość
odstającą na prawo. W rezultacie wartość średnia jest znacznie bardziej
przesunięta na prawo od mody. Wyobraźmy sobie sytuację, że młody
absolwent wyższej uczelni zgłasza się na rozmowę kwalifikacyjną do pracy
w tej firmie. Kandydat dowiaduje się że będzie zarabiać 3500 zł, ale że
średnie miesięczne wynagrodzenie w wynosi 4500 zł. Zatem, powiada
wiceprezes firmy, ma pan przed sobą wspaniałe możliwości awansu i
znaczne wyższego wynagrodzenia. Czego kandydat nie słyszy, to tego, że
około 2/3 pracowników tej firmy zarabia mniej niż 4100 zł. Wysoka
średnia jest konsekwencją bardzo wysokich zarobków kierownictwa.
Kandydat dobrze by uczynił, gdyby zapytał się o inny wskaźnik położenia
zwany medianą.
3.2.1.2 Mediana
Medianą jest środkowa wartość (drugi kwartyl) próby uporządkowanej
niemalejąco, od wartości najmniejszej do wartości największej. Niech
\(x_{(1)}, \dots, x_{(n)}\) będzie
uporządkowana niemalejąco próba naszych danych. Wtedy medianę w próbie
nazywamy następującą wielkość: \[x_{\text
{med }}= \begin{cases}x_{((n+1) / 2)}, & \text { gdy } n \text {
nieparzyste } \\ \frac{1}{2}\left(x_{(n / 2)}+x_{(n / 2+1)}\right),
& \text { gdy } n \text { parzyste. }\end{cases}\]
## [1] 4100
Ciąg dalszy przykładu:
Mediana znacznie lepiej oddaje zarobkowe perspektywy kandydata.
Wszystko dlatego, że mediana jest niewrażliwa na wartości
odstające.
Warto jednak zwrócić uwagę, że jeżeli próba jest mała i kolejne
uporządkowane elementy próby są od siebie dość odległe, to mediana jest
niestabilnym wskaźnikiem położenia. W celu uniknięcia niestabilności
można użyć jeszcze innego wskaźnika położenia jakim jest średnia
ucinana.
Ogólniej, można używać statystyki kwartyli, które dla
uporządkowanego zbioru danych wyznaczają miejsca gdzie dane są dzielone
na poszczególne segmenty. Na przykład:
pierwszy kwartyl (Q1), to wartość poniżej której znajduje się 25%
danych,
trzeci kwartyl (Q3), to wartość poniżej której znajduje się 75%
danych.
3.2.1.3 Średnia
ucinana
Powstaje ona przez obliczenie średniej próby powstałej przez
usunięcie z próby oryginalnej \(k\)
wartości najmniejszych i \(k\) wartości
największych:
\[\bar{x}_{tk}=\frac{1}{n-2k}\sum_{i=k+1}^{n-k}
x_{(i)}.\] Aby obliczyć średnia ucinaną w R, wystarczy użyć
argumentu trim do funkcji mean() w którym
podamy frakcje elementów, o które chcemy uciąć naszą średnią.
## [1] 4168
Podamy przykład jeszcze jednego wskaźnika położenia. Którym warto
posłużyć się na przykład wtedy, gdy skrajne obserwacje charakteryzują
się dużą niepewnością co do ich rzeczywistych wartości.
3.2.1.4 Średnia
Winsorowska
Podobnie jak w przypadku średniej ucinanej, obliczanie średniej
winsorowskiej opiera się na \(n-2k\)
środkowych elementach próby, ale uwzględniając fakt pojawienia się w
próbie \(k\) wartości większych niż
\(x_{(k+1)}\) oraz \(k\) wartości mniejszych niż \(x_{(n-k)}\) przy obliczeniu średniej tak
jakby \(x_{(k+1)}\) oraz \(x_{(n-k)}\) wystąpiły dodatkowo \(k\) razy.
Formalnie, można wyrazić ją za pomocą poniższego wzoru: \[\bar{x}_{w k}=\frac{1}{n}\left[(k+1) \cdot
x_{(k+1)}+\sum_{i=k+2}^{n-k-1} x_{(i)}+(k+1)\cdot x_{(n-k)}\right]
.\] W R należy użyć funkcji winsor.mean() z pakietu
psych.
## [1] 4183.871
3.2.2 Wskaźniki
rozproszenia
Zredukowanie informacji o rozkładzie cech w próbie do wartości
średniej jest zdecydowanie zbyt drastyczne. Potrzebujemy jeszcze
informacji o stopniu rozproszenia próby wokół średniej.
3.2.2.1 Rozstęp
Rozstęp to różnica pomiędzy największą, a najmniejsza
wartością w próbie.
## [1] 9500
3.2.2.2 Wariancja (drugi
moment centralny)
Gdy wariancja jest mała, to oznacza to, że obserwacje są bardziej
skupione wokół centrum. Natomiast, gdy wariancja jest duża to rozstęp
jest duży i obserwacje znajdują się szeroko wokół swojego centrum.
Wariancję (var()) w próbie nazwiemy następującą
wielkość: \[s^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})^2.\]
Pierwiastek z wariancji nazywamy odchyleniem standardowym
(funkcja w R to sd()).
Warto zauważyć, że wariancja jest wskaźnikiem, który nie jest odporny
na wartości odstające.
## [1] 60.32462
## [1] 246.904
par(mfrow=c(2,1))
hist(dane$cisnienie.rozkurczowe, xlim = c(57,178))
hist(dane$cisnienie.skurczowe, xlim = c(57,178))
Na podstawie histogramów również jesteśmy w stanie ocenić, czy dane
mają dużą czy małą wariancję.
3.2.3 Wskaźniki
kształtu
Informują nasz jaki jest kształt rozkładu naszych danych
3.2.3.1 Skośność (Trzeci
moment centralny)
Jest miarą asymetrii. Wyraża się poniższym wzorem \[\gamma=
\frac{n\sum_{i=1}^n(x_i-\bar{x})^3}{(n-1)(n-2)s^3}.\] Może
przyjmować wartości dodatnie, ujemne i zerowe mając wpływ na to jak
wygląda rozkład:
- gdy \(\gamma=0\), to rozkład jest
symetryczny,
- gdy \(\gamma<0\), to rozkład
jest lewo skośny,
- gdy \(\gamma>0\), to rozkład
jest prawo skośny.
Wpływ skośności na kształt rozkładu
Aby obliczyć skośność w R wystarczy użyć funkcji
skewness().
3.2.3.2 Kurtoza
Wyraża się poniższym wzorem:
\[\text{kurtoza}=
\frac{n(n+1)\sum_{i=1}^n(x_i-\bar{x})^4}{(n-1)(n-2)(n-3)s^4}-3\frac{(n-1)^2}{(n-2)(n-3)}.\]
Również może przyjmować wartości dodatnie, ujemne i zerowe. W zależności
od tego wpływa na kształt ogonów rozkładu.
Wpływ kurtozy na rozkład
Jeżeli kurtoza jest ujemna, to wtedy ogony rozkładu stają się
ciężkie, coraz więcej obserwacji przyjmuje wartości skrajne.
Dla rozkładu normalnego kurtoza=3. Dlatego też często stosuje się
skorygowaną kurtozę=kurtoza-3, która dla rozkładu normalnego
wynosi 0. Następnie względem niego porównujemy ciężar ogonów
rozkładu.
Kurtozę można obliczyć za pomocą funkcji kurtosis z
pakietu moments.
3.2.1 Wskaźniki położenia
3.2.1.1 Średnia
Niech \(x_1, \dots, x_n\) będzie próbą o liczności \(n\). Wtedy \[\bar{x}=\frac{1}{n}\sum_{i=1}^n x_i\] nazywamy średnią arytmetyczną wartości cechy w próbie.
Aby obliczyć średnią wystarczy użyć funkcji mean().
Przykład: Rozważmy rozkład miesięcznych zasadniczych wynagrodzeń pracowników z wyższym wykształceniem zatrudnionych w pewnej firmie:
salary <- c(rep(3500,6), rep(4000,8), rep(4100, 7), rep(4500, 4), rep(5000, 3), rep(6000,2), 13000 )
mean(salary) #średnie wynagrodzenie## [1] 4506.452
Zauważmy, że histogram jest prawostronnie skośny i posiada wartość odstającą na prawo. W rezultacie wartość średnia jest znacznie bardziej przesunięta na prawo od mody. Wyobraźmy sobie sytuację, że młody absolwent wyższej uczelni zgłasza się na rozmowę kwalifikacyjną do pracy w tej firmie. Kandydat dowiaduje się że będzie zarabiać 3500 zł, ale że średnie miesięczne wynagrodzenie w wynosi 4500 zł. Zatem, powiada wiceprezes firmy, ma pan przed sobą wspaniałe możliwości awansu i znaczne wyższego wynagrodzenia. Czego kandydat nie słyszy, to tego, że około 2/3 pracowników tej firmy zarabia mniej niż 4100 zł. Wysoka średnia jest konsekwencją bardzo wysokich zarobków kierownictwa. Kandydat dobrze by uczynił, gdyby zapytał się o inny wskaźnik położenia zwany medianą.
3.2.1.2 Mediana
Medianą jest środkowa wartość (drugi kwartyl) próby uporządkowanej niemalejąco, od wartości najmniejszej do wartości największej. Niech \(x_{(1)}, \dots, x_{(n)}\) będzie uporządkowana niemalejąco próba naszych danych. Wtedy medianę w próbie nazywamy następującą wielkość: \[x_{\text {med }}= \begin{cases}x_{((n+1) / 2)}, & \text { gdy } n \text { nieparzyste } \\ \frac{1}{2}\left(x_{(n / 2)}+x_{(n / 2+1)}\right), & \text { gdy } n \text { parzyste. }\end{cases}\]
## [1] 4100Ciąg dalszy przykładu:
Mediana znacznie lepiej oddaje zarobkowe perspektywy kandydata.
Wszystko dlatego, że mediana jest niewrażliwa na wartości odstające.
Warto jednak zwrócić uwagę, że jeżeli próba jest mała i kolejne uporządkowane elementy próby są od siebie dość odległe, to mediana jest niestabilnym wskaźnikiem położenia. W celu uniknięcia niestabilności można użyć jeszcze innego wskaźnika położenia jakim jest średnia ucinana.
Ogólniej, można używać statystyki kwartyli, które dla uporządkowanego zbioru danych wyznaczają miejsca gdzie dane są dzielone na poszczególne segmenty. Na przykład:
pierwszy kwartyl (Q1), to wartość poniżej której znajduje się 25% danych,
trzeci kwartyl (Q3), to wartość poniżej której znajduje się 75% danych.
3.2.1.3 Średnia ucinana
Powstaje ona przez obliczenie średniej próby powstałej przez usunięcie z próby oryginalnej \(k\) wartości najmniejszych i \(k\) wartości największych:
\[\bar{x}_{tk}=\frac{1}{n-2k}\sum_{i=k+1}^{n-k}
x_{(i)}.\] Aby obliczyć średnia ucinaną w R, wystarczy użyć
argumentu trim do funkcji mean() w którym
podamy frakcje elementów, o które chcemy uciąć naszą średnią.
## [1] 4168Podamy przykład jeszcze jednego wskaźnika położenia. Którym warto posłużyć się na przykład wtedy, gdy skrajne obserwacje charakteryzują się dużą niepewnością co do ich rzeczywistych wartości.
3.2.1.4 Średnia Winsorowska
Podobnie jak w przypadku średniej ucinanej, obliczanie średniej winsorowskiej opiera się na \(n-2k\) środkowych elementach próby, ale uwzględniając fakt pojawienia się w próbie \(k\) wartości większych niż \(x_{(k+1)}\) oraz \(k\) wartości mniejszych niż \(x_{(n-k)}\) przy obliczeniu średniej tak jakby \(x_{(k+1)}\) oraz \(x_{(n-k)}\) wystąpiły dodatkowo \(k\) razy.
Formalnie, można wyrazić ją za pomocą poniższego wzoru: \[\bar{x}_{w k}=\frac{1}{n}\left[(k+1) \cdot
x_{(k+1)}+\sum_{i=k+2}^{n-k-1} x_{(i)}+(k+1)\cdot x_{(n-k)}\right]
.\] W R należy użyć funkcji winsor.mean() z pakietu
psych.
## [1] 4183.8713.2.2 Wskaźniki rozproszenia
Zredukowanie informacji o rozkładzie cech w próbie do wartości średniej jest zdecydowanie zbyt drastyczne. Potrzebujemy jeszcze informacji o stopniu rozproszenia próby wokół średniej.
3.2.2.1 Rozstęp
Rozstęp to różnica pomiędzy największą, a najmniejsza wartością w próbie.
## [1] 95003.2.2.2 Wariancja (drugi moment centralny)
Gdy wariancja jest mała, to oznacza to, że obserwacje są bardziej skupione wokół centrum. Natomiast, gdy wariancja jest duża to rozstęp jest duży i obserwacje znajdują się szeroko wokół swojego centrum.
Wariancję (var()) w próbie nazwiemy następującą
wielkość: \[s^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})^2.\]
Pierwiastek z wariancji nazywamy odchyleniem standardowym
(funkcja w R to sd()).
Warto zauważyć, że wariancja jest wskaźnikiem, który nie jest odporny na wartości odstające.
## [1] 60.32462## [1] 246.904par(mfrow=c(2,1))
hist(dane$cisnienie.rozkurczowe, xlim = c(57,178))
hist(dane$cisnienie.skurczowe, xlim = c(57,178))
Na podstawie histogramów również jesteśmy w stanie ocenić, czy dane mają dużą czy małą wariancję.
3.2.3 Wskaźniki kształtu
Informują nasz jaki jest kształt rozkładu naszych danych
3.2.3.1 Skośność (Trzeci moment centralny)
Jest miarą asymetrii. Wyraża się poniższym wzorem \[\gamma= \frac{n\sum_{i=1}^n(x_i-\bar{x})^3}{(n-1)(n-2)s^3}.\] Może przyjmować wartości dodatnie, ujemne i zerowe mając wpływ na to jak wygląda rozkład:
- gdy \(\gamma=0\), to rozkład jest symetryczny,
- gdy \(\gamma<0\), to rozkład jest lewo skośny,
- gdy \(\gamma>0\), to rozkład jest prawo skośny.
Aby obliczyć skośność w R wystarczy użyć funkcji
skewness().
3.2.3.2 Kurtoza
Wyraża się poniższym wzorem:
\[\text{kurtoza}= \frac{n(n+1)\sum_{i=1}^n(x_i-\bar{x})^4}{(n-1)(n-2)(n-3)s^4}-3\frac{(n-1)^2}{(n-2)(n-3)}.\] Również może przyjmować wartości dodatnie, ujemne i zerowe. W zależności od tego wpływa na kształt ogonów rozkładu.
Jeżeli kurtoza jest ujemna, to wtedy ogony rozkładu stają się ciężkie, coraz więcej obserwacji przyjmuje wartości skrajne.
Dla rozkładu normalnego kurtoza=3. Dlatego też często stosuje się skorygowaną kurtozę=kurtoza-3, która dla rozkładu normalnego wynosi 0. Następnie względem niego porównujemy ciężar ogonów rozkładu.
Kurtozę można obliczyć za pomocą funkcji kurtosis z
pakietu moments.
4 Wstępna analiza danych
W większości przypadków zebrane dane wymagają wstępnego przetwarzania zanim zostanie na nich wykonana analiza statystyczna.
4.1 Kwestia brakujących obserwacji
Dość powszechnym zjawiskiem jest brak co najmniej jednej wartości charakteryzującej analizowane próbki. Przyczyną może być np. jakiś błąd w procesie gromadzenia informacji, niemożliwość wykorzystania pewnych form pomiarów lub podczas wypełniania formularza pewne pola mogły pozostać puste.
W R brakujące dane oznaczane są symbolem NA.
Zignorowanie tych wartości może prowadzić do trudności z obliczeniem
podstawowych statystyk lub do błędnych wnioskowań.
Wczytajmy zbiór danych airquality z brakującymi
obserwacjami.
## [1] 153 6W celu sprawdzenia czy w zbiorze danych są i ile jest brakujących
obserwacji możemy użyć funkcji sum(is.na()):
## [1] 44Zatem, posiadamy 44 brakujące obserwacje. Możemy również użyć funkcji
summary(), która również może służyć do wykrycia
brakujących wartości:
## Ozone Solar.R Wind Temp
## Min. : 1.00 Min. : 7.0 Min. : 1.700 Min. :56.00
## 1st Qu.: 18.00 1st Qu.:115.8 1st Qu.: 7.400 1st Qu.:72.00
## Median : 31.50 Median :205.0 Median : 9.700 Median :79.00
## Mean : 42.13 Mean :185.9 Mean : 9.958 Mean :77.88
## 3rd Qu.: 63.25 3rd Qu.:258.8 3rd Qu.:11.500 3rd Qu.:85.00
## Max. :168.00 Max. :334.0 Max. :20.700 Max. :97.00
## NA's :37 NA's :7
## Month Day
## Min. :5.000 Min. : 1.0
## 1st Qu.:6.000 1st Qu.: 8.0
## Median :7.000 Median :16.0
## Mean :6.993 Mean :15.8
## 3rd Qu.:8.000 3rd Qu.:23.0
## Max. :9.000 Max. :31.0
## Jednym z najprostszych sposobów jest usunięcie całych wierszy (lub
kolumn) ze zbioru danych. Można to zrobić za pomocą funkcji
na.omit()
## [1] 111 6Pomysł ten jest zdecydowanie nie najlepszym rozwiązaniem. Usunięcie
kilkudziesięciu procent zebranych wyników, ponieważ brakuje jednego
pomiaru, prowadzi do niepotrzebnego pogorszenia właściwości procedur
statystycznych, zwiększenia wariancji estymatorów, zmniejszenia mocy
testów, itp. Zatem usunięcie wierszy w których pojawiają się wartości
NA ma sens, gdy jest ich mało i pojawiają się one losowo
(Missing at Random). Zamiast usuwać całe wiersze, można w ich miejsce
wstawić sztucznie inny pomiar. Można to zrobić na dwa sposoby:
W miejsce brakującej wartości wstawia się wartość charakterystyczną dla danej zmiennej. Najczęściej medianę, średnią, ewentualnie modę dla danej zmiennej. Bardziej wyrafinowanym sposobem jest dopasowanie modelu regresji i wstawienie wartości wyznaczonej z modelu regresji (o tym w przyszłości).
Wygenerowanie wielu replikacji (powtórzeń) zbioru danych z powtórzeniami. Dla każdego powtórzenia brakujące obserwacje zastępowane są przez obserwacje wygenerowane z rozkładu określonego dla danej zmiennej. Kolejne kroki analizy są wykonywane dla każdej z replikacji zbioru danych, a następnie wyniki poszczególnych kroków są ze sobą zestawiane (więcej na temat tej metody w kolejnych laboratoriach).
Uwaga: W przypadku, gdy analizujemy zbiór danych najlepiej dowiedzieć się o nim jak najwięcej, ponieważ w przypadku, gdy brakujące wartości powstają w inny sposób na przykład dlatego, że wychodzą poza skalę czujnika pomiaru, to wstawianie tam średniej/mediany jest błędne.
Teraz skupimy się na pierwszej metodzie. W celu podstawienia za
brakujące wartości innych wygodnie jest posłużyć się funkcją
inpute() z pakiety Hmisc. Za pomocą argumentu
fun= możemy wybrać, czy podstawiamy średnią
(mean), czy medianę (jest ona domyślnym parametrem tej
funkcji). Jeżeli fun=random, to brakujące wartości będą
zastąpione przez pozostałe wartości z próby, wybrane metodą losowania ze
zwracaniem.
4.2 Zadania - brakujące dane
Załaduj zbiór danych mieszkania z pakietu
Przewodnik. Następnie dla kolumny cena dla indeksów
c(1,2,3,50, 60, 80, 100,105, 110, 120, 180, 181,188, 190, 195,196, 197, 198, 199, 200)
wpisz wartość NA. Następnie, usuń i uzupełnij brakujące
dane poznanymi metodami. Porównaj otrzymany zestaw danych do
oryginalnego. Między innymi oblicz błąd tj. wartość bezwzględną różnicy
średniej i mediany przed usunięciem danych i po zastosowaniu metody
radzenia sobie z wartościami NA. Dla której metody błąd był
największy, a dla której najmniejszy?
5 Skalowanie
Skalowanie zmiennych jest ważnym krokiem w przygotowaniu danych do analizy, ale nie zawsze jest konieczne.
Skalujemy zmienne gdy chcemy użyć ich:
do algorytmów na odległość (\(k\)-means, KNN, SVD),
regresja (liniowa, logistyczna, grzbietowa, LASSO),
analiza składowych głównych (PCA),
sieci neuronowe.
Nie skalujemy zmiennych: gdy chcemy użyć ich:
drzew losowych, lasów losowych, Gradient Boosting Machines (GBM),
algorytmów, które działają na bazie rang lub klasyfikacji.
Pojawiło się tutaj kilka nowych skrótów od nazw algorytmów, które w przyszłość może większość z Was pozna (np. na II stopniu na Statistical Learning, Machine Learning, Wstęp do analizy danych, Deep Learning). Na razie jest to informacja tylko poglądowa. Niestety, nie damy rady omówić ich na tym przedmiocie.
5.1 Normalizacja
Zakresy zmiennych na ogół się różnią (np. z powodu jednostek/zakresów), co może wpływać na to, że niektóre zmienne w większym zakresie będą miały wpływ na wyniki.
Normalizacja min-max skaluje wartości do przedziału \([0,1]\) według poniższego wzoru: \[x_{\text{norm}}^{(i)}=\frac{x^{(i)}-x_{\min}}{x_{\max}-x_{\min}},\] gdzie \(x^{(i)}\) oznacza daną próbkę, \(x_{\max}\), \(x_{\min}\) to odpowiednio maksimum i minimum dla wartości próbki w danej kolumnie.
Aby znormalizować dane w poniższym przykładzie użyjemy funkcji
lapply():
df <- data.frame("Wiek"=c(40, 13, 15, 45),
"Wzrost[cm]"=c(168, 146, 170, 165),
"Waga[kg]"=c(85,27, 42, 67))
# Normalizacja
normalize <- function(x) {
return((x - min(x)) / (max(x) - min(x)))
}
df_normalized <- as.data.frame(lapply(df, normalize))
df_normalized## Wiek Wzrost.cm. Waga.kg.
## 1 0.84375 0.9166667 1.0000000
## 2 0.00000 0.0000000 0.0000000
## 3 0.06250 1.0000000 0.2586207
## 4 1.00000 0.7916667 0.6896552Normalizacja utrudnia interpretację danych, ponieważ są one oderwane od jednostek i skali. W szczególności najmłodsza osoba teraz ma wiek 0, a najstarsza 1.
Normalizacja jest bardzo wrażliwa na wartości odstające, jeżeli chociaż jedna obserwacja znacznie odbiega od pozostałych, to wpłynie to na proces normalizacji.
5.2 Standaryzacja
Standaryzacja skaluje wartości tak, żeby miały średnią \(0\) i odchylenie standardowe \(1\).
Procedurę standaryzacji można zapisać za pomocą poniższego wzoru: \[x_{\text{std}}^{(i)}=\frac{x^{(i)}-\bar{x}^{(i)}}{s^{(i)}},\] gdzie \(\bar{x}^{(i)}\) jest średnia próbek z kolumny \(i\), a \(s^{(i)}\) jest związanym z nią odchyleniem standardowym.
Aby zestandaryzować dane możemy użyć funkcji
scale():
## Wiek Wzrost.cm. Waga.kg.
## 1 0.7078014 0.5214718 1.1531499
## 2 -0.9186358 -1.4737246 -1.0950079
## 3 -0.7981590 0.7028533 -0.5135878
## 4 1.0089934 0.2493996 0.4554457Standaryzacja zachowuje informacje o wartościach odstających, ale nie jest na nie tak podatna. Pokazuje również jak daleko w jakim kierunku od średniej (tj. \(0\)) są obserwacje.
5.3 Zadania
Zadanie 1 Za pomocą funkcji cov()
wykonaj macierz kowariancji dla zbioru danych daneSoc z
pakietu Przewodnik dla danych nie skalowanych i po
standaryzacji. Do macierz kowariancji weź wszystkie zmienne ilościowe.
Przypomnij co to jest macierz kowariancji i jakie ma własności. Porównaj
dwie otrzymacie macierze, czym się różnią i jakie wnioski na ich
podstawie możny wyciągnąć?
6 Zadanie - wstępna analiza danych
W zależności od ostatniej cyfry indeksu wybierz odpowiedni zbiór danych. Zbiory danych znajdują się na stronie przedmiotu. Pod linkami niżej znajdują się opisy zawartości zbiorów danych:
0. Australian Health Service Utilization Data - dane0.csv - dane 0;
1. House Prices in the city of Windsor, Canada - dane1.csv - dane 1;
2. Doctoral Publications - dane2. csv - dane 2;
3. Travel Mode Choice Data - dane3.csv - dane 3;
4. Economic Journals Data Set - dane4.csv - dane 4;
5. Wife Working Hours - dane5.csv - dane 5;
6. Epilepsy Data - dane6.csv - dane 6;
7. Days not Spent at School - dane7.csv - dane 7;
8. Porsche and Jaguar Prices - dane8.csv - dane 8;
9. Overdrawn Checking Account - dane9.csv - dane 9.
Załaduj zbiór danych i wykonaj następujące kroki:
a. Opisz zbiór danych, jeśli jest taka informacja to opisz sposób, w jaki dane zostały pozyskane.
b. Określ które zmienne są ilościowe, a które są jakościowe. Sprawdź, czy są wartości brakujące, jeśli tak to użyj jednej z podanych metod uzupełniania danych.
Uwaga: często zmienne jakościowe są przez język R
rozpoznane jako charakter. Aby R interpretował je jako
zmienne jakościowe musimy wykonać as.factor(zmiennej) i
przypisać ją do “starej” zmiennej.
c. Wykonaj przykładowe wykresy i tabele liczności dla wybranych zmiennych jakościowych i opisz co na ich podstawie możemy powiedzieć o posiadanym zbiorze danych. Nie trzeba wykonywać wykresów dla wszystkich zmiennych jakościowych, możemy wybrać takie, aby móc wyciągnąć jakieś logiczne wnioski.
d. Wykonaj wykresy poznane w tym laboratorium dla wybranej zmiennej ilościowej. Oblicz wskaźniki położenia (średnia, mediana, moda), rozrzutu (odchylenie standardowe, wariancja), skośności i kurtozy. Porównaj obliczone wskaźniki z wyglądem histogramu, opisując, jak wskaźniki te wpływają na kształt i interpretację wykresu.