1 / 26

Analýza dat v R

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. Základné štatistické funkcie. Zpamätajte si I.

judith-lynn
Download Presentation

Analýza dat v R

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  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í!

More Related