NMSA230 - Úvod do programování v R

Zimný semester 2019/2020 | Cvičenie 4 | St 20/11/19



IV. Štatistická analýza v R a reportovanie cez 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ž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ú 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 možne 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 interfacov 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.: Velmi stručný úvod do R. (PDF súbor)
  • De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808)


R knižnica Sweave umožnuje vzájomne prepojenie programov R a LaTeX. Úvodu a základom LaTeX-u 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 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 - 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 akonzistentne implementované. Tieto nedostatky je nutné vždy najprv analyzovať a v prípade potreby ich opravit. Samotné opravy a akékoľvek zásahy do datového súboru musia byť spätne vystopovateľne a v prípade nutnosti taktiež 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).

  • Exploratorná 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ú vhodné 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 a nevyužíva žiadne komplexnejšie pravdepodobnostné/štatistické modely, ale mala by poskytovať informáciu nutnú pre následne rozhodnutie ohľadom modelu, ktorý bude aplikovaný v ďalšom kroku.

  • Štatistická analýza - v určitom zmysle sa jedná o najpodstatnejší krok, ktorý ale viacmenej 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 modely (napr. rôzne pravdepodobnostné rozdelenia, testy hypotéz, regresné modely, a pod.).

  • Záverečný report a interpretácia výsledkov - pre samotneho 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 výsledkov a ich interpretácia. 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 a to tak, aby v prípade potreby bolo možne experiment a analýzu dat z daneho experimentu 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 dokonca existuje knižnica (Sweave), ktorá umožnuje priame prepojenie Rka a TeX-u (viď niekoľko detailov nižšie´).



V nasledujúcej časti budeme uvažovať data z konrétneho experimentu uskutočneného na CZU v Prahe. Data, resp. datové súbory sú vo formáte, ako ich pripravili odborníci 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 zkúmaniu pohybovej aktivity a agresivity sumca v závislosti na 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ť 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čítane popisné charakteristiky vhodne doplňte obrázkami.
  • 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.
  • Pokúste sa navrhnúť nejaký vhodnžý pravdepodobnostný/štatistický model, ktoý by dokázal ponúknuť relevantnú (štatistickú) odpoveď na otázku, či pohybová aktivita sumca závisí, resp. ako závisí na hladine jeho stresových hormónov, resp. na prítomnosti rôzných látok v krvy.
  • Aké problémy s experimentom, jeho návrhom a jeho uskutočnením by ste označili ako za najvážnejšie nedostatky, resp. zásadné nedostatky?



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)

cat("\\begin{table}[ht]\n", file = "tabulka.tex", append = T)
cat("\\centering", file = "tabulka.tex", append = T)
cat("\\begin{tabular}{lcccc}\n", file = "tabulka.tex", append = T)
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 & ", round(mean(agresivita$A_biting[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), " {\\it(", round(sd(agresivita$A_biting[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), ")} &" ,
                           round(mean(agresivita$A_biting[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), " {\\it(", round(sd(agresivita$A_biting[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), ")} &" ,
                          round(mean(agresivita$A_biting[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), " {\\it(", round(sd(agresivita$A_biting[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), ")} &" ,
                          round(mean(agresivita$A_biting[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), " {\\it(", round(sd(agresivita$A_biting[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), ")} \\\\\n", sep = "") , file = "tabulka.tex", append = T)
    
    cat(paste("Naháňanie & ", round(mean(agresivita$A_chasing[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), " {\\it(", round(sd(agresivita$A_chasing[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), ")} &" ,
                           round(mean(agresivita$A_chasing[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), " {\\it(", round(sd(agresivita$A_chasing[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), ")} &" ,
                          round(mean(agresivita$A_chasing[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), " {\\it(", round(sd(agresivita$A_chasing[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), ")} &" ,
                          round(mean(agresivita$A_chasing[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), " {\\it(", round(sd(agresivita$A_chasing[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), ")} \\\\\n", sep = "") , file = "tabulka.tex", append = T)
        
        cat(paste("Bočný útok & ", round(mean(agresivita$A_latdsipl[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), " {\\it(", round(sd(agresivita$A_latdsipl[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), ")} &" ,
                           round(mean(agresivita$A_latdsipl[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), " {\\it(", round(sd(agresivita$A_latdsipl[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), ")} &" ,
                          round(mean(agresivita$A_latdsipl[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), " {\\it(", round(sd(agresivita$A_latdsipl[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), ")} &" ,
                          round(mean(agresivita$A_latdsipl[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), " {\\it(", round(sd(agresivita$A_latdsipl[agresivita$cas3 >= 18 & agresivita$cas3 < 24]), 2), ")} \\\\\n", sep = "") , file = "tabulka.tex", append = T)
            
        cat(paste("Čelný útok & ", round(mean(agresivita$A_frontdispl[agresivita$cas2 >= 0 & agresivita$cas3 < 6]), 2), " {\\it(", round(sd(agresivita$A_frontdispl[agresivita$cas3 >= 0 & agresivita$cas3 < 6]), 2), ")} &" ,
                           round(mean(agresivita$A_frontdispl[agresivita$cas2 >= 6 & agresivita$cas3 < 12]), 2), " {\\it(", round(sd(agresivita$A_frontdispl[agresivita$cas3 >= 6 & agresivita$cas3 < 12]), 2), ")} &" ,
                          round(mean(agresivita$A_frontdispl[agresivita$cas2 >= 12 & agresivita$cas3 < 18]), 2), " {\\it(", round(sd(agresivita$A_frontdispl[agresivita$cas3 >= 12 & agresivita$cas3 < 18]), 2), ")} &" ,
                          round(mean(agresivita$A_frontdispl[agresivita$cas2 >= 18 & agresivita$cas3 < 24]), 2), " {\\it(", round(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)
cat("\\end{table}", 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: 5. cvičenie | St: 04.12.2019)

  • Pomocou knižnice Sweave preneste popisné výberové charakteristiky, ktore ste napočítali k Vášmu nagenerovanému datovému súboru do PDF reportu. Popisné charakteristiky doplňte obrázkami, ktoré ste si pripravili.
  • Do PDF reportu pridajte stručný popis dat (aspoň niekoľko riadkov) a stručnú interpretáciu, resp. popis výberových charakteristík a obrázkov.