Chapter 4 Use case 2: joint_model() with site-level covariates

This second use case uses the same goby data as in use case 1, except this time we will include site-level covariates that affect the sensitivity of eDNA relative to traditional surveys.



4.1 Site-level covariate data

library(eDNAjoint)
data(goby_data)

In addition to count and qPCR data, the goby data includes site-level covariates, which is optional when implementing joint_model(). Here, the data represent salinity, mean time to filter eDNA samples, density of other fish, habitat size, and vegetation presence at each site. Two important notes:

  1. Notice that the continuous covariate data is normalized. This is useful since this data will be used in a linear regression, and it helps algorithm stability for all covariate data to be on the same scale. Similarly, one should use dummy variables for categorical variables (like the “Veg” variable).

  2. The columns in the matrix should be named, since these identifiers will be used when fitting the model.

head(goby_data$site_cov)
##        Salinity Filter_time Other_fishes   Hab_size Veg
## [1,] -0.7114925       -1.17   -0.4738419 -0.2715560   0
## [2,] -0.2109183       -1.24   -0.4738419 -0.2663009   0
## [3,] -1.1602831       -1.29   -0.4738419 -0.2717707   0
## [4,] -0.5561419        0.11    0.5479118 -0.2164312   1
## [5,] -0.9876713       -0.70    0.2437353  4.9981956   1
## [6,]  1.2562818       -0.55   -0.3512823 -0.2934710   0

One way to normalize your covariate data:

\[ \frac{x-\mu}{\sigma} \]

cov_norm <- (cov - mean(cov)) / sd(cov)

For more data formatting guidance, see section 2.1.1.

4.2 Fit the model

Now that we understand our data, let’s fit the joint model. The key arguments of this function include:

  1. data: list of pcr_k, pcr_n, count, and site_cov matrices
  2. cov: character vector of site-level covariates (this model will only include mean eDNA sample filter time and salinity)
  3. family: probability distribution used to model the seine count data. A poisson distribution is chosen here.
  4. p10_priors: Beta distribution parameters for the prior on the probability of false positive eDNA detection, \(p_{10}\). c(1,20) is the default specification. More on this later.
  5. q: logical value indicating the presence of multiple traditional gear types. Here, we’re only using data from one traditional method.

More parameters exist to further customize the MCMC sampling, but we’ll stick with the defaults.

# run the joint model with two covariates
goby_fit_cov1 <- joint_model(data = goby_data,
                             cov = c("Filter_time", "Salinity"),
                             family = "poisson", p10_priors = c(1, 20),
                             q = FALSE, multicore = TRUE)

goby_fit_cov1 is a list containing:

  1. model fit (goby_fit_cov1$model) of the class ‘stanfit’ and can be accessed and interpreted using all functions in the rstan package.
  2. initial values used for each chain in MCMC (goby_fit_cov1$inits)

4.3 Model selection

We previously made a choice to include two site-level covariates. Perhaps we want to test how that model specification compares to a model specification with different site-level covariates.

# fit a new model with one site-level covariate
goby_fit_cov2 <- joint_model(data = goby_data, cov = "Veg",
                             family = "poisson", p10_priors = c(1, 20),
                             q = FALSE, multicore = TRUE)

We can now compare the fit of these model to our data using the joint_select() function, which performs leave-one-out cross validation with functions from the loo package.

# perform model selection
joint_select(model_fits = list(goby_fit_cov1$model, goby_fit_cov2$model))
##        elpd_diff se_diff
## model1   0.0       0.0  
## model2 -35.7      29.4

These results tell us that model1 has a higher Bayesian LOO estimate of the expected log pointwise predictive density (elpd_loo). This means that goby_fit_cov1 is likely a better fit to the data.

You could keep going with this further and include/exclude different covariates, or compare to a null model without covariates.

4.4 Interpret the output

4.4.1 Summarize posterior distributions

Let’s interpret goby_fit_cov1. Use joint_summarize() to see the posterior summaries of the model parameters.

joint_summarize(goby_fit_cov1$model, par = c("p10", "alpha"))
##            mean se_mean    sd   2.5%  97.5%     n_eff Rhat
## p10       0.003   0.000 0.001  0.001  0.007 14055.193    1
## alpha[1]  0.542   0.001 0.099  0.346  0.734  9325.393    1
## alpha[2]  1.020   0.001 0.118  0.786  1.249 10450.422    1
## alpha[3] -0.350   0.001 0.105 -0.553 -0.147 12519.196    1

This summarizes the mean, standard deviation (sd), and quantiles of the posterior estimates of \(p_{10}\) and \(\alpha\), as well as the effective sample size (n_eff) and \(\hat{R}\) (Rhat) for the parameters. More informative about effective sample size and \(\hat{R}\) can be found in the algorithm convergence section, but briefly, the \(\hat{R}\) value is the frequently used statistic for assessing model convergence. We typically want \(\hat{R}\) to be less than 1.05, so it looks like our model converged.

The mean estimated probability of a false positive eDNA detection is 0.001, and the 2.5% and 97.5% quantiles show the bounds of the 95% credibility interval (the equal tailed credibility interval, to be specific.)

The vector \(\alpha\) represents the regression covariates that scales this relationship (see model description for more). alpha[1] corresponds to the intercept of the regression with site-level covariates. alpha[2] corresponds to the regression coefficient associated with Filter_time, and alpha[3] corresponds to the regression coefficient associated with Salinity. Positive regression coefficients indicate an inverse relationship between the covariate and eDNA sensitivity. So here, for example, longer filter time means lower eDNA sensitivity, and higher salinity means higher eDNA sensitivity.

In this example, equation 4 in the model description would be:

\[ \beta_i = \alpha_1 + \alpha_2 \times filtertime_i + \alpha_3 \times salinity_i \]

The parameter \(\beta_i\) represents the site-specific sensitivity of eDNA relative to traditional sampling. We use the index \(i\) to reference the sites. For example, here are the first few \(\beta_i\):

head(joint_summarize(goby_fit_cov1$model, par = "beta"))
##           mean se_mean    sd   2.5%  97.5%     n_eff Rhat
## beta[1] -0.403   0.002 0.167 -0.735 -0.083  9126.907    1
## beta[2] -0.650   0.002 0.189 -1.028 -0.283  8720.117    1
## beta[3] -0.369   0.002 0.176 -0.714 -0.028  9935.716    1
## beta[4]  0.848   0.001 0.107  0.639  1.056 13135.879    1
## beta[5]  0.173   0.001 0.136 -0.097  0.433 10769.079    1
## beta[6] -0.459   0.002 0.222 -0.900 -0.030  9121.241    1

We can also use functions from the bayesplot package to examine the posterior distributions and chain convergence.

First let’s look at the posterior distribution for \(p_{10}\).

library(bayesplot)
# plot posterior distribution, highlighting median and 80% credibility interval
mcmc_areas(as.matrix(goby_fit_cov1$model), pars = "p10", prob = 0.8)

Next let’s look at chain convergence for \(p_{10}\) and \(\mu_{i=1,k=1}\).

# this will plot the MCMC chains for p10 and mu at site 1
mcmc_trace(rstan::extract(goby_fit_cov1$model, permuted = FALSE),
           pars = c("p10", "mu[1,1]"))

These trace plots show that the algorithm has converged. The chains are overlapping and stationary (i.e., are moving around the same mean and have a constant variance). See more about trace plots in the algorithm convergence section.

4.4.2 Effort necessary to detect presence

To further highlight the relative sensitivity of eDNA and traditional sampling, we can use detection_calculate() to find the units of survey effort necessary to detect presence of the species. Here, detecting presence refers to producing at least one true positive eDNA detection or catching at least one individual in a traditional survey.

This function is finding the median number of survey units necessary to detect species presence if the expected catch rate, \(\mu\) is 0.1, 0.5, or 1. The cov_val argument indicates the value of the covariates used for the prediction. Since the covariate data was normalized, c(0, 0) indicates that the prediction is made at the mean Filter_time and Salinity values. For example, this means that all \(\beta_i = \alpha_1\).

detection_calculate(goby_fit_cov1$model, mu = c(0.1, 0.5, 1),
                    cov_val = c(0, 0), probability = 0.9)
##       mu n_traditional n_eDNA
## [1,] 0.1            24     14
## [2,] 0.5             5      4
## [3,] 1.0             3      2

We can see that at the mean covariate values, it takes 14 eDNA samples or 24 seine samples to detect goby presence with 0.9 probability if the expected catch rate is 0.1.

Now let’s perform the same calculation under a condition where the Filter_time covariate value is 0.5 z-scores above the mean. This means that equation 4 in the model description would be:

\[ \beta_i = \alpha_1 + \alpha_2 \times 0.5 + \alpha_3 \times 0 \]

detection_calculate(goby_fit_cov1$model, mu = c(0.1, 0.5, 1),
                    cov_val = c(0.5, 0), probability = 0.9)
##       mu n_traditional n_eDNA
## [1,] 0.1            24     23
## [2,] 0.5             5      5
## [3,] 1.0             3      3

At sites with a longer eDNA sample filter time, it would now take 22 eDNA samples or 24 seine samples to detect goby presence if the expected catch rate is 0.1.

Let’s do the same for salinity. This means that equation 4 in the model description would be:

\[ \beta_i = \alpha_1 + \alpha_2 \times 0 + \alpha_3 \times 0.5 \]

detection_calculate(goby_fit_cov1$model, mu = c(0.1, 0.5, 1),
                    cov_val = c(0, 0.5), probability = 0.9)
##       mu n_traditional n_eDNA
## [1,] 0.1            24     12
## [2,] 0.5             5      3
## [3,] 1.0             3      2

At sites with higher salinity, it would now take 12 eDNA samples or 24 seine samples to detect goby presence if the expected catch rate is 0.1.

We can also plot these comparisons. mu_min and mu_max define the x-axis in the plot.

detection_plot(goby_fit_cov1$model, mu_min = 0.1, mu_max = 1,
               cov_val = c(0, 0), probability = 0.9)

4.4.3 Calculate \(\mu_{critical}\)

The probability of a true positive eDNA detection, \(p_{11}\), is a function of the expected catch rate, \(\mu\). Low values of \(\mu\) correspond to low probability of eDNA detection. Since the probability of a false-positive eDNA detection is non-zero, the probability of a false positive detection may be higher than the probability of a true positive detection at very low values of \(\mu\).

\(\mu_{critical}\) describes the value of \(\mu\) where the probability of a false positive eDNA detection equals the probability of a true positive eDNA detection. This value can be calculated using mu_critical(). Here, we will calculate this value at the mean covariate values.

mu_critical(goby_fit_cov1$model, cov_val = c(0, 0), ci = 0.9)
## $median
## [1] 0.005290696
## 
## $lower_ci
## Highest Density Interval: 1.77e-03
## 
## $upper_ci
## Highest Density Interval: 9.73e-03

This function calculates \(\mu_{critical}\) using the entire posterior distributions of parameters from the model, and ‘HDI’ corresponds to the 90% credibility interval calculated using the highest density interval.

4.5 Initial values

By default, eDNAjoint will provide initial values for parameters estimated by the model, but you can provide your own initial values if you prefer. Here is an example of providing initial values for parameters, mu,p10, and alpha, as an input in joint_model().

# set number of chains
n_chain <- 4

# number of sites
nsites <- dim(goby_data$count)[1]

# initial values should be a list of named lists
inits <- list()
for (i in 1:n_chain) {
  inits[[i]] <- list(
    # length should equal the number of sites for each chain
    mu = stats::runif(nsites, 0.01, 5),
    # length should equal 1 for each chain
    p10 = stats::runif(1, 0.0001, 0.08),
    # length should equal the number of covariates plus 1
    # (to account for intercept in regression)
    alpha = rep(0.1, length(c("Filter_time", "Salinity")) + 1)
  )
}
# now fit the model
fit_inits <- joint_model(data = goby_data, cov = c("Filter_time", "Salinity"),
                         initial_values = inits, multicore = TRUE)
# check to see the initial values that were used
fit_inits$inits
## $chain1
## $chain1$mu_trad
##  [1] 4.8131153 3.9205129 2.3668318 4.8441400 0.3487256 1.5676713 4.7750355 1.2615131 2.2986378 0.8719349
## [11] 3.4325538 1.9360174 3.9451822 0.7391554 1.0671090 4.6332140 0.2975017 2.2343365 2.3768214 0.5212608
## [21] 1.5667693 3.8062090 2.1946372 2.7909638 3.4398058 3.1814049 0.4085425 0.3239513 2.0851465 2.4491118
## [31] 3.0486861 2.6646518 2.0904671 0.5318556 4.7902815 3.1955694 1.6497569 1.3257332 2.4432931
## 
## $chain1$mu
##  [1] 4.8131153 3.9205129 2.3668318 4.8441400 0.3487256 1.5676713 4.7750355 1.2615131 2.2986378 0.8719349
## [11] 3.4325538 1.9360174 3.9451822 0.7391554 1.0671090 4.6332140 0.2975017 2.2343365 2.3768214 0.5212608
## [21] 1.5667693 3.8062090 2.1946372 2.7909638 3.4398058 3.1814049 0.4085425 0.3239513 2.0851465 2.4491118
## [31] 3.0486861 2.6646518 2.0904671 0.5318556 4.7902815 3.1955694 1.6497569 1.3257332 2.4432931
## 
## $chain1$log_p10
## [1] -2.931961
## 
## $chain1$alpha
## [1] 0.1 0.1 0.1
## 
## $chain1$p_dna
## numeric(0)
## 
## $chain1$p11_dna
## numeric(0)
## 
## 
## $chain2
## $chain2$mu_trad
##  [1] 4.2659825 1.9850733 3.1577655 2.7529505 1.3220470 2.1086788 1.2654867 1.1321412 1.4984807 4.6095900
## [11] 2.8755171 1.4442823 3.6007801 1.1281642 0.9067242 0.3457110 0.7240758 1.6459009 1.5099618 3.3451803
## [21] 2.4693550 1.8407941 1.6213274 0.1578897 0.8476885 4.0830482 2.9399878 1.7770375 4.2939067 3.1025464
## [31] 2.4222215 1.2931031 2.1368233 3.2147419 4.0403259 4.0932881 4.9916557 2.9197916 2.7862565
## 
## $chain2$mu
##  [1] 4.2659825 1.9850733 3.1577655 2.7529505 1.3220470 2.1086788 1.2654867 1.1321412 1.4984807 4.6095900
## [11] 2.8755171 1.4442823 3.6007801 1.1281642 0.9067242 0.3457110 0.7240758 1.6459009 1.5099618 3.3451803
## [21] 2.4693550 1.8407941 1.6213274 0.1578897 0.8476885 4.0830482 2.9399878 1.7770375 4.2939067 3.1025464
## [31] 2.4222215 1.2931031 2.1368233 3.2147419 4.0403259 4.0932881 4.9916557 2.9197916 2.7862565
## 
## $chain2$log_p10
## [1] -3.747031
## 
## $chain2$alpha
## [1] 0.1 0.1 0.1
## 
## $chain2$p_dna
## numeric(0)
## 
## $chain2$p11_dna
## numeric(0)
## 
## 
## $chain3
## $chain3$mu_trad
##  [1] 1.2809299 4.0004529 0.5510164 4.9904098 1.3035072 4.3937898 3.1200058 2.5404326 3.1511372 4.7709753
## [11] 2.8539551 3.9966173 3.0224519 3.4739038 1.3061269 0.8342674 1.3104262 2.0631685 1.7564850 3.2157833
## [21] 2.4151580 3.5268270 4.3664013 2.5888662 1.0070164 4.7733519 3.5538149 3.4249273 3.7902915 3.3015182
## [31] 1.3888095 1.6065453 2.7066245 1.2812336 4.8685010 2.2748162 0.3511626 0.5876060 1.0484282
## 
## $chain3$mu
##  [1] 1.2809299 4.0004529 0.5510164 4.9904098 1.3035072 4.3937898 3.1200058 2.5404326 3.1511372 4.7709753
## [11] 2.8539551 3.9966173 3.0224519 3.4739038 1.3061269 0.8342674 1.3104262 2.0631685 1.7564850 3.2157833
## [21] 2.4151580 3.5268270 4.3664013 2.5888662 1.0070164 4.7733519 3.5538149 3.4249273 3.7902915 3.3015182
## [31] 1.3888095 1.6065453 2.7066245 1.2812336 4.8685010 2.2748162 0.3511626 0.5876060 1.0484282
## 
## $chain3$log_p10
## [1] -2.622226
## 
## $chain3$alpha
## [1] 0.1 0.1 0.1
## 
## $chain3$p_dna
## numeric(0)
## 
## $chain3$p11_dna
## numeric(0)
## 
## 
## $chain4
## $chain4$mu_trad
##  [1] 2.3728858 1.8963169 0.6446509 0.3862684 1.0285025 3.0203035 2.6307789 1.1526760 3.3557736 3.5787119
## [11] 0.5883730 0.6293207 2.9198424 4.7622327 2.2674280 1.9950383 4.7033466 2.6880021 1.4285339 4.6340902
## [21] 4.3317243 2.3806921 0.9219941 2.1517698 3.2230118 4.3802262 3.8195084 3.5773529 2.3940610 4.8488464
## [31] 0.5748097 1.1642955 4.0765834 3.3803395 3.0418844 0.9352231 3.0276323 4.7296121 0.2322767
## 
## $chain4$mu
##  [1] 2.3728858 1.8963169 0.6446509 0.3862684 1.0285025 3.0203035 2.6307789 1.1526760 3.3557736 3.5787119
## [11] 0.5883730 0.6293207 2.9198424 4.7622327 2.2674280 1.9950383 4.7033466 2.6880021 1.4285339 4.6340902
## [21] 4.3317243 2.3806921 0.9219941 2.1517698 3.2230118 4.3802262 3.8195084 3.5773529 2.3940610 4.8488464
## [31] 0.5748097 1.1642955 4.0765834 3.3803395 3.0418844 0.9352231 3.0276323 4.7296121 0.2322767
## 
## $chain4$log_p10
## [1] -3.201502
## 
## $chain4$alpha
## [1] 0.1 0.1 0.1
## 
## $chain4$p_dna
## numeric(0)
## 
## $chain4$p11_dna
## numeric(0)