NMSA230 - Úvod do programování v R

Zimný semester 2021/2022 | Cvičenie 4 | St 01/12/21



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ž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 naučiť sa niektore základy vhodné pre reportovanie výsledkov štatistickej analýzy pomocou R knižnice Sweave. Táto knižnica umožnuje 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 - 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).

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