####################################################################################### #################### Cvičení 4 - Základy statistiky ################################### ########## Martina Litschmannová, Adéla Vrtková, Lenka Přibylová ###################### ################# Katedra aplikované matematiky, FEI, VŠB-TUO ######################### ####################################################################################### # Zobrazuje-li se vám skript s chybným kódováním, využijte příkaz File / Reopen with Encoding # Náplň je inspirována online kurzy na https://www.datacamp.com/ ## Obsah skriptu # 1. Popisná statistika kvantitavní proměnné # 1.1 Preprocesing # 1.2 Výpočet a interpretace základních charakteristik (základní R) # 1.3 Práce s chybějícími údaji (NA) # 1.4 Zaokrouhlování základních charakteristik # 1.5 Příprava tabulky základních charakteristik (dplyr) # 2. Vizualizace kvantitavní proměnné # 2.1 Krabicový graf (základní R) # 2.2 Krabicový graf (ggplot2) # 2.3 Histogram (základní R) # 2.4 Histogram (ggplot2) # 3. Identifikace a ev. odstranění odlehlých pozorování # Instalace a načtení balíčků # install.packages("readxl") # pro načtení xlsx souborů # install.packages("dplyr") # pro efektivnější manipulaci s datovým souborem # install.packages("moments") # pro výpočet šikmosti a špičatosti # install.packages("ggplot2") # pro hezčí vizualizace library(readxl) library(dplyr) library(moments) library(ggplot2) # Úkol: ## Analyzujme měsíční životní náklady respondentů našeho šetření (proměnná naklady). ############################################################################################### ## 1. Popisná statistika kvantitativní proměnné ############################################### #///////////////////////////////////////////////////////////////////////////////////// # 1.1 Preprocesing # Veškerá analýza je prezentována na výstupech z dotazníkového šetření určeného # studentům Základů statistiky - viz http://am-nas.vsb.cz/lit40/ZS/dotaznik.pdf. # Importujeme do R data uložena v souboru dotaznik.xlsx - viz http://am-nas.vsb.cz/lit40/ZS/dotaznik.xlsx # Datová soubor si uložte do pracovního adresáře a před vlastním exportem do R # proveďte potřebný preprocessing, popřípadě využijte data, v nichž již první úpravy # byly provedeny - viz http://am-nas.vsb.cz/lit40/ZS/Excel/dotaznik_pp.xlsx. # Nastavení pracovního adresáře getwd() setwd("C:/Users/Martina/Výuka/Základy statistiky") # V upravte odkaz na vlastní pracovní adresář dotaznik = read_excel("dotaznik_pp.xlsx")# Název souboru případně upravte. # Zjednodušme pojmenování proměnných. colnames(dotaznik)=c("ID","prace","koleje","naklady","prijmy","plat","pocitac_studium","soc_site", "os","leg_software","ikb_znalost","soc_inz","soc_inz_techniky","sexting","darknet", "hesla","piti","sport","pocitac_hry","student_typ","hry","odpocivarna","serialy", "prokrastinace","hobby","bez_stresu","bez_problemu") # Analýzu budeme provádět pouze pro respondenty, u nichž máme úplné záznamy s výjimkou odpovědí na otázku # „S jakými technikami sociálního inženýrství jste se setkal(a)? (soc_inz_techniky)". # V minulém cvičení jsme se naučili, jak tato data vyfiltrovat a uložit do datového rámce dotaznik2: dotaznik2 = dotaznik %>% filter(complete.cases(select(dotaznik,-soc_inz_techniky))) # Pozn.: V celém skriptu budeme pracovat s datovým rámcem dotaznik2. # Pouze v úvodu do problematiky bude na práci s datovým rámcem dotaznik # prezentována práce s proměnnou obsahující chybějící údaje (dále NA). # Nejdříve ověřme klasifikaci (typ) analyzované proměnné. class(dotaznik2$naklady) # Pokud proměnná není klasifikována jako typ numeric, musíme ji překódovat. # K nesprávné klasifikaci proměnné dochází zejména v případech chybně nastaveného oddělovače des. míst, # překlepech v hodnotách proměnné (hodnota je vyhodnocena jako textový řetězec), ... dotaznik2$naklady=as.numeric(dotaznik2$naklady) # --------------------------------------------------------------------------------------------- #///////////////////////////////////////////////////////////////////////////////////// # 1.2 Výpočet základních charakteristik (základní R) # Nejdříve nastavíme požadovaný formát vypočtených číselných charakteristik options(scipen = 100, # 100 (desetinné číslo), -100 (vědecká notace) digits = 10) # max. počet platných cifer # Začněme výpočtem základních měr polohy min(dotaznik2$naklady) # minimum quantile(dotaznik2$naklady,probs = 0.25) # dolní kvartil quantile(dotaznik2$naklady,probs = 0.5) # medián mean(dotaznik2$naklady) # průměr quantile(dotaznik2$naklady,probs = 0.75) # horní kvartil max(dotaznik2$naklady) # maximum # Základní míry polohy lze získat rovněž pomocí funkce summary. summary(dotaznik2$naklady) # Úkol A ------------------------------------------------------------------------------ # Zjištěné výsledky komentujte a rozhodněte se o dalším přístupu k analýze. # --- # Zamysleme se nad dalším postupem. Jak si vysvětlit průměrné mesíční náklady studenta # ve výši 1 Kč? A co 400 Kč? Jak se k těmto zjištěním dále postavit? # Nemůže mít na uváděnou výši nákladů vliv nepřesně formulována otázka? # Zjistěme kolik studentů uvádí průměrné měsíční náklady nižší než # 1 000 Kč, tj. méně než 35 Kč/den. dotaznik2 %>% filter(naklady<1000) # 9 z 93 respondentů uvádí, že jejich průměrné měsíční náklady jsou nižší než 1 000 Kč. # Je zřejmé, že naše data k této otázce nejsou dostatečně relevantní a nelze je # smysluplně vyhodnotit. ## Analyzujme představy respondentů našeho šetření ## o měsíčním nástupním platu (čistého, proměnná plat). class(dotaznik2$plat) # kontrola klasifikace proměnné plat summary(dotaznik2$plat) length(dotaznik2$plat) # počet respondentů # Úkol B ------------------------------------------------------------------------------ # Zjištěné výsledky komentujte. # -------------------------------------------------------------------------------------- # Zjistěme kolik studentů uvedlo očekávaný měsíční plat vyšší než # 100 000 Kč. (Částka 100 000 Kč byla stanovena pouze expertním odhadem, # tímto se snažíme pouze o první náhled na pozorovaná data.) dotaznik2 %>% filter(plat>100000) # Pouze jeden student (z 93) uvedl očekávaný měsíční plat vyšší než 100 000 Kč. # Na první pohled není v datech vidět výraznější "problémy", pokusíme se # v analýze pokračovat... # V tuto chvíli bychom se v praxi určitě zabývali tím, zda v analýze ponechat záznam # studenta, který očekává měsíční plat po ukončení studia ve výši 500 000,- Kč. # My se k tomuto problému vrátíme až na závěr skriptu v rámci práce # s odlehlými hodnotami. # Určeme základní míry variability očekávaných platů. var(dotaznik2$plat) # rozptyl (Připomeňme, že jednotka je kvadrátem jednotky měřené veličiny!) sd(dotaznik2$plat) # směrodatná odchylka 100*sd(dotaznik2$plat)/mean(dotaznik2$plat) # variační koeficient (%) # Úkol C ------------------------------------------------------------------------------ # Komentujte variabilitu očekávaných platů. # Využijte Čebyševovu nerovnost a variační koeficient. # Zamysleme se čím je vysoká variabilita očekávaných platů způsobena. # -------------------------------------------------------------------------------------- # Určeme základní míry šikmosti a špičatosti. # Funkce pro výpočet standardizované šikmosti a standardizované špičatosti # nejsou v základním R implementovány, je nutno instalovat některý z balíčků, # který tyto funkce obsahuje (např. moments). # POZOR! Funkce pod názvem skewness (šikmost) a kurtosis (špičatost) jsou součástí mnoha balíčků. # U funkce kurtosis() si musíme dát pozor na to, zda normálnímu rozdělení odpovídá # hodnota 3 (špičatost) nebo hodnota 0 (standardizovaná špičatost). Zapamatujme si, že v balíčku # moments se pod názvem kurtosis ukrývá funkce pro výpočet špičatosti, tj. špičatosti normálního # rozdělení odpovídá hodnota 3. skewness(dotaznik2$plat) # standardizovaná šikmost (hodnoty mezi -2 a 2 odpovídají přibližně symetrickému rozdělení) # Použití funkce kurtosis()z balíčku moments si vynutíme tím, že před název funkce uvedeme # název balíčku a "::". moments::kurtosis(dotaznik2$plat)-3 # standardizovaná špičatost (hodnoty mezi -2 a 2 odpovídají přibližně špičatosti normálního rozdělení) # Úkol D ------------------------------------------------------------------------------ # Komentujte šikmost a špičatost očekávaných platů. # Zkuste si představit, jak bude vypadat histogram očekávaných platů. # -------------------------------------------------------------------------------------- #///////////////////////////////////////////////////////////////////////////////////// # 1.3 Práce s chybějícími údaji (NA) ## Určeme průměrný měsíční příjem respondentů našeho šetření (prijmy). ## Využijme maximální možnou informaci, tj. vycházejte z původních neúplných záznamů (dotaznik). mean(dotaznik$prijmy) # Výstup NA získáme většinou tehdy, když někteří respondenti nezadali odpověď, # tj. máme neúplná data - data obsahující chybějící údaje (obvykle značeno NA). # Ve výpočtu je pak nutno dát toto najevo: mean(na.omit(dotaznik$prijmy)) # Do výpočtu průměru vstupují pouze úplné údaje dané proměnné. # nebo mean(dotaznik$prijmy,na.rm = TRUE) # Funkce pro výpočet základních charakteristik mají parametr na.rm. # Nastavíme-li na.rm = TRUE, bude funkce pracovat pouze s úplnými údaji dané proměnné. #///////////////////////////////////////////////////////////////////////////////////// # 1.4 Zaokrouhlování základních charakteristik # Chceme-li číselné charakteristiky publikovat, měli bychom je "rozumně" zaokrouhlit. # Pravidla pro zaokrouhlování číselných charakteristik najdete na # http://am-nas.vsb.cz/lit40/PRASTA/zaokrouhlovani.pdf (str. 2-3, dále manuál). ## Rozvažme si, jak zaokrouhlit základní číselné charakteristiky pro proměnnou plat. ## Při výpočtu charakteristik ošetřeme možnou přítomnost NA. (Zvykejme si :-)) # Nejdříve zjistíme rozsah výběru. length(na.omit(dotaznik2$plat)) # Zjistíme směrodatnou odchylku a rozhodneme se o jejím zaokrouhlení. sd(dotaznik2$plat, na.rm = TRUE) # Dle manuálu (Tabulka 1) bychom měli směrodatnou odchylku zaokrouhlit nahoru na 3 platné cifry, # tj. v našem případě na stovky, tj. na 50 600 Kč (funkce pro zaokrouhlení nahoru v R není). # Dále budeme tolerovat i zaokrouhlení směr. odchylky dle aritmetického zaokrouhlení na daný počet # platných cifer. round(sd(dotaznik2$plat, na.rm = TRUE),digits = -2)# zaokrouhlení na stovky # Vypočtené míry polohy (dolní kvartil, medián, průměr, horní kvartil) # zaokrouhlíme na stejný řád jako směrodatnou odchylku (v našem případě na stovky). # Minimum a maximum nezaokrouhlujeme (jde přímo o vybrané hodnoty analyzované proměnné). # Šikmost a špičatost zaokrouhlujeme na desetiny. # Variační koeficient (%) zaokrouhlujeme na desetiny. #///////////////////////////////////////////////////////////////////////////////////// # 1.5 Příprava tabulky základních charakteristik kvantitativní proměnné (dplyr) # Úkol E ------------------------------------------------------------------------------ # Vytvořte tabulku základních charakteristik pro očekávané platy respondentů. # Využijte funkce balíčku dplyr. # Počítejte i s tím, že by se v datech mohly objevit NA hodnoty. # Tabulku exportujte do pracovního adresáře. # Při publikování výsledků nezapomeňte charakteristiky správně zaokrouhlit # viz Manuál pro zaokrouhlování: http://am-nas.vsb.cz/lit40/PRASTA/zaokrouhlovani.pdf # -------------------------------------------------------------------------------------- # Je zřejmé, že výpis kódu pro tvorbu tabulky základních chrakteristik je poměrně dlouhý. # Vytvořme proto funkci, která nám dále přípravu tabulky základních charakteristik pro kvantitativní proměnné usnadní. #------------------------------------------------------------------------------------ # Poznámka (jak vytvořit funkci). # Každá funkce obsahuje název (např. ff), parametry (zde - a,b - ale může jich být více či méně podle potřeby) # a tzv. tělo, které se uvádí do {}. Uveďme si nejprve # příklad jednoduché funkce s názvem ff, která zobrazí součet dvou čísel. # Definování funkce ff=function(a,b){ a b paste("a+b =", a+b) } # Výstup funkce pro zadané parametry ff(a = 3,b = 4) #------------------------------------------------------------------------------------- # O programování funkcí využívajících ve svém těle funkce z balíčků dplyr (komplikovanější zadání # proměnné z daného datového rámce jako argumentu funkce) si můžete přečíst # na https://dplyr.tidyverse.org/articles/programming.html. souhrn_num = function(df,x){ # df ... daný datový rámec # x ... proměnná z datového rámce df library(dplyr) library(moments) xx = enquo(x) souhrn = df %>% summarise(rozsah_vyberu = length(na.omit(!!xx)), pocet_NA = sum(!complete.cases(!!xx)), minimum = min(!!xx,na.rm = TRUE), Q1 = quantile(!!xx,na.rm = TRUE, probs = 0.25), Q2 = quantile(!!xx,na.rm = TRUE, probs = 0.5), prumer = mean(!!xx,na.rm = TRUE), Q3 = quantile(!!xx,na.rm = TRUE, probs = 0.75), maximum = max(!!xx,na.rm = TRUE), smer_odchylka = sd(!!xx,na.rm = TRUE), var_koeficient = 100*smer_odchylka/prumer, sikmost = skewness(!!xx,na.rm = TRUE), spicatost = moments::kurtosis(!!xx,na.rm = TRUE)-3) colnames(souhrn) = c("rozsah výběru","počet chyb. hodnot", "minimum","dolní kvartil","medián","průměr","horní kvartil","maximum", "směr. odchylka","var. koeficient (%)", "stand. šikmost", "stand. špičatost") return(t(souhrn)) # definování výstupu funkce } souhrn_num(df = dotaznik2,x = plat) # výpočet tabulky číselných charakteristik pomocí funkce souhrn_num # Funkci souhrn_num si můžete překopírovat do samostatného skriptu a vždycky, když ji budete potřebovat, tak si jej # otevřete, funkci "projedete", tzn. uložíte do Environment, a můžete s ní dále pracovat. # Úkol F ------------------------------------------------------------------------------ # Pomocí funkce souhrn_num vytvořte tabulku základních charakteristik # pro průměrné měsíční příjmy respondentů našeho šetření (prijmy). Výsledky komentujte. # -------------------------------------------------------------------------------------- ################################################################################################### # 2. Vizualizace kvantitativní proměnné ########################################################### #///////////////////////////////////////////////////////////////////////////////////// # 2.1. Krabicový graf (základní R) # Pomoci krabicového grafu vizualizujte očekávané platy respondentů (proměnná plat) # Nastavení vzhledu grafického okna - nutné pouze pro základní R (v ggplotu řešeno jinak) par(mfrow = c(1,1), # rozdělení grafického okna - 1 řádek, 1 sloupec mar = c(4,4,1,2)) # okraje - c(dole, vlevo, nahoře, vpravo), viz # http://rstudio-pubs-static.s3.amazonaws.com/315576_85cccd774c29428ba46969316cbc76c0.html#mar boxplot(dotaznik2$plat) # základní krabicový graf # Úprava grafického výstupu boxplot(dotaznik2$plat, main="", xlab="Očekávaný měsíční plat (Kč)", ylab="", ylim=c(0,1.1*max(dotaznik2$plat)), col="grey", # barva výplně "krabice" boxwex = 0.5,# relativní změna šířky "krabice" vzhledem k původnímu nastavení horizontal = TRUE) # horizontální poloha "krabice" # Doplnění bodu reprezentujícího průměr (low level funkce - vykresluje do stávajícího grafu) points(mean(dotaznik2$plat,na.rm=TRUE),1, # souřadnice bodu pch=3, # typ bodu (možno volit křížek, kolečko, ...) col="red") # barva bodu # Úkol G ------------------------------------------------------------------------------ # Pomocí krabicového grafu vizualizujte průměrné měsíční příjmy respondentů (prijmy). # Krabicový graf umístěte vertikálně a doplňte jej o bod reprezentující průměr # (pozor na pořadí souřadnic). boxplot(dotaznik2$prijmy, main="", xlab="", ylab="Průměrný měsíční příjem (Kč)", ylim=c(0,1.1*max(dotaznik2$prijmy))) # Doplnění bodu reprezentujícího průměr points(1,mean(dotaznik2$prijmy,na.rm=TRUE), # souřadnice bodu pch=3, col="red") # -------------------------------------------------------------------------------------- #///////////////////////////////////////////////////////////////////////////////////// # 2.2. Krabicový graf (ggplot2) # Pomocí krabicového grafu vytvořeného pomocí balíčku ggplot2 vizualizujte # očekávané platy respondentů (proměnná plat) # Základní verze ggplot(dotaznik2,aes(x="",y=plat))+ # estetika geom_boxplot() # specifikace výstupu - krabicový graf # Upravujeme design grafu... ggplot(dotaznik2,aes(x="",y=plat))+ stat_boxplot(geom = "errorbar", width = 0.1) + # ohraničení fousů (musí být před geom_boxplot) geom_boxplot()+ coord_flip()+ # horizontální orientace labs(x="", y="Očekávaný měsíční plat (Kč)", title="")+ theme_bw()+ # přidání bodu reprezentujícího průměr geom_point(aes(x="",y=mean(dotaznik2$plat,na.rm=T)), # souřadnice (POZOR! Jde o součást funkce ggplot, neměníme pořadí souřadnic!) color="red", shape=3)+ theme(plot.margin = unit(c(1,1,1,1), "cm")) # nastavení okrajů (nahoře, vpravo, dole, vlevo) - pokud potřebujeme... # Úkol H ------------------------------------------------------------------------------ # Pomocí krabicového grafu vytvořeného pomocí balíčku ggplot2 vizualizujte průměrné měsíční příjmy respondentů (prijmy). # Krabicový graf umístěte vertikálně a doplňte jej o bod reprezentující průměr # (pozor na pořadí souřadnic). "Tahák" k balíčku ggplot2 najdete na https://rstudio.com/wp-content/uploads/2015/03/ggplot2-cheatsheet.pdf # --------------------------------------------------------------------------------------------- #///////////////////////////////////////////////////////////////////////////////////// # 2.3. Histogram (základní R) # Pomoci histogramu vizualizujte očekávané platy respondentů (proměnná plat) hist(dotaznik2$plat) # Základní histogram hist(dotaznik2$plat,breaks=20) # Vyzkoušejte si, jak různá nastavení parametru breaks ovlivní vzhled grafu. # Prohlédněte si nápovědu k funkci hist() a posuďte význam argumentu breaks. # Zadáváte-li argument breaks jako číslo, udává počet tříd (sloupečků), do nichž je kvantitativní proměnná tříděna. # Implicitní nastavení parametru breaks odpovídá počtu tříd navrženému Sturgesovým pravidlem. # Upravujeme design histogramu... hist(dotaznik2$plat, breaks = 20, main = "", xlab = "Očekávaný měsíční plat (Kč)", ylab = "počet respondentů", ylim = c(0,90), xlim = c(0,500000), col = "grey", # barva výplně sloupců border = "black", # barva ohraničení sloupců labels = TRUE) # popisky (počet respondentů v jednotlivých třídách) # Kvantitativní znaky se velmi často vizualizují pomocí tzv. hustoty pravděpodobnosti. # Zjednodušeně řečeno lze říci, že hustota pravděpodobnosti je křivka, jejíž tvar odpovídá # "obálce" histogramu. Tato funkce je z hlediska osy y normována tak, že plocha pod ní je rovna 1. # Vizualizujme očekávané platy respondentů (proměnná plat) pomocí hustoty pravděpodobnosti. # Význam hustoty pravděpodobnosti demonstrujme tím, že její graf doplníme o příslušně normovaný histogram. # Nejdříve musíme vykreslit histogram a pomocí parametru freq nastavit jeho normování tak, # aby plocha vymezena sloupci histogramu byla rovna 1. hist(dotaznik2$plat, breaks = 20, main = "", xlab ="Očekávaný měsíční plat (Kč)", ylab="Hustota pravděpodobnosti", ylim = c(0,0.00005), xlim = c(0,500000), col="grey", border="black", freq=FALSE) # normování histogramu (změna měřítka na ose y) # Následně přidáme graf tzv. empirické hustoty pravděpodobnosti (odhadnuté z naměřených dat). lines (density(dotaznik2$plat), col = "blue", lwd = 2)# tloušťka čáry # Graf empirické hustoty pravděpodobnosti lze použít pro empirické posouzení, # zda lze analyzovanou proměnnou modelovat tzv. normálním rozdělením. # Graf hustota pravděpodobnosti normálního rozdělení je znám pod názvem Gaussova křivka. # Pro definování normálního rozdělení je nutno znát jeho parametry, jimiž jsou střední hodnota # a směrodatná odchylka. Chceme-li naměřená data aproximovat normálním rozdělením, # pak střední hodnotu odhadujeme průměrem naměřených dat a směrodatnou odchylku (výběrovou) směr. odchylkou # naměřených dat. # Doplňme náš graf o graf hustoty pravděpodobnosti normálního rozdělení, kterým bychom rozdělení # analyzovaného znaku mohli modelovat. xfit=seq(min(dotaznik2$plat, na.rm=TRUE), max(dotaznik2$plat, na.rm=TRUE), length = 100) # hustota pravděpodobnosti normálního rozdělení, které by mohlo sloužit jako model pro naše data yfit=dnorm(xfit, mean=mean(dotaznik2$plat, na.rm=TRUE), sd=sd(dotaznik2$plat, na.rm=TRUE)) lines(xfit, yfit, col="red", lwd=2) # Logika: # xfit - podle našich zjištěných hodnot si vygenerujeme posloupnost od minimálního zjištěného # platu do maximálního zjištěného platu # yfit - kdyby naše data pocházela z normálního rozdělení, jak by vypadaly y-souřadnice # Z této vizualizace se zdá, že očekávaný měsíční plat po ukončení studia # není "rozumné" modelovat normálním rozdělením. # Říkáme, že se zdá, že data nejsou výběrem z normálního rozdělení. # Úkol I ------------------------------------------------------------------------------ # Pomocí hustoty pravděpodobnosti vizualizujte průměrné měsíční příjmy respondentů (prijmy). # Význam hustoty pravděpodobnosti demonstrujme tím, že její graf doplníme o příslušně normovaný histogram. # Graf doplňte rovněž o hustotu p-sti normálního rozdělení, které by mohlo sloužit # k modelování průměrných měsíčních příjmů studentů 1. ročníku FEI. # --------------------------------------------------------------------------------------------- #///////////////////////////////////////////////////////////////////////////////////// # 2.4. Histogram (ggplot2). # Pomoci histogramu připraveného za pomocí balíčku ggplot2 vizualizujte očekávané platy respondentů (proměnná plat) ggplot(dotaznik2,aes(x=plat))+ geom_histogram(bins=20) # bins = počet sloupců # Úpravujeme design grafu... ggplot(dotaznik2,aes(x=plat))+ geom_histogram(bins = 20, color="black", fill="grey")+ labs(x="Očekávaný měsíční plat (Kč)", y="počet respondentů", title="")+ theme_bw() # Pomocí hustoty pravděpodobnosti vizualizujme očekávaný měsíční plat respondentů (plat). # Význam hustoty pravděpodobnosti demonstrujme tím, že její graf doplníme o příslušně normovaný histogram. # Graf doplňme rovněž o hustotu p-sti normálního rozdělení, které by mohlo sloužit # k modelování měsíčních platů, které studenti 1. ročníku FEI očekávají po ukončení studia. a = ggplot(dotaznik2,aes(x=plat))+ geom_histogram(aes(y=..density..), # normování histogramu (změna měřítka na ose y) bins = 20, color="black", fill="grey")+ geom_density(color="blue")+ # empirická hustota p-sti labs(x="Očekávaný měsíční plat (Kč)", y="Hustota pravděpodobnosti", title="")+ theme_bw() a # Vytvoříme pomocný datový rámec obsahující dva sloupce odpovídající souřadnicím pro generování # hustoty pravděpodobnosti normálního rozdělení, kterým bychom rozdělení # analyzovaného znaku mohli modelovat. pom = data.frame(xfit=seq(min(dotaznik2$plat, na.rm=TRUE), max(dotaznik2$plat, na.rm=TRUE), length = 100)) pom$yfit=dnorm(pom$xfit, mean=mean(dotaznik2$plat, na.rm=TRUE), sd=sd(dotaznik2$plat, na.rm=TRUE)) # Dříve vygenerovaný graf doplníme o hustotu p-sti normálního rozdělení, které by mohlo sloužit # k modelování očekávaných měsíčních platů. a + geom_line(data = pom,aes(x = xfit,y = yfit),color = "red") # POZOR! Tvar histogramu v základním R a v ggplotu může být lehce jiný! Důvodem je lehce odlišný # přístup při "tvorbě sloupců". # --------------------------------------------------------------------------------------------- #///////////////////////////////////////////////////////////////////////////////////// # 3. Identifikace odlehlých pozorování (outliers) a jejich ev. odstranění # Zjistěme, zda některé pozorované hodnoty očekávaných platů nelze považovat za odlehlá pozorování. # První možností je zjistit, zda některé pozorované hodnoty očekávaných platů neleží mimo vnitřní hradby. # Určíme vnitřní hradby. IQR(dotaznik2$plat,na.rm = TRUE) # mezikvartilové rozpětí (mohli bychom určit i "ručně" jako rozdíl horního a dolního kvartilu) dolni_mez = quantile(dotaznik2$plat,probs = 0.25,na.rm = TRUE)-1.5*IQR(dotaznik2$plat,na.rm = TRUE) horni_mez = quantile(dotaznik2$plat,probs = 0.75,na.rm = TRUE)+1.5*IQR(dotaznik2$plat,na.rm = TRUE) # Vyfiltrujeme pozorované hodnoty očekávaných platů, které leží mimo vnitřní hradby, tj. jsou odlehlými pozorováními. dotaznik2 %>% select(plat) %>% filter(plathorni_mez) # Je zřejmé, že 6 respondentů uvedlo v našem šetření očekávaný plat, # který se "velmi" lišil od údajů, které uvedli ostatní respondenti. # Tyto hodnoty identifikované metodou vnitřních hradeb nazýváme odlehlými pozorováními. # Jde přesně o ty hodnoty, které jsou v krabicovém grafu reprezentovány body. boxplot(dotaznik2$plat) # Podívejte se znovu na výpis odlehlých pozorování # a rozvažte si, proč na grafu vidíme pouze 3 body zatímco # ve výpisu je uvedeno 6 hodnot. # Druhou možností, jak získat výpis odlehlých pozorování je využití funkce boxplot. boxplot(dotaznik2$plat,plot = FALSE) # Tímto způsobem získáme výpis několika složek funkce boxplot(). # $stats uvádí uspořádanou pětici čísel definující: # konec dolního vousu,dolní kvartil, medián, horní kvartil, konec horního vousu # $n udává rozsah výběru # $out ... výpis odlehlých hodnot (angl. outliers) # Ostatní atributy této funkce komentovat nebudeme. # Odlehlá pozorování tedy můžeme snadno získat funkcí boxplot(dotaznik2$plat)$out # Rozhodnutí o tom, jak se zjištěním odlehlých pozorování naložíme je # na rozhodnutí analytika, resp. na výsledku jeho konzultace se zadavatelem studie. # V každém případě bychom měli ve výzkumné zprávě uvést, # že jsme odlehlá pozorování identifikovali a jak jsme s nimi naložili. # Dále si ukážeme, jak postupovat, když se odlehlá pozorování rozhodneme # z dalšího zpracování vyloučit. # Nikdy bychom si neměli "zničit" původní data, # proto si vytvoříme pomocnou proměnnou,v níž budou odlehlá pozorování vymazána, # tj. nahrazená NA. outliers = boxplot(dotaznik2$plat)$out dotaznik2 = dotaznik2 %>% mutate(plat_out = ifelse(plat %in% outliers,NA,plat)) # Funkce ifelse() má stejnou strukturu jako funkce KDYŽ() v MS Excel. # Podívejme se jak jsou rozloženy očekávané měsíční platy poté # co jsme se rozhodli do zpracování odlehlá pozorování nezařadit. boxplot(dotaznik2$plat_out) hist(dotaznik2$plat_out) # Poznámka 1: Před publikováním bychom si s designem grafů měli "pořádně pohrát" # - nejlépe s využitím balíčku ggplot2. # Poznámka 2: Kdybychom na krabicovém grafu viděli znovu "nějaká" odlehlá pozorování, # další odstranění už bychom nedělali. Nebyla by to už odlehlá pozorování z původních dat. # Poznámka 3: Nejjednoduším způsobem posouzení odlehlých hodnot je individuální posouzení, # neboli jako expert se rozhodnu, že taková a taková pozorování budou vyloučena, protože # k tomu mám odborný vědecký důvod. Vše opět zapíšu a vysvětlím ve výzkumné zprávě. # Úkol J ------------------------------------------------------------------------------ # Identifikujte, zda pozorované průměrné měsíční příjmy respondentů (prijmy) neobsahují odlehlá pozorování. # Pokud ano, vylučte je z dalšího zpracování a srovnejte číselné charakteristiky původní proměnné # a proměnné, z níž jste odstranili odlehlá pozorování. # Obdobné srovnání proveďte s grafickými výstupy. # --------------------------------------------------------------------------------------------- # Prohlédněte si jak velký vliv mělo odstranění odlehlých pozorování na jednotlivé charakteristiky. # Všimněte si, že tímto způsobem jsme se nezbavili problému s tím, že někteří respondenti uváděli nulové příjmy... # Rozhodnutí, jak s tímto zjištěním naložit je na nás. Jen o tom nesmíme # zapomenout učinit záznam ve výzkumné zprávě... ###################################################################################### #### Shrnutí aneb Jak postupovat při analýze kvantitativní proměnné ################## ###################################################################################### # 1.) Ověříme, zda jsou data klasifikována jako typ numeric. # 2.) Vizuálně posoudíme data dle krabicového grafu (odlehlá pozorování, rozpětí pozorovaných dat, ...) # 3.) Nastavíme formát číselných charakteristik (options()) # 4.) Pokud uznáme, že je vhodné data "vyčistit" (opravit či odstranit odlehlá pozorování, nesmyslné údaje apod.), # učiníme tak a NEZAPOMENEME o tom učinit záznam do výzkumné zprávy. # 5.) Posoudíme základní číselné charakteristiky. # 6.) Připravíme grafické výstupy (krabicový graf, histogram, # popř. histogram doplněný o empirickou hustotu a hustotu normálního rozdělení, # pomocí něhož bychom teoreticky mohli analyzovanou proměnnou modelovat). # 7.) Na základě histogramu, šikmosti a špičatosti se pokusíme posoudit možnost modelování # analyzovaných dat normálním rozdělením. # 8.) Výsledky naší analýzy shrneme do výzkumné zprávy. Při prezentaci číselných # charakteristik NEZAPOMENEME zaokrouhlit číselné charakteristiky podle # pravidel zaokrouhlování (viz manuál). Zaokrouhlené číselné charakteristiky # uvádíme nejen v tabulkách, ale i v komentářích (v textu).