User Tools

Site Tools


en:history:2015-03-20-anadatr

2015-03-20 PCA, NMDS, env. variables in unconstrained ordination, Bonferroni correction of P-values

anadatr20032015.R
### PCA on environmental data
 
chem <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/chemistry.txt', row.names = 1)
chem <- chem[,-15]
chem <- scale (chem, center = T, scale = T)
 
PCA <- rda (chem)
biplot (PCA, display = 'species')
 
evplot (PCA$CA$eig)
 
scores (PCA, display = 'si', choices = 1)
 
### NMDS on compositional data
vltava.spe <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/vltava-spe.txt', row.names = 1)
NMDS <- metaMDS (vltava.spe)
ordiplot (NMDS, display = 'sp', type = 't')
stressplot (NMDS)
 
### Solution of Exercise 3 from PCoA & NMDS
tikus.spe <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/tikus-spe.txt', row.names = 1)
tikus.env <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/tikus-env.txt')
 
# select the years 81,83 and 85 only
tikus.spe.sel <- tikus.spe[tikus.env$Year %in% c('1981', '1983', '1985'),]
tikus.env.sel <- tikus.env[tikus.env$Year %in% c('1981', '1983', '1985'),]
 
year <- as.numeric (as.factor (tikus.env.sel$Year))
tikus.spe.sel <- log1p (tikus.spe.sel)
 
library (vegan)
nmds.bray <- metaMDS (tikus.spe.sel)
nmds.eucl <- metaMDS (tikus.spe.sel, dist = 'eucl')
 
par (mfrow = c(1,2))
ordiplot (nmds.bray, display = 'si', type = 'n')
points (nmds.bray, col = year, pch = year)
 
ordiplot (nmds.eucl, display = 'si', type = 'n')
points (nmds.eucl, col = year, pch = year)
 
### Passively projected environmental variables 
 
vltava.spe <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/vltava-spe.txt', row.names = 1)
vltava.env <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/vltava-env.txt')
 
vltava.spe.hell <- decostand (log1p (vltava.spe), 'hellinger')
PCA <- rda (vltava.spe.hell)
ordiplot (PCA, display = 'si')
#biplot (PCA)
 
names (vltava.env)
vltava.env.sel <- vltava.env[,c("ASPSW", "SOILDPT", "pH.H")]
 
ef <- envfit (PCA, vltava.env.sel)
plot (ef)
ef$vectors$pvals
ef.adj <- ef
ef.adj$vectors$pvals <- p.adjust (ef$vectors$pvals, method = 'bonf')
en/history/2015-03-20-anadatr.txt · Last modified: 2018/03/30 23:04 (external edit)