Skip to contents

This function uses the full posterior distributions of parameters estimated by jointModel() to calculate mu_critical, or the expected catch rate at which the probabilities of a false positive eDNA detection and true positive eDNA detection are equal. See more examples in the Package Vignette.

Usage

muCritical(modelfit, cov.val = NULL, ci = 0.9)

Arguments

modelfit

An object of class stanfit

cov.val

A numeric vector indicating the values of site-level covariates to use for prediction. Default is NULL.

ci

Credible interval calculated using highest density interval (HDI). Default is 0.9 (i.e., 90% credible interval).

Value

A list with median mu_critical and lower and upper bounds on the credible interval. If multiple gear types are used, a table of mu_critical and lower and upper credible interval bounds is returned with one column for each gear type.

Note

Before fitting the model, this function checks to ensure that the function is possible given the inputs. These checks include:

  • Input model fit is an object of class 'stanfit'.

  • Input credible interval is a univariate numeric value greater than 0 and less than 1.

  • Input model fit contains p10 parameter.

  • If model fit contains alpha, cov.val must be provided.

  • Input cov.val is numeric.

  • Input cov.val is the same length as the number of estimated covariates.

  • Input model fit has converged (i.e. no divergent transitions after warm-up).

If any of these checks fail, the function returns an error message.

Examples

# \donttest{
# Ex. 1: Calculating mu_critical with site-level covariates

# Load data
data(gobyData)

# Fit a model including 'Filter_time' and 'Salinity' site-level covariates
fit.cov <- jointModel(data = gobyData, cov = c('Filter_time','Salinity'),
                      family = "poisson", p10priors = c(1,20), q = FALSE,
                      multicore = FALSE)
#> 
#> SAMPLING FOR MODEL 'joint_binary_cov_pois' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 4.9e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.49 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 3000 [  0%]  (Warmup)
#> Chain 1: Iteration:  500 / 3000 [ 16%]  (Warmup)
#> Chain 1: Iteration:  501 / 3000 [ 16%]  (Sampling)
#> Chain 1: Iteration: 1000 / 3000 [ 33%]  (Sampling)
#> Chain 1: Iteration: 1500 / 3000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 2000 / 3000 [ 66%]  (Sampling)
#> Chain 1: Iteration: 2500 / 3000 [ 83%]  (Sampling)
#> Chain 1: Iteration: 3000 / 3000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.734 seconds (Warm-up)
#> Chain 1:                1.708 seconds (Sampling)
#> Chain 1:                2.442 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'joint_binary_cov_pois' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 4.4e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 3000 [  0%]  (Warmup)
#> Chain 2: Iteration:  500 / 3000 [ 16%]  (Warmup)
#> Chain 2: Iteration:  501 / 3000 [ 16%]  (Sampling)
#> Chain 2: Iteration: 1000 / 3000 [ 33%]  (Sampling)
#> Chain 2: Iteration: 1500 / 3000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 2000 / 3000 [ 66%]  (Sampling)
#> Chain 2: Iteration: 2500 / 3000 [ 83%]  (Sampling)
#> Chain 2: Iteration: 3000 / 3000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 6.703 seconds (Warm-up)
#> Chain 2:                1.735 seconds (Sampling)
#> Chain 2:                8.438 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'joint_binary_cov_pois' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 4.6e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.46 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 3000 [  0%]  (Warmup)
#> Chain 3: Iteration:  500 / 3000 [ 16%]  (Warmup)
#> Chain 3: Iteration:  501 / 3000 [ 16%]  (Sampling)
#> Chain 3: Iteration: 1000 / 3000 [ 33%]  (Sampling)
#> Chain 3: Iteration: 1500 / 3000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 2000 / 3000 [ 66%]  (Sampling)
#> Chain 3: Iteration: 2500 / 3000 [ 83%]  (Sampling)
#> Chain 3: Iteration: 3000 / 3000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 0.751 seconds (Warm-up)
#> Chain 3:                1.754 seconds (Sampling)
#> Chain 3:                2.505 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'joint_binary_cov_pois' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 4.8e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.48 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 3000 [  0%]  (Warmup)
#> Chain 4: Iteration:  500 / 3000 [ 16%]  (Warmup)
#> Chain 4: Iteration:  501 / 3000 [ 16%]  (Sampling)
#> Chain 4: Iteration: 1000 / 3000 [ 33%]  (Sampling)
#> Chain 4: Iteration: 1500 / 3000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 2000 / 3000 [ 66%]  (Sampling)
#> Chain 4: Iteration: 2500 / 3000 [ 83%]  (Sampling)
#> Chain 4: Iteration: 3000 / 3000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 0.667 seconds (Warm-up)
#> Chain 4:                1.701 seconds (Sampling)
#> Chain 4:                2.368 seconds (Total)
#> Chain 4: 
#> Refer to the eDNAjoint guide for visualization tips:  https://ednajoint.netlify.app/tips#visualization-tips 

# Calculate mu_critical at the mean covariate values (covariates are
# standardized, so mean = 0)
muCritical(fit.cov$model, cov.val = c(0,0), ci = 0.9)
#> $median
#> [1] 0.005256258
#> 
#> $lower_ci
#> Highest Density Interval: 1.80e-03
#> 
#> $upper_ci
#> Highest Density Interval: 9.79e-03
#> 

# Calculate mu_critical at habitat size 0.5 z-scores greater than the mean
muCritical(fit.cov$model, cov.val = c(0,0.5), ci = 0.9)
#> $median
#> [1] 0.004406221
#> 
#> $lower_ci
#> Highest Density Interval: 1.39e-03
#> 
#> $upper_ci
#> Highest Density Interval: 8.16e-03
#> 

# Ex. 2: Calculating mu_critical with multiple traditional gear types

# Load data
data(greencrabData)

# Fit a model with no site-level covariates
fit.q <- jointModel(data = greencrabData, cov = NULL, family = "negbin",
                    p10priors = c(1,20), q = TRUE, multicore = FALSE)
#> 
#> SAMPLING FOR MODEL 'joint_binary_catchability_negbin' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000286 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.86 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 3000 [  0%]  (Warmup)
#> Chain 1: Iteration:  500 / 3000 [ 16%]  (Warmup)
#> Chain 1: Iteration:  501 / 3000 [ 16%]  (Sampling)
#> Chain 1: Iteration: 1000 / 3000 [ 33%]  (Sampling)
#> Chain 1: Iteration: 1500 / 3000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 2000 / 3000 [ 66%]  (Sampling)
#> Chain 1: Iteration: 2500 / 3000 [ 83%]  (Sampling)
#> Chain 1: Iteration: 3000 / 3000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 3.672 seconds (Warm-up)
#> Chain 1:                11.653 seconds (Sampling)
#> Chain 1:                15.325 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'joint_binary_catchability_negbin' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 0.000334 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.34 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 3000 [  0%]  (Warmup)
#> Chain 2: Iteration:  500 / 3000 [ 16%]  (Warmup)
#> Chain 2: Iteration:  501 / 3000 [ 16%]  (Sampling)
#> Chain 2: Iteration: 1000 / 3000 [ 33%]  (Sampling)
#> Chain 2: Iteration: 1500 / 3000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 2000 / 3000 [ 66%]  (Sampling)
#> Chain 2: Iteration: 2500 / 3000 [ 83%]  (Sampling)
#> Chain 2: Iteration: 3000 / 3000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 3.319 seconds (Warm-up)
#> Chain 2:                11.503 seconds (Sampling)
#> Chain 2:                14.822 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'joint_binary_catchability_negbin' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 0.000327 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.27 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 3000 [  0%]  (Warmup)
#> Chain 3: Iteration:  500 / 3000 [ 16%]  (Warmup)
#> Chain 3: Iteration:  501 / 3000 [ 16%]  (Sampling)
#> Chain 3: Iteration: 1000 / 3000 [ 33%]  (Sampling)
#> Chain 3: Iteration: 1500 / 3000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 2000 / 3000 [ 66%]  (Sampling)
#> Chain 3: Iteration: 2500 / 3000 [ 83%]  (Sampling)
#> Chain 3: Iteration: 3000 / 3000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 3.415 seconds (Warm-up)
#> Chain 3:                10.946 seconds (Sampling)
#> Chain 3:                14.361 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'joint_binary_catchability_negbin' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 0.000338 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 3.38 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 3000 [  0%]  (Warmup)
#> Chain 4: Iteration:  500 / 3000 [ 16%]  (Warmup)
#> Chain 4: Iteration:  501 / 3000 [ 16%]  (Sampling)
#> Chain 4: Iteration: 1000 / 3000 [ 33%]  (Sampling)
#> Chain 4: Iteration: 1500 / 3000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 2000 / 3000 [ 66%]  (Sampling)
#> Chain 4: Iteration: 2500 / 3000 [ 83%]  (Sampling)
#> Chain 4: Iteration: 3000 / 3000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 3.446 seconds (Warm-up)
#> Chain 4:                8.753 seconds (Sampling)
#> Chain 4:                12.199 seconds (Total)
#> Chain 4: 
#> Refer to the eDNAjoint guide for visualization tips:  https://ednajoint.netlify.app/tips#visualization-tips 

# Calculate mu_critical
muCritical(fit.q$model, cov.val = NULL, ci = 0.9)
#>               gear_1      gear_2
#> median   0.059070936 0.046617182
#> lower_ci 0.009265172 0.006740918
#> upper_ci 0.135199672 0.104594860
# }