1+ # ' Calculate expected sample coverage C_hat
2+ # '
3+ # ' Returns expected sample coverage of a sample `x` for a smaller than observed
4+ # ' sample size `m` (Chao & Jost, 2012). This code was copied from INEXT's internal
5+ # ' function \code{iNEXT::Chat.Ind} (Hsieh et al 2016).
6+ # '
7+ # ' @param x integer vector (species abundances)
8+ # ' @param m integer a number of individuals that is smaller than observed total
9+ # ' community abundance.
10+ # '
11+ # ' @return a numeric value that is the expected coverage.
12+ # '
13+ # ' @references
14+ # ' Chao, A., and L. Jost. 2012. Coverage-based rarefaction and extrapolation:
15+ # ' standardizing samples by completeness rather than size. Ecology 93:2533–2547.
16+ # '
17+ # ' Anne Chao, Nicholas J. Gotelli, T. C. Hsieh, Elizabeth L. Sander, K. H. Ma,
18+ # ' Robert K. Colwell, and Aaron M. Ellison 2014. Rarefaction and extrapolation
19+ # ' with Hill numbers: a framework for sampling and estimation in species
20+ # ' diversity studies. Ecological Monographs 84:45-67.
21+ # '
22+ # ' T. C. Hsieh, K. H. Ma and Anne Chao. 2024.
23+ # ' iNEXT: iNterpolation and EXTrapolation for
24+ # ' species diversity. R package version 3.0.1
25+ # ' URL: http://chao.stat.nthu.edu.tw/wordpress/software-download/.
26+ # '
27+ # '
28+ # ' @export
29+ # '
30+ # ' @examples
31+ # ' data(inv_comm)
32+ # ' # What is the expected coverage at a sample size of 50 at the gamma scale?
33+ # ' Chat(colSums(inv_comm), 50)
34+ Chat <- function (x , m )
35+ {
36+ x <- x [x > 0 ]
37+ n <- sum(x )
38+ f1 <- sum(x == 1 )
39+ f2 <- sum(x == 2 )
40+ f0.hat <- ifelse(f2 == 0 , (n - 1 ) / n * f1 * (f1 - 1 ) / 2 , (n -
41+ 1 ) / n * f1 ^
42+ 2 / 2 / f2 )
43+ A <- ifelse(f1 > 0 , n * f0.hat / (n * f0.hat + f1 ), 1 )
44+ Sub <- function (m ) {
45+ if (m < n ) {
46+ xx <- x [(n - x ) > = m ]
47+ out <- 1 - sum(xx / n * exp(
48+ lgamma(n - xx + 1 ) - lgamma(n -
49+ xx - m + 1 ) - lgamma(n ) + lgamma(n - m )
50+ ))
51+ }
52+ if (m == n )
53+ out <- 1 - f1 / n * A
54+ if (m > n )
55+ out <- 1 - f1 / n * A ^ (m - n + 1 )
56+ out
57+ }
58+ sapply(m , Sub )
59+ }
60+
61+ # ' Number of individuals corresponding to a desired coverage (inverse C_hat)
62+ # '
63+ # ' If you wanted to resample a vector to a certain expected sample coverage, how
64+ # ' many individuals would you have to draw? This is C_hat solved for the number
65+ # ' of individuals. This code is a modification of INEXT's internal function
66+ # ' `invChat.Ind` (Hsieh et al 2016).
67+ # '
68+ # ' @param x integer vector (species abundances)
69+ # ' @param C coverage value between 0 and 1
70+ # '
71+ # ' @return a numeric value which is the number of individuals for a given
72+ # ' level of coverage \code{C}.
73+ # ' @references
74+ # ' Chao, A., and L. Jost. 2012. Coverage-based rarefaction and extrapolation:
75+ # ' standardizing samples by completeness rather than size. Ecology 93:2533–2547.
76+ # '
77+ # ' Anne Chao, Nicholas J. Gotelli, T. C. Hsieh, Elizabeth L. Sander, K. H. Ma,
78+ # ' Robert K. Colwell, and Aaron M. Ellison 2014. Rarefaction and extrapolation
79+ # ' with Hill numbers: a framework for sampling and estimation in species
80+ # ' diversity studies. Ecological Monographs 84:45-67.
81+ # '
82+ # ' T. C. Hsieh, K. H. Ma and Anne Chao. 2024.
83+ # ' iNEXT: iNterpolation and EXTrapolation for
84+ # ' species diversity. R package version 3.0.1
85+ # ' URL: http://chao.stat.nthu.edu.tw/wordpress/software-download/.
86+ # ' @seealso \code{\link{calc_S_C}}
87+ # ' @export
88+ # ' @importFrom stats optimize
89+ # ' @examples
90+ # ' data(inv_comm)
91+ # ' # What sample size corresponds to an expected sample coverage of 55%?
92+ # ' invChat(colSums(inv_comm), 0.55)
93+ # '
94+ invChat <- function (x , C )
95+ {
96+ m <- NULL
97+ n <- sum(x )
98+ refC <- Chat(x , n )
99+ f <- function (m , C )
100+ abs(Chat(x , m ) - C )
101+ # for interpolation
102+ if (refC > C ) {
103+ opt <- stats :: optimize(f ,
104+ C = C ,
105+ lower = 0 ,
106+ upper = sum(x ))
107+ mm <- opt $ minimum
108+ }
109+ # for extrapolation
110+ if (refC < = C ) {
111+ f1 <- sum(x == 1 )
112+ f2 <- sum(x == 2 )
113+ if (f1 > 0 & f2 > 0 ) {
114+ A <- (n - 1 ) * f1 / ((n - 1 ) * f1 + 2 * f2 )
115+ }
116+ if (f1 > 1 & f2 == 0 ) {
117+ A <- (n - 1 ) * (f1 - 1 ) / ((n - 1 ) * (f1 - 1 ) + 2 )
118+ }
119+ if (f1 == 1 & f2 == 0 ) {
120+ A <- 1
121+ }
122+ if (f1 == 0 & f2 == 0 ) {
123+ A <- 1
124+ }
125+ mm <- (log(n / f1 ) + log(1 - C )) / log(A ) - 1
126+ mm <- n + mm
127+
128+ }
129+ if (mm > 2 * n )
130+ warning(
131+ " The maximum size of the extrapolation exceeds double reference sample size, the results for q = 0 may be subject to large prediction bias."
132+ )
133+ return (mm )
134+ }
135+
136+
137+ # ' Calculate species richness for a given coverage level.
138+ # '
139+ # ' This function uses coverage-based rarefaction to compute species richness.
140+ # ' Specifically, the metric is computed as the
141+ # '
142+ # ' @param x a site by species matrix or a species abundance distribution
143+ # ' @param C_target target coverage between 0 and 1 (default is NULL). If not
144+ # ' provided then target coverage is computed by \code{\link{calc_C_target}}
145+ # ' @param extrapolate logical. Defaults to TRUE in which case richness is
146+ # ' extrapolated to sample sizes larger than observed in the dataset.
147+ # ' @param interrupt logical. Should the function throw an error when \code{C_target}
148+ # ' exceeds the maximum recommendable coverage?
149+ # '
150+ # ' @returns numeric value which is the species richness at a specific level of
151+ # ' coverage.
152+ # ' @references
153+ # ' Chao, A., and L. Jost. 2012. Coverage-based rarefaction and extrapolation:
154+ # ' standardizing samples by completeness rather than size. Ecology 93:2533–2547.
155+ # '
156+ # ' Anne Chao, Nicholas J. Gotelli, T. C. Hsieh, Elizabeth L. Sander, K. H. Ma,
157+ # ' Robert K. Colwell, and Aaron M. Ellison 2014. Rarefaction and extrapolation
158+ # ' with Hill numbers: a framework for sampling and estimation in species
159+ # ' diversity studies. Ecological Monographs 84:45-67.
160+ # '
161+ # ' T. C. Hsieh, K. H. Ma and Anne Chao. 2024.
162+ # ' iNEXT: iNterpolation and EXTrapolation for
163+ # ' species diversity. R package version 3.0.1
164+ # ' URL: http://chao.stat.nthu.edu.tw/wordpress/software-download/.
165+ # '
166+ # ' @seealso \code{\link{invChat}}
167+ # ' @export
168+ # '
169+ # ' @examples
170+ # ' data(tank_comm)
171+ # ' # What is species richness for a coverage value of 60%?
172+ # ' calc_S_C(tank_comm, C_target = 0.6)
173+ calc_S_C <- function (x ,
174+ C_target = NULL ,
175+ extrapolate = TRUE ,
176+ interrupt = TRUE ) {
177+ x <- as.matrix(x )
178+ if (any(dim(x ) == 1 ))
179+ sad <- as.numeric(x )
180+ else
181+ sad <- colSums(x )
182+ if (is.null(C_target ))
183+ C_target <- calc_C_target(x )
184+ N <- round(invChat(sad , C_target ))
185+ C_max = calc_C_target(x , factor = ifelse(extrapolate , 2 , 1 ))
186+ if (C_target > C_max & interrupt ) {
187+ if (extrapolate ) {
188+ stop(
189+ paste0(
190+ " Coverage exceeds the maximum possible value recommendable for extrapolation (i.e. C_target = "
191+ round(C_max , 4 ),
192+ " ). Reduce the value of C_target."
193+ )
194+ )
195+ } else {
196+ stop(
197+ paste0(
198+ " Coverage exceeds the maximum possible value for interpolation (i.e. C_target = "
199+ round(C_max , 4 ),
200+ " ). Use extrapolation or reduce the value of C_target."
201+ )
202+ )
203+ }
204+ }
205+ if (N > 1 ) {
206+ S_C = rarefaction(x = sad ,
207+ method = " IBR"
208+ effort = N ,
209+ extrapolate = extrapolate ,
210+ quiet_mode = TRUE )
211+ } else {
212+ S_C = NA
213+ }
214+ attr(S_C , " C" = C_target
215+ attr(S_C , " N" = N
216+ return (S_C )
217+ }
218+
219+ # ' Calculate the recommended target coverage value for the computation of beta_C
220+ # '
221+ # ' Returns the estimated gamma-scale coverage that corresponds to the largest
222+ # ' allowable sample size (i.e. the smallest observed sample size at the alpha
223+ # ' scale multiplied by an extrapolation factor). The default (factor = 2) allows
224+ # ' for extrapolation up to 2 times the observed sample size of the smallest
225+ # ' alpha sample. For factor= 1, only interpolation is applied. Factors larger
226+ # ' than 2 are not recommended.
227+ # '
228+ # ' @param x a site by species abundance matrix
229+ # ' @param factor numeric. A multiplier for how much larger than total community
230+ # ' abundance to extrapolate to. Defaults to 2.
231+ # '
232+ # ' @return numeric value
233+ # ' @export
234+ # '
235+ # ' @examples
236+ # ' data(tank_comm)
237+ # '
238+ # ' # What is the largest possible C that I can use to calculate beta_C
239+ # ' calc_C_target(tank_comm)
240+ calc_C_target <- function (x , factor = 2 ) {
241+ x <- as.matrix(x )
242+ if (any(dim(x ) == 1 )) {
243+ n <- factor * sum(x )
244+ C_target <- Chat(x , n )
245+ }
246+ else {
247+ n <- min(factor * rowSums(x ))
248+ C_target <- Chat(colSums(x ), n )
249+ }
250+ return (C_target )
251+ }
0 commit comments