Analýza dat v R

1. Analýza dat v R Viacrozmerné štatistické metódy v R, Analýza reálneho príkladu Eva Budinská Ivana Ihnatová 3.-5. máj 2013 Brno

2. Základné štatistické funkcie

3. Zpamätajte si I. • VŽDY si ukladajte všetky skripty. Najlepší spôsob práce je písať príkazy v textovom súbore a vykonávať ich pomocou CTRL+R v R konzole. • VŽDY si robte poznámky k skriptom (za pár týždňov či mesiacov si nebudete pamätať čo ste tým mysleli). • Udržujte skripty prehľadné, ako v názvoch tak v štruktúre, môže sa stať, že po Vás bude niekto analýzu preberať, musí sa v tom vyznať. • Ak používate niektoré skripty často, spravte z nich FUNKCIE s viacerými argumentami ktoré umožňujú flexibilitu: • Vaše skripty budú kratšie a prehľadnejšie • Minimalizujete chyby • Funkciu môžete aplikovať na akýkoľvek podobný dataset bez nutnosti prepisovania premenných • Používajte jednoznačné a pochopiteľné názvy funkcií!

4. Zpamätajte si II. V R-ku pracujú desaťtisíce analytikov po celom svete už dlhé roky, preto: • Ak sa Vám zdá v R niečo príliš komplikované alebo zdĺhavé, je veľmi pravdepodobné, že: • existuje omnoho jednoduchší spôsob, ale vy ho nepoznáte, lebo postrádate správnu funkciu • určite už pred Vami niekto riešil rovnaký problém • a preto určite nájdete odpoveď na webe alebo v helpe! • Help je najlepší spôsob ako pochopiť funkciu, VŽDY sú uvedené príklady. V R-ku help naozaj pomáha!

5. Kapitola 5 Viacrozmerné štatistické metódy v R

6. Dátové súbory • 3 dátové tabuľky vytvorené z dát dostupných na stránke ČSU a wikipédii. • crops.csv – hektárové výnosy 14 plodín na 14 miestach • places.csv – popis miest • plants.csv – popis plodín • Crops crop15 crop16 crop17 place2 5.58 4.08 5.04 place3 5.23 3.63 4.58 place4 4.74 3.19 4.39 • Places Namegroupsize place2 Hl. m. Praha Czech 496 place3 StředočeskýCzech 11016 place4 JihočeskýCzech 10056 • Plants Name type obiloviny crop15 pšenice ozimní obiloviny Y crop16 pšenice jarní obiloviny Y crop17 ječmen ozimní obiloviny Y

7. Načítanie dát • crops<-read.table("crops.csv", header=TRUE, sep=";", row.names=1) • places<-read.table("places.csv", header=TRUE, sep=";", row.names=1) • plants<-read.table("plants.csv", header=TRUE, sep=";", row.names=1) • Kontrola načítania: • dim(crops) • dim(places) • dim(plants) • crops[1:5,1:5] • places[1:5,1:5] • plants[1:5,1:5]

8. Testovanie hypotéz • t.test(x, y = NULL, alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, ...) • x,y - numericky vektor • alternative - typ alternatívnej hypotézy • mu – skutočná priemerná hodnota (rozdiel priemerov pri dvojvyb. teste) • paired – párové / nepárové usporiadanie • var.equal - predpoklad rovnosti rozptylov • conf.level - hladina významnosti • Výstup – špecifická trieda htest • Hodnota testovej štatistiky dostupná cez $statistic, p-hodnota $p-value • Dostupnéfunkcie pre rýchlejšie výpočty v genomike: rowttests (z genefilter)

9. Testovaniehypotéz - príklad • informácia o výnose plodín je uložená v súbore crops • plodiny sú uložené ako stĺpce, preto sa na ne budeme odvolávať pomocou crops[,i] • informácia o príslušnosti miest ku skupine je v súbore places • uložme si informáciu o skupine do premennejgroup group = places$group • Funkcia preT-test: t.test() • Možný spôsob vstupných hodnôt • t.test(hodnoty v skupine 1, hodnoty v skupine 2) • t.test(hodnoty~skupina) Konkrétne pre otestovanie i-tej plodiny (stĺpca) • t.test(crops[group=="Czech",i], crops[group=="Morav",i]) • t.test(crops[,i]~group)

10. Testovaniehypotéz - príklad i=1 vysledok = t.test(crops[,i]~group) vysledok #zobrazime výsledok testu Ako vyextrahujeme p-hodnotu??? summary(vysledok) #ukáže nám všetky elementy výsledného objektu Je medzi nimi i hodnota p.value. Získame ju teda odvolaním sa na ňu cez $p.value: vysledok$p.value [1] 0.32249 Nakreslime si obrázok, aby sme videli rozdiely medzi skupinami boxplot(crops[,i]~group)

11. Výsledný boxplot

12. Testovaniehypotéz- príklad • Chceme spočítať T-test pre každú plodinu. Aplikujeme teda funkciu apply na súbor crops po riadkoch, a z t.testu vyberieme len p-hodnoty: vsetky.ttest = apply(crops, 1, FUN=function(x) t.test(x~ group)$p.value) • Chceme zistiť, výnosy ktorých plodín sú štatisticky významné, a teda ktoré p-hodnoty sú menšie ako 0.05 vyznamne.t.test = vsetky.ttest<0.05 #podmienkou vytvoríme vektor false/true Úloha: Ako zistíme, koľko ich je? Aké iné (aspoň 3) štatistické testy R ponúka? Nápoveda: balík stats

13. Mnohonásobné testovanie hypotéz • Pri testovaní veľkého počtu hypotéz je potrebná korekcia p-hodnôt. • 10 000 hypotéz na hladine významnosti 5% -> 500 falošne pozitívných výsledkov • p.adjust(p, method = p.adjust.methods, n = length(p) • p –vektor p-hodnôt • method – napr. Bonferoniho korekcia "bonferroni", Benjamini & Hochberg "BH" alebo "fdr" Príklad: p.adj<- p.adjust(vsetky.ttest, "bonferroni")

14. Zhluková analýza • Hierarchické zhlukovanie: • Z balíka cluster: • hclust(d, method = "complete", members=NULL) • d – matica vzdialenosti medzi objektami • method – metóda: "ward", "single", "complete", "average", "mcquitty", "median"alebo"centroid" • members – počet pozorovaní pri zhlukovaní zhlukov • Výstupom je objekt triedy hclust, pre ktorý existujú zvláštne metódy pre print, plot atď.

15. Zhlukovanie - príklad • Chceme vedieť, ktoré zo skúmaných plodín tvoria spoločné zhluky • Takisto chceme vedieť, či skúmané plodiny dokážu rozdeliť miesta podľa definovanej skupiny • Urobíme preto hierarchické zhlukovanie jednak plodín a jednak miest.

16. Metrika vzdialeností I. • Ak chceme zhlukovať, najprv musíme definovať maticu vzdialeností d, ktorá je základným vstupom do funkcie hclust. Použijeme dve: • Euklideovskú vzdialenosť – funkcia dist() dist(x,method) parametre: x – dátová tabuľka, objekty medzi ktorými počítame vzdialenosť musia byť v riadkoch method – metrika vzdialenosti euclid.plant = dist(t(crops), method="euclidean")# vzdialenosť medzi plodinami euclid.place = dist(crops, method="euclidean")# vzdialenosť medzi miestami

17. Metrika vzdialeností II. • Metrika založená na Spearmanovej korelácii – cor() cor(x, method) – vypočíta koreláciu medzi stĺpcami matice x – matica method – "pearson" alebo "spearman" cor.plant = cor(crops, method="spearman")#počítame koreláciu medzi plodinami Potrebujeme ale vzdialenosť!!! cor.plant.dist = as.dist(1-cor.plant)# prepočítame korelácie na maticu vzdialeností cor.place = cor(t(crops), method="spearman") #počítame koreláciu medzi miestami cor.place.dist = as.dist(1-cor.place) # prepočítame na vzdialenosť

18. Hierarchické zhlukovanie • Zhlukujeme plodiny algoritmom average linking na euklideovskej vzdialenosti: hclust.plant.eucl.avg = hclust(euclid.plant, method="average") a na vzdialenosti odvodenej od korelácie hclust.plant.spear.avg = hclust(cor.plant.dist, method="average") • Podobne zhlukujeme miesta: na euklideovskej vzdialenosti: hclust.place.eucl.avg = hclust(euclid.place, method="average") a na vzdialenosti odvodenej od korelácie hclust.place.spear.avg = hclust(cor.place.dist, method="average")

19. Dendrogram • Máme teda 4 objekty hierarchického zhlukovania hclust.plant.eucl.avg #plodiny na euklideovskej vzd. hclust.plant.spear.avg #plodiny na koreláciách hclust.place.eucl.avg #miesta na euklideovskej vzd. hclust.place.spear.avg #miesta na koreláciách • Ak chceme výsledky zobraziť, musíme vytvoriť dendrogramy dendr.hclust.plant.eucl.avg = as.dendrogram(hclust.plant.eucl.avg, labels=group) dendr.hclust.plant.spear.avg = as.dendrogram(hclust.plant.spear.avg) dendr.hclust.place.eucl.avg = as.dendrogram(hclust.place.eucl.avg) dendr.hclust.place.spear.avg = as.dendrogram(hclust.place.spear.avg)

20. Dendrogram - zobrazenie • Vykreslíme dendrogramy do jedného grafu pre plodiny: par(mfrow=c(2,1)) plot(dendr.hclust.plant.eucl.avg) plot(dendr.hclust.plant.spear.avg) savePlot("dendrogram.plants.png", type="png") a pre miesta windows(14,14) par(mfrow=c(2,1), cex=0.8) plot(dendr.hclust.place.eucl.avg) plot(dendr.hclust.place.spear.avg) savePlot("dendrogram.place.png", type="png")

21. Analýza hlavných komponent (PCA) • princomp – R-mode PCA, vlastné čísla, viac pozorovaní ako parametrov • prcomp – Q-mode PCA, singulárny rozklad • pr1<-prcomp(plants[,6:9]) • pr2<-princomp(plants[,6:9]) • biplot(pr1) • biplot(pr2)

22. Heatmapa • Dostupných niekoľko funkcií: heatmap, heatmap.2 (z balíka gplots), • heatmap(x, Rowv, Colv, ColSideColors, labRow, cexRow, col, …) • Vybrané argumenty: • x – numerická matica • Rowv, Colv – dendrogram pre riadky alebo stĺpce • ColSideColors – farebné označenia stĺpcov (podobne pre riadky RowSideColors) • labRow – textové popisky riadkov (pre stĺpce labCol) • cexRow – relatívna veľkosť písma pre popisky riadkov (pre stĺpce cexCol) • col – farebná škála

23. Heatmapa - príklad • Nakreslíme heatmapu – častý výstup analýzy microarrayí Nastavíme farebnú škálu zelenej cez čiernu po červenú, kde zelená budú málo výnosné plodiny a červená veľmi výnosné plodiny. Čierna farba znamená priemerný výnos. farby = colorRampPalette(c("green", "black","red"))(15) Crops<-as.matrix(crops) heatmap(Crops) heatmap(Crops[p.adj<0.05,], col=farby) pridáme farby skupín (parameter RowSideColors): heatmap(Crops, col=farby, RowSideColors=as.character(as.numeric(group)))

24. Heatmapa a vlastné zhlukovanie • Ak chceme heatmapu na základe nášho vlastného zhlukovania, musíme definovať parametre Rowv (pre dendrogram riadkov - miest) a Colv (pre dendrogram stĺpcov-plodín) • Pre miesta použijeme dendrogram zo zhlukovania na spearmanovej korelácii dendr.hclust.place.spear.avg • Pre plodiny použijeme dendrogram zo zhlukovania na euklideovskej vzdialenosti dendr.hclust.plant.eucl.avg heatmap(Crops, col=farby, RowSideColors=as.character(as.numeric(group)), Colv= dendr.hclust.plant.spear.avg, Rowv= dendr.hclust.place.eucl.avg) Zadáme dendrogram z korelácií i u plodín: heatmap(Crops, col=farby, RowSideColors=as.character(as.numeric(group)), Colv= dendr.hclust.plant.spear.avg, Rowv=dendr.hclust.place.spear.avg)

25. Upravujeme heatmapu • Chceme vidieť názvy plodín (plants$Name), zadefinujeme preto parameter labCol a zväčšíme písmo v stĺpcoch – parameter cexCol: heatmap(Crops, col=farby, RowSideColors=as.character(as.numeric(group)), Colv= dendr.hclust.plant.spear.avg, Rowv= dendr.hclust.place.spear.avg, labCol=plants$Name, cexCol=1.3) • Ak nechceme zhlukovanie riadkov alebo stĺpcov, zadáme hodnotu NA do parametrov Colv alebo Rowv heatmap(Crops, col=farby, RowSideColors=as.character(as.numeric(group)), Colv= dendr.hclust.plant.spear.avg, Rowv=NA, labCol=plants$Name, cexCol=1.3)

26. Upravujeme heatmapu II. Heatmapu uložíme napríklad príkazom png v definovanej veľkosti: png("heatmapa.png", height=1200, width=600) heatmap(Crops, col=farby, RowSideColors=as.character(as.numeric(group)), Colv= dendr.hclust.plant.spear.avg, Rowv= dendr.hclust.place.spear.avg, labCol=plants$Name, cexCol=1.3) dev.off() Pretože heatmapa tejto funkcie je zadefinovaná cez layout a nie je tu možné meniť ho, nie je možné ani jej zväčšovanie. Ak sa Vám nechce prerábať zdrojový kód, skúste funkciu heatmap.2 z gplots, tá už má viac možností!