library(ggplot2)
te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f",
clay = "#b5534e", gold = "#cda23f", line = "#dad9ca",
ink = "#16241d", paper = "#f5f4ee")
theme_te <- function() {
theme_minimal(base_size = 12) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = "#e7e6dc"),
plot.background = element_rect(fill = "#f5f4ee", colour = NA),
panel.background = element_rect(fill = "#f5f4ee", colour = NA),
plot.title = element_text(face = "bold", colour = te_pal$ink),
axis.title = element_text(colour = "#2c3a31"))
}Sampling design for eDNA surveys
R
eDNA
occupancy
survey design
ecology tutorial
ggplot2
Split an eDNA budget between sites, water samples and PCR replicates in R, with exact asymptotic precision for occupancy and a measured optimal design.
Corrected 15 October 2026: the post explained the shift towards more sites for a rare species by the binomial term of the variance being small; it is the false-absence term that is small, and the binomial term, which only more sites can shrink, carries most of the variance, as the post now says.
Corrected 7 October 2026: the post said that a survey maximising the chance of at least one detection should use fewer sites sampled harder; in this model that objective favours more sites sampled more lightly, and the post now says so. It also said that at a low sample positivity the best design needed ten bottles per site; ten is the most the grid allows, and the text now compares with the best design on the grid.
Updated 25 September 2026: the post assumed a clean assay. A new section, When the assay is not clean, puts replicate-level contamination back into the design calculation: the two-replicate optimum loses precision once contamination is possible, and three replicates hedge across the rates tried.
An eDNA survey has three numbers in it and one budget to cover all three: how many sites to visit, how many water samples to filter at each site, how many PCR replicates to run from each sample. They are not interchangeable. A site visit costs travel and a day, a sample costs a filter and an extraction, a replicate costs one well on a plate, and each level buys a different kind of certainty.
The two-level version of this problem is already solved in How many visits? Occupancy survey design, which is the prerequisite here: a fixed budget of site visits, exact asymptotic standard errors from the Fisher information, and an answer for how many visits a site needs. This tutorial does the three-level version with a different unit cost at each level, computes the precision every affordable design delivers, and then spends the money.
The likelihood, rebuilt
A site is occupied with probability psi. Inside an occupied site, each water sample carries detectable DNA with probability theta. Inside such a sample, each PCR replicate comes up positive with probability p. Nothing else is ever positive, so this post assumes the contamination problem away and works with a clean assay; a later section puts contamination back.
The data from one site are the positive-replicate counts, one count per sample, each between zero and K. The distribution of a single count under occupancy is a zero-inflated binomial: with probability theta the sample holds DNA and the count is binomial, otherwise the count is zero for certain.
sample_probs <- function(theta, p, n_rep) {
counts <- 0:n_rep
theta * dbinom(counts, n_rep, p) + (1 - theta) * (counts == 0)
}
psi_true <- 0.4
theta_true <- 0.5
p_true <- 0.7
cat("design parameters: psi", psi_true, " theta", theta_true, " p", p_true, "\n")design parameters: psi 0.4 theta 0.5 p 0.7 cat("per-sample outcome probabilities at two replicates:",
round(sample_probs(theta_true, p_true, 2), 4), "\n")per-sample outcome probabilities at two replicates: 0.545 0.21 0.245 cat("per-sample detection probability, one to four replicates:",
round(theta_true * (1 - (1 - p_true)^(1:4)), 4), "\n")per-sample detection probability, one to four replicates: 0.35 0.455 0.4865 0.496 At two replicates a sample from an occupied site already tests positive with probability 0.455 out of a ceiling of 0.5, and a third replicate raises that only to 0.4865. Those numbers are the whole argument of this post in miniature, and the rest of it is about what else the same money could have bought.
Expected information, without simulating
The asymptotic standard error of psi comes from the expected Fisher information, and for this model it can be written down rather than simulated. Given occupancy, the J samples of a site are independent draws from the outcome vector of length K + 1 computed above, so a site is a multinomial draw over K + 1 categories. The log-likelihood of any site with at least one positive replicate is linear in those category counts, which means every second moment of the score reduces to a sum over K + 1 terms instead of an enumeration of outcomes. The all-zero site is the only awkward case, because it is the one outcome that mixes the occupied and the unoccupied branch, so it gets its own term.
info_site <- function(psi, theta, p, n_samp, n_rep) {
counts <- 0:n_rep
bin <- dbinom(counts, n_rep, p)
qc <- theta * bin + (1 - theta) * (counts == 0)
dq_theta <- bin - (counts == 0)
dq_p <- theta * bin * (counts / p - (n_rep - counts) / (1 - p))
ratio_theta <- dq_theta / qc
ratio_p <- dq_p / qc
full <- matrix(0, 3, 3)
full[1, 1] <- 1 / psi^2
full[2, 2] <- n_samp * sum(qc * ratio_theta^2)
full[3, 3] <- n_samp * sum(qc * ratio_p^2)
full[2, 3] <- full[3, 2] <- n_samp * sum(qc * ratio_theta * ratio_p)
score_zero <- c(1 / psi, n_samp * ratio_theta[1], n_samp * ratio_p[1])
prob_zero <- psi * qc[1]^n_samp + 1 - psi
grad_zero <- c(qc[1]^n_samp - 1,
psi * n_samp * qc[1]^(n_samp - 1) * dq_theta[1],
psi * n_samp * qc[1]^(n_samp - 1) * dq_p[1])
psi * (full - qc[1]^n_samp * outer(score_zero, score_zero)) +
outer(grad_zero, grad_zero) / prob_zero
}
se_psi <- function(psi, theta, p, n_sites, n_samp, n_rep) {
info <- n_sites * info_site(psi, theta, p, n_samp, n_rep)
if (1 / kappa(info) < 1e-10) return(NA_real_)
sqrt(solve(info)[1, 1])
}
cat("reciprocal condition number of the information matrix\n")reciprocal condition number of the information matrixcat(" one sample per site:",
format(1 / kappa(info_site(psi_true, theta_true, p_true, 1, 3)), digits = 3),
" one replicate per sample:",
format(1 / kappa(info_site(psi_true, theta_true, p_true, 4, 1)), digits = 3),
"\n") one sample per site: 1.35e-17 one replicate per sample: 6.08e-17 cat(" four samples, two replicates:",
round(1 / kappa(info_site(psi_true, theta_true, p_true, 4, 2)), 4), "\n") four samples, two replicates: 0.3457 Two structural facts arrive before any design work. With one sample per site the data depend on psi and theta only through their product; with one replicate per sample the same happens to theta and p. In both cases the information matrix is singular to machine precision, so no amount of budget rescues the design and the grid below starts at two samples and two replicates.
Now the check that makes the formula usable. Simulate data at a design, fit the model by maximum likelihood, and compare the spread of the estimates against what the formula promised.
sim_counts <- function(psi, theta, p, n_sites, n_samp, n_rep) {
occupied <- rbinom(n_sites, 1, psi)
has_dna <- rbinom(n_sites * n_samp, 1, theta * rep(occupied, n_samp))
matrix(rbinom(n_sites * n_samp, n_rep, p * has_dna), n_sites, n_samp)
}
neg_loglik <- function(par, ymat, n_rep) {
qc <- sample_probs(plogis(par[2]), plogis(par[3]), n_rep)
n_samp <- ncol(ymat)
per_site <- rowSums(matrix(log(qc)[ymat + 1], nrow(ymat), n_samp))
found <- rowSums(ymat) > 0
-(sum(log(plogis(par[1])) + per_site[found]) +
sum(!found) * log(plogis(par[1]) * qc[1]^n_samp + 1 - plogis(par[1])))
}
fit_psi <- function(ymat, n_rep) {
opt <- optim(c(0, 0, 0), neg_loglik, ymat = ymat, n_rep = n_rep,
method = "Nelder-Mead",
control = list(reltol = 1e-10, maxit = 800))
c(plogis(opt$par[1]), opt$convergence)
}
set.seed(908)
n_fits <- 300
validate <- NULL
for (cfg in list(c(58, 4, 2), c(52, 4, 3))) {
draws <- t(replicate(n_fits, fit_psi(
sim_counts(psi_true, theta_true, p_true, cfg[1], cfg[2], cfg[3]), cfg[3])))
validate <- rbind(validate, data.frame(
sites = cfg[1], samples = cfg[2], replicates = cfg[3],
mean_psi_hat = round(mean(draws[, 1]), 4),
simulated_sd = round(sd(draws[, 1]), 4),
asymptotic_se = round(se_psi(psi_true, theta_true, p_true,
cfg[1], cfg[2], cfg[3]), 4),
failures = sum(draws[, 2] != 0)))
}
validate$ratio <- round(validate$simulated_sd / validate$asymptotic_se, 4)
cat("maximum likelihood fits per design:", n_fits, "\n")maximum likelihood fits per design: 300 print(validate, row.names = FALSE) sites samples replicates mean_psi_hat simulated_sd asymptotic_se failures
58 4 2 0.4117 0.0716 0.0716 0
52 4 3 0.4092 0.0759 0.0736 0
ratio
1.0000
1.0312Three hundred maximum likelihood fits at each of two designs. At 58 sites with four samples and two replicates the formula says 0.0716 and the simulation gives 0.0716, a ratio of 1. At 52 sites with three replicates the formula says 0.0736 and the simulation gives 0.0759, a ratio of 1.0312. With 300 fits the Monte Carlo scatter on a ratio of standard deviations runs to a few per cent either way, so a gap of that size is not evidence against the asymptotic formula at this sample size. Both designs recover psi close to the true 0.4, and on this run both land slightly high (0.4117 and 0.4092). The likelihood is skewed and the estimator is not unbiased in finite samples, so a small offset in either direction is what 300 fits will show.
Spending a budget
The cost model has to be written down explicitly because everything below depends on it. A site visit costs 10, a water sample costs 4, a PCR replicate costs 1, in the same arbitrary units, and the total spend for S sites is S times the sum of the site cost and J samples each carrying its own replicates. Only the ratios matter, so substitute your own: if your sites are a boat charter apart and your plates are cheap, the site cost rises and the answer moves.
cost_site <- 10
cost_sample <- 4
cost_replicate <- 1
budget <- 2000
design_cost <- function(n_sites, n_samp, n_rep) {
n_sites * (cost_site + n_samp * (cost_sample + n_rep * cost_replicate))
}
design_grid <- function(psi, theta, p, budget, samp_max = 10, rep_max = 9) {
out <- NULL
for (j in 2:samp_max) {
for (k in 2:rep_max) {
n_sites <- floor(budget / design_cost(1, j, k))
if (n_sites < 2) next
out <- rbind(out, data.frame(samples = j, replicates = k,
sites = n_sites,
spent = design_cost(n_sites, j, k),
se = se_psi(psi, theta, p, n_sites, j, k)))
}
}
out <- out[order(out$se), ]
out$excess_pct <- round(100 * (out$se / out$se[1] - 1), 2)
out
}
base_grid <- design_grid(psi_true, theta_true, p_true, budget)
top_designs <- head(base_grid, 6)
top_designs$se <- round(top_designs$se, 4)
cat("unit costs: site", cost_site, " sample", cost_sample,
" replicate", cost_replicate, "\n")unit costs: site 10 sample 4 replicate 1 cat("budget", budget, " affordable designs on the grid:", nrow(base_grid), "\n")budget 2000 affordable designs on the grid: 72 print(top_designs, row.names = FALSE) samples replicates sites spent se excess_pct
4 2 58 1972 0.0716 0.00
5 2 50 2000 0.0729 1.82
4 3 52 1976 0.0736 2.82
3 2 71 1988 0.0741 3.49
3 3 64 1984 0.0743 3.74
5 3 44 1980 0.0766 6.97The winner is 58 sites, four samples per site and two replicates per sample, spending 1972 of the 2000 available for a standard error of 0.0716. The five designs behind it are worse by 1.82, 2.82, 3.49, 3.74 and 6.97 per cent, which is the result worth carrying away: the optimum is flat. Anything in the region of three to five samples and two to three replicates gives essentially the same precision, and a survey that lands there gives up at most 6.97 per cent by not solving this problem exactly, on a clean assay; the two-replicate end of that region is what contamination breaks (When the assay is not clean).
best_row <- base_grid[1, ]
ggplot(base_grid, aes(replicates, samples, fill = se)) +
geom_tile(colour = "#f5f4ee") +
geom_point(data = best_row, shape = 21, size = 6, stroke = 1.2,
fill = NA, colour = te_pal$clay) +
scale_fill_gradient(low = te_pal$paper, high = te_pal$forest,
name = "SE of psi") +
scale_x_continuous(breaks = 2:9) +
scale_y_continuous(breaks = 2:10) +
labs(title = "A flat optimum with a steep wall behind it",
subtitle = "sites set by the remaining budget; pale is better precision",
x = "PCR replicates per sample", y = "water samples per site") +
theme_te() +
theme(panel.grid = element_blank())
The tile plot shows why the folklore survives. There is a broad pale plateau and then a steep climb once replicates or samples eat the money that sites needed.
Replicates saturate
The saturation is arithmetic, not statistical. On a clean assay a replicate can only raise the chance of catching a sample that already holds DNA, and that chance is one minus the failure probability raised to the number of replicates, so it approaches its ceiling geometrically. Once it is close, extra wells buy almost nothing (with contamination they also separate contaminated samples from real ones, which is a later section’s point) while still costing money that could have been a site.
best_at_k <- function(psi, theta, p, budget, n_rep, samp_max = 12) {
vals <- sapply(2:samp_max, function(j) {
n_sites <- floor(budget / design_cost(1, j, n_rep))
if (n_sites < 2) return(NA_real_)
se_psi(psi, theta, p, n_sites, j, n_rep)
})
j_best <- (2:samp_max)[which.min(vals)]
data.frame(replicates = n_rep, samples = j_best,
sites = floor(budget / design_cost(1, j_best, n_rep)),
se = min(vals, na.rm = TRUE))
}
sat <- NULL
for (th in c(0.2, 0.8)) {
for (pp in c(0.4, 0.7)) {
rows <- do.call(rbind, lapply(2:9, function(k)
best_at_k(psi_true, th, pp, budget, k)))
rows$theta <- th
rows$p <- pp
sat <- rbind(sat, rows)
}
}
for (th in c(0.2, 0.8)) {
for (pp in c(0.4, 0.7)) {
sub_rows <- sat[abs(sat$theta - th) < 1e-9 & abs(sat$p - pp) < 1e-9, ]
k_best <- sub_rows$replicates[which.min(sub_rows$se)]
cat("theta", th, " p", pp, " SE at two to nine replicates:",
round(sub_rows$se, 5), "\n best at", k_best,
"replicates; nine replicates cost an extra",
round(100 * (sub_rows$se[sub_rows$replicates == 9] /
min(sub_rows$se) - 1), 1), "percent\n")
}
}theta 0.2 p 0.4 SE at two to nine replicates: 0.13004 0.1251 0.12599 0.12855 0.13254 0.13836 0.14444 0.14911
best at 3 replicates; nine replicates cost an extra 19.2 percent
theta 0.2 p 0.7 SE at two to nine replicates: 0.10913 0.11236 0.11842 0.12438 0.12951 0.13648 0.14349 0.14852
best at 2 replicates; nine replicates cost an extra 36.1 percent
theta 0.8 p 0.4 SE at two to nine replicates: 0.06815 0.06497 0.06449 0.06397 0.06499 0.06632 0.06816 0.06975
best at 5 replicates; nine replicates cost an extra 9 percent
theta 0.8 p 0.7 SE at two to nine replicates: 0.05733 0.05729 0.05928 0.06116 0.06339 0.06538 0.0676 0.06941
best at 3 replicates; nine replicates cost an extra 21.2 percentWith a weak assay (p = 0.4) the curve has a real minimum: three replicates at low sample positivity, five at high positivity. With a good assay (p = 0.7) the minimum has moved down to two or three replicates and the curve climbs from there. Running nine replicates instead of the best number costs from 9 to 36.1 per cent of precision depending on the row, all of it spent on wells that were nearly certain to agree with each other.
sat$theta_label <- factor(ifelse(sat$theta < 0.5,
"low sample positivity (theta = 0.2)",
"high sample positivity (theta = 0.8)"),
levels = c("low sample positivity (theta = 0.2)",
"high sample positivity (theta = 0.8)"))
sat$p_label <- ifelse(sat$p < 0.55, "weak assay (p = 0.4)",
"good assay (p = 0.7)")
minima <- do.call(rbind, lapply(split(sat, list(sat$theta, sat$p)),
function(z) z[which.min(z$se), ]))
ggplot(sat, aes(replicates, se, colour = p_label)) +
geom_line(linewidth = 0.9) +
geom_point(data = minima, size = 3) +
facet_wrap(~theta_label, scales = "free_y") +
scale_colour_manual(values = c(`weak assay (p = 0.4)` = te_pal$clay,
`good assay (p = 0.7)` = te_pal$forest),
name = NULL) +
scale_x_continuous(breaks = 2:9) +
labs(title = "Past a handful of replicates the money is wasted",
subtitle = "points mark the best number of replicates in each series",
x = "PCR replicates per sample", y = "SE of psi") +
theme_te() +
theme(legend.position = "top")
What moves the optimum
Sample positivity and occupancy both move the answer, and they move it in different directions.
move_tab <- NULL
for (th in c(0.2, 0.8)) {
top <- design_grid(psi_true, th, p_true, budget)[1, ]
move_tab <- rbind(move_tab, data.frame(varied = "theta", value = th,
samples = top$samples,
replicates = top$replicates,
sites = top$sites,
se = round(top$se, 4)))
}
for (ps in c(0.1, 0.8)) {
top <- design_grid(ps, theta_true, p_true, budget)[1, ]
move_tab <- rbind(move_tab, data.frame(varied = "psi", value = ps,
samples = top$samples,
replicates = top$replicates,
sites = top$sites,
se = round(top$se, 4)))
}
move_tab$site_detection <- round(
1 - (1 - ifelse(move_tab$varied == "theta", move_tab$value, theta_true) *
(1 - (1 - p_true)^move_tab$replicates))^move_tab$samples, 4)
print(move_tab, row.names = FALSE) varied value samples replicates sites se site_detection
theta 0.2 9 2 31 0.1091 0.8360
theta 0.8 2 3 83 0.0573 0.9509
psi 0.1 3 2 71 0.0424 0.8381
psi 0.8 5 2 50 0.0650 0.9519Sample positivity dominates. At theta = 0.2 the optimum is 31 sites with nine samples each and two replicates; at theta = 0.8 it is 83 sites with two samples each and three replicates. When most water samples from an occupied site contain nothing, the only way to learn that the site is occupied is to take more bottles, and the plate cannot help: a replicate cannot find DNA that is not in the tube. This is the sharpest practical difference between eDNA design and the two-level visit problem, where there is no such level to be starved.
Occupancy moves the optimum much less, and towards more effort per site as the species gets commoner: three samples and 71 sites at psi = 0.1, five samples and 50 sites at psi = 0.8, with the per-site detection probability rising from 0.8381 to 0.9519. The reason is visible in the two contributions to the variance of psi. One is the binomial spread of a proportion, which shrinks as psi leaves 0.5, and the other is the cost of not knowing whether a site with no positives was empty or merely missed. For a rare species most sites are genuinely empty, the second term is small and most of the variance is the binomial term, which only more sites can shrink, so extra sites are the better buy; for a common species almost every all-zero site is a candidate false absence and resolving it is what the money is for. The same logic drives the two-level result in the visits tutorial, which is reassuring, because the two calculations share no code.
Designing at a guessed parameter
The optimum is optimal for the parameters you assumed, and those parameters are what the survey exists to estimate. So the honest question is not what the optimum is but how much a wrong guess costs. Fix the design chosen at psi = 0.4 and theta = 0.5, then evaluate it at other truths and compare with the design that would have been chosen there.
excess_curve <- function(vary, values) {
do.call(rbind, lapply(values, function(v) {
ps <- if (vary == "psi") v else psi_true
th <- if (vary == "theta") v else theta_true
opt <- design_grid(ps, th, p_true, budget)[1, ]
fixed <- se_psi(ps, th, p_true, best_row$sites, best_row$samples,
best_row$replicates)
data.frame(vary = vary, value = v, opt_samples = opt$samples,
opt_replicates = opt$replicates,
excess_pct = round(100 * (fixed / opt$se - 1), 2))
}))
}
sens <- rbind(excess_curve("psi", seq(0.05, 0.9, by = 0.05)),
excess_curve("theta", seq(0.15, 0.95, by = 0.05)))
is_psi <- sens$vary == "psi"
cat("worst extra SE from a wrong occupancy guess:",
max(sens$excess_pct[is_psi]), "percent; over psi up to 0.6 only",
max(sens$excess_pct[is_psi & sens$value < 0.61]), "percent\n")worst extra SE from a wrong occupancy guess: 17.23 percent; over psi up to 0.6 only 0.52 percentcat("worst extra SE from a wrong positivity guess:",
max(sens$excess_pct[!is_psi]), "percent\n")worst extra SE from a wrong positivity guess: 65.25 percentprint(sens[round(sens$value, 2) %in% c(0.05, 0.15, 0.6, 0.9), ],
row.names = FALSE) vary value opt_samples opt_replicates excess_pct
psi 0.05 3 2 0.52
psi 0.15 4 2 0.00
psi 0.60 5 2 0.44
psi 0.90 7 2 17.23
theta 0.15 10 2 65.25
theta 0.60 3 2 2.52
theta 0.90 2 2 19.90Guessing occupancy wrong is close to free over the range a survey is usually in. Up to psi = 0.6 the design chosen at psi = 0.4 loses at most 0.52 per cent of precision, and its worst case anywhere in 0.05 to 0.9 is 17.23 per cent, at the far end where the species occupies nine sites in ten. Guessing sample positivity wrong is expensive: at a true theta of 0.15 the design chosen at 0.5 gives a standard error 65.25 per cent larger than the best design on the grid, because it took four bottles where that design takes ten, the most the grid allows, so a grid that allowed more bottles per site could only widen the gap.
That asymmetry is the useful design advice in this post. A pilot study that only pins down rough occupancy has bought little; a pilot that measures how often a bottle from a known-occupied site tests positive has bought most of what the design needs; if the assay may be contaminated, the pilot should also run blanks, because whether the assay is clean moves the best replicate count (below).
sens$panel <- factor(ifelse(sens$vary == "psi",
"occupancy psi (design assumed 0.4)",
"sample positivity theta (design assumed 0.5)"),
levels = c("occupancy psi (design assumed 0.4)",
"sample positivity theta (design assumed 0.5)"))
ggplot(sens, aes(value, excess_pct)) +
geom_hline(yintercept = 0, colour = te_pal$line) +
geom_line(colour = te_pal$forest, linewidth = 0.9) +
geom_point(colour = te_pal$forest, size = 1.8) +
facet_wrap(~panel, scales = "free_x") +
labs(title = "A wrong positivity guess costs far more than a wrong psi",
subtitle = "extra SE of psi against the best design available at that truth",
x = "true value of the parameter",
y = "extra SE (percent)") +
theme_te()
Scoring the three-replicate rule
Field protocols routinely recommend three PCR replicates per sample, and the tile plot suggests the rule is not bad. Measure it: for each pair of assumed theta and p, compare the best design overall with the best design that is forced to three replicates.
folklore <- NULL
for (th in c(0.2, 0.35, 0.5, 0.8)) {
for (pp in c(0.4, 0.55, 0.7, 0.9)) {
gr <- design_grid(psi_true, th, pp, budget)
forced <- gr[gr$replicates == 3, ][1, ]
folklore <- rbind(folklore, data.frame(
theta = th, p = pp, best_replicates = gr$replicates[1],
best_samples = gr$samples[1], best_se = round(gr$se[1], 4),
rule_samples = forced$samples, rule_se = round(forced$se, 4),
loss_pct = forced$excess_pct))
}
}
print(folklore, row.names = FALSE) theta p best_replicates best_samples best_se rule_samples rule_se loss_pct
0.20 0.40 3 10 0.1251 10 0.1251 0.00
0.20 0.55 3 10 0.1157 10 0.1157 0.00
0.20 0.70 2 9 0.1091 10 0.1124 2.96
0.20 0.90 2 9 0.1045 10 0.1112 6.41
0.35 0.40 3 6 0.0962 6 0.0962 0.00
0.35 0.55 2 7 0.0896 6 0.0896 0.05
0.35 0.70 2 6 0.0842 5 0.0866 2.90
0.35 0.90 2 5 0.0805 5 0.0855 6.12
0.50 0.40 4 4 0.0812 5 0.0814 0.23
0.50 0.55 3 4 0.0756 4 0.0756 0.00
0.50 0.70 2 4 0.0716 4 0.0736 2.82
0.50 0.90 2 4 0.0693 4 0.0729 5.29
0.80 0.40 5 2 0.0640 3 0.0650 1.57
0.80 0.55 3 2 0.0598 2 0.0598 0.00
0.80 0.70 3 2 0.0573 2 0.0573 0.00
0.80 0.90 2 2 0.0545 2 0.0565 3.68cat("worst loss from the three-replicate rule:", max(folklore$loss_pct),
"percent, median loss:", median(folklore$loss_pct), "percent\n")worst loss from the three-replicate rule: 6.41 percent, median loss: 0.9 percentcat("cells where the rule is exactly optimal:",
sum(folklore$loss_pct < 0.005), "of", nrow(folklore), "\n")cells where the rule is exactly optimal: 6 of 16 The rule is close to optimal over most of this table. Its worst cell loses 6.41 per cent of precision, the median loss is 0.9 per cent, and it is exactly optimal in six of the sixteen cells. It mostly fails towards a good assay: at p = 0.9 the third replicate is nearly redundant and the loss runs from 3.68 to 6.41 per cent depending on theta. At the weak end (p = 0.4) it can miss the other way, three replicates being too few, but that costs at most 1.57 per cent. Against the flatness of the optimum and the uncertainty in theta, a fixed rule of three replicates costs less than a bad guess at sample positivity does, which is a fair defence of the folklore.
When the assay is not clean
Everything so far assumed that a positive replicate means target DNA. Put contamination back: a replicate from a DNA-free sample, and so every replicate from an unoccupied site, now amplifies with probability p10, while a replicate from a sample with DNA amplifies with probability p11, which keeps the value p = 0.7 used above. This is the four-parameter mixture of False positives in eDNA surveys. That post scores the replicate-count rule and measures what blanks buy, both at one fixed survey design; it does not ask where the budget optimum goes. For two-level occupancy surveys, Clement (2016) searched numerically for the designs that minimise the mean square error of occupancy when false positives occur, and in the model with two observation states found that false positives raise the recommended number of repeat surveys; what follows is the three-level, fixed-budget version of that question at this post’s costs, in a model where no positive is certain. Here p10 is estimated along with psi, theta and p11, at the same costs, budget and grid, for contamination rates of 0.01 and 0.05.
The shortcut that made the clean information matrix cheap is gone. An unoccupied site can now produce any pattern of counts, so every site outcome mixes the two branches, not only the all-zero one, and the expected information has to be summed over all possible site outcomes. The J samples of a site are exchangeable, so an outcome is just the number of samples at each count from zero to K, and there are choose(J + K, K) of them: 92378 at ten samples and nine replicates, few enough to enumerate. The score of each outcome is analytic. Two checks come first: the enumerated probabilities must sum to one, and with p10 fixed at zero and treated as known the new code must return the se_psi of the clean model.
site_outcomes <- function(n_samp, n_rep) {
bars <- combn(n_samp + n_rep, n_rep)
t(diff(rbind(0, bars, n_samp + n_rep + 1)) - 1)
}
info_mix <- function(psi, theta, p11, p10, n_samp, n_rep, p10_known = FALSE,
outcomes = site_outcomes(n_samp, n_rep)) {
counts <- 0:n_rep
bin11 <- dbinom(counts, n_rep, p11)
bin10 <- dbinom(counts, n_rep, p10)
qa <- theta * bin11 + (1 - theta) * bin10
score10 <- if (p10 > 0) {
counts / p10 - (n_rep - counts) / (1 - p10)
} else {
0 * counts
}
dlog_qa <- cbind((bin11 - bin10) / qa,
theta * bin11 * (counts / p11 - (n_rep - counts) / (1 - p11)) / qa,
(1 - theta) * bin10 * score10 / qa)
log_occ <- drop(outcomes %*% log(qa))
log_empty <- drop(outcomes %*% ifelse(bin10 > 0, log(bin10), -1e300))
top_log <- pmax(log_occ, log_empty)
e_occ <- exp(log_occ - top_log)
e_empty <- exp(log_empty - top_log)
site_lik <- psi * e_occ + (1 - psi) * e_empty
log_mult <- lgamma(n_samp + 1) - rowSums(lgamma(outcomes + 1))
prob <- exp(log_mult + top_log + log(site_lik))
w_occ <- psi * e_occ / site_lik
w_empty <- (1 - psi) * e_empty / site_lik
sum_occ <- outcomes %*% dlog_qa
scores <- cbind(psi = (e_occ - e_empty) / site_lik,
theta = w_occ * sum_occ[, 1],
p11 = w_occ * sum_occ[, 2],
p10 = w_occ * sum_occ[, 3] +
w_empty * drop(outcomes %*% score10))
if (p10_known) scores <- scores[, 1:3]
info <- crossprod(scores * sqrt(prob))
attr(info, "total_prob") <- sum(prob)
info
}
p10_values <- c(0.01, 0.05)
mix_grid <- NULL
enum_error <- 0
clean_gap <- 0
for (j in 2:10) {
for (k in 2:9) {
outc <- site_outcomes(j, k)
n_st <- floor(budget / design_cost(1, j, k))
info0 <- n_st * info_mix(psi_true, theta_true, p_true, 0, j, k, TRUE, outc)
enum_error <- max(enum_error, abs(attr(info0, "total_prob") - 1))
clean_gap <- max(clean_gap, abs(sqrt(solve(info0)[1, 1]) -
se_psi(psi_true, theta_true, p_true,
n_st, j, k)))
for (p10 in p10_values) {
info4 <- n_st * info_mix(psi_true, theta_true, p_true, p10, j, k,
FALSE, outc)
enum_error <- max(enum_error, abs(attr(info4, "total_prob") - 1))
vc4 <- solve(info4)
mix_grid <- rbind(mix_grid, data.frame(
p10 = p10, samples = j, replicates = k, sites = n_st,
se = sqrt(vc4[1, 1]),
se_p10_known = sqrt(solve(info4[1:3, 1:3])[1, 1]),
se_of_p10 = sqrt(vc4[4, 4]),
rcond = 1 / kappa(info4)))
}
}
}
cat("largest error in the enumerated probabilities:",
formatC(enum_error, format = "e", digits = 2), "\n")largest error in the enumerated probabilities: 2.44e-15 cat("largest gap to se_psi with p10 known to be zero:",
formatC(clean_gap, format = "e", digits = 2), "\n")largest gap to se_psi with p10 known to be zero: 9.71e-17 mix_best <- NULL
mix_top <- list()
for (p10 in p10_values) {
gr10 <- mix_grid[mix_grid$p10 == p10, ]
gr10 <- gr10[order(gr10$se), ]
gr10$rank <- seq_len(nrow(gr10))
gr10$excess_pct <- 100 * (gr10$se / gr10$se[1] - 1)
clean_cell <- gr10[gr10$samples == best_row$samples &
gr10$replicates == best_row$replicates, ]
mix_best <- rbind(mix_best, cbind(design = "contamination-aware", gr10[1, ]),
cbind(design = "clean optimum", clean_cell))
mix_top[[as.character(p10)]] <- gr10
cat("p10", p10, " smallest reciprocal condition number on the grid:",
format(min(gr10$rcond), digits = 3), "\n")
print(head(gr10[, c("samples", "replicates", "sites", "se",
"excess_pct")], 4), row.names = FALSE, digits = 4)
}p10 0.01 smallest reciprocal condition number on the grid: 0.000354
samples replicates sites se excess_pct
4 4 47 0.07944 0.0000
4 3 52 0.08013 0.8624
5 3 44 0.08070 1.5780
3 4 58 0.08101 1.9701
p10 0.05 smallest reciprocal condition number on the grid: 0.00148
samples replicates sites se excess_pct
4 4 47 0.08556 0.000000
4 5 43 0.08557 0.005267
5 4 40 0.08628 0.834055
4 6 40 0.08651 1.104948print(mix_best[, c("design", "p10", "sites", "samples", "replicates", "se",
"se_p10_known", "se_of_p10", "rank", "excess_pct")],
row.names = FALSE, digits = 4) design p10 sites samples replicates se se_p10_known
contamination-aware 0.01 47 4 4 0.07944 0.07923
clean optimum 0.01 58 4 2 0.09314 0.07910
contamination-aware 0.05 47 4 4 0.08556 0.08445
clean optimum 0.05 58 4 2 0.10991 0.09230
se_of_p10 rank excess_pct
0.00486 1 0.00
0.01120 31 17.24
0.01068 1 0.00
0.02084 55 28.46aware_row <- function(p10) {
mix_best[mix_best$p10 == p10 & mix_best$design != "clean optimum", ]
}
stopifnot(all(aware_row(0.01)[, c("sites", "samples", "replicates")] ==
aware_row(0.05)[, c("sites", "samples", "replicates")]))
clean_row <- function(p10) {
mix_best[mix_best$p10 == p10 & mix_best$design == "clean optimum", ]
}
insurance_pct <- base_grid$excess_pct[
base_grid$samples == aware_row(0.05)$samples &
base_grid$replicates == aware_row(0.05)$replicates]
cat("the contamination-aware design on a clean assay: extra SE",
insurance_pct, "percent\n")the contamination-aware design on a clean assay: extra SE 7.42 percenthedge <- merge(base_grid[, c("samples", "replicates", "sites", "excess_pct")],
mix_top[["0.01"]][, c("samples", "replicates", "excess_pct")],
by = c("samples", "replicates"), suffixes = c("_clean", "_01"))
hedge <- merge(hedge, mix_top[["0.05"]][, c("samples", "replicates",
"excess_pct")],
by = c("samples", "replicates"))
names(hedge)[names(hedge) == "excess_pct"] <- "excess_pct_05"
hedge$worst <- pmax(hedge$excess_pct_clean, hedge$excess_pct_01,
hedge$excess_pct_05)
hedge <- hedge[order(hedge$worst), ]
stopifnot(hedge$replicates[1] == 3)
cat("smallest worst-case extra SE over clean, 0.01 and 0.05:",
with(hedge[1, ], sprintf("%d x %d x %d, %.2f percent", sites, samples,
replicates, worst)), "\n")smallest worst-case extra SE over clean, 0.01 and 0.05: 52 x 4 x 3, 2.82 percent print(head(hedge, 3), row.names = FALSE, digits = 3) samples replicates sites excess_pct_clean excess_pct_01 excess_pct_05 worst
4 3 52 2.82 0.862 2.69 2.82
5 3 44 6.97 1.578 1.92 6.97
4 4 47 7.42 0.000 0.00 7.42Both checks pass to floating-point rounding, and the largest gaps are printed above: the enumerated probabilities sum to one, and with p10 known to be zero the new code agrees with se_psi on all 72 designs. Every design on the grid also identifies all four parameters from the survey alone: the smallest reciprocal condition number of the information matrix is 0.00035 at p10 = 0.01, nowhere near the \(10^{-10}\) at which se_psi gives up. So in this model, where contamination strikes each replicate independently, blanks are not needed to make it estimable, and none are added to the design here. That is a property of the model rather than of real assays. Contamination that arrives a whole sample at a time (a dirty bottle, one bad plate) can give counts that look like DNA, no replicate pattern separates them, and then only blanks or sites of known status identify the rate, which is the point False positives in eDNA surveys ends on; Guillera-Arroita et al (2017) show that false positives at two sampling stages need at least two sources of extra information. Everything below is the optimum for independent, replicate-level contamination.
Identified is not the same as well estimated, and that is where the optimum moves. At both contamination rates the best design is 47 sites, 4 samples and 4 replicates, against 58 x 4 x 2 for the clean assay: 2 more replicates per sample, paid for with 11 fewer sites. The clean optimum falls to rank 31 of 72 at p10 = 0.01 and rank 55 at p10 = 0.05, and its standard error of psi is 17.24 and 28.46 per cent larger than the contamination-aware optimum’s (0.0931 against 0.0794, and 0.1099 against 0.0856). These are asymptotic values, and the simulation below checks them. The new optimum is as flat as the old one: at p10 = 0.05 the runner-up, 43 x 4 x 5, is 0.005 per cent behind, and the four best designs use 4 to 6 replicates.
The se_of_p10 column says why. With two replicates, one positive out of two is what a sample with DNA gives 42 per cent of the time, and it is also what a contaminated sample that amplifies at all gives 99.5 per cent of the time at p10 = 0.01 and 97.4 per cent at p10 = 0.05 (the chance of one positive given at least one, 2(1 - p10) / (2 - p10)), so the counts carry little about p10: at p10 = 0.01 its standard error on the clean design is 0.0112, larger than the rate itself. With four replicates a sample with DNA shows at most one positive only 8.4 per cent of the time, the two kinds of sample separate (the separation that False positives in eDNA surveys describes, here priced in sites), and the standard error of p10 falls to 0.0049.
Fixing p10 at its true value, which is what an unlimited supply of blanks sharing the field samples’ contamination would do, separates the two sources of that loss. At p10 = 0.01 the clean design with p10 known has a standard error of 0.0791, against 0.0794 for the contamination-aware optimum with p10 estimated, so the whole of that 17.24 per cent is the price of learning p10 from two replicates. At p10 = 0.05 the clean design with p10 known is still 7.88 per cent worse, so part of the loss remains even with p10 known: two replicates separate a contaminated sample from one with DNA less well than four do, and no number of blanks changes that. On the four-replicate design the survey learns its own contamination rate nearly for free: estimating p10 instead of knowing it costs 0.27 per cent at p10 = 0.01 and 1.31 per cent at p10 = 0.05.
An exact formula for a mixture deserves a simulation check, at the four designs the comparison rests on: the clean and the contamination-aware optimum at each contamination rate. The fits keep the ordering constraint of the neighbour post, p11 = p10 + (1 - p10) * gap, which rules out label switching, and start from three points, keeping the best likelihood, because the likelihood surface can still stop an optimiser short of its best point.
sim_mix <- function(psi, theta, p11, p10, n_sites, n_samp, n_rep) {
occupied <- rbinom(n_sites, 1, psi)
has_dna <- rbinom(n_sites * n_samp, 1, theta * rep(occupied, n_samp))
matrix(rbinom(n_sites * n_samp, n_rep, ifelse(has_dna == 1, p11, p10)),
n_sites, n_samp)
}
nll_mix <- function(par, ymat, n_rep) {
psi <- plogis(par[1])
theta <- plogis(par[2])
p10 <- plogis(par[3])
p11 <- p10 + (1 - p10) * plogis(par[4])
l_occ <- rowSums(log(theta * dbinom(ymat, n_rep, p11) +
(1 - theta) * dbinom(ymat, n_rep, p10)))
l_empty <- rowSums(dbinom(ymat, n_rep, p10, log = TRUE))
top_log <- pmax(l_occ, l_empty)
ll <- sum(top_log + log(psi * exp(l_occ - top_log) +
(1 - psi) * exp(l_empty - top_log)))
if (!is.finite(ll)) return(1e10)
-ll
}
mix_starts <- list(c(0, 0, -3, 0.5), c(-1, 1, -4, 1), c(1, -1, -2, 0))
fit_mix <- function(ymat, n_rep) {
fits <- lapply(mix_starts, function(st) optim(st, nll_mix, ymat = ymat,
n_rep = n_rep, method = "BFGS", control = list(maxit = 500)))
best_fit <- fits[[which.min(sapply(fits, function(f) f$value))]]
c(plogis(best_fit$par[1]), plogis(best_fit$par[3]), best_fit$convergence)
}
sd_mcse <- function(x) {
s <- sd(x)
sqrt((mean((x - mean(x))^4) - s^4) / length(x)) / (2 * s)
}
set.seed(7331)
n_mix_fits <- 600
sim_check <- NULL
sim_est <- list()
for (i in seq_len(nrow(mix_best))) {
d <- mix_best[i, ]
est <- t(replicate(n_mix_fits, fit_mix(
sim_mix(psi_true, theta_true, p_true, d$p10, d$sites, d$samples,
d$replicates), d$replicates)))
sim_est[[paste(d$p10, d$design)]] <- est
sim_check <- rbind(sim_check, data.frame(
p10 = d$p10, design = d$design,
cell = paste(d$sites, d$samples, d$replicates, sep = " x "),
mean_psi_hat = mean(est[, 1]), simulated_sd = sd(est[, 1]),
sd_mcse = sd_mcse(est[, 1]), asymptotic_se = d$se,
p10_hat_near_zero = mean(est[, 2] < 0.001),
failures = sum(est[, 3] != 0)))
}
sim_check$ratio <- sim_check$simulated_sd / sim_check$asymptotic_se
sim_check$ratio_mcse <- sim_check$sd_mcse / sim_check$asymptotic_se
cat("maximum likelihood fits per design:", n_mix_fits, "\n")maximum likelihood fits per design: 600 print(sim_check, row.names = FALSE, digits = 4) p10 design cell mean_psi_hat simulated_sd sd_mcse
0.01 contamination-aware 47 x 4 x 4 0.4029 0.08082 0.002328
0.01 clean optimum 58 x 4 x 2 0.3994 0.08682 0.002347
0.05 contamination-aware 47 x 4 x 4 0.4074 0.08469 0.002731
0.05 clean optimum 58 x 4 x 2 0.4037 0.11163 0.003672
asymptotic_se p10_hat_near_zero failures ratio ratio_mcse
0.07944 0.0150 0 1.0173 0.02930
0.09314 0.2333 0 0.9321 0.02520
0.08556 0.0000 0 0.9898 0.03192
0.10991 0.0150 0 1.0156 0.03341sim_loss <- NULL
for (p10 in p10_values) {
at_rate <- sim_check[sim_check$p10 == p10, ]
aw <- at_rate[at_rate$design != "clean optimum", ]
cl <- at_rate[at_rate$design == "clean optimum", ]
ratio_sd <- cl$simulated_sd / aw$simulated_sd
sim_loss <- rbind(sim_loss, data.frame(
p10 = p10,
asymptotic_loss_pct = clean_row(p10)$excess_pct,
simulated_loss_pct = 100 * (ratio_sd - 1),
mcse_pct = 100 * ratio_sd * sqrt((cl$sd_mcse / cl$simulated_sd)^2 +
(aw$sd_mcse / aw$simulated_sd)^2)))
}
print(sim_loss, row.names = FALSE, digits = 3) p10 asymptotic_loss_pct simulated_loss_pct mcse_pct
0.01 17.2 7.42 4.24
0.05 28.5 31.81 6.07sim_row <- function(p10, dsg) {
sim_check[sim_check$p10 == p10 & sim_check$design == dsg, ]
}
ratio_z <- (sim_check$ratio - 1) / sim_check$ratio_mcse
names(ratio_z) <- paste(sim_check$p10, sim_check$design)
loss_gap_z <- (sim_loss$simulated_loss_pct - sim_loss$asymptotic_loss_pct) /
sim_loss$mcse_pct
print(round(ratio_z, 2))0.01 contamination-aware 0.01 clean optimum 0.05 contamination-aware
0.59 -2.69 -0.32
0.05 clean optimum
0.47 cat("simulated minus asymptotic loss, in Monte Carlo standard errors:",
round(loss_gap_z, 2), "\n")simulated minus asymptotic loss, in Monte Carlo standard errors: -2.31 0.55 est_edge <- sim_est[["0.01 clean optimum"]]
on_edge <- est_edge[, 2] < 0.001
cat("clean optimum at p10 = 0.01: correlation of psi and p10 estimates",
round(cor(est_edge[, 1], est_edge[, 2]), 3), "\n")clean optimum at p10 = 0.01: correlation of psi and p10 estimates -0.427 cat(" mean psi estimate with fitted p10 below 0.001:",
round(mean(est_edge[on_edge, 1]), 4), " otherwise:",
round(mean(est_edge[!on_edge, 1]), 4), "\n") mean psi estimate with fitted p10 below 0.001: 0.4579 otherwise: 0.3816 No fit failed to converge. The ratios of simulated to asymptotic standard deviation are 1.017 and 0.990 on the four-replicate optimum at p10 = 0.01 and 0.05, and 1.016 on the clean optimum at p10 = 0.05, each within 0.59 Monte Carlo standard errors (about 0.03 each) of one, so at those three designs the formula holds. The fourth is the exception: on the clean optimum at p10 = 0.01 the ratio is 0.932, 2.7 standard errors below one, so psi varies less across surveys than the formula says. That is the design on which the formula already predicted trouble with p10, and the boundary shows: its fitted contamination rate is below 0.001 in 23 per cent of surveys, which is what a standard error larger than the rate looks like in practice, against 1.5 per cent on the four-replicate design. A rate that cannot go below zero cannot swing as far as the normal approximation behind the formula allows, and psi is tied to it: the two estimates have a correlation of -0.43, and the surveys with p10 on zero have a mean psi estimate of 0.458 against 0.382 in the rest. The surveys that would have wanted a negative contamination rate are the high-psi ones, and the likely reason for the narrower spread is that the boundary holds them back.
The loss measured directly, as a ratio of two simulated standard deviations, is noisier than either ratio alone. At p10 = 0.05 the clean optimum’s spread of psi is 31.8 per cent wider than the contamination-aware optimum’s (Monte Carlo standard error 6.1), against 28.5 per cent from the formula, 0.6 standard errors apart: at that rate the formula’s loss stands. At p10 = 0.01 it is 7.4 per cent (Monte Carlo standard error 4.2) against 17.2, 2.3 standard errors short of the formula, which is the boundary effect above acting on the two-replicate design. So at p10 = 0.01 and 58 sites the asymptotic 17.24 per cent overstates the loss. The simulated loss is still positive, 1.75 standard errors above zero from 600 fits per design, so the direction probably holds, and the size is smaller than the formula says.
clean_curve <- do.call(rbind, lapply(2:9, function(k) {
sub_k <- base_grid[base_grid$replicates == k, ]
data.frame(replicates = k, se = min(sub_k$se))
}))
clean_curve$assay <- "clean assay (p10 = 0, known)"
dirty_curve <- do.call(rbind, lapply(
split(mix_grid, list(mix_grid$p10, mix_grid$replicates)),
function(z) data.frame(replicates = z$replicates[1], se = min(z$se),
p10 = z$p10[1])))
dirty_curve$assay <- sprintf("p10 = %.2f, estimated", dirty_curve$p10)
curves <- rbind(clean_curve, dirty_curve[, c("replicates", "se", "assay")])
curves$assay <- factor(curves$assay, levels = c("clean assay (p10 = 0, known)",
"p10 = 0.01, estimated",
"p10 = 0.05, estimated"))
curve_min <- do.call(rbind, lapply(split(curves, curves$assay),
function(z) z[which.min(z$se), ]))
ggplot(curves, aes(replicates, se, colour = assay)) +
geom_line(linewidth = 0.9) +
geom_point(size = 1.6) +
geom_point(data = curve_min, size = 4, shape = 21, fill = NA, stroke = 1.2) +
scale_colour_manual(values = c(te_pal$forest, te_pal$gold, te_pal$clay),
name = NULL) +
scale_x_continuous(breaks = 2:9) +
labs(title = "Contamination moves the best replicate count",
subtitle = "best SE of psi at each replicate count, samples and sites chosen freely",
x = "PCR replicates per sample", y = "SE of psi") +
theme_te() +
theme(legend.position = "top")
For a clean assay the best achievable standard error climbs from two replicates onwards, which is the saturation argument of the earlier section. Contamination lifts the left end of the curve, because two replicates separate a contaminated sample from one with DNA poorly, and the minimum moves to four. From there the three curves climb together and close up by nine replicates, where the components are well separated and every extra well costs the same sites whatever p10 is.
So the clean optimum does not survive contamination at this cost ratio. Insurance is cheap, and the cheapest policy is not the four-replicate optimum. Run on an assay that turns out to be clean and analysed with the clean model, the four-replicate design gives up 7.42 per cent against the two-replicate optimum. The three-replicate design 52 x 4 x 3 is never more than 2.82 per cent behind the best design, on the asymptotic standard errors, for any of the three assays tried here (clean, p10 = 0.01 and p10 = 0.05), which is the three-replicate rule of the previous section, defended on a new ground. If the contamination rate of the assay is unknown, it is a parameter to guess at before a design is tuned to it, like sample positivity; three replicates hedge across the rates tried here.
What the optimum does not know
Three limits, and none of them is small.
The design depends on parameters the survey is meant to estimate. The sensitivity measured above bounds the damage for psi and for theta separately at this cost ratio, but a design optimised at a guess is not a design with a guarantee, and the loss from a wrong theta was large enough to change what a pilot study should measure.
The cost model is a caricature. Travel is not linear in the number of sites: a route that visits 58 sites is cheaper per site than one that visits five, and a remote site can cost more than a whole plate. Replicates come in steps because a plate has a fixed number of wells, so the grid above has affordable designs on it that no lab would run. Substituting realistic costs is easy; the point is that the answer is a function of them, and reporting an optimum without the cost ratios is not reproducible.
Minimising the standard error of psi is a choice, not the only objective. A survey that only needs to know whether a species is present anywhere in a catchment should maximise the probability of at least one detection, which in this model pushes towards more sites than the standard-error optimum, each sampled more lightly, because every extra site is another chance of reaching an occupied one, and a survey comparing two areas needs the standard error of a difference instead. Each objective has its own optimum on the same grid, which is cheap to re-score once the information matrix is in hand.
Where to go next
The next tutorial in this cluster turns the machinery on a finished analysis rather than a planned one: what to check before believing a three-level eDNA occupancy fit, how to tell a boundary estimate from a real one, and which of the diagnostics from the single-season occupancy world still apply when the detection process has two nested layers.
References
Guillera-Arroita G, Lahoz-Monfort JJ 2012 Methods in Ecology and Evolution 3(5):860-869 (10.1111/j.2041-210X.2012.00225.x)
MacKenzie DI, Royle JA 2005 Journal of Applied Ecology 42(6):1105-1114 (10.1111/j.1365-2664.2005.01098.x)
Bailey LL, Hines JE, Nichols JD, MacKenzie DI 2007 Ecological Applications 17(1):281-290 (10.1890/1051-0761(2007)017[0281:SDTIOS]2.0.CO;2)
Dorazio RM, Erickson RA 2018 Molecular Ecology Resources 18(2):368-380 (10.1111/1755-0998.12735)
Clement MJ 2016 Methods in Ecology and Evolution 7(12):1538-1547 (10.1111/2041-210X.12617)
Guillera-Arroita G, Lahoz-Monfort JJ, van Rooyen AR, Weeks AR, Tingley R 2017 Methods in Ecology and Evolution 8(9):1081-1091 (10.1111/2041-210X.12743)