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")
}}