eval_util_R: 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
)

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.

Arguments

settings: A data frame that specifies a set of conditions under which utility is evaluated. It must include columns named J, K, and N, which specify the number of sites, number of replicates per site, and sequencing depth per replicate, respectively. J, K, and N must be numeric vectors greater than 0. When J and K contain decimal values, they are discarded and treated as integers. Additional columns are ignored, but may be included.
fit: An occumbFit object.
psi: Sample values of the site occupancy probabilities of species stored in a matrix with sample $\times$ species dimensions or an array with sample $\times$ species $\times$ site dimensions.
theta: Sample values of sequence capture probabilities of species stored in a matrix with sample $\times$ species dimensions or an array with sample $\times$ species $\times$ site dimensions.
phi: Sample values of sequence relative dominance of species stored in a matrix with sample $\times$ species dimensions or an array with sample $\times$ species $\times$ site dimensions.
N_rep: Controls the sample size for the Monte Carlo integration. The integral is evaluated using a total of N_sample * N_rep random samples, where N_sample is the maximum size of the MCMC sample in the fit argument and the parameter sample in the psi, theta, and phi arguments.
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.

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. tools:::Rd_expr_doi("10.1111/2041-210X.13732")

Examples

Run this code

# \donttest{
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))))
# }

Run the code above in your browser using DataLab