User Tools

Site Tools


en:rarefaction_examples_rscript
hp.abund <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/hp.abund.txt')
hp.incid <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/hp.incid.txt')
hp.elev <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/hp.elev.txt', row.names = 1) 
head (hp.abund)
 
 
# We will use package iNEXT (interpolation-extrapolation) developed by the lab of prof. Anne Chao (Tsinghua University, Hsinchu; http://chao.stat.nthu.edu.tw/wordpress/software_download/inext-online/)
# install.packages ('iNEXT')
library (iNEXT)
 
 
DataInfo (hp.abund, datatype = 'abundance')
DataInfo (hp.incid, datatype = 'incidence')
 
D_abund <- iNEXT (hp.abund, datatype = 'abundance')
plot (D_abund)
 
D_abund_232 <- iNEXT (hp.abund, datatype = 'abundance', endpoint = 232)
plot (D_abund_232)
 
est_D_abund <- estimateD (hp.abund, datatype = 'abundance', conf = NULL)
est_D_abund
 
D_individuals <- est_D_abund$`q = 0`
 
plot (D_abund, type = 3)
plot (D_abund, type = 3, xlim = c(.95, 1))
 
est_D_abund_coverage <- estimateD (hp.abund, datatype = 'abundance', base = 'coverage', conf = NULL)
est_D_abund_coverage
D_coverage <- est_D_abund_coverage$`q = 0`
D_area <- DataInfo (hp.abund, datatype = 'abundance')$S.obs
 
D_est <- cbind (D_area, D_individuals, D_coverage)
rownames (D_est) <- D_area <- DataInfo (hp.abund, datatype = 'abundance')$site
 
barplot ( (D_est), beside = T, legend.text = T, col = heat.colors (7), xlab = 'Sample standardisation', ylab = 'Species richness', names.arg = c('sample area', '# of individuals', 'sample coverage'))
 
matplot (scale (D_est), type = 'b', axes = F, xlab = 'Locality', ylab = 'Standardized species richness')
axis (1, at = 1:7, labels = rownames (D_est))
axis (2)
box ()
legend ('topright', title = 'Diversity standardised by:', legend = c('sample area', '# individuals', 'sample coverage'), pch = as.character (1:3), col = 1:3, lty = 1:3)
 
 
# The function iNEXT calculates rarefaction curves (can even extrapolate it beyond the observed number of samples/individuals)
D_incid <- iNEXT (hp.incid, datatype = 'incidence_freq')  # this is for incidence-based data in hp.incid
plot (D_incid)
 
D_abund <- iNEXT (hp.abund, datatype = 'abundance', endpoint = 464)  # and this for abundance-based data for hp.abund
D_abund$DataInfo
plot (D_abund)
 
D_abund_est <- estimateD (hp.abund, datatype = 'abundance', conf = NULL)   # function estimateD returns estimated diversities from given dataset (in the default setting, the estimate is done to the lowest number of samples/individuals in the dataset, here 232 individuals as in locality FT)
plot (D_abund_est[D_abund_est$order == 0, 'qD'] ~ elev$Elevation)  # draw pattern of diversity along elevation
 
# calculate rarefaction for completeness/coverage
D_abund <- iNEXT (hp.abund, datatype = 'abundance')
plot (D_abund, type = 2, ylim = c(0.95, 1))
plot (D_abund, type = 3)
abline (v = 0.9658)
 
D_abund_est_coverage <- estimateD (hp.abund, datatype = 'abundance', base = 'coverage', conf = NULL)
D_abund_est_coverage
 
D_area <- D_abund$DataInfo$S.obs
D_indi <- D_abund_est$`q = 0`
D_cove <- D_abund_est_coverage$`q = 0`
 
Div <- cbind (D_area, D_indi, D_cove)
rownames (Div) <- D_abund$DataInfo$site
Div
matplot (t(scale (Div)), type = 'b')
en/rarefaction_examples_rscript.txt · Last modified: 2018/06/08 12:14 by David Zelený