NMSA230 - Úvod do programování v R Zimný semester 2022/2023 | Cvičenie 1 | St 05/10/22 • zdrojový Rmd súbor I. Štatistický program R (úvod a základy) Základným cieľom predmetu NMSA230 je oboznámiť sa so štatistickým programom R (R Development Core Team, 2020), ktorý je jedným z najčastejšie používaných (vývojových) softwarových nástrojov určených k spracovaniu, príprave a následnej štatistickej analýze dat. V yužitie a uplatnenie nájde ale aj pri reportovaní výsledkov, vytváraní prezentácii, alebo plne automatických analytických webových stránok. Jedná sa o GNU projekt, podobný programu S a primárne je navrhnutý pre štatisticku analýzu dat. Program R (dostupný pod GNU GPL licenciou) je k dispozícii k stiahnutiu (free of charge) na adrese https://www.r-project.org K dispozícii sú distribúcie s podporou pre OS Windows, Linux, aj Macintosh. Základnú inštaláciu programu R je možne jednoducho rozšíriť pomocou dodatočných knižníc (balíčkov), 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). Jednotlivé R knižnice sú ale tvorené samotnými užívateľmi programu R a ich správne fungovanie nie je garantované – je preto namieste určitá opatrnosť a hlavne aktívne premýšľanie (a dostatočná kontrola) pri ich používaní. V prípade nutnosti použíť niektorý z dodatočných balíčkov (nad rámec štandardnej inštalácie), bude konkrétny príkaz nutný k inštalácii explicitne uvedený v R skripte. K dispozícii sú aj rôzne nástroje a knižnice pre interakciu s iným softwarom (optimalizačný toolbox mosek, databazove systémy ako SQL, developerské nástroje ako C++ a podobne). Pre užívateľov programu R sú k dispozícii aj user-friendly grafické rozhrania, ktoré je možné dodatočne nainštalovať a umožňujú (v určitých smeroch) jednoduchšiu a prehľadnejšiu prácu. Najznámejší a pravdepodobne aj jeden z najlepších R interface je RStudio. K dispozícii je opäť pre operačný systém Windows, Linux, aj Macintosh. Užitočné materiály pre prácu so štatistickým programom 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.: Stručný úvod do R. (PDF súbor) Scott, T.: An Introduction to R (PDF súbor) De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808) Štatisticky software R bude počas výuky predmetu NMSA230 priamo používaný a jeho používanie je vyžadovane aj od študentov. V prípade štandardnej prezenčnej výuky môžu študenti využívať počítače, ktoré sú k dispozícii v posluchárni K4 (program R a RStudio sú nauinštalované na všetkých počítačoch), ale v prípade distančnej online výuky je nutné, aby každý študent mal program R nainštalovaný a plne funkčný. V prípade potreby doinštalovať dodatočnú knižnicu je potrebné, aby bol počítač pripojený na internet. Výuka počas semestru bude prebiehať formou predpripravených R skriptov s rôznymi praktickými ukážkami, jednoduchými príkladmi a úlohami, ktoré majú primárne slúžiť ako inštruktáž a ilustrácia fungovania programu R (jednak konkrétna implementácia jednotlivých príkazov, ale aj práca s datami a príprava podkladov na tvorbu vyzkumnej správy). Jednotlivé ukážky a príklady sú pripravené formou copy/paste zdrojového R kódu, ktorý stačí nakopírovať do príkazového riadku (R konzoly) a spustiť pomocou tlačítka ENTER. Na rozdiel od pasívneho kopírovania uvedených zdrojových kódov do Rka sa ale od študentov očakáva aj aktívne zapojenie sa formou porozumenia a pochopenia daného príkazu, prípravy vlastného zdrojového kódu a samostatného riešenia zadaných úloh. V prípade potreby/záujimu je možné využiť podrobnejšieho “sprievodcu” na cestu základmi Rka: Hrátky s R (autor: doc. Arnošt Komárek) – potrebné sú ale dátové súbory: auta2004.dat, auta2004.csv a auta2004.xls) 1. Základné matematické operácie v R Najzákladnejšou a najjednoduchšou úlohou programu R je využitie R konzoly pre účely jednoduchej kalkulačky. Vyskúšajste do R konzoly potupne zadať nasledujúce riadky a zakaždým príkaz potvrdiť klávesou ENTER. 2 + (4 * 5)^2 ## za znakom '#' nasleduje komentar, ktory R ignoruje pri vyhodnoteni ignoruje 1:5 * 2 sqrt(2) + sin(pi/2) Jednotlivé príkazy v programe R sú intuitivné (pozri napr. goniometrické funkcie pomocou príkazu ?Trig) a získané výsledky je možné uložiť do vhodne zvoleného objektu, na ktorý sa da následne odkazovať a plnohodnotne daný objekt pri práci ďalej využívať. v1 <- 2 + (4 * 5)^2 ## za '#' nasleduje komentar, ktory R ignoruje v2 <- 1:5 * 2 v3 <- seq(1,5, length = 3) v4 <- rep(4,3) Jednotlivé výsledky je možné spolu kombinovať (aplikovať na ne vhodné matematické oprácie, alebo využiť iné nematematické manipulácie) a opäť uložiť ako nový výsledok (resp. prepísať pôvodný objekt). v1 + v2 v5 <- v3 - 2 v5 <- v5 * 2 Momentálne je v prostredí R uložených niekoľko objektov. Zoznam týchto objektov, resp. výsledkov, získame pomocou príkazu ls(): ls() V prípade potreby môžeme objekty vymažeme buď selektívne, alebo hromadne pomocou jedného z príkazov rm(v1) ## vymaze objekt v1 rm(list = ls()) ## vymaze vsetky aktivne objekty Prácu v programe R je vždy dobré začínať “s čistým štítom”, teda príkazom rm(list = ls()). Názvy jednotlivých objektov by mali byť intuitivné a zároveň dostatočne stručne. Umožní to rýchlejšie vytváranie (písanie) zdrojového kódu a tiež sa tým minimalizuje rizíko chyby (resp. nežiaducej zámeny medzi existujúcimu objektami). Užitočné Program R je objektový nástroj a zvláda aj matematické operacie, ktoré štandardne nie sú definované. Napr. súčet dvoch vektorov, ktoré majú rôznu dĺžku, operácie s maticami, ktoré majú rôzne typy/dimenzie. Je podstatné správne pochopiť, ako program R v takýchto prípadoch funguje. Zabráni sa možným problémom v budúcnosti a v mnohých prípadoch to zase na druhu stranu može zjednodušiť výpočet a ušetriť čas. Pokúste sa odhaliť princíp fungovania takýchto operácii. Napríklad, porovnajte nasleduúce dva príkazy a príslušné výsledky: c(1, 2, 3, 4) * c(1, 2, 3) 1:6 + 1:3 Štandardné operácie v programe R získame pomocou štandardných znamienok +, -, *, \, a symbolu ^ pro mocniny; Komplexnejšie matematické výrazy môžeme formulovať pomocou zátvoriek. V Rku slúžia k tomuto účelu pouze okrúhle zátvorky (). Hranaté zátvorky [] a zátvorky {} maju v programe R vlastnú, špecifickú úlohu; Ku každému príkazu v R je k dispozícii návod, ako príkaz správne funguje - help sa zavolá zadaním samotného príkazu, ktorému predchádza otazník: napr. ?c zobrazí help pre funkciu c(), analogicky ?seq zobrazí help k funkcii seq(). Help obsahuje popis samotnej funkcie, vysvetlenie implementácie a tiež niekoľko ilustratívnych príkladov a odkazov na dalšie príbuzne funkcie, alebo funkcie podobného charakteru. Podrobný návod je k dispozícii aj na internete a v prípade problému stačí použiť google a riešenie sa určite ľahko nájde; 2. Práca s vektormi a maticami V prípade, že máme v objektoch uložené nejaké hodnoty (vektor, maticu, pole, list, atď.), je možné sa v Rku odkazovať aj na jednotlivé elementy v daných objektoch, používať ich v dalších výpočtoch, alebo s nimi manipulovať. K tomuto účelu práve slúžia hranaté zátvorky []. v3[1] ### zobrazi prvu zlozku vektora v3 v3[-c(2,3)] ### rovnaky vystup, vektor v3 bez druhej a tretej zlozky A to isté platí aj pre matice (a analogicky aj pre pole). m1 <- rbind(v3, v3 * 2, v3 * 3) m1[1:2,1:2] ## elementy z prveho a druheho riadku a stlpca matice m1 Príkazy je možné ľubovoľne kombinovať a vytvárať nové, zložitejšie objekty. m2 <- rbind(cbind(m1[1:2,1:2], c(0,0)), seq(10,11, length = 3)) l1 <- list(m1, m2) ## list, ktory obsahuje dve rozne matice, m1 a m2 l1[[2]][1,] # vypise na obrazovku prvy riadok matice m2 Objekty rôzneho typu môžu byť pohromade uložené napr. v objekte typu list() li <- list(v3, m1) ## list, ktory obsahuje vektor v3 a maticu m1 Na jednotlivé elementy v liste sa odkazuje následovne: li[[1]] ## zobrazi prvy element listu, vektor v3 li[[2]] ## zobrazi druhy element listu, maticu m1 Vždy je vhodné voliť kompaktnú a čo najjednoduchšiu formu zápisu. Pre porovnanie, všetky nasledujúce príkazy vytvoria vo výsledku jednotkovú maticu typu 4x4, ale najkompaktnejši a najstručnejší je pouze posledný zápis. Vyžaduje ale znalosť príkazu diag() a jeho spávnu implementáciu. m1 <- rbind(c(1,0,0,0), c(0,1,0,0), c(0,0,1,0), c(0,0,0,1)) m2 <- cbind(c(1,0,0,0), c(0,1,0,0), c(0,0,1,0), c(0,0,0,1)) m3 <- matrix(rep(0,16), nrow = 4) diag(m3) <- 1 m4 <- NULL for (i in 1:4){ m4 <- rbind(m4, c(rep(0, i - 1), 1, rep(0, 4 - i))) } m5 <- diag(1, 4, 4) ### resp. s vyuzitim defaultneho nastaveni by stacil prikaz `diag(4)` Užitočné V programe R existuje niekoľko prikazov, ktoré vytvoria vektor: napr. c(), seq(), rep(), replicate() a pod.; Podobne existuje niekoľko príkazov, ktoré vytvoria maticu: napr. cbind(), rbind(), matrix(), diag(), upper.tri(), lower.tri(), data.frame() a pod.; Užitočné príkazy pre vytváranie a pracovanie s vektormi a maticami sú as.numeric(), as.matrix(), array(). Pre prácu s objektom typu list() sú užitočné príkazy unlist(), as.list(), pairlist() a pod. Pre maticové násobenie je nutné použiť operátor %*%, ktorý výnasobi dve matice (ak majú správne rozmery). Vhodnosť rozmerov sa jednoducho skontroluje pomocou príkazu dim() (užitočne je poznať tiež príkazy nrow(), ncol(), alebo dimnames()). Inverznú maticu lze získať príkazom solve(), vlastne čísla/vektory matice pomocou príkazu eigen(), rozklady pomocou príkazov eiv(), svd(), qr(), alebo chol(). Štandardne binárne operátory +, -, *, /, ^ fungujú ‘po zložkách’ - toto ‘fungovanie’ je ale R-špecifické - porovnajte nasledujúce výstupy:`` 1:4 + 1 ### sucet po zlozkach 1:4 * 1:2 ### nasobenie po zlozkach v paroch (dlzka jedneho vektora je celociselny nasobok dlzky druheho) 1:4 / 1:4 ### nasobenie v paroch (rovnake dlzky oboch vektorov) 1:5 + 1:2 ## Error - preco? matrix(1:9, nrow = 3) + 1:3 ### ako je ziskany vysledok? matrix(1:9, nrow = 3, byrow = 3) + 1:3 ### ako je ziskany vysledok teraz? matrix(1:9, nrow = 3) + 1:4 ### Error - preco? 3. Náhodné data - generátory náhodných čísel Program R je primárne štatistický nástroj - umožňuje preto aj prácu s náhodnymi hodnotami. Implementované sú rôzne generátory (pseudo) náhodných čísel. Pre rôzne pravdepodobnostné rozdelenia sú k dispozícii príslušne generátory, ktoré simulujú náhodné hodnoty z daného pravdepodobnostného rozdelenia. Diskrétne rozdelenia: rbinom(), rpois(), rgeom(), a ďalšie; Spojité rozdelenia: runif(), rexp(), rnorm(), a ďalšie; Pre správne použitie geneátorov je potrebné preštudovať prislušný návod - napr. ?rnorm zobrazí help pre generátor náhodných čísel s normálnym rozdelením. Je vhodné si tiež všimnúť, že existujú podobné príkazy, označené ako pnorm(), dnorm() a rnorm(), ktoré postupne značia distribučnú funkciu daného rozdelenia, hustotu a kvantilovú funkciu. Užitočné Je nutné dôsledne si preštudovať implementáciu jednotlivých generátorov a funkcií, ktoré súvisia s pravdepodobnostnými rozdeleniami - je totiž pomerne časté, že parametre, ktoré špecifikujú nejaké rozdelenie, majú iný význam. Napríklad: rnorm(10, 0, 2) ## simuluje 10 pozorovani z N(0, 4) - sigma^2 = 4 rexp(10, 4) ## simuluje 10 pozorovani z Exp(lambda = 4), kde EX = 1/4 V prípade použítia nejakého generátora náhodných čísel je vhodné predom nastaviť tzv. seed. Umožní to následnu replikáciu generovania s garanciou rovnakých výsledkov. set.seed(1234) Výsledok použitia generátora náhodných čisel môže byť napríklad vektor, alebo matica (každý stĺpec ma stejný počet hodnôt), alebo list (umožňuje rôzny počet hodnôt v jednotlivých elemntoch). set.seed(1234) g1 <- rnorm(20, 160, 20) ## napr. vyska 20 nahodnych lidi s prumernou vyskou 160 cm a rozptylom 400 g2 <- cbind(g1, rexp(20, 1), rbinom(20, size = 1, prob = 1/2)) ## matica g3 <- list(g1, rexp(10, 1), rbinom(50, size = 1, prob = 1/2)) ## list Generátory náhodných čísel sú častym nástrojom používaným v pravdepodobnosti a štatistike, pre overenie teoretických vlastnosti rôznych modelov. Jednoduchý príklad na ilustráciu: Pre náhodný výber \(X_1, \dots, X_n\) z normálneho rozdelenia \(N(\mu, \sigma^2)\) vieme, že 95 % interval spoľahlivosti pre neznámy parameter \(\mu\) je \[ \Big(\overline{X}_n - u_{1 - \alpha/2} \frac{\sigma}{\sqrt{n}}, \overline{X}_n + u_{1 - \alpha/2} \frac{\sigma}{\sqrt{n}}\Big), \] kde \(u_{1 - \alpha/2} = 1.959964\) (príkaz qnorm(0.975)) je kvantil normálneho rozdelenia pre \(\alpha = 0.05\) (t.j., 95% interval spoľahlivosti). Tento interval teda pokrýva neznámu hodnotu parametru \(\mu \in \mathbb{R}\) s pravdepodobnosťou 0.95. Pomocou generátora náhodných čísel ověříme, že to tak naozaj je. Inými slovami, ak budeme experiment opakovať/generovať dostatočne veľakrát, mali by sme sledovať cca 95 % prípadov, kedy skonštruovaný interval pokryje neznámu hodnotu parametru \(\mu \in \mathbb{R}\) a cca 5 % prípadov, keď sa tak nestane. mu <- 0 ### parameter strednej hodnoty sigma <- 1 ### parameter smerodatnej chyby n <- 1000 ### rozsah nahodneho vyberu B <- 100 ### pocet opakovani experimentu plot(0,0, xlim = c(0,B), ylim = c(mu - 8 * sigma/sqrt(n), mu + 8 * sigma/sqrt(n)), xlab = "Opakovanie", ylab = "Interval spolahlivosti", pch = "") abline(mu, 0, col = "grey", lwd = 2) set.seed(1234) empirickePokryti <- 0 for (opakovanie in 1:B){ vyber <- rnorm(n, mu, sigma) interval <- c(mean(vyber) - qnorm(0.975) * sigma/sqrt(n), mean(vyber) + qnorm(0.975) * sigma/sqrt(n)) if (data.table::between(mu, interval[1], interval[2])){ empirickePokryti <- empirickePokryti + 1 lines(interval ~ c(opakovanie,opakovanie), col = "green") } else { lines(interval ~ c(opakovanie,opakovanie), col = "red") } } text(0, mu + 8 * sigma/sqrt(n), paste("Empirické pokrytie neznámej strednej hodnoty: ", round(100 * empirickePokryti/B, digits = 2), "%", sep = ""), col = "red", pos = 4) Bez nutnosti porozumieť samotnému zdrojovému kódu, pokúste sa aspoň intuitivne vytušiť princíp jeho fungovania. Pomocou nastavenia rôzných hodnôt pre parametere \(\mu \in \mathbb{R}\) a \(\sigma > 0\) a \(n, B \in \mathbb{N}\) experiment zopakujťe a porovnajte výsledky. Pokúste sa vytvoriť analogický experiment k dokázaniu nejakej inej vlastnosti (napr. konzistencia odhadu, a podobne). 4. Reálne data a vlastné datové súbory V programe R je k dispozícii niekoľko vzorových ilustračných datových súborov. Ich zoznam je možne zobraziť zavolaním príkazu data(). Konkrétny datovy súbor je následne zobrazený zavolaním mena prislušného súboru z príkazoveho riadku v R - napr. Orange, alebo mtcars. Vo väčšíne prípadov sa jedná o objekt typu data.frame - resp. špeciálny typ matice (pozri help ?data.frame). Na rozdiel od klasickej matice umožňuje type data.frame()použiť rôzne datové typy (integer, boolean, character, etc.) v jednotlivých stĺpcoch - premenných. Je tak možné vytvoriť datový súbor, ktorý bude obsahovať premenné rôznych typov - spojité veličiny, disrétne hodnoty, ale aj faktory, alebo štandardný popis vo forme textu. stlpec1 <- seq(1:5) stlpec2 <- c(rep("male", 3), rep("female", 2)) stlpec3 <- c("a", "b", "c", "d", "e") ### nie je mozne vytvorit maticu cbind(stlpec1, stlpec2, stlpec3) ### ale je mozne vytvorit data.frame data.frame(stlpec1, stlpec2, stlpec3) S ľubovolným datovým súborom sa následne pracuje ako so štandardnou maticou - môžeme sa odkazovať ne jednotlivé elementy, celé stĺpce (väčšinou jednotlivé premenné), alebo konkrétne řádky (väčšinou jednotlivé pozorovania). Nasledujúce príkazy dávajú stejné výstupy, spôsob implementácie príkazu je ale zakaždým iný. ´ head(Orange) ### zaujima nas pouze promena circumference (treti sloupec) attach(Orange) circumference Orange$circumference Orange[,3] K načítaniu vlastných datových súborov slúži niekoľko príkazov. V programe R je technicky možné načítať vpodstate akýkoľvek datovy súbor (txt, csv, xls, ale aj jpg, png, mp3 a mnoho ďalších). Najčastejšie používané sú read.table() a read.csv(), ktoré vyžaduju súbor uložený vo formáte txt resp. csv. Dodatočné parametre je možné využiť pri volani príkazu - tieto parametre upresňujú format súboru a spôsob načítania do R. Viac podrobnosti v príslušnom helpe ?read.table prípadne ?read.csv. Analogicky je možné stejne súbory pomocou programu R uložiť a zapísať konkrétny subor na disk. Užitočné Program R umožňuje priamo načítať aj online zdroje. Potrebné je iba špecifikovať link a použiť vhodny príkaz. Napríklad: test.data <- read.table("http://www.karlin.mff.cuni.cz/~maciak/NMSA407/NMSA407-1617-HW1.txt", header=T) Použijte google a help v programe R a pokúste sa načítať do prostredia R napr. obrázok (jpg súbor). Podívajte sa, čo ste do Rka načítali (resp. aký objekt je v Rku po načítaní uložený) 5. Práca s datami a s textom (regular expressions) Program R je prispôsobený aj na manipuláciu s datami a prácu s textom, čo môže byť užitočné hlavne v prípadoch, keď je nutné nejaký dátový súbor predpripraviť k samotnej analýze. Program R pracuje s tzv. regular expressions, čím sa rozsah môžnosti používania programu R na prácu textu výrazne rožiruje. Zaujímavý návod, rôzné ukážky a praktické rady, ako pracovať s textom v Rku, je k dispozícii na stránke http://gastonsanchez.com/. Reálne data, ktoré dostane väčšinou štatistik k spracovaniu/vyhodnoteniu, obsahujú často rôzne chyby, preklepy, alebo nesprávne uvedené hodnoty. Základnym predpokladom pre korektne spracovanie dat je nikdy nezasahovať do samotného datového súboru. Priame zásahy (napr. opravy v excel súbore) sú väčšínou nezaznamenané a spätne nevystopovateľné a v prípade následného updatu datového súboru často nie je možné spraviť rovnaké zmeny (úpravy) opakovane. Akákoľvek manipulácia s datami by mala striktne prebiehať vramci samotného R scriptu, ktorý v prvej časti pripravi data k analýze a v druhej časti prevedie samotnu analýzu. Nasledujúci vektor kóduje pouze dve rôzne pohlavia, program R ale rozlíšuje 13 rôznych hodnôt (predpokladame, že 0 značí muža a 1 označuje ženu). sex <- c("M", "m", "f", "Female", "Zena", "Z", "male", "Male", "Muz", "muz", "Zena", "0", "1") table(sex) Takýto vektor je nutné korektne upraviť. Namiesto manuálneho zásahu je výhodnejšie využiť súbor R príkazov, ktoré zmenu prevedu automaticky. sexNew <- sex sexNew[grep("[fFZ1]", sex)] <- "F" sexNew[grep("[Mm0]", sexNew)] <- "M" Nový vektor sexNew naďalej kóduje pouze dva rôzne pohlavia, ale tentokrát aj samotný program R rozumie, že sa jedna pouze o dve rôzne kategórie. Podobný problém, ktorý vyžaduje prípravu dat, je napríklad vek respondentov uvádzaný vo forme dátumu narodenia. age <- c("02/1990", "02/1980", "06/2000", "05/1976", "01/1991") age <- 2021 - as.numeric(unlist(strsplit(age, split = "/"))[2 * (1:length(age))]) Užitočné Zaujímavé príkazy na prácu s textom sú napr.: paste(), strsplit(), grep(), sub(), match(), pmatch() a mnoho jiných. Pre prácu s časovými hodnotami je výhodné poznať funkcie date(), as.Date() a pod.. Pomocou helpu k jednotlivých Rkovým funkciám sa dá jednoducho zistiť, čo dané funkcie robia a ako sú implementované. V helpe sú k dispozícii aj jednoduhé ilustratívne príklady. Domáca úloha (Deadline: 2. seminár | St: 19.10.2022 V programe R si pripravte jednoduchý skript, ktorý vytvorí jednoduchý datovy súbor. Data budú obsahovať aspoň 30 pozorovaní (riadkov) a aspoň 6 premenných. Premyslíte si jednoduchý (hypotetický) experiment a k tomuto experimentu v Rku pripravte (nasimulujte) data. Aspoň jedná veličina musí byť faktorová s aspoň dvoma rôznymi úrovňami. Aspoň dve veličiny musia byť generované náhodným generátorom pre spojité rozdelnie. Aspoň dve veličiny musia byť generované náhodným generátorom pre diskrétne rozdelenie. Aspoň jedna veličina bude vytvorená nejakým vhodným využitím predchádzajúcich veličín (manipuláciou s predchadzajúcimu R objektmi). Data uložte ako data.frame() objekt a výsledny objekt si pomocou príkazu write.table() uložte na disk k ďalšiemu použitiu.