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 tej bibliotece
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 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 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"))
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" ))
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") Domyślnie funkcja pie() przy macierzach nieliczność nie
podpisuje wycinków koła. Należy nazwy przekazać do funkcji za pomocą
argumentu labels.
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
2.3.1 Zadanie - zmienne jakościowe + wykresy
set.seed(111111) # W miejsce "jedynek" wpisz swój numer indeksu
dane_j <- data.frame(
Płeć = as.factor(sample(c("Mężczyzna", "Kobieta"), 100, replace = TRUE)),
Rok = as.numeric(sample(c("1", "2", "3", "4"), 100, replace = TRUE)),
Matematyka = as.factor(sample(c("Lubi M", "Nie lubi M"), 100, replace = TRUE)),
Programowanie = as.factor(sample(c("Lubi Pr", "Nie lubi Pr"), 100, replace = TRUE)),
Sztuka = as.factor(sample(c("Lubi Sz", "Nie lubi Sz"), 100, replace = TRUE))
)
# Podgląd danych
head(dane_j)
summary(dane_j)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 (podstawowa grafika).
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 (podstawowa grafika).
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.
W wyniku powinniśmy dostać coś podobnego do:
##
## Lubi Pr Nie lubi Pr
## 60 40##
## Lubi Pr Nie lubi Pr
## 0.6 0.4##
## Lubi Sz Nie lubi Sz
## 44 56##
## Lubi Sz Nie lubi Sz
## 0.44 0.56
##
## Lubi Pr Nie lubi Pr
## Kobieta 25 18
## Mężczyzna 35 22##
## Lubi Pr Nie lubi Pr
## Kobieta 0.25 0.18
## Mężczyzna 0.35 0.22
2.3.2 Zadanie -wykresy w
ggplot2
Wykonaj wykresy z powyższego zadania w wykorzystaniem biblioteki
ggplot2.
W wyniku powinniśmy dostać (mniej więcej) wykresy: 
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.
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 Scota -
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\). 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:
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:
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
stan cywilny jest wyróżniony za pomocą różnych 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:
3.1.6 Zadanie - wykresy dla danych jakościowych
Załaduj zbiór danych diamonds z biblioteki
ggplot. 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.1.6.1 Wskazówki
Co powinniśmy otrzymać (wykorzystanie różnych bibliotek nie jest konieczne)
## carat cut color clarity depth
## Min. :0.2000 Fair : 1610 D: 6775 SI1 :13065 Min. :43.00
## 1st Qu.:0.4000 Good : 4906 E: 9797 VS2 :12258 1st Qu.:61.00
## Median :0.7000 Very Good:12082 F: 9542 SI2 : 9194 Median :61.80
## Mean :0.7979 Premium :13791 G:11292 VS1 : 8171 Mean :61.75
## 3rd Qu.:1.0400 Ideal :21551 H: 8304 VVS2 : 5066 3rd Qu.:62.50
## Max. :5.0100 I: 5422 VVS1 : 3655 Max. :79.00
## J: 2808 (Other): 2531
## table price x y
## Min. :43.00 Min. : 326 Min. : 0.000 Min. : 0.000
## 1st Qu.:56.00 1st Qu.: 950 1st Qu.: 4.710 1st Qu.: 4.720
## Median :57.00 Median : 2401 Median : 5.700 Median : 5.710
## Mean :57.46 Mean : 3933 Mean : 5.731 Mean : 5.735
## 3rd Qu.:59.00 3rd Qu.: 5324 3rd Qu.: 6.540 3rd Qu.: 6.540
## Max. :95.00 Max. :18823 Max. :10.740 Max. :58.900
##
## z
## Min. : 0.000
## 1st Qu.: 2.910
## Median : 3.530
## Mean : 3.539
## 3rd Qu.: 4.040
## Max. :31.800
## 
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] 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 labolatoriach).
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ą
impute() 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 Zadanie - 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 Zadanie - macierz kowariancji i standaryzacja
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 (tutaj linki do opisów, same dane na stronie):
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.