User Tools

Site Tools


en:rarefaction_examples_rscript

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

en:rarefaction_examples_rscript [2018/06/08 12:14] (current)
David Zelený created
Line 1: Line 1:
 +<code rsplus>
 +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'​)
 +
 +</​code>​
en/rarefaction_examples_rscript.txt · Last modified: 2018/06/08 12:14 by David Zelený