An analysis of the case-control study without group testing would require a direct comparison of the frequency of antibody-specific reactions between cases and controls across the large number of antibodies, as depicted in Fig. 1. These frequencies are usually based on thresholding a quantitative serological assay or can be directly assessed with a qualitative assay that is inherently dichotomous. The case-control analysis with 15,000 antibodies would require researchers to analyze 15,000 multiplied by the total study sample size in number of assays. The large number of assays required for a sufficiently powered study would make this approach infeasible. We propose a group testing strategy to substantially reduce the number of required assays (tests) without sacrificing much power.
Our inferential goal is to test for case-control differences for each antibody where we control the point-wise error rate (e.g., each antibody-specific comparison between cases and controls has a type I error rate of \(\alpha\)). We recognize the number of antibodies is large and that we would expect an average number of false discoveries of \(\alpha\) multiplied by the number of antibodies.
We propose a two-phase design where in the first phase we screen antibodies using group testing and only proceed to a second stage when there is a good indication of an effect. We describe the procedure as follows.
Phase 1
In phase 1, we use group testing to estimate the prevalence of antibodies, compare the prevalence estimates between cases and controls, and use this comparison to select antibodies. We split the observations by case-control status into groups of equal size. We then test each group for each antibody. We estimate the case and control prevalences for each antibody using the Burrows estimator [7, 8] and references within. This estimator is given by
$$\widehat p=1-\left(1-\frac x{n+v}\right)^\frac1k\;with\;v=\frac{k-1}{2k}$$
where x is the number of positive groups, k is the group size, and n is the number of groups.
Although prevalence can be estimated using maximum-likelihood (MLE), this estimator will be biased. An alternative estimator was proposed by Burrows that eliminates most of the bias. In addition, Burrows showed empirically that his estimator not only improves on the bias but yields a smaller mean-square error (MSE) than the MLE for all values of p considered (p ≤ 0.5) [7].
We compute a two-sided z-test for each antibody to evaluate evidence for a case-control difference,
$${Z}_{antibody}=\frac{{\widehat{p}}_{case}-{\widehat{p}}_{control}}{\sqrt{\widehat{var}\left({\widehat{p}}_{case}\right)+\widehat{var}\left({\widehat{p}}_{control}\right)}}$$
where \({\widehat{p}}_{case}\) and \({\widehat{p}}_{control}\) are Burrow’s estimators for the cases and controls, respectively. The variances of these estimators are computed as
$$\widehat{var}\left(\widehat p\right)=\frac{\left(1-\theta\right)\left(1-\widehat p\right)^2}{k^2}\ast\left(\frac1{n\theta}+\frac{2\left(1-\theta\right)v^2}{(n{\theta)}^2}\right)-\left(\frac{v\left(1-2v\right)\left(1-\theta\right)\left(1-\widehat p\right)\left(1+\theta\right)\left(1-v\right)}{6n^2\theta^2}\right)^2with\;\theta:\left(1-\widehat p\right)^\text{k}$$
We then use the two-sided p-value from the calculated z-statistic to determine whether there is enough evidence of a difference for that antibody to advance to phase 2 individual testing. If the p-value is less than the phase 1 cutoff (c1), we conduct individual testing; if it is greater, we assume there is no effect.
Phase 2
For those antibodies that proceed to Phase 2, we conduct a Fisher’s exact test and conclude there is a case-control difference if the resulting p-value is less than cutoff c2. The type I error rate of the final test is a function of both c1 and c2. Therefore, given c1, we need to determine c2 to control the final type I error rate at the nominal \(\alpha\) level.
Calibration of c2
We use a Monte-Carlo approach to compute c2 as a function of the antibody prevalence, by applying the two-phase design to data that was generated under the null hypothesis of no case-control effects, with 10,000 realizations for each prevalence value.
Figure 2 illustrates how to choose the p-value used in phase 2 testing to achieve a final \(\alpha\) level test. The figure shows the observed p-values in phase 2 testing under the null distribution, for an example antibody prevalence of 0.20. The \(1-\alpha\) percentile of the resulting phase 2 p-values determines the cutoff value c2. Rather than applying the Monte-Carlo procedure for each of the large number of antibodies (e.g., 15,000), we evaluate c2 as a function of prevalence by partitioning prevalence in units ranging from 0 to 1 by steps of size 0.01 (this requires only performing 100 Monte-Carlo simulations). Figure 3 shows the Monte-Carlo p-value cutoffs (c2) as a function of prevalence, and these values were used for phase 2 testing. Noting that the resulting curve was not continuous, we also applied Lowess smoothing in order to construct a continuous curve of c2 as a function of prevalence. However, smoothing the curve showed little differences in testing characteristics relative to simply interpolating between discrete sequence values so we used the non-continuous values for simplicity. We identified statistically significant antibody effects by comparing the p-value from the Fisher’s exact test to c2.
Standard group testing approach
An alternative to the two-phase design described above uses group testing to reconstruct the complete data. In this design, group testing is applied to the entire dataset in the following manner: for groups that are negative, we assume all individuals in that group are negative and for groups that are positive, we retest individual samples to reconstruct the individual data on which standing Fisher’s exact tests can be applied.
Simulation
We compare the proposed two-phase design with both a standard case-control and group testing design in terms of expected numbers of tests and statistical power. We generate data of 15,000 antibodies for 500 cases and 500 controls. The probability for a particular antibody j for individual i is given by
$${p}_{antibod{y}_{ij}}={\Phi }({\alpha }_{0j}+{\alpha }_{1j}{y}_{i}+{b}_{i})$$
$$\mathrm{with}\;{\mathrm\alpha}_{0\mathrm j}\overset{}{\sim\mathrm N}(-3,0.5)\;;\;\alpha_{1[1-200]}=0.73\;and\;\alpha_{1[201-\text{15,000}]}=0;$$
$$y_i=\left\{\begin{array}{cc}1&case\\0&control\end{array}\right.;b_i\overset{}{\sim N}(0,2)\;;\;\Phi\;\text{is the cumulative distribution function of the N}(0,1)\;\text{distribution}$$
With the resulting antibody probabilities (\(mean\left({p}_{antibod{y}_{ij}}\right)=0.1\)), we generate 15,000 antibody outcomes for each individual using a binomial distribution. The random effect \({b}_{i}\) incorporates an exchangeable correlation structure between antibody responses on the same individual.
The 0.73 in the above equation reflects the case-control differences on the probit scale for the first 200 antibodies. The remaining 14,800 antibodies have no case-control differences. In the following simulations we evaluate power based on the first 200 antibodies and type I error from the remaining antibodies.