Clustered and Unclustered Group Testing for Biosecurity

Group-testing is an important element of biosecurity operations, designed to reduce the risk of introducing exotic pests and pathogens with imported agricultural products. Groups of units, such as seeds, are selected from a consignment, and tested for contamination, with a positive or negative test returned for each group. These schemes are usually designed such that the probability of detecting contamination is high assuming random mixing and a somewhat arbitrary design prevalence. We pro-pose supplementing this approach with an assessment of the distribution of the number of contaminated units conditional on testing results. We develop beta-binomial models allowing for between-consignment variability in contamination levels, with a further layer of nesting to allow for possible clustering within the groups for testing. The latent beta distributions can be considered as priors and chosen based on expert judgement, or estimated from historical test results. We show that the parameter controlling within-group clustering is, unsurprisingly, eﬀectively non-identiﬁable. It can be handled by sensitivity analysis, but we demonstrate theoretically and empirically that the probability of a consignment with contamination evading detection is almost perfectly robust to mis-speciﬁcation of this clustering. We apply the new models to large cucur-bit seed lots imported into Australia where they provide important new insights for biosecurity regulation.


INTRODUCTION
A fundamental challenge for the protection of agricultural productivity and trade in many countries is to confirm the absence of exotic pests and pathogens in arriving or departing consignments. The volume and movement of agricultural goods is increasing rapidly (Chapman et al., 2017) but, although biosecurity regulations require inspections, costs and practicalities limit the extent of sampling and testing.
Group testing is a technique commonly applied to the inspection of consignments (e.g., shipments) of units (or items), where units may be prawns, grain, cut flowers (Hepworth, 1996), perishable goods or seeds (Constable et al., 2018). Units are selected and tested collectively in groups, where group size may be determined by the efficacy of the testing process (e.g., prawns), or by the cost of testing and other practical considerations (e.g., seeds).
The test may be non-destructive (e.g., cut flowers) or destructive (e.g., seeds), and the test result is either positive or negative for each group. Typically, the number of contaminated items in a group that tests positive remains unknown, with the detected contamination either actioned or the consignment rejected for importation -both at considerable cost. We focus on biosecurity, but forms of group testing are also used in other contexts including Covid-19 management (Lokuge et al., 2021) and environmental DNA surveys (Furlan et al., 2016).
Effective group-sampling schemes must take account of the potential heterogeneity of contamination within consignments, whether this be due to the clustering of contaminants within consignments of goods, or clustering related to country or region of origin or other spatial variation. However, because the extent of clustering within groups may be unidentifiable, it is common to assume that groups are random samples of units. The number of groups tested is often small and zero contamination is detected in the great majority of cases, adding to the challenge of analysing test results and assessing risk.
A key concept is the leakage (or slippage), which is the number of contaminated units that are imported. If there are any contaminated units in the sampled groups, the leakage is zero, because importation is then blocked. If there is contamination in the consignment but none in the inspected sample, then leakage is positive.
Group testing regimes are often designed by specifying a maximum tolerable or plausible prevalence of contaminated seeds (p max ). The number of groups to be sampled is calculated such that the probability that contamination occurs in the sample is at some high level if the prevalence equals p max . Leakage is therefore unlikely if p max obtains. If groups are random samples, and the number of contaminated seeds in the consignment is treated as fixed, then the number of seeds in the sample is hypergeometrically distributed (e.g., IPPC 2008).
The binomial distribution is sometimes used as a simpler approximation, where algebraic results are possible, provided the sampling fraction is small. Overdispersed extensions of binomial distributions have also been used. Trouvé and Robinson (2020) consider a twolevel model with normal random effects to allow for between-consignment heterogeneity.
They conduct sensitivity analysis on the consignment-level variance component, rather than estimating this parameter, and analytical formulas are not available as they involve integration over the latent variable. Beta-binomial distributions have also been suggested in non-group-testing contexts (Venette et al., 2002;IPPC, 2008;Trouvé et al., 2022). This paper was motivated by a specific biosecurity problem faced by an operational division of the Australian government. Australia imports numerous consignments of seeds for a variety of agricultural commodities (e.g., Constable et al. 2018;Dall et al. 2019). Biosecurity regulations require that consignments are tested, and are rejected for importation if contamination is detected. Because import volumes are increasing and testing is expensive and destructive, there is an imperative to reduce costs for importers, but also to understand the likely consequences for risk with any changes to the regulations.
Motivated by these challenges, we answer the following questions about group sampling and testing: • Can we obtain distributions for the leakage conditional on group test results, both when groups are clustered and when they are randomly formed? Existing literature has focussed on the probability of detecting contamination assuming p max rather than inferring leakage.
• Can we estimate model parameters reflecting between-group and between-consignment heterogeneity using historical group testing results?
• How sensitive are estimates of model parameters and leakage rates to assumptions regarding clustering within groups? Section 2 defines our notation and summarises selected models in the literature. Section 3 describes our proposed model when groups are assumed to be randomly formed. We use Beta distributions to reflect heterogeneity or uncertainty about the prevalence in a particular consignment. Our model is similar to that of Trouvé et al. (2022) but allows for group sampling rather than full observation of the number of contaminated units. New approximations are derived for the distribution of contaminated units in a consignment given group testing results. The distribution of the leakage and the worst case for expected leakage are also derived.
Section 4 generalizes this model to allow for clustering of contamination within groups. A heuristic approximation reveals that the parameter controlling clustering is not identifiable, and so must be assumed a-priori. We show that if the model is estimated from group testing outcome data, estimates of the probability of leakage from a given shipment are approximately correct even if the clustering parameter is mis-specified. This is a novel and important finding, because it means that the possibility of clustering can be ignored when the primary concern is with the number of consignments with undetected contamination.
The number of contaminated units entering the country, however, is sensitive to clustering.
Section 5 is a case study, showing how our new framework would lead to more informed biosecurity decisions in a particular example. A simulation study motivated by the case study is used to validate estimation and inference methods. Section 6 is a discussion.

Problem Definition
We consider consignments (e.g., shipments, consignments, or subdivisions of consignments), containing units (e.g., seeds) which are grouped into equally sized groups. Typically, there are many consignments of a particular kind imported into a country annually, and we refer to ongoing importation of a particular type of product in a particular form as a pathway.
Each consignment may or may not contain some contaminated units, and the focus is on inferring the risk of importation of contaminated units, using models which may be assumed a priori, or estimated from historical data on multiple consignments' inspection results.
Each consignment is assumed to contain N units, which are grouped into B equally sized groups, each containingN = N/B units. (We assume that N is divisible byN , which would be a negligible approximation in practice.) A simple random sample without replacement of b groups is selected, and tested. We let X ij be the number of contaminated units in group j of consignment i, j = 1, . . . , B, so that X ij may take on values 0, 1, . . . ,N . We denote the indexes of the sampled groups as 1, . . . , b, and the non-sampled groups are b + 1, . . . , B. The total number of units in selected groups is n = bN .
The testing procedure is such that we do not directly observe X ij even for sampled groups. Instead, all that is available for sampled groups j is Y ij which equals 1 if there is any contamination in the group and 0 otherwise, i.e., Y ij = I (X ij > 0) where I(.) is the indicator function. We write T Xi = B j=1 X ij and t xi = b j=1 X ij for the population and sample number of contaminated units, respectively, in consignment i. We let T Y i = B j=1 Y ij and t yi = b j=1 Y ij be the population and sample number of contaminated groups respectively. Error in measurement of X ij or Y ij is sometimes allowed for in the literature (Cowling et al., 1999;Liu et al., 2012), although only with randomly formed groups. We assume throughout that testing is perfect, to simplify the exposition given the other innovations made, and because this is realistic in the case study in Section 5.
We define the leakage (sometimes referred to as slippage) to equal L i = T Xi I (t yi = 0), which is the number of contaminated units imported in consignment i, reflecting the fact that consignments are blocked if any contamination is detected.
An important aim in practice, which is not usually considered explicitly, is to infer T Xi for a particular consignment, using the observed data t yi , particularly the case when t yi = 0.
The leakage L i is also a key quantity. We consider the expected leakage E(L i ) and the probability of leakage P(L i > 0) under plausible statistical models. These quantities provide important new measures of the ongoing risk from an importation pathway.

Review of Models assuming Random Group Formation
Suppose that groups are simple random samples without replacement of units, constrained to be non-overlapping. The sampling distribution of t yi , treating T Xi as a fixed parameter, can then be obtained without any further modelling assumptions. It is clear in this case that t xi follows a hypergeometric distribution with parameters N , T Xi and n (e.g., McSorley and Littel 1993). The distribution of t yi is less obvious but has been derived by Theobald and Davie (2014) who considered group testing with random groups. Interest has often focussed on the probability that t yi = 0, which is: This equation can be inverted numerically to find the sample size n such that the right hand side of (1) is sufficiently low given a postulated value T Xi = p max N (e.g., Lane et al. 2018;Constable et al. 2018). When n ≪ N , the distribution of t xi given T xi is approximately bin (n, T Xi /N ) (Brunk et al., 1968). The required sample size is then n ≈ log P (t yi = 0) / log (1 − T Xi /N ) (e.g., McSorley and Littel 1993;Lane et al. 2018). More generally, the symmetries of the hypergeometric distribution allow for four distinct binomial approximations, with performance dependent on conditions (Lieberman and Owen, 1961). Hepworth (1996Hepworth ( , 2019 also considers binomial models for group sampling, assuming random group formation, and focusing on exact inference and bias, respectively. Cowling et al. (1999) considers inference about the prevalence from unclustered group sampling, but does not consider leakage.

Review of Models assuming Non-Random Group Formation
It is possible that contamination tends to co-occur in the same group or groups. This will result in X ij and t xi being overdispersed relative to the hypergeometric or binomial distributions discussed under the previous heading. McArdle (1990) discusses the possibility of clustering in sampling of plots for the presence of rare species. IPPC (2008) comment that clustering of contamination within consignments is possible, and suggest the beta-binomial distribution to model it. Venette et al. (2002) also suggest the beta-binomial distribution to model clustering of t xi , noting that the beta-binomial distribution is approximated by the negative binomial when the contamination rate is low (citing Madden et al. 1996). None of these papers discussing beta-binomial models consider group sampling.

Review of Models allowing for Consignment-Level Heterogeneity
Trouvé and Robinson (2020) consider the situation where the only data available is whether any detections are made from a whole consignment. They suggest an overdispersed model for t xi relative to the binomial distribution, with overdispersion controlled by a parameter α. This results in higher probability that t yi = 0 when there is overdispersion. They use data from multiple consignments with different sample sizes per consignment, assuming a common value of α and T Xi /N , enabling estimation of both T Xi /N and α. Simulation of different values of α shows that the bias from ignoring overdispersion can be substantial. Trouvé et al. (2022) model consignment-level heterogeneity by assuming that a latent propensity for each consignment follows a Beta distribution, equivalent to the model we will use in the next section. However, they assume that the counts t xi of contaminated units are observed, and do not consider group testing. They build on this basic model to infer the consequences if observed consignments are filtered based on testing in the source country.

Model Definition
We propose the following model: for consignments i and groups j = 1, . . . , B, where α and β are shape parameters. When we are considering a single consignment of interest, as is often the focus in biosecurity processes, and there are no suitable data to estimate α and β, then (3) is a prior distribution representing the state of knowledge or uncertainty about the prevalence of contamination.
If historical data from a set of consignments are available, then α and β can be estimated from these data by maximum likelihood.
Note that (2) and (3) imply that T Xi ∼ BB (N, α, β), where BB denotes the beta-binomial Under group testing, we do not know t xi , but instead only observe t yi for one or more consignments i. Under the above model, Y ij are independent and identically distributed Bernoulli random variables conditional on p i , with probability given by As φ i is a function of p i , it is also a random variable. To emphasise this functional relationship, we will write φ i = φ (p i ) when convenient.

Likelihood for Data on t yi
It is clear that t yi |p i ∼ bin (b, φ i ) and therefore If observations of t yi for a sample of m consignments are available, then the likelihood is given by the product of (5) over i = 1, . . . , m. Numerical integration is required.

Distribution of Contamination Levels Conditional on Group Test Results
Given t yi for a particular consignment i, the distribution of p i conditional on a group test result t yi follows directly: For the special case t yi = 0, the right hand side of (7) is proportional to a Beta density, with This is a standard result in binomial models with beta priors.
The distribution of T Xi conditional on an observed value of t yi would often be of interest.
Section S1 of the Supplementary Materials describes an algorithm for sampling from this conditional distribution. An analytical formula is also available (see Supplementary Section S2) but it is quite complex, so we now derive a new simple approximation which should work well when p i is close to zero with high probability, as would often be the case in biosecurity where breaches are rare. Applying Taylor's Theorem, (Note that (9) can also be obtained by assuming that there is at most one contaminated unit within each group.) We directly obtain the following approximation to the posterior of T Xi (which is exact when t yi = 0):

Properties of the Leakage
We now consider the distribution of the leakage. We write f BB (k; n, α, β) to denote the probability function of a beta-binomial random variable with parameters n, α and β, evaluated at k. Assuming (2) and (3), the probability function of L i is defined by: where B(., .) is the Beta function. So L i is a zero-inflated beta-binomial (N − n, α, β) distribution. (See Section S3 of the supporting information for the derivation of (11).) We also consider E (L i ) and derive an upper bound for it. Note that Using the law of iterated expectations and the properties of the Beta distribution, we immediately obtain To maximise E (L i ) with respect to α and β, note first that E (L i |p i ) in (12) is maximised by p i = p max where p max = (n + 1) −1 by elementary calculus. It follows that with equality obtaining if and only if p i is equal to p max with probability of 1. This is the limiting case of the beta prior in (3) when both α and β tend to infinity with α/ (α + β) = 1/(n + 1) (or equivalently, with α = nβ). Therefore, the worst case prior for p i is the distribution where p i = (n + 1) −1 with probability 1. Lane et al. (2018, p41) also found that p max = 1/(n + 1) in a similar setup, but using a hypergeometric model. However, their leakage rate stated in (E.1) appears to be higher than the one derived here by a factor of ((N − n)/N ) 2 , possibly because an upper bound was derived rather than E(L i ) itself, with a further approximation to a combinatorial term between their equations (D.5) and (E.1).
It is notable the worst case is not simply when p i is high (or has high prior expected value). This is because when p i is high, it is likely that at least one sampled unit will be contaminated, in which case the consignment would be rejected and no leakage would occur.
The worst case is when p i is small enough that there may well be no contamination in the sample, but large enough that there are substantial contaminated units in the consignment.
The worst case of p i = p max could be used for planning purposes, although it may be too pessimistic. However, a prior distribution which has mean much lower than 1/(n + 1) should be used with caution, as results will sensitively depend on this prior if no contamination is detected in the sample.

Model Definition
If contamination tends to strike in the same group, then the binomial model (2) is not appropriate, and we would expect to see over-dispersion relative to this model. This section develops nested Beta models to model variation within and between groups. We assume the following prior and model: The hyper-parameter θ > 0 controls the level of clustering; when θ = ∞, p ij = p i and (16) is the same as (2). In the other extreme, when θ ↓ 0, the Beta distribution in (15) becomes Bernoulli with parameter p i (see Supplementary Materials Section S4).
We treat θ as a known hyper-parameter and we propose conducting sensitivity analysis over a range of plausible values of θ. We take this approach because, even if observations of t yi are available for multiple consignments, the model appears to be practically unidentifiable from observations of t yi if θ is unknown. This is not surprising given that all we observe is whether group totals X ij are zero or not, so identifying clustering within groups is unlikely to be feasible. The heuristic in Section 4.4 confirms this, but also shows that the estimated probability of leakage occurring in a consignment is insensitive to the assumed value of θ.
Conditional on p i , y ij are independent with where By elementary manipulations of beta functions, (18) can be re-expressed as: which is more computationally convenient. Note that φ i is a function of p i , which we will express as φ i = φ (p i ) when convenient. An equivalent form is Since θ (θ + k) −1 ≤ 1, it is clear that the right hand side of (20) has maximum value of φ i = 1 − (1 − p i )N , and the maximum is achieved as θ → ∞. This corresponds to the model in Section 3 and is the best case for φ i (i.e., greatest chance of detecting contamination). It is also clear that the worst possibility for φ i is θ ↓ 0, in which case φ i = p i .

Likelihood and Posterior
Note that t yi |p i ∼ bin (b, φ i ) and therefore where φ(p i ) is defined by (18). If observations of t yi for a sample of consignments are available, then the likelihood is given by the product of (21) over consignments i, with numerical integration required in practice.
The posterior of p i for a given consignment conditional on its value of t yi is The integral of the right hand side of (22) with respect to p i over [0, 1] does not have a simple form, but can be calculated numerically. The posterior of T Xi could also be obtained, using the fact that X ij |p i ∼ BB N , θp i , θ (1 − p i ) are independent.

Leakage for the Nested Beta Model
The distribution of the leakage is much more complex than in Section 3.4 when groups were unclustered. We focus on the probability that the leakage is zero and the expected leakage: These quantities can be obtained by numerical integration using the distribution of p i in (14) and the definition of φ as a function of p in (18).
Given that θ is unlikely to be identifiable in practice, we consider the worst case value of θ for the expected leakage. Combining (20) with (24), we obtain By inspection, it is obvious that the right hand side of (26) is maximised by θ = 0, in which case it becomes Because θ = 0 maximises E (L i |p i ) for every p i , it must also maximise E (L i ). By elementary maximisation, the right hand side of (27) is maximised with respect to p i by 1/(b + 1).
Using arguments similar to those made in Section 3.4, it follows that the worst case for E (L i ) with respect to θ, α and β is θ = 0 and α, β → ∞ with β = bα.

Heuristic Approximation of the Model
We now consider a heuristic approximations for the distribution of t yi , to show that θ is indeed not identifiable in practice, and to determine the effect of assuming an incorrect value of θ.
We consider the case where p i is small with high probability, which is often the case in practice (including in the case study in Section 5). Dropping all terms involving powers of p i greater than one, equation (20) leads to: where g (θ) = N −1 k=0 θ (θ + k) −1 . Note that g(0) = 1 and g(∞) =N . The quantity g (θ) has a real world interpretation: the approximate expected number of contaminated units per contaminated group is easily found to beN /g (θ), which is an intuitive measure of clustering within groups.
We further note that, for α ≪ β, the Beta (α, β) distribution is approximated by p i a ∼ Gamma (α, β −1 ) (shape-scale parameterisation) (Johnson et al., 1995), and therefore We have that t yi |p i ∼ bin (b, φ i ). If φ i ≪ 1 (again, a reasonable assumption in most biosecurity applications), then Putting together (29) with (30), and making use of the fact that Poisson-gamma mixtures are negative binomially distributed, we have that: using the mean-shape parameterisation of the negative-binomial distribution.
We conclude from (31) that it is possible to identify α from a dataset of observations of t yi , but it is not possible to identify both θ and β. The best that can be done is to estimate the product g (θ) β −1 . Therefore it is not possible to identify the marginal prevalence If we assume θ = ∞ but in reality θ = 0, then we will under-estimate µ by a factor of g (∞) /g (0) =N .
This sounds like a hopeless situation if we believe that there may be clustering of contamination within groups, but things are not quite as bad as they seem. Suppose that we incorrectly assume θ = θ 1 , when the true value is θ = θ 0 . We fit the multi-level beta model to data on t yi . Equation (31) means that we will then obtain a reasonable estimate of α, but we will mis-estimate β withβ ≈ βg (θ 1 ) /g (θ 0 ). If we use the fitted model to infer the distribution of φ i , our estimated distribution will then be Gamma α,β −1 g (θ 1 ) which is approximately the same as Gamma (α, β −1 g (θ 0 )), which is the correct distribution. So, even though ourβ may be well out, our estimated distribution of φ i will be about right.
If we use the fitted model to estimate the probability of positive leakage using (23), we see that this calculation only involves the distribution of φ i . We should therefore obtain reasonable estimates of the leakage probability even if we badly mis-state the value of θ.
However, if we are more interested in the per-unit expected prevalance (α/ (α + β)) or the expected leakage in (25) (which depends on the distribution of p i as well as on φ i ), then we would expect these results to be sensitive to mis-specification of θ.
The approximation in (31) also suggests a simple method for fitting model (14) and (16), rather than full maximum likelihood which requires numerical integration. Instead, a negative binomial model can be fitted to data on t yi , and (31) can be used to obtain estimates of α and β, for any assumed value of θ.

Background and Data
Cucurbit fruits such as melons, cucumbers and zucchinis are grown and traded in large quantities. Agricultural production of cucurbit fruits derives from cucurbit seeds, that are themselves commonly produced, processed and transported in bulk across multiple national boundaries, through highly structured supply chains (Bonny, 2017).
Cucurbit plants are susceptible to infection by many viruses, some of which can be transmitted in seeds from infected plants to emerging seedlings (DAWR, 2017). In order to manage biosecurity risks associated with imported seeds, Australia sets broad phytosanitary entry conditions, and in some cases, also requires testing for one or more seed-borne pests, such as the tobamovirus cucumber green mottle mosaic virus (CGMMV) (DAWR, 2017). In such circumstances, tested seed lots must be certified to be free from the detectable presence of the specified pest(s) before they are permitted entry. CGMMV has been repeatedly detected in seed lots intended for entry at the Australian border (Constable et al., 2018).
We refer to "infection" of cucurbit seeds by one or more viruses, which may be regarded as synonymous with "contamination" defined in Sections 1 -4 of the paper.
We focus on large seed lots (defined as more than 2kg for cucumber and melon seeds, which is approximately 80,000 seeds; DAWR 2017). A sample of b = 94 groups is taken from each large lot, with each group containingN = 100 seeds. We assume a total of B = 800 groups in each lot. Each sampled group is tested using the highly specific serology-based enzyme-linked immunosorbent assay (ELISA) protocol, resulting in a presence or absence finding. It is typically assumed that groups are random samples of seeds (or equivalently, that infections are spread at random across groups), although this may not always be perfectly satisfied, so that it is important to understand the sensitivity to this assumption.

Estimating Parameters from Historical Data
Test results from 2016 seed imports are used for this case study, including cucumber, melon, pumpkin, squash and watermelon seeds. No infection is recorded in 68 of the 72 large seed lots imported. From the other 4 large seed lots, the numbers of groups testing positive (i.e., t yi ) are 4, 4, 7 and 21. These results are a subset of the data discussed by Constable et al. (2018). In the terminology of Sections 1 -4, seed lots are consignments.
We fitted the nested beta model defined by (14) where f β is the beta density of p, and q k are the k/(K + 1) quantiles of this beta distribution where we set K = 9, 999. The R function qbeta used to calculate q k also gives warnings of numerical imprecision in some cases, but the approximation is still very precise in this example. The likelihood is maximised with respect to log(α) and log(β) using the optim function in R with default settings. Confidence intervals are constructed by profile likelihood as the simulation study summarised in Section 5.3 finds that Wald confidence intervals perform poorly here. The likelihood is calculated over a grid of approximately 2.1 × 10 5 values of α and β, and the confidence intervals are the range of values of each parameter of interest where the deviance is less than the relevant quantile of the χ 2 1 distribution.
Maximum likelihood estimators and profile confidence intervals are calculated for α, β, µ = E(p i ) = α/(α + β), the expected leakage E(L i ) from (13) and the probability of non-zero leakage P(L i > 0) from (11). The latter two quantities are obtained from (25) and (23), respectively, with numerical integration required in general. Table 1 shows the fitted models for the two assumed values of θ. We firstly discuss the results when it is assumed that there is no clustering of infections within groups (θ = ∞).
The confidence intervals are wide, unsurprisingly given that infection was only detected in four out of 72 lots. Nevertheless, the modelling provides some important insights. We see thatP(L i > 0) = 0.035, implying (on average) 2.5 seed lots with undetected infection in a year's worth of importing (about 70 lots). Even the lower end of the confidence interval (0.0100) for P(L i > 0) suggests leakage from 0.7 seed lots annually. This gives a useful new perspective on the risk from this pathway. The expected number of infected seeds admitted into the country is estimated atÊ(L i ) = 0.14 per lot or about 10 per year. The estimated first shape parameter is very small atα = 0.016. This parameter must be positive and the beta distribution has a mode at zero when α < 1. A value as low as 0.016 means a very sharp peak at zero, with high probability of a value close to zero. This reflects the data where the great majority of lots have zero detections out of 9400 seeds. Table 1 also shows results for the most extreme possible level of clustering of infection within groups (θ = 0). As predicted in the discussion in Section 4.4, the estimates of α and P(L i > 0) are virtually unchanged (compared to the unclustered case where θ = ∞), while the estimate of E(L i ) has increased by a factor ofN = 100. The estimate of β has also increased by a factor in the ballpark of 100. This shows that, provided data is available to estimate the model parameters, estimates of α and P(L i > 0) are insensitive to the assumption of random mixing, consistent with the theoretical approximation discussed in Section 4.4. If the proportion of lots admitted with infection (P(L i > 0)) is of primary concern, the random mixing model can be used with confidence that violations of this assumption will have little effect. But if the number of infected seeds in those lots is more important, then wrongly assuming random mixing will mislead.  Seed lots are blocked from importation if any infections are detected, so one question of interest is the distribution of T Xi conditional on t yi = 0. This calculation could be made by assuming any particular value of the clustering parameter θ, but we consider the unclustered case θ = ∞, which corresponds to the model in Section 3.1. Figure 2 shows the estimated distribution of T Xi unconditionally and conditional on t yi = 0, assuming model (2) -(3).
The figure shows that the conditional probability that T Xi = 0 is much higher at 0.97 than the unconditional probability of 0.91. Moreover the conditional probability of positive values of T Xi decreases rapidly with T Xi , particularly for T Xi ≥ 5), whereas the corresponding decrease in the unconditional probabilities is much more gentle. This demonstrates that the testing process greatly reduces the risk of any leakage at all, and virtually rules out the risk of more than about 5 infected seeds being imported from any given consignment.

Simulation Study for the Estimation of the Beta Distribution Parameters
The dataset in the case study is sparse, with positive detections in only 4 of 72 consignments.
A simulation study was conducted to evaluate estimation and inference of model parameters in this situation. The parameters α and β are estimated by maximum likelihood assuming model (2) and (3); this model assumes no clustering of contamination within groups. Nominal 95% Wald confidence intervals are calculated using the inverse of the observed information.
We simulated from a range of values of α spread around the estimate of α from the case study. Clustering is either non-existent (θ = ∞) (so that the generating model matched the fitted model), moderate (θ = 0.37) or perfect (θ = 0). The "moderate" value of θ was based on E(L) being roughly halfway between the two asymptotes in Figure 1(e). The values of µ (and hence β) were chosen such that the proportion of groups testing positive roughly matched the proportion in the case study. 1000 datasets were simulated for each scenario.
Results are tabulated in Section S5 of the Supplementary Information. The simulation found that estimates were highly variable reflecting the sparsity of detections. Nevertheless, it was confirmed thatα did estimated reasonably well for all values of α and all three values of θ, but as expected,μ did very poorly if θ was 0 or 0.37. Profile confidence intervals performed well whereas Wald intervals were sub-standard in some cases, in the sense that their non-coverage rates were far from the nominal value of 0.05.

DISCUSSION
We propose two and three-level models utilising latent Beta random variables to model heterogeneity between consignments and between groups. Using these models, we are able to obtain distributions for the leakage and for the number of infected units given testing results, via approximate or exact analytic results for random groups, and numerical calculation for non-random groups. The three-level model includes a parameter controlling the clustering within groups that cannot be identified from typically available data; however sensitivity analysis can be used. We show heuristically that the probability that there is non-zero leakage from a consignment is robust to mis-specification of this clustering. Empirical and simulation results back this up.
The novelty of the framework developed here is that it allows inference on ongoing leakage risk from a pathway. This adds significant value to the more standard approach where sampling is designed to achieve a high probability of detection assuming a somewhat arbitrary design prevalence with no allowance for heterogeneity between groups or consignments.
Group-sampling strategies are commonly applied in current biosecurity processes for many classes of imports in many countries, as well as in environmental sampling and public health.
This work provides new perspectives on the design criteria and associated data analysis required to manage risk.

CONFLICT OF INTEREST STATEMENT
The authors have no conflicts of interest to declare.  (14) -(16) based on cucurbit surveillance results. The value of θ is held fixed at either ∞ (infection unclustered within groups, which corresponds to model (2) -(3)) or 0 (infection perfectly clustered within groups). Confidence intervals are based on profile likelihoods. The expected leakage E(L i ) and the probability of leakage P(L i > 0) are the functions of model parameters α and β given by (25) and (23) (14) - (16) based on cucurbit surveillance results, using the approximate model (31). The value of θ is held fixed at either ∞ (infection unclustered within groups, which corresponds to model (2) -(3)) or 0 (infection perfectly clustered within groups). Confidence intervals are based on profile likelihoods. The expected leakage E(L i ) and the probability of leakage P(L i > 0) are the functions of model parameters α and β given by (25) and (23), respectively. Figure 1 Estimates of various parameters, and deviance, vs assumed values of θ Estimated distribution of the number of infected seeds (T Xi ) in a consignment, based on maximum likelihood estimates of α and β where θ = ∞ is assumed (i.e., no clustering within groups). Both the unconditional distribution and the conditional distribution given no positive test results (t yi = 0) are shown.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. supplementarymaterial13092022.pdf