NMSA230 - Softwarové prostředky pro matematiku a stochastiku Zimný semester 2017/2018 | Cvičenie 2 | Ut 22/10/18 II. Štatistický software R (štatistika a grafika) Všeobecné informáie Na cvičeni budeme pracovať s programom R (Team Development Core, 2016), ktorý je dostupný pod GNU GPL licenciou (free of charge) na nasledujúcej adrese https://www.r-project.org. K dispozícii sú distribúcie s priamou podporou pre OS Windows, Linux aj Macintosh. Základnú inštaláciu programu R je možne rozšíriť o ďalšie dodatočné knižníce (packages), ktoré sú k dispozícii na rôznych online repozitároch (zoznam hlavných repozitárov je na adrese https://cran.r-project.org/mirrors.html. Na inštaláciu baličku v programe R slúži príkaz install.packages("***"), kde namiesto *** je potrebné uviesť správne meno príslušného balíčka. Okrem samotného programu R a rozširujúcich knižníc sú k dispozícii aj rôzne užívateľské rozhrania (user interfaces) , ktoré je možne dodatočne nainštalovať a umožňujú (v určitých smeroch) jednoduchšiu a prehľadnejšiu prácu s Rkom. Najznámejší a pravdepodobne aj jeden z najlepších R interface je RStudio (defaultne nainštalované na každom počítači v tejto posluchárni). Program druhého cvičenia NMSA 230 Oboznámenie sa so zkladnými funkciami a nástrojmi, ktoré su v programe R určené k vytváraniu rôznych grafických výstupov (grafy, boxploty, diagramy). Práca s reálnymi datami: vytvorenie jednoduchých popisných štatistík, ich interpretácia a následne zobrazeie pomocou vhodných graficých výstupov. Práca s niektorými komplexnejšími príkazmi, pomocou ktorých je možné vytvárať automatizované R skripty (napr. for, while, if, list.files(), which(), match(), a mnoho ďalších) Podrobnejšieho “sprievodca” pre jednoduchú prácu s R programom: Hrátky s R (autor: doc. Arnošt Komárek). (k Hrátkam s R sú potrebné datove súbory: auta2004.dat, auta2004.csv a auta2004.xls) Užitočné materiály pre prácu so štatistickým softwarom R Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (PDF súbor) Komárek, A.: Základy práce s R. (PDF súbor) Kulich, M.: Velmi stručný úvod do R. (PDF súbor) De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808) 1. Základné grafické nástroje v programe R V programe R je k dispozícii veľký výber rôznych funkcii a príkazov určených k vytváraniu variabilných grafických výstupov. Dalšíe dodatočné príkazy a funkcie je možne získať pri rozšírení základnej nainštalovaním niektorých odatočných balíčkov (packages) - tie sú väčšinou špecificky zamerané na konkrétnu oblasť, balíček často poskytuje komplexnú sadu nástrojov na riešenie určitého špecifického problému. Okrem toho je samozrejme možné vytvárať si vlastne príkazy a funkcie a defaultne grafické nástroje v Rku v mnohom zmeniť, vylepšiť a zdokonaliť. Grafické funkcie v programe R možeme podľa spôsobu fungovania začleniť do troch základných kategórii: Funkcie, ktoré vytvárajú nové grafické rozhranie Pri volaní ktorejkoľvek fukcie z tejto kategórie automaticky dôjde k vytvoreniu nového grafického okna (graph window) v ktorom sa vykresli samotný výsledný obrázok(graf). Toto okno je ponechané ako aktívne a pomocou ďalších funkcii z druhej skupiny (viď nižšie) je možné dopĺňať ďalšie objekty do již existujúceho obrázku, prípadne meniť rôzne nastavenia a vlastnosti. Základné príkazy (napr.): plot(), dotchart(), hist(), barplot(), pie(), boxplot(), pairs(), atď.; Funkcie, ktoré dopĺňajú již vytvorený graf Do druhej kategórie patria funkcie a príkazy, ktoré fungujú len v prípade, že bol pred ich použitím zavolaný niektorá funkcia z prvej kategórie a vytvorené okno je stále aktívne. Pomocou príkazov a funkcii z tejto kategórie je možné dopĺňať ďalšie (často netriviálne) objekty do již existujúceho grafu (grafického okna - graph window). Základné príkazy (napr.): points(), lines(), abline(), text(), legend(), atď.; Funkcie určené na manipuláciu s grafmi Do tejto kategórie patria funkcie a príkazy, ktoré primárne nie su schopné vytvárať obrázok, ani dopĺňať objekty do již existujúceho obrázku, ale služia na jednoduchú manipuláciu s vytvorenými grafickými oknami (napr. uloženie obr., vytvorenie niekoľkýh grafov v jednom okne, vytvorenie nového okna, a pod.) Základné príkazy (napr.): dev.new(), X11(), dev.off(), postscript(), pdf(), png(), par(), layout(), atď.; Pomocou R helpu zistite, k čomu jednotlivé príkazy a funkcie slúžia, ako sú implementované, aké dodatočné parametre je možné pri volaní funkcie využiť a každú z uvedených funkcii použijte aspoň v jednom grafe, ktorý samostatne vytvoríte. Niekoľko názorných príkladov, ktoré využívajú niektoré z uvedených funkcii, sú uvedené nižšie. Všimnite si aj úlohu dodatočných parametrov, ktoré sa pri volaní konkrétnej funkcie spomínaju. Príklad 1 | datový súbor ‘mtcars’ Datový súbor ‘mtcars’, ktorý je štandardne k dispozícii v klasickej inštalácii programu R obsahuje údaje o 32 rôznych automobiloch predávaných v USA (typ auta, motoru, spotreba, atď.). x <- mtcars[order(mtcars$mpg),] # sort by mpg x$cyl <- factor(x$cyl) # it must be a factor x$color[x$cyl==4] <- "red" x$color[x$cyl==6] <- "blue" x$color[x$cyl==8] <- "darkgreen" dotchart(x$mpg,labels=row.names(x),cex=.7,groups= x$cyl, main="Gas Milage for Car Models\ngrouped by cylinder", xlab="Miles Per Gallon", gcolor="black", color=x$color) Alternatívný náhľad dostaneme prostredníctvom matice scatterplotov, pomocou príkazu pairs(): pairs(mtcars[,c(1,3,4,6)], col=mtcars$cyl, pch = mtcars$am) Sofistikovanejšie verzie aj s odhadom hustoty a regresnej závislosti možeme konštruovať pomocou knižnice library(car) príkazom scatterplot.matrix(): panel.hist <- function(x, ...) { usr <- par("usr"); on.exit(par(usr)) par(usr = c(usr[1:2], 0, 1.5) ) h <- hist(x, plot = FALSE) breaks <- h$breaks; nB <- length(breaks) y <- h$counts; y <- y/max(y) rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...) } panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) { usr <- par("usr"); on.exit(par(usr)) par(usr = c(0, 1, 0, 1)) r <- abs(cor(x, y)) txt <- format(c(r, 0.123456789), digits = digits)[1] txt <- paste0(prefix, txt) if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt) text(0.5, 0.5, txt, cex = 2) } pairs(~mpg+disp+drat+wt + cyl, data=mtcars, lower.panel=panel.smooth, upper.panel = panel.cor, diag.panel = panel.hist, pch = 21) Existujú samozrejme mnohé iné možnosti a ďalšie grafiké nástroje v Rku, pomocou ktorých vytvoriť ‘case-specific’ obrázky a grafy. Hlavný dôraz je pritom kladený na vizuálnu stránku, graf by mal dostatočne názorne zobrazovať podstatnú štruktúru datového súboru. V prípade menšieho počtu premenných je celkom zaujímavou možnosťou aj zobrazenie marginálnych hustot pre nejaké dve konkrétne veličíny (z dôvodu obmedzeného vnímania viac-rozmerných grafov sa väčšinou nepoužíva viac, ako dve veličiny súčastne): library(MASS) library("RColorBrewer") colors <- brewer.pal(5, "Reds") z <- kde2d(mtcars[,1], mtcars[,4], n=5) plot(mtcars[,c(1,4)], xlab="Miles per Galon", ylab="Horse Power", pch=21, bg = "yellow") contour(z, drawlabels=FALSE, nlevels=4, col=colors, add=TRUE) Knižnica MASS je potrebná k spočítaniu jadrového odhadu dvojrozmernej hustoty (príkaz kde2d()) a knižnica RColorBrewer je potrebná na vytvorenie farebnej palety, podľa ktorej sú vyfarbené kontúry. Boxploty Predstavu o štruktúre datového súboru získame (aj keď len v určitej obmedzenej miere) aj pomocou boxplotov boxplot(mpg ~ cyl, col = brewer.pal(3, "Blues"), xlab = "number of cylinders", ylab = "miles per galon", main= "", data = mtcars) prípadne verzia, ktorá zohľadňuje dve rôzne kategorické premenné (v tomto konkrétnom prípade spotrebu udávanu v míľach prejdených na jeden galón a identifikátor automatickej prevodovky: Automatic = 0 a Manual = 1). boxplot(mpg ~ cyl + am, col = c(brewer.pal(3, "Greens"), brewer.pal(3, "Blues")), xlab = "No. of cylinders (vs.) Auto Transmission Indicator", data = mtcars) Alebo dokonca zložitejší prípad zobrazujúci spotrebu v závislosti na dvoch rôznych faktoroch a 95% konfidenčným intervalom pre skutočnú hodnotu médianu (s voľbou parametru notched = TRUE). boxplot(mpg ~ cyl * am, col = c(brewer.pal(3, "Greens"), brewer.pal(3, "Blues")), xlab = "Number of Cylinders", ylab = "Miles per Galon | Automatic Transmission", main= "", notch = T, data = mtcars) Alebo celkom sofistikovany ‘histogram’, ktorý inkorporuje aj odhad hustoty (pomocou príkazu vioplot() z knižnice ‘vioplot’). library(vioplot) x1 <- mtcars$mpg[mtcars$cyl==4] x2 <- mtcars$mpg[mtcars$cyl==6] x3 <- mtcars$mpg[mtcars$cyl==8] vioplot(x1, x2, x3, names=c("4 cyl", "6 cyl", "8 cyl"), col="gold") title("Violin Plots of Miles Per Gallon") S využitím dodatočnej knižnice MASS je možné zavolať funkciu parcoord() ktorá je užitočna, keď chceme vizuálne posúdiť mieru korelácie medzi jednotlivými sledovanými premennými. Ako by ste nasledujúci graf interpretovali? Resp. aké závery je možné z tohto grafu vyvodiť? library(MASS) attach(mtcars) colorVector <- rep("brown", dim(mtcars)[1]) colorVector[am == 0] <- "brown1" colorVector[cyl == 6] <- "blue4" colorVector[cyl == 6 & am == 0] <- "blue1" colorVector[cyl == 4] <- "darkgreen" colorVector[cyl == 4 & am == 0] <- "green" parcoord(mtcars, col = colorVector, lty = gear) Príklad 2 | datový súbor ‘Orange’ V druhom príklade použijeme datový súbor ‘Orange’, s ktorým sme už pracovali. V datasete sú zaznamenané informácie o piatich druhoch pomarančovníkoch - zaznamenaváme vek stromu v dňoch a priemer kmeňa v cm. Najprv využijeme štandardný príkaz plot() a následne sa pokusíme graf vylepšiť pomocou dodatočných parametrov a dalších funkcii a príkazov z druhej kategórie. plot(Orange[,"circumference"] ~ Orange[,"age"]) A teraz vylepšená (snáď krajšia) verzia toho istého grafu: plot(0,0, xlim = c(0, 1800), ylim = c(0, 250), pch = "", xlab = "age [days]", ylab = "circumference [cm]", main = "Orange Trees") for (i in 1:5){ points(Orange[Orange[,1] == i,3] ~ Orange[Orange[,1] == i,2], pch = 21, bg = i) lines(Orange[Orange[,1] == i,3] ~ Orange[Orange[,1] == i,2], col = i) } legend(0, 250, legend = c("Tree type #1", "Tree type #2", "Tree type #3", "Tree type #4", "Tree type #5"), col = "black", bg = "lightgray" , pch = 21, lty = c(1,1,1,1,1), pt.bg = 1:5) Príklad 3 | multiple plot output Pomocou príkazov par() a layout() je možné kombinovať niekoľko grafov/obrázkov v jednom grafickom okne. par(mfrow = c(1,2)) hist(Orange[,"circumference"], col = "lightblue", xlab = "Circumference [cm]", ylab = "Frequency", freq = F, main = "") lines(density(Orange[,3]), col = "red", lwd = 2, lty = 2) barplot(Orange[,3], col = sort(rep(1:5,7)), xlab = "Measurements by Trees", ylab = "Circumference [cm]") Užitočné Program R poskytuje obrovskú škálu nástrojov a funkcii k vytváraniu a modifikovaniu grafickych výstpov. Ako názorná ukážka možnosti, ktoré sú k dispozícii, môže poslúžiť jednoduché demo, ktoré je defaultne k dispozícii v Rku: príkaz demo("graphics"). Analogickym spôsobom funguje demo ku všetkym knižniciam, ktoré su v Rku nainštalované (príkaz demo() a následne demo("***"), kde namiesto *** je potrebné dopĺniť názov príslušného balíčku - package. Každá zo spomínaných funkcii je implementována s dodatočnými parametrami, ktore menia jej defaultne nastavenie a ovplyvňujú finálny vysledok výsledného grafu (viď help k príslušnej funkcii). Je dôležité tieto parametre poznať a náležite ich pri práci aj využívať. 2. Použitie na reálnych datach Z datovych súborov, ktoré sú defaultne k dispozícii v štandardnej inštalácii programu R (príkaz data()) si vyberte jeden datový súbor a pomocou grafických nástrojov zostrojte niekoľko obrázkov, ktoré nejakým spôsobom vypovedajú o nameraných hodnotách. Pomocov dvojice príkazov pdf() a dev.off() získame výsledný grafický výstup vo foprmáte PDF (ktorý sa automaticky vytvorí pri volaní funkcie pdf() a ukonči/zprístupni volaním funkcie dev.off() na záver). pdf("obrazok1.pdf", width = 6, height = 4) par(mfrow = c(1,2)) hist(Orange[,"circumference"], col = "lightblue", xlab = "Circumference [cm]", ylab = "Frequency", freq = F, main = "") lines(density(Orange[,3]), col = "red", lwd = 2, lty = 2) barplot(Orange[,3], col = sort(rep(1:5,7)), xlab = "Measurements by Trees", ylab = "Circumference [cm]") dev.off() Pre prácu s LaTeXom je užitočné poznať aj príkaz na vytvorenie obrázku/grafu vo vektorovej grafike, ktorý je všeobecne kvalitnejší. Slúži k tomu príkaz postscript(). Dodatočné grafické parametre sú totožné, ako pri vytváraní výstupov vo formáte PDF, teda príkaze pdf(): postscript("obrazok1.eps", width = 6, height = 4) par(mfrow = c(1,2)) hist(Orange[,"circumference"], col = "lightblue", xlab = "Circumference [cm]", ylab = "Frequency", freq = F, main = "") lines(density(Orange[,3]), col = "red", lwd = 2, lty = 2) barplot(Orange[,3], col = sort(rep(1:5,7)), xlab = "Measurements by Trees", ylab = "Circumference [cm]") dev.off() Užitočné Ako už bolo uvedené, v programe R je k dispozícii celá škála rôznych balíčkov a knižníc, z ktorých su mnohé venované aj grafike. Podrobnosti sú napr. na webovej stranke http://www.r-graph-gallery.com. Vôčšína z uvedených príkladov vyžaduje inštaláciu potrebnej knižnice (package) - inštalácia pomocou príkazu install.packages() a následna inicializácia knižnice pomocou príkazu library(). 3. Knižnica ggplot2 Alternatívu ku klasickým grafickým nástrojom, ktoré su k dispozícii v programe R, ponúka dodatočna knižnica ggplot2, ktorú naištalujete a inicializujete štandardnými príkazmi install.packages("ggplot2") library(ggplot2) Základnou funkciou v tejto knižnici je funkcia qplot(), ktorá na rozdiel od štandardnej funkcie plot(), ponúka ďaleko väčšie možnosti pri vytváraní grafov a obrázkov (pozri help ?qplot()). V datovom súbore mtcars najprv upravíme niektoré premenné - faktorizujeme ich a priradíme príslušne názvý jednotlivých úrovni. mtcars$gear <- factor(mtcars$gear,levels=c(3,4,5), labels=c("3gears","4gears","5gears")) mtcars$am <- factor(mtcars$am,levels=c(0,1), labels=c("Automatic","Manual")) mtcars$cyl <- factor(mtcars$cyl,levels=c(4,6,8), labels=c("4cyl","6cyl","8cyl")) Následne už možeme použiť funkciu qplot() a vykresliť odhad hustoty pre spotrebu (miles per gallon) v závislosti na počte prevodových stupňov (premenná gear). qplot(mpg, data=mtcars, geom="density", fill=gear, alpha=I(.5), main="Distribution of Gas Milage", xlab="Miles Per Gallon", ylab="Density") Pomocou tej istej funkcie qplot(), ale s inou voľbou parametru geom = "...", je možne celkovo meniť vzhľad a typ výsledného grafu. qplot(gear, mpg, data=mtcars, geom=c("boxplot", "jitter"), fill=gear, main="Mileage by Gear Number", xlab="", ylab="Miles per Gallon") Prípadne vykresliť regresnú závislosť spotreby na váhe automobilu (aj s príslušnými pásmi spoľahlivosti pre regresnú přímku). qplot(wt, mpg, data=mtcars, geom=c("point", "smooth"), method="lm", formula=y~x, color=cyl, main="Regression of MPG on Weight", xlab="Weight", ylab="Miles per Gallon") Užitočné Knižnica ggplot2 funguje na princípe postupných vrstiev (layers), v ktorých sa vykresľuje výsledný graf. Jednotlivé vrstvy sa prídávajú pomocou vhodných príkazov, ako napr. geom_points(), geom_abline(), geom_area() geom_bar(), geom_contour(), geom_polygon(), atď. (zoznam napr. na stránke http://sape.inf.usi.ch/quick-reference/ggplot2/geom). Jednoduchý príklad s využitím niekoľkých vrstiev: dataSum <- plyr::ddply(mtcars, "gear", plyr::summarise, mean = mean(mpg), sd = sd(mpg)) ggplot() + geom_point(data = mtcars, aes(x = gear, y = mpg)) + geom_point(data = dataSum, aes(x = gear, y = mean), colour = 'red', size = 3) + geom_errorbar(data = dataSum, aes(x = gear, y = mean, ymin = mean - sd, ymax = mean + sd), colour = 'red', width = 0.4) + ggtitle("Car consumption summary given the number of gears") Domáca úloha (Deadline: 3. cvičenie | Ut: 06.11.2018) Použijte stejný datový súbor, ktorý ste použili pri prvej domácej úlohe (výpočet základných popisných charakteritik) a k napočítaným popisným štatistikam vytvorte vhodný graf (obrárok), na ktorom bude možné spočítané hodnoty vizuálne odčítať a v kontexte súboru dat interpretovať. Skript, v ktorom počítate základne popisné charakteristiky v datach, doplňte o časti kódu, ktorý bude vytvárať príslušne obrázky. Vytvorené grafické výstupy náležite okomentujte a finálne grafy uložte ako eps, alebo pdf súbory. Výsledný R skript si opäť uložte (napr. ako ‘Rskript.R’, alebo plain txt ‘Rskript.txt’). V prípade, že si skript uložíte ako ‘R script’ (súbor ‘menoSuboru.R’), je možné tento súbor priamo načítať do programu R pomocou príkazu source("menoSuboru.R"). Po načítaní súboru programautomaticky spracuje predpripravené príkazy.