install.packages(vegan) library(vegan) getwd() setwd("c:/DANKA HLAVNY ADRESAR/Vyuka/Bi5980 Biodiverzita/Biodiverzita vyuka 2015") sneky <- read.csv("06_2015_DataSneky_CaseStudy.csv",sep = ";" ,dec = ",", header = T, row.names = 1) # nazvy riadkov su povodny prvy stlpec sneky head(sneky) tail(sneky) snekyDruh <- sneky[ ,1:56] # nacteni casti tabulky s druhy na lokalitach jako snekyDruh snekyDruh snekyEnv <- sneky[ ,57:67] # nacteni casti tabulky s vlastnostmi lokalit jako snekyEnv snekyEnv snekyDruh.bray <- vegdist(snekyDruh, "bray", diag = T, upper = T) snekyDruh.bray hclust(snekyDruh.bray, method = "complete") snekyDruh.bray.complete <- hclust(snekyDruh.bray, method = "complete") plot(snekyDruh.bray.complete, main = "Bray-Curtis, Complete linkage") plot(snekyDruh.bray.complete, labels = rownames(snekyDruh.bray), hang = -1, main = "Bray-Curtis, Complete linkage") plot(snekyDruh.bray.complete, labels = rownames(snekyDruh.bray), hang = -1, axes = FALSE, main = "Bray-Curtis, Complete linkage") plot(snekyDruh.bray.complete, labels = rownames(snekyDruh.bray), hang = -1, axes = FALSE, ann = FALSE) ----------------------------------- oznaceni shluku v dendrogramu ------------------------------------------------------------------- plot (snekyDruh.bray.complete, labels = rownames(snekyDruh.bray), hang = -1, axes = FALSE, main = 'Complete linkage') rect.hclust (snekyDruh.bray.complete, k = 2) rect.hclust (snekyDruh.bray.complete, k = 3, border = "blue") rect.hclust (snekyDruh.bray.complete, k = 4, border = "green") rect.hclust (snekyDruh.bray.complete, k = 5) prislusnost.5shluky <- cutree (snekyDruh.bray.complete, k = 5) # vrati prislusnost objektu ke shlukum prislusnost.5shluky ---------------------------------------------------- PCA ---------------------------------------------------------------------------- PCA je použitelná pro homogenní druhová data; nebo po správné transformaci dat - (doporučuje se Hellingerova transformace; je to odmocnina z relativní početnosti, tj. z dominance; tj. stejné hodnoty dostaneme i z matice abundancí i z matice dominancí) snekyDruh.hellinger <- decostand(snekyDruh, "hellinger") snekyDruh.hellinger ?rda # rda (library vegan) - ak nie sú špecifikované environmentální proměnné, tak funkce "rda" z balíku "vegan" udělá unconstrained ordination (PCA) PCA.snekyDruh.hellinger <- rda(snekyDruh.hellinger) PCA.snekyDruh.hellinger # celkova inerce a vlastni hodnoty prvnich 8 os summary(rda(snekyDruh.hellinger)) # vsechny vysledky plot(rda(snekyDruh.hellinger)) # alebo: plot(PCA.snekyDruh.hellinger) ... druhy i lokality v jednom ordinacnim grafu, druhy jako body biplot (PCA.snekyDruh.hellinger, display = 'species') # při použití funkce "ordiplot" se vykreslí druhy i lokality jako centroidy biplot (PCA.snekyDruh.hellinger, display = 'sites') par(mfrow=c(1,2)) # druhy a lokality nakreslime vedle sebe do jednoho okna biplot (PCA.snekyDruh.hellinger, display = 'species') biplot (PCA.snekyDruh.hellinger, display = 'sites') source ('http://www.davidzeleny.net/anadat-r/doku.php/en:numecolr:cleanplot.pca?do=export_code&codeblock=0') # funkce cleanplot.pca pro dva grafy s ruznym skalovanim cleanplot.pca (PCA.snekyDruh.hellinger) ------ jak se rozhodnou o poctu hlavnich komponent, ktere budeme interpretovat? -------------------------------------------------- source ("http://www.davidzeleny.net/anadat-r/doku.php/en:numecolr:evplot?do=export_code&codeblock=0") # definujeme funkci k "evplot" ev <- PCA.snekyDruh.hellinger$CA$eig # select the data frame with eigenvalues of particular axes: evplot (ev) # spocita dulezitost os a vykresly do sloupcoveho grafu install.packages(BiodiversityR) library(BiodiversityR) sig <- PCAsignificance (PCA.snekyDruh.hellinger, axes = 25) # spocita vyznam uvedeneho poctu hlavnich komponent sig # procento vysvetleneho rozptylu muzeme vykreslit do sloupcoveho grafu barplot (sig[c('percentage of variance', 'broken-stick percentage'), ], beside = T, xlab = 'PCA axis', ylab = 'explained variation [%]', col = c('grey', 'black')) legend ('topright', legend = c('PCA', 'broken-stick'), pt.bg = c('grey', 'black'), pch = 22) --------------------------- CA ------------------------------------------------------------------------------------------------- CA.snekyDruh <- cca(snekyDruh) CA.snekyDruh # celkova inerce a vlastni hodnoty prvnich 8 os summary(cca(snekyDruh)) # vsechny vysledky plot(cca(snekyDruh)) # druhy i lokality v jednom ordinacnim grafu plot(CA.snekyDruh) # to iste ordiplot(CA.snekyDruh) # funkce ordiplot nakresli pouze body bez labels, muzeme je tam ale dopsat, # pripadne ordiplot (CA.snekyDruh, type = 't') napise texty namiesto bodov ordiplot(CA.snekyDruh, choices = c(1, 2)) # choices - osi, ktere se nakresli CAfig <- ordiplot(CA.snekyDruh, type = "none") # nebudou nakresleny body pro druhy a lokality points(CAfig, "sites", pch=21, col="red", bg="yellow") # nakresli body pro lokality (ne popisky) text(CAfig, "species", col="blue", cex=0.9) # nakresli popisky pro druhy (ne body, popisky jsou na pozici bodu) # chceme nakreslit dva obrazky - jeden pouze s druhy, druhy pouze s lokalitami par(mfrow=c(1,2)) CAfig <- ordiplot(CA.snekyDruh, type = "none") text(CAfig, "species", col="blue") CAfig <- ordiplot(CA.snekyDruh, type = "none") text(CAfig, "sites", col="red", cex=0.7) # cex=0.7 zmensena velikost fontu --------------------------- DCA ------------------------------------------------------------------------------------------------ DCA.snekyDruh <- decorana(snekyDruh) # kdyz je to potrebne, muzeme pouzit logaritmickou transformaci, tj. log1p(x) - vypocita log(1+x) DCA.snekyDruh # ve vysledku nedostaneme celkovou inerci, tu musime ziskat z CA. DCA nam slouzi k zmereni delky gradientu. ordiplot(DCA.snekyDruh) ordiplot(DCA.snekyDruh, type = "t") DCAfig <- ordiplot(DCA.snekyDruh, type = "none") text(DCAfig, "species", col="blue") DCAfig <- ordiplot(DCA.snekyDruh, type = "none") text(DCAfig, "sites", col="red") # chceme nakreslit dva obrazky - jeden pouze s druhy, druhy pouze s lokalitami par(mfrow=c(1,2)) DCAfig <- ordiplot(DCA.snekyDruh, type = "none") text(DCAfig, "species", col="blue", cex=0.7) # prida popisky druhu/lokalit DCAfig <- ordiplot(DCA.snekyDruh, type = "none") text(DCAfig, "sites", col="red", cex=0.7) # anebo jednoduchseji par(mfrow=c(1,2)) ordiplot (DCA.snekyDruh, display = 'species', type = 't') ordiplot (DCA.snekyDruh, display = 'sites', type = 't') # znazornit pouze druhy a to tak, aby se nazvy neprekryvaly ordiplot (DCA.snekyDruh, display = 'species', type = 'n') # nakresli prazdny ordinacni diagram orditorp (DCA.snekyDruh, display = 'species') # orditorp prida labels do existujiciho ordinacniho diagramu tak ze jsou citelne # muzeme urcit, ktere labels budou videt (parametr select) --------------------------- NMDS ----------------------------------------------------------------------------------------------- snekyDruh.bray <- vegdist(snekyDruh, "bray", diag = T, upper = T) snekyDruh.bray ?metaMDS NMDS.snekyDruh <- metaMDS(snekyDruh, distance = "bray", k = 2) NMDS.snekyDruh ordiplot (NMDS.snekyDruh, type = 't') NMDS.snekyDruh <- metaMDS(snekyDruh, distance = "bray", k = 2, autotransform = FALSE) ordiplot (NMDS.snekyDruh, type = 't') par (mfrow = c(1,2)) stressplot (NMDS.snekyDruh) Cvykresli Shepardov diagram plot (NMDS.snekyDruh, display = 'sites', type = 't', main = 'Goodness of fit') # vykresli lokality ako text, tj. labels points (NMDS.snekyDruh, display = 'sites', cex = goodness (NMDS.snekyDruh)*200) # prida body s velkostou odrazajucou goodness-of-fit (cim vacsi tym horsie) plot (NMDS.snekyDruh, display = 'species', type = 't') orditorp (NMDS.snekyDruh, display = 'sites') --------------------------- suplementary variables - k vysvetleni vyznamu os v neprime ordinaci --------------------------------- # mame PCA, obrazek s lokalitami; funkce envfit pouzije k fitovani suplementary variables mnohonasobni regresi; vyznamnost se testuje permutacnim testem biplot (PCA.snekyDruh.hellinger, display = 'sites') PCAsuplementary <- envfit (PCA.snekyDruh.hellinger, snekyEnv [, c('Typ', 'Area_odhad', 'Coverage', 'Altitude_opr')], perm = 999) PCAsuplementary plot(PCAsuplementary, cex=0.7) # cex - mensi velikost pisma # mame CA, obrazek s lokalitami; ordiplot (CA.snekyDruh, display = 'sites', type = 't') CAsuplementary <- envfit (CA.snekyDruh, snekyEnv [, c('Typ', 'Area_odhad', 'Coverage', 'Altitude_opr')], perm = 999) CAsuplementary plot(CAsuplementary, cex=0.7) # mame DCA, obrazek s lokalitami; DCAsuplementary <- envfit (DCA.snekyDruh, snekyEnv [, c('Typ', 'Area_odhad', 'Coverage', 'Altitude_opr')], perm = 999) DCAsuplementary plot(DCAsuplementary, cex=0.7) -------------------------------------------------- CCA -------------------------------------------------------------------------- CCA.snekyDruh <- cca(snekyDruh, snekyEnv[ , c(2,9,10)]) CCA.snekyDruh ordiplot (CCA.snekyDruh, type = 't') ordiplot (CCA.snekyDruh, type = 't', choices=c(1,4)) ---------------------------------------------------------------------------------------------------------------------------------