| eval_util_R {occumb} | R Documentation |
Expected utility for regional species diversity assessments.
Description
eval_util_R() evaluates the expected utility of a regional
species diversity assessment using Monte Carlo integration.
Usage
eval_util_R(
settings,
fit = NULL,
psi = NULL,
theta = NULL,
phi = NULL,
N_rep = 1,
cores = 1L
)
Arguments
settings |
A data frame that specifies a set of conditions under which
utility is evaluated. It must include columns named |
fit |
An |
psi |
Sample values of the site occupancy probabilities of species
stored in a matrix with sample |
theta |
Sample values of sequence capture probabilities of species
stored in a matrix with sample |
phi |
Sample values of sequence relative dominance of species stored in
a matrix with sample |
N_rep |
Controls the sample size for the Monte Carlo integration.
The integral is evaluated using a total of |
cores |
The number of cores to use for parallelization. |
Details
The utility of a regional species diversity assessment can be defined as
the number of species expected to be detected in the region of interest
(Fukaya et al. 2022). eval_util_R() evaluates this utility for the region
modeled in the occumbFit object for the combination of J, K, and N
values specified in the conditions argument.
Such evaluations can be used to balance J, K, and N to maximize the
utility under a constant budget (possible combinations of J, K, and N
under a specified budget and cost values are easily obtained using
list_cond_R(); see the example below).
It is also possible to examine how the utility varies with different J,
K, and N values without setting a budget level, which may be useful in determining
the satisfactory levels of J, K, and N from a purely technical point of
view.
The expected utility is defined as the expected value of the conditional utility in the form:
U(J, K, N \mid \boldsymbol{r}, \boldsymbol{u}) = \sum_{i = 1}^{I}\left\{1 - \prod_{j = 1}^{J}\prod_{k = 1}^{K}\left(1 - \frac{u_{ijk}r_{ijk}}{\sum_{m = 1}^{I}u_{mjk}r_{mjk}} \right)^N \right\}
where u_{ijk} is a latent indicator variable representing
the inclusion of the sequence of species i in replicate k
at site j, and r_{ijk} is a latent variable that
is proportional to the relative frequency of the sequence of species
i, conditional on its presence in replicate k at site
j (Fukaya et al. 2022).
Expectations are taken with respect to the posterior (or possibly prior)
predictive distributions of \boldsymbol{r} = \{r_{ijk}\} and
\boldsymbol{u} = \{u_{ijk}\}, which are evaluated numerically using
Monte Carlo integration. The predictive distributions of
\boldsymbol{r} and \boldsymbol{u} depend on the model
parameters \psi, \theta, and \phi values.
Their posterior (or prior) distribution is specified by supplying an
occumbFit object containing their posterior samples via the fit argument,
or by supplying a matrix or array of posterior (or prior) samples of
parameter values via the psi, theta, and phi arguments. Higher
approximation accuracy can be obtained by increasing the value of N_rep.
The eval_util_R() function can be executed by supplying the fit argument
without specifying the psi, theta, and phi arguments, by supplying the
three psi, theta, and phi arguments without the fit argument, or by
supplying the fit argument and any or all of the psi, theta, and phi
arguments. If the psi, theta, or phi arguments are specified in addition
to the fit, the parameter values given in these arguments are preferentially
used to evaluate the expected utility. If the sample sizes differed among
parameters, parameters with smaller sample sizes are resampled with
replacements to align the sample sizes across parameters.
The expected utility is evaluated assuming homogeneity of replicates, in the
sense that \theta and \phi, the model parameters
associated with the species detection process, are constant across
replicates within a site. For this reason, eval_util_R() does not accept
replicate-specific \theta and \phi. If the
occumbFit object supplied in the fit argument has a replicate-specific
parameter, the parameter samples to be used in the utility evaluation must be
provided explicitly via the theta or phi arguments.
If the parameters are modeled as a function of site covariates in the fit
object, or if the psi, theta, and/or phi arguments have site dimensions,
the expected utility is evaluated to account for the site heterogeneity of
the parameters. To incorporate site heterogeneity, the
parameter values for each J site are determined by selecting site-specific
parameter values in the fit, or those supplied in psi, theta, and phi
via random sampling with replacement. Thus, expected utility is
evaluated by assuming a set of supplied parameter values as a statistical
population of site-specific parameters.
The Monte Carlo integration is executed in parallel on multiple CPU cores, where
the cores argument controls the degree of parallelization.
Value
A data frame with a column named Utility in which the estimates of
the expected utility are stored. This is obtained by adding the Utility column
to the data frame provided in the settings argument.
References
K. Fukaya, N. I. Kondo, S. S. Matsuzaki and T. Kadoya (2022) Multispecies site occupancy modelling and study design for spatially replicated environmental DNA metabarcoding. Methods in Ecology and Evolution 13:183–193. doi:10.1111/2041-210X.13732
Examples
set.seed(1)
# Generate a random dataset (20 species * 2 sites * 2 reps)
I <- 20 # Number of species
J <- 2 # Number of sites
K <- 2 # Number of replicates
data <- occumbData(
y = array(sample.int(I * J * K), dim = c(I, J, K)))
# Fitting a null model
fit <- occumb(data = data)
## Estimate expected utility
# Arbitrary J, K, and N values
(util1 <- eval_util_R(expand.grid(J = 1:3, K = 1:3, N = c(1E3, 1E4, 1E5)),
fit))
# J, K, and N values under specified budget and cost
(util2 <- eval_util_R(list_cond_R(budget = 50000,
lambda1 = 0.01,
lambda2 = 5000,
lambda3 = 5000),
fit))
# K values restricted
(util3 <- eval_util_R(list_cond_R(budget = 50000,
lambda1 = 0.01,
lambda2 = 5000,
lambda3 = 5000,
K = 1:5),
fit))
# J and K values restricted
(util4 <- eval_util_R(list_cond_R(budget = 50000,
lambda1 = 0.01,
lambda2 = 5000,
lambda3 = 5000,
J = 1:3, K = 1:5),
fit))
# theta and phi values supplied
(util5 <- eval_util_R(list_cond_R(budget = 50000,
lambda1 = 0.01,
lambda2 = 5000,
lambda3 = 5000,
J = 1:3, K = 1:5),
fit,
theta = array(0.5, dim = c(4000, I, J)),
phi = array(1, dim = c(4000, I, J))))
# psi, theta, and phi values, but no fit object supplied
(util6 <- eval_util_R(list_cond_R(budget = 50000,
lambda1 = 0.01,
lambda2 = 5000,
lambda3 = 5000,
J = 1:3, K = 1:5),
fit = NULL,
psi = array(0.9, dim = c(4000, I, J)),
theta = array(0.9, dim = c(4000, I, J)),
phi = array(1, dim = c(4000, I, J))))