NMSA230 - Úvod do programování v R Zimný semester 2022/2023 | Cvičenie 4 | St 16/11/22 • zdrojový Rmd súbor IV. Štatistická analýza v R a reportovanie pomocou Sweave 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 priamou podporou pre OS Windows, Linux aj Macintosh. Základnú inštaláciu programu R je možné jednoducho rozšíriť pomocou dodatočných knižníc (balíčkov), ktoré sú k dispozícii v online repozitároch (zoznam hlavných repozitárov je na adrese https://cran.r-project.org/mirrors.html). Jednotlivé R knižnice sú tvorené samotnými užívateľmi softwaru R a ich správne fungovanie nie je garantované - je preto namieste určitá opatrnosť a hlavne aktívne premýšľanie pri ich používaní. Pre užívateľov programu R sú k dispozícii aj rôzne grafické rozhrania, ktoré je po inštalácii 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. 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.: 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) Cieľom tohto semináru je ukázať si niektoré užitočné nástroje (v programe R) vhodné pre reportovanie výsledkov štatistickej analýzy – hlavne sa bude jednať o knižnicu Sweave. Táto knižnica umožňuje vzájomne prepojenie programov R a LaTeX. Čo sa týka samotného programu LaTeX-u, tomu je na MFF venovaný samostatný seminár. Užitočné informácie je možné násť aj v nasledujúcich dokumentoch, alebo na internete. Užitočné materiály pre prácu s LaTeX-om Tobias Oetiker a kol.: Ne příliš stručný úvod do systému LaTeX 2ε. (PDF súbor) Kolektív AF UPOL: Drsný úvod do LaTeX-u. (PDF súbor) Warbrick, J.: Essential LaTeX ++ (PDF súbor) Olšák, P.: Typografický systém TeX, Konvoj, 1995. Knuth, D.: Computer and Typesetting Series, Vol. A: The TeX Book, Addison Wesley, 1986. Knuth, D.: Computer and Typesetting Series, Vol. B: TeX The Program, Addison Wesley, 1986. 1. Štatistická analýza v programe R Na začiatku každej štatistickej analýzy je vždy nejaký konkrétny datový súbor a s ním spojená základná (výzkumná/vedecká) otázka, na ktorú sa pomocou dat a ich analýzy snažíme najsť odpoveď. Každá štatistická analýza konkrétneho datového súboru vyžaduje niekoľko na seba naväzujích krokov. Celkovým výsledkom týchto niekoľkých krokov je potom formulácia odpovede na predom kladenú otázku, resp. dokument, ktorý výsledky analýzy dostatočne podrobne popisuje, podáva plnohodnotnú odpoveď na kladenú otázku a jednotlivé časti náležite interpretuje spôsobom, ktorý je zrozumiteľný a pochopiteľný aj pre neodborníka – t.j., nematematika/neštatistika. Celkový proces štatistickej analýzy sa dá rozdeliť do niekoľkých samostatných krokov: Príprava datového súboru – je viacmenej nereálne očakávať, že datový súbor, ktorý dostane štatistík k spracovaniu, je ideálne koncipovaný, bez zásadných chýb, preklepov a nedostatkov, resp. jednotlivé hodnoty v datach budú vždy korektne a konzistentne implementované. Tieto nedostatky je nutné najprv analyzovať a v prípade potreby ich opraviť. Samotné opravy a akékoľvek zásahy do datového súboru musia byť korektne zaznamenané, spätne vystopovateľne a v prípade nutnosti taktiež plnohodnotne reprodukovateľné. Výsledkom prvého kroku je datový súbor, ktorý je vhodne pripravený k následnej štatistickej analýze (od správneho formátu nutného k načítaniu softwarom, až po vhodnú reprezentáciu jednotlivých pozorovaní a premenných). Exploratívná analýza – jedná sa o prvý krok štatistickej analýzy dát, ktorý ale výchádza výhradne z popisných (výberových) charakteristík spočítaných z pozorovaní, ktoré sú k dispozícii. Súčasťou exploratívnej analýzy sú textové, tabuľkové, aj grafiké výstupy, ktoré sa vzájomne logicky dopĺňajú a poskytujú ucelený prehľad o vnútornej štruktúre datového súboru, resp. časti relevantnej pre kladenú otázku. Exploratívna analýza nevýchádza z komplexných pravdepodobnostných/štatistických modelov, ale mala by poskytovať informáciu nutnú pre následne rozhodnutie sa ohľadom modelu, ktorý bude následne korektne aplikovaný v ďalšom kroku. Štatistická analýza – v určitom zmysle sa jedná o najpodstatnejší krok, ktorý ale nie je možné uskutočniť korektne, ak nebol vhodne pripravený datový súbor v prvom kroku a ak nemáme dostatočný prehľad o datách a ich vnútornej štruktúre, získaný z exploratívnej analýzy v druhom kroku. Štatistická analýza sa opiera o rôzne pravdepodobnostné a štatistické modely (napr. rôzne pravdepodobnostné rozdelenia, testy hypotéz, regresné modely, a pod.). Záverečný report a interpretácia výsledkov – pre samotného zadávateľa/zákazníka je najpodstatnejšia práve posledná, štvrtá časť, z ktorej sa dozvie odpoveď na otázku položenú v úvode. Jedná sa o reportovanie získaných štatistických výsledkov a ich interpretáciu. Súčasťou je väčšinou aj popis metodológie a celkového postupu pri analýze dat. Idealne je spracovanie v písomnej forme. Popis použitej metodológie musí byť dostatočne podrobný a presný a mal by poskytovať ucelenú informáciu o celkovom experimente, podkladových datach aj ich analýze. V prípade potreby by malo byť možne pomocou tohto popisu experiment a analýzu dat plnohodnotne zreplikovať. Oproti tomu popis výsledkov a ich interpretácia musí byť zmysluplná aj pre človeka, ktorý nie je odborníkom v štatistike/matematike a vhodne by mala byť doplnená o podstatné grafy a obrázky, ktoré použitý model vizualizujú a pomáhajú celkovému pochopeniu. V prípade vypracovania štatistickej analýzy v programe R je ideálnym kandidátom na vytvorenie reportu TeX. V programe R existuje knižnica (Sweave), ktorá umožnuje priame prepojenie Rka a TeX-u (viď niekoľko príkladov nižšie). V nasledujúcej časti budeme uvažovať data z konrétneho experimentu, ktorý sa uskutočnil na CZU v Prahe. Data, resp. datové súbory sú vo formáte, ako ich pripravili na CZU. Cieľom je data spracovať a pripraviť k analýze - zameriame sa pouze na exploratívnu časť. Niektoré zistenia a výsledky budeme následne pomocou knižnice Sweave a TeX-u reportovať priamo do výsledného PDF dokumentu. Príklad 1 | pohybová aktivita a stres u sumcov Ako názorný a pomerné jedoduchý príklad (jednoduchý hlavne z hľadiska nutnej prípravy dat) poslúži experiment, ktorý bol navrhnutý k zisťovaniu pohybovej aktivity a agresivity sumca v závislosti na jeho zafarbení, hladine stresových hormónov a krvného rozboru. Potrebné data sú postupne uložené v troch datových súboroch, ktoré je potrebné načítať samostatne: agresivita <- read.csv("http://www.karlin.mff.cuni.cz/~maciak/NMSA230/agresivita.csv", header=T, sep = " ", dec = ".", stringsAsFactors = TRUE) stres <- read.csv("http://www.karlin.mff.cuni.cz/~maciak/NMSA230/stresoveHormony.csv", header=T, sep = " ", dec = ".", stringsAsFactors = TRUE) krv <- read.csv("http://www.karlin.mff.cuni.cz/~maciak/NMSA230/krvneTesty.csv", header=T, sep = " ", dec = ".", stringsAsFactors = TRUE) Zaujíma nás hlavne celková pohybová aktivita z hľadiska súčtu niektorých konkrétnych premenných a to hlavne v štyroch segmentoch dňa, konkrétne (napr. priemerná) aktivita v čase medzi 0:00 až 6:00, potom 6:00 až 12:00, 12:00 až 18:00 a nakoniec segment 18:00 až 24:00. Informácia o čase je ale v jednotlivých datových súboroch zaznamenaná rôzne. V prvom rade budeme sledovať aktivitu formou binárnej premennej (došlo k nejakej agresívnej aktivite, alebo nedošlo). Za agresívnu aktivitu budeme považovať napr. kusnutie iného jedinca (bitting), naháňanie iného jedinca (chasing), útok z boku (lateral) a útok z predu (frontal). agresivita$total <- agresivita$A_biting + agresivita$A_chasing + agresivita$A_latdsipl + agresivita$A_frontdispl agresivita$total01 <- agresivita$total agresivita$total01[agresivita$total > 1] <- 1 Následne je potrebné zmeniť čas zaznamenaný v jednotlivých datových säboroch tak, aby bol priamo porovnateľný medzi súbormi. agresivita$cas3 <- sapply(strsplit(as.character(agresivita$cas2),":"), function(x) { x <- as.numeric(x) if (length(x) == 2){return(x[1]+x[2]/60)} if (length(x) == 3){x[1]+x[2]/60 + x[3]/3600} } ) krv$cas2 <- sapply(strsplit(as.character(krv$cas),":"), function(x) { x <- as.numeric(x) if (length(x) == 2){return(x[1]+x[2]/60)} if (length(x) == 3){x[1]+x[2]/60 + x[3]/3600} } ) Následne môžeme informáciu o čase segmentovať podla požadovaných segmentov vrámci daného dňa. agresivita$segment <-1*(agresivita$cas3 > 0 & agresivita$cas3 <=6) + 2*(agresivita$cas3 > 6 & agresivita$cas3 <= 12) + 3*(agresivita$cas3 > 12 & agresivita$cas3 <= 18) + 4*(agresivita$cas3 > 18 & agresivita$cas3 <=24) stres$segment <- 1*(stres$cas == "0:00" | stres$cas == "4:00") + 2 * (stres$cas == "6:50" | stres$cas == "8:00") + 3 * (stres$cas == "12:00" | stres$cas == "16:00") + 4 * (stres$cas == "18:50" | stres$cas == "20:00") krv$segment <- 1*(krv$cas== "0:00" | krv$cas == "4:00") + 2 * (krv$cas == "6:50" | krv$cas == "8:00") + 3 * (krv$cas == "12:00" | krv$cas == "16:00") + 4 * (krv$cas == "18:50" | krv$cas == "20:00") S takto upravenými datovými súbormi sa už da rozumne pracovať a dokážeme dokonca získať úrčitu grafickú vizualizáciu štruktúry v datach. Zameriame sa pouze na tri pohybové aktivity (zdroj, pokus a total) a pouze tri stresové ukazovatele (SOD v nameraný v ziabrach, SOD zaznamenaný v mozgu a CAT zaznamenaný opäť v žiabrach). vars <- c("zdroj.W_L", "W_L_pokus", "total") varStress <- c("SODzabra", "SODmozek", "CATzabra") par(mfrow = c(3,3)) for (c1 in vars){ for (c2 in varStress){ par(mar = c(5, 4, 4, 4) + 0.3) plot(jitter(agresivita[,c1]) ~ I(jitter(agresivita[,"segment"], factor = 0.5) - 0.1), data = agresivita, pch = 21, bg = "lightblue", cex = 0.8, xlim = c(0.5, 4.5), xlab = "Time segments", ylab = "") Means <- aggregate(agresivita[,c1], by = list(Segment = agresivita$segment), mean)[,2] lines(Means ~ c(1,2,3,4), col = "blue", lwd = 2) mtext(c1 ,side=2,line=3, col = "black") axis(2,col="black",col.axis="black") par(new=TRUE) Means2 <- aggregate(stres[,c2], by = list(Segment = stres$segment), mean)[,2] plot(Means2 ~ c(1,2,3,4), axes = FALSE, bty = "n", xlab = "", ylab = "", type = "l", col ="red", ylim = c(min(Means2), max(Means2)), lwd = 2) axis(4,col="red" ,col.axis= "red") mtext(c2,side=4,line=3, col = "red") }} Samostatne Použijte uvedené datové súbory a pokúste sa o nasledujúce: Spočítajte niektoré ďalšie popisné charakteristiky, ktoré bude možne využiť k následnému získaniu informácie a vzájomnej súvislosti medzi pohybovou aktivitou sumca a jeho stresovými a krvnými testami. Spočítané popisné charakteristiky vhodne doplňte obrázkami – vizualizujte hlavne faktory, ktoré pomôžu zodpovedať otázku, či agresivita sumcov závisí na ich sfarbení, stresových enzýmoch a jednotlivých hodnotách krvných testov. Pokúste sa vlastnými slovami výsledky exploratívnej analýzy interpretovať a vyjadriť svoje očakávania, ako by mohla vyzerať, resp. dopadnǔť štatistické analýza – teda aplikácia nejakého konkrétneho pravdepodobnostného/štatistického modelu. Užitočné V prípade, že používate pre zápis český, alebo slovenský text (tzn. znaky s diakritikou), je dôležité vedieť, že program R defaultne nevyexportuje obrázky (pomocou príkazu pdf()) do PDF s tzv. embedovanými fontami (embedded fonts). To môže mať za následok nesprávne vykreslenie obrázkov a hlavne príslušných popiskov. Jedna z možnosti, ako pomocou programu R vyexportovať obrázky tak, aby boli fonty embedované, je použiť príkaz cairo_pdf(). Viac informácii je možné nájsť napr. na tejto stránke: https://hansjoerg.me/2018/02/15/font-embedding-for-latex-and-r-users/ 2. Reportovanie pomocou R knižnice Sweave Sweave je balíček pre štatistický program R, ktorý umožňuje integráciu Rkového kódu priamo v LaTeX-ovom kóde. Umožňuje tak vytvárať dynamické reporty, ktoré je možne jednoducho a priamočiaro updatovať v prípade, že dôjde k zmene v Rkových výpočtoch, alebo podkladovom datovom súbore. Jednoduchý návod na používanie Sweave je na webovej stránke https://support.rstudio.com/ Povedane stručne a veľmi jednoducho, Sweave je knižnica a tiež funkcia v programe R (teda Sweave()), ktorá integruje Rko a LaTeX. Základom je vytvorenie klasického tex súboru, ktorý je ale uložený (predbežne, kým ho program R zkompiluje) vo formáte *.Rnw. Do tohto suboru je možné okrem štandardného tex kódu vkládať aj vhodne označený Rkový kód. Následne volaním funkcie Sweave("nazov_suboru.Rnw") program R zkompiluje celý dokument, Rkový kód vyhodnotí a na dané miesto doplní výsledok získaný spracovaním daného Rkového kódu. Podkladový tex kód R program ignoruje a výstupom je *.tex súbor, s ktorým je ďalej možné pracovať ako s klasickými tex súborom. V prípade, že použijeme k práci user-friendly interface (asi jedna z najlepších možnosti je RStudio) https://www.rstudio.com tak je možné okamžite vytvárať PDF súbory priamo z Rnw súborov a to pouze stlačením príslučného tlačítka v menu RStudia (program R a LaTeX prebehnú na pozadí a výsledkom bude finálny PDF report). Podrobnejší návod na prácu so Sweave napr. tu: http://gosset.wharton.upenn.edu/ alebo v klasickom helpe: ´ help("Sweave", package="utils") Pre ilustráciu, ako vkládať Rkový kód do Rnw súboru, slúži táto malá ukážka: Download Rnw súbor (Sweave | UTF 8 kódovanie) a výsledný PDF súbor Download PDF súbor Následná kompilácia a teda aj automatické vytvorenie finálneho PDF súboru sa dosiahne pomocou stlačenia tlačítka Compile PDF (v interface RStudio) - viď obrázok nižšie. Samostatne Stiahnite si podkladový Rnw súbor a otvorte ho v programe RStudio. Súbor najrpv zkompilujte (pomocou tlačítka Compile PDF v menu) a následne súbor doplňte o vlastnú čast R-kového aj LaTeX-ového zdrojového kódu. Pridajte aspoň jeden obrázok a jednu tabuľku (s príslušnym popiskom tak, aby bolo možné sa na obrázok a tabuľku v texte následne odkazovať) a finálny Rnw súbor opäť skompilujte a vytvorte výsledný PDF súbor. Podívajte sa na podkladovú tex súbor, ktorý vznikol použitim príkazu Sweave() z R knižnice Sweave. 3. R to LaTeX Niekedy je užitočné nevytvárať pomocou programu R celý zdrojový kód pre LaTeX, ale vytvoriť pouze určitú menšiu časť (napr. tabuľky), ktorú do LaTeX-ového zdrojového kódu pouze prídáme pomocou príkazu \input{}. V programe R existuje celá řada rôznych balíčkov, ktoré umožňuju vytvárať tabuľky, časti textu, prípadne celé PDF subory. Väčšinou sú ale v niečom vždy príliš obmedzujúce a ponúkaju len určitú škálu možných nastavení a adaptácii, ktoré nemusia vždy stačiť, prípadne vyhovovať naším požiadavkám. V takom prípade je výhodne naučiť sa pracovať s R-kovou funkciou cat(), ktorá dokáže (okrem mnohých iných vecí) vytvoriť aj podkladový LaTeX-ový súbor (tex súbor) priamo ušitý na mieru, presne podľa naších preferencii. Konkrétne fungovanie a spôsob implementácie zistite pomocou helpu: ?cat Nasledujúci R-kový kod vyprodukuje LaTeX-ový súbor (tex súbor) tabulka.tex, ktorý následne môžeme vložiť do LaTeX-ového zdrojového súboru. agresivita$cas2 <- as.numeric(agresivita$cas2) f <- function(x, d = 2){### stejny pocet 'd' desetinnych mist return(format(round(x, d), nsmall = d)) } cat("\\begin{tabular}{lcccc}\n", file = "tabulka.tex", append = F) cat("\\hline\n", file = "tabulka.tex", append = T) cat("& \\textbf{00:00 - 06:00} & \\textbf{06:00 -- 12:00} & \\textbf{12:00 -- 18:00} & \\textbf{18:00 -- 24:00}\\\\\n", file = "tabulka.tex", append = T) cat("\\hline\n", file = "tabulka.tex", append = T) cat(paste("Pohryzenie & ", f(mean(agresivita$A_biting[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 3), " {\\it(", f(sd(agresivita$A_biting[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), ")} &" , f(mean(agresivita$A_biting[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 3), " {\\it(", f(sd(agresivita$A_biting[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), ")} &" , f(mean(agresivita$A_biting[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 3), " {\\it(", f(sd(agresivita$A_biting[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), ")} &" , f(mean(agresivita$A_biting[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 3), " {\\it(", f(sd(agresivita$A_biting[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), ")} \\\\\n", sep = "") , file = "tabulka.tex", append = T) cat(paste("Naháňanie & ", f(mean(agresivita$A_chasing[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 3), " {\\it(", f(sd(agresivita$A_chasing[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), ")} &" , f(mean(agresivita$A_chasing[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 3), " {\\it(", f(sd(agresivita$A_chasing[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), ")} &" , f(mean(agresivita$A_chasing[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 3), " {\\it(", f(sd(agresivita$A_chasing[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), ")} &" , f(mean(agresivita$A_chasing[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 3), " {\\it(", f(sd(agresivita$A_chasing[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), ")} \\\\\n", sep = "") , file = "tabulka.tex", append = T) cat(paste("Bočný útok & ", f(mean(agresivita$A_latdsipl[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), " {\\it(", f(sd(agresivita$A_latdsipl[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), ")} &" , f(mean(agresivita$A_latdsipl[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), " {\\it(", f(sd(agresivita$A_latdsipl[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), ")} &" , f(mean(agresivita$A_latdsipl[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), " {\\it(", f(sd(agresivita$A_latdsipl[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), ")} &" , f(mean(agresivita$A_latdsipl[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), " {\\it(", f(sd(agresivita$A_latdsipl[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), ")} \\\\\n", sep = "") , file = "tabulka.tex", append = T) cat(paste("Čelný útok & ", f(mean(agresivita$A_frontdispl[agresivita$cas2 >= 0 & agresivita$cas3 < 6]), 2), " {\\it(", f(sd(agresivita$A_frontdispl[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), ")} &" , f(mean(agresivita$A_frontdispl[agresivita$cas2 >= 6 & agresivita$cas3 < 12]), 2), " {\\it(", f(sd(agresivita$A_frontdispl[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), ")} &" , f(mean(agresivita$A_frontdispl[agresivita$cas2 >= 12 & agresivita$cas3 < 18]), 2), " {\\it(", f(sd(agresivita$A_frontdispl[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), ")} &" , f(mean(agresivita$A_frontdispl[agresivita$cas2 >= 18 & agresivita$cas3 < 24]), 2), " {\\it(", f(sd(agresivita$A_frontdispl[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), ")} \\\\", sep = "") , file = "tabulka.tex", append = T) cat("\\hline\n", file = "tabulka.tex", append = T) cat("\\end{tabular}\n", file = "tabulka.tex", append = T) Samostatne Použijte R-kovú funkciu cat(). a pomocou nej vytvorte jednoduchu tabuľku (súbor tabulka.tex), ktorá bude obsahovať základné popisné charakteristiky k nejakému vhodne zvolenému dátovému súboru (napr. ten, ktorý ste si na začiatku semestra zvolili). Pre správne fungovanie príkazu cat() je nutné zvoliť hodnotu dodatočného parametru ako append = T. Vyskúšajte, ako by fungoval príkaz bez tohto dodatočného parametru (resp. s jeho defaultným nastavením append = FALSE). Domáca úloha (Deadline: 6. cvičenie | St: 16.12.2021) Pomocou knižnice Sweave preneste popisné výberové charakteristiky, ktoré ste napočítali k Vášmu nagenerovanému datovému súboru do PDF reportu. Popisné charakteristiky doplňte obrázkami, ktoré ste vytvorili. Doplňte zmysluplný text (napr. popis vygenerovaných dat), stručnú interpretáciu popisných charakteristík, prípadne obrázkov. Výsledný Rnw súbor zkompilujte pomocou knižnice Sweave a vytvorte PDF dokument.