NMSA230 - Úvod do programování v R

Zimný semester 2023/2024 | Cvičenie 5 | Čt 21/12/23

Program R (dostupný pod GNU GPL licenciou) je k dispozícii k stiahnutiu (free of charge) na adrese

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.: 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 piatého semináru NMSA 230 je práca s reálnymi datami, názorná ukážka jednoduchých simulácii a tiež R knižnica Knitr (ktorá slúži na prípravu HTML Markdown-u, ako je napríklad táto HTML stránka).

V. Analýza dat pomocou knižnice Knitr (HTML)

Tématicky tento seminár naväzuje na minulý seminár, ktorého cieľom bolo ukázať a vyskúšať si vzájomne prepojenie a fungovanie programu R s programom LaTeX (pričom toto vzájomné prepojenie a spolupráca oboch programov je užitočná hlavne v prípade vytvárania rôznych PDF reportov s výsledkami štatistickej analýzy).

Analogickým spôsobom je prepojený aj program R so štandardným zdrojovým kódom HTML stránok – a slúži k tomu R knižnica knitr.

1. Knižnica ‘Knitr’

Knitr (GNU General Public License) je knižnica pre štatistický program R, ktorá umožňuje integráciu R-kového zdrojového kódu
do HTML jazyka (v určitom zmysle analogickým spôsobom, ako knižnica Sweave umožňovala integráciu programu R a LaTeXu). Výslednym produktom knižnice Knitr je súbor html (ale využitie knižnice je všeobecnejšie a umožňuje aj kompiláciu do PDF dokumentu, prípadne do .doc súboru), ktorý je možné prezerať pomocou internetového prehliadača. Takto bola vytvorená napríklad aj táto html stránka (zdrojový kód stránky je v tomto súbore).

Inštalácia a inicializácia knižnice pomocou štandardných príkazov:

install.packages("knitr")
library("knitr")

U niektorých grafických interface (napr. RStudio) môže byť už tento balíček nainštalovaný štandardne a volanie príslušnej funkcie je redukované iba na stlačenie príslušného tlačítka Knit v hlavnom menu. Vstupným suborom je podkladový súbor typu Rmd, ktorý obsahuje jednak zdrojový kód pre html stránku (prípadne php zdrojovým kódom) a tiež vhodne označené príkazy pre program R. Pri kompilácii podkladového súboru sú R-kové príkazy najprv vyhodnotené a výsledky sú dopĺnené do zdrojového html kódu.

Niekoľko príkladov využitia knižnice ‘Knitr’ v programe R v podobe podkladových Rmd súborov z predchádzajúcich seminárov (niektoré z nich ale môžu pri kompilácii vyžadovať súbory, ktoré nejsou automatickou súčasťou podkladového Rmd súboru - v takom prípade napr. postači príslušnú časť zdrojového kódu zmazať). Podkladové Rmd súbory je možné otvoriť v programe RStudio (UTF8 kódovanie – v menu programu využiť možnosť Reopen with encoding) a následne skompilovať stlačením tlačítka “Knit” na hlavnej lište (na hlavnom paneli).

podkladový súbor z prvého cvičenia predmetu NMSA230: Rmd súbor (kódovanie UTF8);
podkladový súbor z druhého cvičenia predmetu NMSA230: Rmd súbor (kódovanie UTF8);
podkladový súbor z tretieho cvičenia predmetu NMSA230: Rmd súbor (kódovanie UTF8);
podkladový súbor zo štvrtého cvičenia predmetu NMSA230: Rmd súbor (kódovanie UTF8);
podkladový súbor z piateho cvičenia predmetu NMSA230: Rmd súbor (kódovanie UTF8);

Samostatne

Na rozdiel od knižnice Sweave je zdrojový kód programu R pre knižnicu Knitr oddelený pomocou párových znakov

pričom v úvode je opäť dovolené používať dodatočné premené (napr. eval = F, čo znamená, že príslušná časť R-kového zdrojového kódu bude pri kompilácii ignorovaná).
Využijte podkladové Rmd súbory uvedené vyššie, otvorte ich v programe RStudio a príslušný Rmd súbor zkompilujte (pomocou tlačítka Knit HTML v menu).
Podkladový Rmd súbor doplňte o vlastnú časť R-kového zdrojového kódu. Pridajte minimálne jeden obrázok a finálny Rmd súbor opäť skompilujte a vytvorte výsledný HTML súbor.

2. Štatistická analýza a jednoduché štatistické simulácie

Nasledujúci príkaz načíta reálne data, ktoré zaznamenávajú rôzne informácie o tehotných pacientkách a následných pôrodoch na gynekologickom oddelení v Liptovskom Mikuláši. Data sú v pôvodnom tvare/formáte, ako ich poskytli samotní lekári.

data <- read.csv("http://www.karlin.mff.cuni.cz/~maciak/NMSA230/gynData.csv", header = T, dec = ".", sep = "")

Data obsahujú niekoľko rôznych premenných z ktorých väčšina je momentálne pre naše účely nepodstatná. V prevažnej väčšine sa ale jedná o rôzne lekárske kritéria, výsledky konkrétnych tehotenských testov, čí prípadné výskyty rôznych komplikácii (informácia v tvare pozitívny výsledok vs. negatívný výsledok) a tiež hodnoty testov na tehotensku cukrovku (veličiny glik1 až glik4). Doplnené sú niektoré základné popisné charakteristiky pacientky (vek, výška, hmotnosť, BMI), popisné charakteristiky narodeného dieťaťa (pohlavie, obvod hlavy a hrudníka, dĺžka, váha, systolický a diastolický krvný tlak a ďalšie).

names(data)

##  [1] "nove"         "povodne"      "glik1"        "glik2"        "glik3"       
##  [6] "glik4"        "Vek"          "Parita"       "GestVek"      "PrekoncHm"   
## [11] "PorHm"        "Vyska"        "BMIinic"      "BMIfinal"     "weightGain"  
## [16] "pohlavie"     "obvodhrud"    "obvodHlava"   "dlzka"        "sTK"         
## [21] "dTK"          "Apgar1"       "Apgar5"       "Apgar10"      "cisarskyREX" 
## [26] "cisarskyVEX"  "DYSTOKIO"     "sposobPorodu" "PEE"          "hospit"      
## [31] "velkyPlod"    "malyPlod"     "hmotnost"     "HYPERBILIRUb" "Bilirubin"   
## [36] "ATB"          "POLYCYTEMIA"  "Hematokrit"   "HYPOGLYKEMIA" "Glykemia"    
## [41] "PORANENIE"    "TRANSFUZIE"

Podstatné ale je prehliadnúť si data a uvedomiť si niektoré chyby (napr. hodnota pz je najskôr preklep a jedná sa o pozitívnu hodnotu testu; podobne hodnota nn asi najskôr značí negatívny výsledok):

table(data$cisarskyREX)

## 
## neg  pn poz 
## 491   1 213

table(data$cisarskyVEX)

## 
##     neg poz 
##   1 693  11

table(data$hospit)

## 
## neg  nn poz 
## 569   1 135

Tieto chyby opravíme pomerne jednoducho, nakoľko sa jedná o jednoduché korekcie. V zložitejších prípadoch je ale vhodné pamätať na možnosti, ktoré v programe R máme vďaka jeho schopnosti pracovať s tzv. regularnými výrazmi (napr. príkazy grep(), grepl(), sub(), gsub(), regexpr()`, atď.).

data$cisarskyREX[which(data$cisarskyREX == "pn")] <- "poz"
data$cisarskyREX <- factor(data$cisarskyREX, levels = c("neg", "poz"))

data$cisarskyVEX[which(data$cisarskyVEX == "")] <- "neg"
data$cisarskyVEX <- factor(data$cisarskyVEX, levels = c("neg", "poz"))

data$hospit[which(data$hospit == "nn")] <- "neg"
data$hospit <- factor(data$hospit, levels = c("neg", "poz"))

Z určitých expertných dôvodov je pre lekárov doležite rozlíšiť staré a nové kritéria (prvé dva stĺpce v datach). Zavedieme preto novú premennú, ktorá bude súhrnne popisovať výsledky oboch kritérii:

data$kriteria <- rep("neg", dim(data)[1])
data$kriteria[data$nove == "poz" & data$povodne == "neg"] <- "posNove"
data$kriteria[data$nove == "neg" & data$povodne == "poz"] <- "posPovd"
data$kriteria[data$nove == "poz" & data$povodne == "poz"] <- "posBoth"
data$kriteria <- as.factor(data$kriteria)

table(data$kriteria, data$nove)

##          
##           neg poz
##   neg     611   0
##   posBoth   0  18
##   posNove   0  64
##   posPovd  12   0

table(data$kriteria, data$povodne)

##          
##           neg poz
##   neg     611   0
##   posBoth   0  18
##   posNove  64   0
##   posPovd   0  12

Takto nejak by mohol v praxi vyzerať tzv. ``data processing step’’, ktorý je väčšinou nutnou súčasťou pri každom štatistickom spracovaní a vyhodnotení dat. Úpravy v datach by sa mali vždy robiť pomocou samotného programu a príslušného skritpu, ktorý všetky zásahy v datach zaznamenáva prostredníctvom zdrojového kódu (spätná traktabilita/vystopovatelnosť všetkých zásahov, ktore sa s datami robili).

Samostatne

Podívajte sa na data a pokúste sa spočítať niektoré základné popisné charakteristiky. Tieto charakteristiky doplňte o vhodne zvolené obrázky, ktoré budú popisné charakteristiky vizuálne reprezentovať.
Pokúste sa pomocou popisných charakteristík povedať niečo podstatné o starých a nových kritériach vzhledem na výskyt rôznych komplikácii (rozlišujte pacientky v závislosti na výsledkoch starých a nových kriterii, resp. na novo-zavedenej premennej kriteria).
Pomocou popisných charakteristík novorodenca sa pokúste ponúknuť odpoveď na otázku, či je rozdiel medzi novonarodeným chlapcom a novonarodenou holkou.

V nasledujúcej časti sa zameriame pouze na hmotnosť novonarodených deti (premenná hmotnost). Histogram so zaznamenanými hodnotami vyzera takto (pre lepšiu vizuálizáciu je súčasťou histogramu aj neparametrický odhad hustoty a hustota normálneho rozdelenia so strednou hodnotou rovnou výberovému priemeru a rozptylom rovným výberovému rozptylu):

hist(data$hmotnost, xlab= "Hmotnosť novorodenca", freq = F, ylab= "Relatívny výskyt", main = "", col = "lightblue", breaks = 20)
lines(density(data$hmotnost), col = "red", lwd = 2)
lines(dnorm(seq(1500, 5000, length = 10000), mean(data$hmotnost), sd(data$hmotnost)) ~ seq(1500, 5000, length = 10000), col = "blue", lty = 2)

Červenou krivkou je zobrazený neparametrický odhad hustoty, ktorý svojim tvarom možno pripomína hustotu normálneho rozdelenia s nejakým vhodným parametrom strednej hodnoty a rozptylu.
Modrou krivkou je zobrazená teoretická hustota normálneho rozdelenia, kde jako stredná hodnota a rozptyl sú použite empirické odhady získane z dat.
```
mean(data$hmotnost)
```
```
## [1] 3414.851
```
```
var(data$hmotnost)
```
```
## [1] 237520.5
```

Formálne by sme mohli napríklad otestovať, čí namerané hodnoty výšky novorodencov pochádzajú z normálneho rozdelenia (nulová hypotéza), alebo pochádzajú z nejakého iného rozdelenia (alternatívna hypotéza). Na testovanie takto definovanej nulovej hypotézy máme v teórii pravdepodobnosti a v štatistike niekoľko rôznych testov, ktoré ale nie sú ekvivalentné.

Kolmogorov-Smirnovov test – v programe R implementovány v príkaze ks.test();

ks.test(data$hmotnost,"pnorm",mean=mean(data$hmotnost),sd=sd(data$hmotnost),exact=FALSE)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  data$hmotnost
## D = 0.048552, p-value = 0.07203
## alternative hypothesis: two-sided

Shapiro-Wilkův test – v programe R implementovány v príkaze shapiro.test();

shapiro.test(data$hmotnost)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$hmotnost
## W = 0.99016, p-value = 0.0001164

Jarque-Bera test – v programe R implementovány v knižnici tseries v príkaze jarque.bera.test();

library(tseries)
jarque.bera.test(data$hmotnost)

## 
##  Jarque Bera Test
## 
## data:  data$hmotnost
## X-squared = 21.75, df = 2, p-value = 1.893e-05

$\boldsymbol{\chi}^2$ test dobrej zhody – k testovaniu normality by bolo možné využiť napr. aj $\chi^2$ test dobrej zhody (v rôznej jeho komplexnosti a zložitosti). Pre výberovy priemer váho novorodenca máme $\overline{X}_n \approx 3414.8$
```
mean(data$hmotnost)
```
```
## [1] 3414.851
```
a pre výberovú smerodatnú chybu $S_n \approx 487.3$
```
sd(data$hmotnost)
```
```
## [1] 487.3607
```
Za predpokladu normality by teda štandardizované hodnoty $(X_i - \overline{X}_n)/S_n$, pre $i = 1, \dots, n$ mali mať štandardné normálne rozddelenie $N(0,1)$. Z pravdepodobnostnej teórie vieme, že $P[X < -1] \approx 0.1586$ a tiež $P[X > 1] \approx 0.1586$, kde $X \sim N(0,1)$. Inými slovami, v intervale $[-1,1]$ (t.j. $\pm \sigma$ je približne $68.28~\%$ celkovej pravdepodobnostnej masy). To znamená, že za platnosti nulovej hypotézy o normalite by približne $!5.86~\%$ z hodnôt $(X_i - \overline{X}_n)/S_n$ malo byť menších (rovných), než hodnota $-1$, približne $34.14~\%$ hodnôt by malo byť väčších ako $-1$, ale menších (rovných) ako nula, rovnako cca $34.14~\%$ hodnôt by malo byť väčších ako nula, ale menších (rovných) ako 1 a konečne, zostávajúcich $15.86~\%$ hodnôt by malo byť väčších ako 1.
```
I1 <- sum((data$hmotnost - mean(data$hmotnost))/sd(data$hmotnost) <= -1)
I2 <- sum((data$hmotnost - mean(data$hmotnost))/sd(data$hmotnost) > -1 & (data$hmotnost - mean(data$hmotnost))/sd(data$hmotnost) <= 0 )
I3 <- sum((data$hmotnost - mean(data$hmotnost))/sd(data$hmotnost) > 0 & (data$hmotnost - mean(data$hmotnost))/sd(data$hmotnost) <= 1 )
I4 <- sum((data$hmotnost - mean(data$hmotnost))/sd(data$hmotnost) > 1)
```
Za platnosti nulovej hypotézy o normalite by vektor $(I_1, I_2, I_3, I_4)^\top$ mal mať multinomické rozdelenie s vektorom parametrov (proporcii) $\boldsymbol{p}_0 = (0.1586, 0.3414, 0.3414, 0.1586)^\top$. Pomocou testu dobrej zhody teda môžeme otestovať nulovú hypotézu, že parameter $\boldsymbol{p} = (p_1, p_2, p_3, P_4)^\top$ (t.j., skutočné pravdepodobnosti v multinomickom rozdelení) je rovný hodnotam v $\boldsymbol{p}_0$.
```
p0 <- c(0.1586, 0.3414, 0.3414, 0.1586)
chisq.test(c(I1, I2, I3, I4), p0)
```
```
## 
##  Pearson's Chi-squared test
## 
## data:  c(I1, I2, I3, I4) and p0
## X-squared = 4, df = 3, p-value = 0.2615
```

Rozhodnutie o nulovej a alternatívnej hypotéze je rôzne na základe rôznych testov. Ako si vybrať ten správny test? K odpovedi nám pomôže mala simulačná štúdia, ktorá sa bude zameriavať na sílu jednotlivých testov (to je pravdepodobnosť, že test zamietne nulovú hypotézu, ak nulová hypotéza neplatí). Toto sa práve simuluje následujúci príklad, v ktorom testujeme (použitím rôznych testov), čí náhodný výber pochádza z rovnomerného rozdelenia, pričom vďaka simuláciam vieme, že výber z normálneho rozdelenia nepochádza – pretože data simulujeme z $\chi^2$ rozdelenia.

N <- c(10, 50, 100, 200, 500, 1000, 2000, 5000)

power <- NULL
set.seed(1234)
for (n in N){
  p1 <- p2 <- p3 <- p4 <- NULL
for (i in 1:100){
  
  x <- rchisq(n, df = 50) ### pozorovania simulujeme z chi^2 rozdelenia
  
  I1 <- sum((x - mean(x))/sd(x) <= -1)
  I2 <- sum((x - mean(x))/sd(x) > -1 & (x - mean(x))/sd(x) <= 0 )
  I3 <- sum((x - mean(x))/sd(x) > 0 & (x - mean(x))/sd(x) <= 1 )
  I4 <- sum((x - mean(x))/sd(x) > 1)
  
  p1 <- c(p1, as.numeric(ks.test(x,"pnorm",mean=mean(x),sd=sd(x),exact=FALSE)$p.value < 0.05))
  p2 <- c(p2, as.numeric(shapiro.test(x)$p.value < 0.05))
  p3 <- c(p3, as.numeric(jarque.bera.test(x)$p.value < 0.05))
  p4 <- c(p4, as.numeric(chisq.test(c(I1, I2, I3, I4), p0)$p.value < 0.05))
}
power <- rbind(power, c(n, mean(p1), mean(p2), mean(p3), mean(p4)))
}

plot(power[,2] ~ log10(power[,1]), col = "red", type = "l", pch = 21, bg = "red", xlab = "Logarimus rozsahu náhodného výberu", ylab = "Empirická síla testu", ylim = c(0,1))
lines(power[,3] ~ log10(power[,1]), col = "blue")
lines(power[,4] ~ log10(power[,1]), col = "green")
lines(power[,5] ~ log10(power[,1]), col = "violet")

points(power[,2] ~ log10(power[,1]), pch = 21, bg = "red")
points(power[,3] ~ log10(power[,1]), pch = 21, bg = "blue")
points(power[,4] ~ log10(power[,1]), pch = 21, bg = "green")
points(power[,5] ~ log10(power[,1]), pch = 21, bg = "violet")

abline(1,0, col = "black", lty = 3)

legend(1, 1, legend = c("Kolmogorov-Smirnov test", "Shapiro-Wilk test", "Jarque-Beta test", "Chi-kvadrat test"), col = c("red", "blue", "green", "violet"), lty= rep(1,3))

S9la testu (empirická) by mala pre rastúci rozsah náhodného výberu – na ose x — konvergovať k jednotke… inými slovami, štatistický test je konzistentný a pravdepodobnosťou idúcou k jednej sa rozhodne nulovú hypotézu zamietnúť, ak je táto hypotéza nepravdivá.

Z tohto pohľadu je najsilnejším testom práve Shapiro-Wilkov test. Je potrebné ale pamätať na to, že testy boli spočítane pri konkrétnej alternatíve - skutočné rozdelenie bolo $\chi^2$ rozdelenie s 50 stupňami voľnosti. Jak sa ale porovnanie v zmysle síly testu zmení, ak zmeníme postup a generovať budeme data z $t$ rozdelenia s ťažkými chvostmi (t.j. nízke stupne voľnosti)?

N <- c(10, 50, 100, 200, 500, 1000, 2000, 5000)
tdf <- c(2, 5, 10, 50)
set.seed(1234)

PWR <- list()
for (d in 1:length(tdf)){
  power <- NULL
for (n in N){
  p1 <- p2 <- p3 <- NULL
for (i in 1:100){
  x <- rt(n, df = tdf[d])
  
  p1 <- c(p1, as.numeric(ks.test(x,"pnorm",mean=mean(x),sd=sd(x),exact=FALSE)$p.value < 0.05))
  p2 <- c(p2, as.numeric(shapiro.test(x)$p.value < 0.05))
  p3 <- c(p3, as.numeric(jarque.bera.test(x)$p.value < 0.05))
}
power <- rbind(power, c(n, mean(p1), mean(p2), mean(p3), mean(p4)))
}
PWR[[d]] <- power
}

par(mfrow = c(2,2))
for (d in 1:length(tdf)){
  power <- PWR[[d]]
  plot(power[,4] ~ log10(power[,1]), col = "red", type = "l", pch = 21, bg = "red", xlab = "Logarimus rozsahu náhodného výberu", ylab = "Empirická síla   testu", main = paste("t rozdelenie (df = ", tdf[d],")", sep = ""), ylim = c(0,1))
  points(power[,4] ~ log10(power[,1]), pch = 21, bg = "red")
  points(power[,3] ~ log10(power[,1]), pch = 21, bg = "blue")
  points(power[,2] ~ log10(power[,1]), pch = 21, bg = "green")
  
  lines(power[,3] ~ log10(power[,1]), col = "blue")
  lines(power[,2] ~ log10(power[,1]), col = "green")
  abline(1,0, col = "black", lty = 3)
  
  legend(1, 1, legend = c("Shapiro-Wilk test", "Jarque-Beta test", "Kolmogorov-Smirnov test"), col = c("blue", "red", "green"), lty= rep(1,3), cex = 0.5)
}

Pri danej alternatíve - studentovo $t$ rozdelenie s 5 stupňami voľnosti sa zdá, že najsilnejším testom je Jarque-Bera test.

Samostatne

Dokážete vysvetliť sledované správanie sa síly jednotlivých testov?
Použijte podobnú simulačnú štúdiu pre iné alternatívy a pre iné štatistické testy.
Rozmyslite si podobnú simulačnú štúdiu nižšie, ktorej cieľom by bolo porovnávať dosiahnutú/empirickú hladinu.
Dokážete vysvetliť sledované rozdiely v jednotlivých obrázkoch?

Nasledúca časť kódu ukazuje simululácie vzhľadom k hladine testu – to je dosiahnutá (empirická) pravdepodobnosť chyby prvého druhu – t.j. pravdepodobnosť, že nulovú hypotézu zamietneme, ak nulová hypotéza platí.

Preto v následujúcom príklade generujeme data z normálneho rozdelenia – t.j. za platnosti nulovej hypotézy a pri zvolenej teoretickej hladine $\alpha = 0.05$ sledujeme, ako často sa test zmýli – tzv., že zamietne nulovú hypotézu (p-hodnota testu bude menšia, než $\alpha$). Správny test by empiricky držať zvolenú teoretickú hladinu $\lpha \in (0,1)$.

N <- c(10, 50, 100, 200, 500, 1000, 2000, 5000)

set.seed(1234)
power <- NULL
for (n in N){
  p1 <- p2 <- p3 <- NULL
for (i in 1:100){
  x <- rnorm(n)
  p1 <- c(p1, as.numeric(ks.test(x,"pnorm",mean=0,sd=1,exact=FALSE)$p.value < 0.05))
  p2 <- c(p2, as.numeric(shapiro.test(x)$p.value < 0.05))
  p3 <- c(p3, as.numeric(jarque.bera.test(x)$p.value < 0.05))
}
power <- rbind(power, c(n, mean(p1), mean(p2), mean(p3)))
}


par(mfrow = c(1,1))
  plot(power[,4] ~ log10(power[,1]), col = "red", type = "l", pch = 21, bg = "red", xlab = "Logarimus rozsahu náhodného výberu", ylab = "Dosiahnutá hladina testu", main = "", ylim = c(0,0.4))
  points(power[,4] ~ log10(power[,1]), pch = 21, bg = "red")
  points(power[,3] ~ log10(power[,1]), pch = 21, bg = "blue")
  points(power[,2] ~ log10(power[,1]), pch = 21, bg = "green")
  
  lines(power[,3] ~ log10(power[,1]), col = "blue")
  lines(power[,2] ~ log10(power[,1]), col = "green")
  abline(0.05,0, col = "black", lty = 2, lwd =2)
  
  legend(1, 0.4, legend = c("Shapiro-Wilk test", "Jarque-Beta test", "Kolmogorov-Smirnov test"), col = c("blue", "red", "green"), lty= rep(1,3))

3. Program R: Čo sa ešte oplatí vedieť, alebo aspoň tušiť?

Niekoľko (snáď) zaujímavých námetov na samostatnú aktivitu…

A) Vlastné funkcie a súbory typu .R a .RData

Program R umožňuje vytváranie rôznych užívateľských funkcií. Vytvorené funkcie je možné jednoducho uložiť do datového súboru typu .R a pri opätovnom načítaní súboru pomocou R príkazu source(), sú funkcie opäť k dispozícii. Nasledujúcim príkazom načítame do programu R vlastné preddefinované funkcie:

rm(list = ls())
source("http://www.karlin.mff.cuni.cz/~maciak/NMSA230/mojeFunkce.R")

Vpodstate sa jedna k štandardný R skript, ktorý je uložený pomocou príkazu save().

Súbor obsahuje jednu preddefinovanú užívateľskú funkciu mojaFunkce1(), ktora je definovaná následovne:

mojaFunkce1 <- function(x){
    m <- mean(x)
    v <- var(x)
    return(paste("Sample mean and sample variance: ", round(m, digits = 2), " (", round(var(x), digits = 3), ")", sep = ""))
}

Štandardným volaním funkcie pomocou jej ména (t.j. mojaFunkce1()) funkciu aplikujeme na konkrétny numerický vektor:

mojeFunkce1(rnorm(1000))

V podkladovom súbore je možné okrem samotných funkcii ukládať aj iné Rkové objekty - napr. výsledky parciálných výpočtov, simulácii a podobne. Nejedná sa ale o samotné objekty vo vorme hotových výsledkov, ale o podkladový Rkový zdrojový kód, ktorý tieto výsledky pri načítaní súboru vyhodnotením zdrojového kódu vytvorí.
Načitaním podkladového .R súboru teda zakaždým dôjde k samotnému výpočtu a vyhodnoteniu všetkých príkazov, ktoré sú v súbore definované. Tento postup samozrejme nie je ideálny v prípade, keď sú výpočty časovo náročné a skôr by sme ocenili pouze jednorázové vyhodnotenie a následne odkazovanie sa na hotové výsledky - v programe R k tomu slúži práve datový typ, R súbor typu .RData.

Nasledujúci príkaz načíta do programu R dva objekty - náhodné výbery o rozsahu 10000, avšak samotné generovanie náhodých výberov už neprebieha, dôjde pouze k načítaniu výsledných hodnôt dvoch náhodných výberov. Samotné generovanie prebehlo vrámci R skriptu, ktorý načítaním súboru nemáme už k dispozícii.

load(url("http://www.karlin.mff.cuni.cz/~maciak/NMSA230/hotoveVysledky.RData"))

Náhodné výbery su (v skutočnosti) vygenerované zo štandardného normálneho rozdelenia a z exponenciálneho rozdelenia so strednou hodnotou 1. Nastavenie seedu bolo zakždým set.seed(1). Empiricky môžeme tvrdenie jednoducho overiť:

set.seed(1)
newSample1 <- rnorm(10000)
sum(newSample1 == sample1)

## [1] 10000

Samostatne

V programe RStudio otvorte ako R skript súbor mojeFunkce.R, ktorý stiahnete tu a vytvorte niekoľko ďalších vlastných funkcii. Analogicky môžete využiť a iný vhodný textový editor. Výsledny súbor opäť uložte, načitajte v programe R a vytvorené funkcie aplikujte.
Vytvorte jednoduchú sumulačnú štúdiu, ktorej parciálne výsledk budete postupne ukládať do súboru typu .RData. Je potrebné ukladanie súboru vhodne implementovať, aby nedochádzalo k prepisovaniu ukladaných výsledkov.
Porovnajte fungovanie R príkazov source() a load().

Užítočné

Dôležité a často používané užívaťeľské funkcie je výhodné uložiť v samostatnom súbore typu .R a podkladový súbor načítať pomocou príkazu source(). Do tohto súboru by mal užívateľ ideálne zasahovať iba vo výnimočnych prípadoch, keď je nutné konkrétnu funkciu pozmeniť, alebo nadefinovať novú.
Postup práce a priebeh (štatistickej) analýzy dat vo forme postupnosti konkrétnych príkazov a riadkov zdrojového kódu je opäť vhodné uložit ako .R súbor, ale z dôvodu častého zasahovania do tohto súboru je doporučené, aby sa jednalo o iný súbor, než ten, v ktorom sú uložené základné uživateľské funkcie.
Jednotlivé riadky zdrojového kódu je užitočne náležite (a dostatočne podrobne) komentovať a k názvom jednotlivých premenných a objektov používať intuitivné značenie.
```
### komentar v programe R, resp. v subore typu .R 
mojaFunkce1(x)
```

B) Datové súbory typu .xls, .csv, alebo .txt

V praxi sa štatistík stretáva najčastejšie s datami ktoré sú reprezentované číselne, väčšinou pomocou tabuľky, ktorá je uložená buť v Exceli, alebo v textových súboroch typu .csv a .txt.

Vyššie spomínaný datový súbor .RData je ideálný pre ukladanie samotných objektov vytvorených v programe R, resp. niektoré vzorové datové súbory určené k výukovým účelom bývajú taktiež uloženém v súbore typu .RData. Tieto objekty nemusia nutne reprezentovať štandardný súbor vo forme tabuľky s hodnotami. Pre uloženie konkrétnych dat vo forme tabuľky sú možno praktickejšie klasické textové súbory typu .xls, .txt, alebo .csv. Pre manipuláciu so súbormi slúžia funkcie read.table() (príkaz z R knižnice ‘gdata’), read,csv() a read.table().

Help k príkazom:

?read.table
?write.table

K uloženiu datovych súborov (ideálne vo forme tabuľky, resp. data.frame() slúžia v programe R príkazy write.xlsx() (z knižnice ‘xlsx’), write.csv() alebo write.table(). K dispozícii je ale množstvo iných knižníc, ktoré ponúkajú alternatívne príkazy buď pre načítanie dat z Excel súboru, alebo dokonca načítanie informácie zo súborov iných typov (napr. obrázky, videa, zvykové záznamy a pod.).

Samostatne

Vytvorte nejaké konkrétne data, alebo využijte datové súbory ktoré sú predinštalované v programe R (príkaz data()) a konkrétnu tabuľku s jednotlivými premennými a pozorovaniami uložte a načítajte ako datové typy .txt, alebo .csv.
Pomocou vhodných dodatočných parametrov sa pokuste formu uložených dat pozmeniť a následne opätovne načítať tak, aby nedošlo k zmene, alebo strate pôvodnej informacie (napr. dodatočné parametre ‘sep’, ‘dec’, ‘skip’, ‘na.string’, ‘as.is’, …).
Pre manipuláciu s datami a hlavne veľkými datovými súbormi sú v programe R užitočné knižnice
```
install.packages("dplyr")
install.packages("plyr")
```
Viac informácii o oboch knižniciach a niekoľko vzorových príkladov napríklad tu: https://www.rdocumentation.org/packages/dplyr/versions/0.7.8 alebo tu https://seananderson.ca/2013/12/01/plyr/.
Užítočne je hlavne podívať sa na funkcie ako aggregate(), apply(), sapply(), lapply(), alebo mapply().
Nainštalujte si knižnicu jpeg (viac podrobnosti napr. tu) a pomocou príkazu readJPEG() načítajte do programu R nejaky obrázok/fotografiu. Podívajte sa, aký objekt ste do programu načítali a pokúste sa tento objekt nejakým sposobom zmeniť. Výsledný/modifikovaný objekt pomocou príkazu writeJPEG() opäť uložte a na obrázok sa podívajte pomocou programu na prezeranie obrázkov.

C). Čo ešte na záver?

Program R využívame hlavne ako nástroj na prípravu a následnú analýzu dat. Program R je možné taktiež využiť k príprave reportu, ktorý k štatistickej analýze väčšinou automaticky patrí. Takýto report by mal obsahovať úvod a zmysluplný popis samotného experimentu, z ktorého data pochádzajú. Taktiež musí dostatočne podrobný popis samotných dat. Z odborného (matematického a štatistického) hľadiska musí obsahovať aj dostatočne presné informácie o použitých metódach pri štatistickom vyhodnotení. Popis musí byť natoľko podrobný, aby s pomocou popisu bolo možné celkový experiment a šštatistickú analýzu replikovať. Závery analýzy ale musia byť interpretované v reči pôvodných dat a to spôsobom, ktorému aj nematematik a neštatistik pochopí. Pri výtvárani takéhoto reportu je užitočné využiť program LaTeX, balíček Sweave, alebo R funkciu cat().

Nie vždy sme ale nútení pripraviť report v PDF formáte. Na druhej strane by ale malo vždy platiť, že samotné zdrojové kódy,defaultné označenie premenných, alebo výstupy z Rkových funkcii by sa v reporte objavovať nemali. To platí aj rôzne tabuľky (napr. tabuľka popisných charakteristík, ako výstup funkcie summary()). V takomto prípade môže byť užitočné poznať R knižnicu ‘formattable’ - viac na https://www.littlemissdata.com/blog/prettytables.

Možnosti programu R sú vpodstate neobmedzené, zaujímavý pohľad na možné nadstavy a rozšírenia ponúka webová stránka https://awesome-r.com/.

Samostatne

Pre účely zápočtu máte pripravený R skript, ktorý obsahuje generovanie dat, jednoduché popisné charakteristiky, vytvorenie obrázkov, aj následne reportovanie pomocou knižnice Sweave. Vyberte si niektorú konkrétnu časť zdrojového kódu a pokúste sa v nejakom zmysle fungovanie programu v danej časti kódu zefektivniť.
Program R ponúka množstvo rôznych nástrojov (napr. rôzne R knižnice) na vytváranie grafických výstupov – obrázkov. Pokúste sa alespoň jeden obrázok vypracovať nejakým alternatívnym spôsobom – t.j., napr. s použitím iného R-kového baličku (knižnice).