NMSA230 - Úvod do programování v R

Zimný semester 2021/2022 | Cvičenie 5 | St 15/12/21

V. Knižnice Knitr & Shiny

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 datami, jednoduché simulácie a R knižnice Knitr (na prípravu HTML Markdownow, ako je napríklad táto stránka) a Shiny (na prípravu interaktívnych webových aplikácii).

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, ktorý je možné prezerať pomocou internetového prehliadača. Takto bola vytvorená napríklad aj táto stránka.

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 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 v pre knižnicu Knitr oddelený pomocou párových znakov

pričom v úvode je opäť dovolené používať dodatočné premenné určené (ako 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 samotni 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íkaz grep()).

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, alebo podobný data processing step je väčšinou nutnou súčasťou pri š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.

Samostatne

Podívajte sa na data a skuste 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 o starých a nových kritériach vzhledem na výskyt rôznych komplikácii (rozlišujete pacientky v závislosti na výsledkoch starých a nových kriterii, resp. na novo-zavedenej premenne 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

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 <- NULL
for (i in 1:100){
  
  x <- rchisq(n, df = 50) ### pozorovania simulujeme z chi^2 rozdelenia
  
  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)))
}

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

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)))
}
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ť sledovane správanie sa síly jednotlivých testov?
Použijte podobnú simulačnú štúdiu pre iné alternatívy a pre iné testy.
Rozmyslite si podobnu simulačnu štúdiu nižšie, ktorej cieľom je 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. Knižnica ‘shiny’

Shiny je (open source) knižnica pre štatistický program R (inštalácia pomocou štandardného príkazu install.packages("shiny")), ktorá umožňuje vytvárať hyperaktívne webové aplikácie, pomocou ktorých dokáže koncový úžívateľ jednoducho meniť nastavenia a kontrolovať priebeh výpočtu/analýzy nejakého konkrétneho procesu (napr. štatistickej analýzy datového súboru).

Knižnica shiny umožňuje jednoducho premeňiť R-kový zdrojový kód na interaktívnu webovú aplikáciu (vyžaduje k tomu ale tzv. shiny server, ktorý beží na hosťovskom počítači spolu s programom R a umožňuje tak webovej aplikácii na klientskom počítači plnohodnotne fungovať - vyhodnocovať R-kové príkazy v závislosti na klientských modifikáciach, ktoré uskutočnuje priamo vo svojom internetovom prehliadači).

Inšpirácia a niekoľko príkladov vytvorených pomocou knižnice shiny: http://shiny.rstudio.com/gallery/

Pre účely tohto cvičenia nebude potrebne inštalovať shiny server na hosťovskom počítači. Niekoľko príkladov si ukážeme pouze lokálne a ako hosťovský server nám bude slúžiť priamo počítač, ktorý ma každý k dispozícii a na ktorom pracuje.
Táto čast slúži pouze jako ilustrácia rôznych možnosti, ktoré má užívateľ programu R k dispozícii.

K vytvoreniu funkčnej shiny aplikácie je potrebných niekoľko krokov:

nainštalovať a načítať knižnicu shiny:
```
library("shiny")
```
vytvoriť (HTML) stránku pomocou R-kového príkazu fluidPage(), ktorá bude slúžiť ako vstup pre program R a zároveň bude poskytovať potrebné zdrojové informácie pre program R.
Pomocou tejto stránky užívateľ kontroluje aplikáciu a modifikuje výpočty/analýzu;
naprogramovať príslušný R-kový kód, ktorý beží na serveri shiny a ktorý pracuje so vstupnými informáciami poskytnutými webovou stránkou vytvorenou pomocou fluidPage() a po spracovaní a vyhodnotení R-kového zdrojového kódu odošle príslušný výstup späť na stránku. Požadovaný výstup sa uživateľovi zobrazí priamo na webovej stránke;
spustit samotnú shiny aplikáciu pomocou príkazu shinyApp();

Najjednoduchší príklad funkčnej shiny aplikácie (ktorá ale v zásade nič nerobí), je vytvorená pomocou nasledujúcej časti R-kového kódu:

library("shiny")

ui <- fluidPage("Hello Word"
   ### su za zadavaju vstupne informacie (inputs) a zobrazuju sa vystupy (outputs)
)

server <- function(input, output){
  ### R-kovy zdrojovy kod, ktory vyuziva list input (vstupne informacie) a vytvara z toho list vystupov (output)
}

### spustenie samotnej aplikacie (otvori webovy prehliadac s pozadovanou strankov vytvorenou vo fluidPage())
shinyApp(ui = ui, server = server)

Obe časti - webovú stránku ui, aj shiny server server - lze následne doplniť o dodatočné informácie, HTML kód, rôzne preddefinované moduly na zadávanie jednotlivých príkazov, či konkrétnych dat do programu R, a tiež moduly na vizualizáciu výstupov (tabuľky, obrázky, grafy a pod.).

Doležité

Pri spustení R-kového kódu uvedeného vyššie, dôjde k otvoreniu webového prehliadača, v ktorom sa otvorí prázdna stránka s uvítaním “Hello World”. Príkazový riadok v R-ku (konzola v RStudio) zostáva aktívna a nie je možné zadavať prostredníctvom konzoly dalšie príkazy. R-ko očakáva, že príkazy bude dostavať prostredníctvom otvoreného a aktivného webového prehliadača. Keďže nateraz nemáme v prehliadači nič vytvorené (príkaz fluidPage() je prázdny), musíme R-ko ukončiť manuálne, buď kliknutím na tlačítko STOP, alebo stlačením kombinácie kláves CTRL + C.

Všetky hodnoty, ktoré potrebujeme R-ku predať ako vstupné hodnoty, sú automaticky uložené v liste (R object) s názvom input a na jednotlivé zložky listu sa odkazujeme pomocou ID nálepky pri vstupe. Analogicky, všetky výstupy, ktoré server pripraví, aby boli dostupné prostredníctvom klientskej webovej stránky, sa ukladajú do listu (R object) s názvom output a na jednotlivé zložky sa opať odkazujeme pomocou príslušných ID nálepiek.

Nasledujúci príklad umožnuje vykresliť histogram pre náhodný výber z normálneho rozdelenia $N(\mu, \sigma^2)$ o rozsahu $n \in \mathbb{N}$, kde všetky tri parametre ($\mu$, $\sigma^2 > 0$ aj $n \in \mathbb{N}$) môže nastaviť užívateľ prostredníctvom webového prehliadača.
Histogram a odhad hustoty sa interaktvívne vykreslí priamo v prehliadači, hneď po zadaní vstupných hodnot.

library("shiny")

ui <- fluidPage(
  
  title = "Simple Example: Simulating data from normal distribution",
  
  sidebarPanel(
    sliderInput(inputId = "mean", label = "Choose Mean Value", min = -100, max = 100, value = 0),
    sliderInput(inputId = "var", label = "Choose Variance Value", min = 0, max = 10, value = 1),
    numericInput("number", label = "Number of Observations", value = 10),
    numericInput("breaks", label = "Number of Breaks", value = 10)
  ),
  
  mainPanel(
    plotOutput(outputId = "density")
  )
)

server <- function(input, output){
  output$density <- renderPlot({sample <- rnorm(input$number, input$mean, sqrt(input$var))
                                hist(sample, breaks = input$breaks, col = "blue", freq = F)
                                lines(density(sample), col = "red", lwd = 2)
                                })
}

shinyApp(ui = ui, server = server)

Ak ste príkazy zadali správne, vo webovom prehliadači by ste mali vidieť niečo podobné:

Pomocou nástrojov na zadávanie vstupných hodnôt na ľavej strane je možné modifikovať nastavenie generátoru náhodných čísel. Následne sa v programe R vygeneruje nový/iný náhodný výber a príslušne sa prekresli histogram a odhad hustoty.

Príkazový riadok v RStudio zostáva naďalej aktívny a nie je možné zadavať príkazy jinak, než pomocou otvoreného webového prehliadača - program R totiž očakáva vstupné hodnoty, ktoré budu zadané pouze prostredníctvom nástrojov na webovej stránke. V prípade, že chceme R sekciu ukončiť, je nutné stlačiť tlačítko STOP alebo kombináciu kláves CTRL + C (prípadne zatvoriť okno webového prehliadača).

Doležité

Všimnite si, akým sposobom funguje zadávanie vstupných hodnôt a akým spôsobom sú tieto hodnoty následne predané programu R na sever shiny

V knižnici shiny je niekoľko rôzných nástrojov na zadávanie vstupných hodnôt - viď napr. http://shiny.rstudio.com/tutorial/lesson3/;
Každý príkaz, ktorý umožnuje zadávať vstupné hodnoty, musí mať priradené jednoznačné ID - pozri prvý argument v príkazoch vyššie. S týmto ID je potom vstupná hodnota predaná v liste input do R-ka: napr. vstupná hodnota pre voľbu strednej hodnoty je uložena v input$mean;
Fiktívny server - v našom prípade užívateľský počítač - využíva hodnoty uložené v liste input a v sérii R-kových príkazov pripraví príslušný výstup/výstupy, ktoré analogicky uloži do listu output;
V druhej časti príkazu fluidPage() umiestňujeme na webovú stránku pripravené výstupy, ktoré predal server po spracovaní programom R;
Analogicky ako v prípade zadávania vstupných hodnôt, aj pre vykresľovanie výstupov je v knižnici shiny k dispozícii niekoľko užitočných nástrojov - viď napr. http://shiny.rstudio.com/tutorial/lesson4/;

Samostatne

Využijte ďalšie možnosti, ktoré ponúka knižnica shiny a pokuste sa vytvorit interaktívnu aplikáciu, pomocou ktorej budete analyzovať nejaky datový súbor, ktorý je k dispozícii v programe R, avšak analýza bude závisieť na uživateľských nastaveniach (napr. výber premenných, voľba podsúboru dat, voľba hladiny spoľahlivosti a pod.).
Pokuste sa využiť rôzne možnosti zadávania vstupných informácii a tiež rôzne možnosti prezentovania výstupu.
Využijte rôzne varianty k peknému a účelnému usporiadaniu nástrojov na zadávanie vstupných informacii a panelov na vykresľovanie spočítaných výstupov.

Komplexnejšia ukážka s využitim viacerých panelov pre rôzne varianty výstupov je napr. nižšie:

library("shiny")

ui <- fluidPage(
  
  title = "Simple Example",
  
  sidebarPanel(
    sliderInput(inputId = "mean", label = "Choose Mean Value", min = -100, max = 100, value = 0),
    sliderInput(inputId = "var", label = "Choose Variance Value", min = 0, max = 10, value = 1),
    numericInput("number", label = "Number of Observations", value = 10),
    numericInput("breaks", label = "Number of Breaks", value = 10)
  ),
  
  mainPanel(
    tabsetPanel(
      tabPanel("Plot",  plotOutput(outputId = "density")), 
      tabPanel("Summary", verbatimTextOutput("summary")),
      tabPanel("Regression", plotOutput(outputId = "regression")),
      tabPanel("Regression Summary", verbatimTextOutput("regsum"))
    )
  )
)

server <- function(input, output){
  
  output$density <- renderPlot({set.seed(1234)
                                data <- data.frame(1:input$number, rnorm(input$number, input$mean, sqrt(input$var)))
                                hist(data[,2], breaks = input$breaks, col = "blue", freq = F)
                                lines(density(data[,2]), col = "red", lwd = 2)
                                })
  
  output$summary <- renderPrint({set.seed(1234)
                                 data <- data.frame(1:input$number, rnorm(input$number, input$mean, sqrt(input$var)))
                                 summary(data[,2])})
  
  output$regression <- renderPlot({set.seed(1234)
                                   data <- data.frame(1:input$number, rnorm(input$number, input$mean, sqrt(input$var)))
                                   plot(data[,2] ~ data[,1])
                                   model <- lm(data[,2] ~ data[,1], data = data, xlab = "Observation Number", ylab = "Observation Value")
                                   abline(model$coeff[1], model$coeff[2], col = "red", lwd = 2)})
  
  output$regsum <- renderPrint({set.seed(1234)
                                data <- data.frame(1:input$number, rnorm(input$number, input$mean, sqrt(input$var)))
                                plot(data[,2] ~ data[,1])
                                model <- lm(data[,2] ~ data[,1], data = data, xlab = "Observation Number", ylab = "Observation Value")
                                summary(model)})
}

shinyApp(ui = ui, server = server)

Vyskúšajte uvedený skript a pokúste sa ho podľa vlastných nápadov vhodne modifikovať. Vygenerovaná webová stránka by mala vyzerať nejak takto:

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

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/.

Požiadavky na zápočet

Nutnou a zároveň postačujúcou podmienkou na získanie zápočtu (ak má študent/študentka predmet riadne zapísaný v aktuálnom školskom roku v systéme SIS) je vypracovanie parciálných úloh zadaných v priebehu semestra a včasné odovzdanie celkového PDF reportu vrámci stanoveného deadlinu.
Parciálne úlohy postupne zadávané na seminároch č.1 až 4 zahrňovali vygenerovanie vlastného datového súboru, vytvorenie niektorých dôležitých popisných charakteristík s vhodným grafickým doplnením a celkové spracovanie pomocou R knižnice ‘Sweave’.
Finálný súbor PDF je potrebné zaslať emailom

na adresu maciak[AT]karlin.mff.cuni.cz najneskôr do konca skúškového obdobia, t.j. do 13.02.2022..