Single-index logistic model for high-dimensional group testing data

1.   Introduction

Group testing, or pooled testing, was first introduced by Dorfman ^[1] to identify syphilis infections among U.S. Army personnel during World War II. This approach involves combining specimens (e.g., blood, plasma, urine, swabs) from multiple individuals and conducting a single test to check for infection. According to Dorfman's procedure, if the combined sample tests negative, all individuals in this sample can be confirmed disease-free. Conversely, a positive result necessitates further testing to identify the affected individuals. This strategy gained prominence during the COVID-19 pandemic ^[2,3,4] and has been applied to detect various infectious diseases, including HIV ^[5,6], chlamydia and gonorrhea ^[7], influenza ^[8], and the Zika virus ^[9]. The primary motivation for pooled testing lies in its economic efficiency; for instance, the State Hygienic Laboratory at the University of Iowa saved approximately fanxiexian_myfh3.1 million over five years by employing a modified Dorfman protocol for testing chlamydia and gonorrhea among residents of Iowa ^[10,11].

Despite its cost-effectiveness, group testing poses significant challenges for statistical analysis due to the absence of individual response data ^[12]. However, advancements in digital technology have provided access to rich covariate information, including demographic data, electronic health records, genomic data, lifestyle data, physiological monitoring data, imaging data, and environmental variables ^[13]. Integrating these covariates into various statistical models for group testing has been shown to enhance accuracy and robustness, as evidenced by studies from Mokalled et al. ^[14], Huang and Warasi ^[15], Haber et al. ^[16]. This integration leads to improved estimations of individual risk probabilities, thereby reducing the number of tests required and overall costs.

In managing covariates, single-index models offer advantages, such as less restrictive assumptions, good interpretability, and adaptability to high-dimensional data ^[17]. For high-dimensional single-index models, Radchenko ^[18] proposed a novel estimation method based on $L_1$ regularization, extending it to generalized linear models. Elmezouar et al. ^[19] developed a functional single index expectile model with a nonparametric estimator to address spatial dependency in financial data, showing strong consistency and practical applicability. Chen and Samworth ^[20] explored generalized additive models, deriving non-parametric estimators for each additive component by maximizing the likelihood function, and adapted this approach to generalized additive index models. Kereta et al. ^[21] employed a k-nearest neighbor estimator, enhanced by geodesic metrics, to extend local linear regression for single-index models. However, research on generalized semi-parametric single-index models in high-dimensional contexts remains limited, particularly in group testing applications, which are still underexplored.

Most current integrations of covariate information with group testing are developed based on parametric regression models. For example, Wang et al. ^[22] introduced a comprehensive binary regression framework, while McMahan et al. ^[11] developed a Bayesian regression framework. Gregory et al. ^[23] adopted an adaptive elastic net method, which remains effective as data dimensionality increases. Ko et al. ^[24] compared commonly used group testing procedures with group lasso regarding true positive selection in high-dimensional genomic data analysis. Furthermore, nonparametric regression methods have gained traction for applying covariates in group testing. Delaigle and Hall ^[25] proposed a nonparametric method for estimating conditional probabilities and testing specificity and sensitivity, addressing the unique dilution effects and complex data structures inherent in group testing. Self et al. ^[26] introduced a Bayesian generalized additive regression method to tackle dilution effects further, while Yuan et al. ^[12] developed a semiparametric monotone regression model using the expectation-maximization (EM) algorithm to navigate the complexities of group testing data. Zuo et al. ^[27] proposed a more flexible generalized nonparametric additive model, utilizing B-splines and group lasso methods for model estimation in high-dimensional data.

This article proposes a generalized single-index group testing model aimed at enhancing flexibility in addressing various nonlinear models and facilitating the selection of important variables. Given the absence of individual disease testing results in group testing data, the EM algorithm is employed to perform the necessary calculations for the model. B-spline functions are utilized to approximate the nonlinear unknown smooth functions, with model parameters estimated by maximizing the likelihood function. In modern group testing, a substantial amount of individual covariate information is typically collected during sample testing. Consequently, a penalty term is incorporated into the likelihood function, promoting the construction of a sparse model and enabling effective variable selection. We apply the method to four group testing strategies: master pool, Dorfman, halving, and array. The method is evaluated using both simulated and real data.

The remaining sections are organized as follows. Section 2 introduces our model with B-spline approximation, detailing the corresponding algorithm employing the EM algorithm. Section 3 elaborates on the E-step in the EM algorithm, facilitating the acceleration of our algorithm's convergence. Sections 4 and 5 present comprehensive simulations and real data application, demonstrating the method's robust performance. Finally, we conclude our findings and provide some discussion in Section 6.

2.   Primary results and methodological advancements

2.1.   Logistic single-index model for high-dimensional covariates

Consider a dataset comprising $n$ individuals. For each $i \in \{1, 2, \dots, n\}$, let the true disease status of the $i$-th individual be denoted by $\tilde{Y}_{i} \in \{0, 1\}$, where $\tilde{Y}_{i} = 1$ indicates disease presence, and $\tilde{Y}_{i} = 0$ indicates absence. Additionally, the dataset includes covariate information for each individual, represented as $\mathit{\boldsymbol{X}}_i = (X_{i1}, \ldots, X_{iq_n})^T \in \mathbb{R}^{q_n}$, where $\mathbb{R}^{q_n}$ denotes a $q_n$-dimensional real vector space. We assume the number of covariates $q_n$ is high-dimensional.

Let the risk probability for the $i$-th individual be defined as $p_i = \operatorname{Pr}(\tilde{Y}_i = 1 \mid \mathit{\boldsymbol{X}}_i)$, where $i \in \{1, 2, \dots, n\}$. In many cases, the influence of covariates may be nonlinear; imposing linearity can result in inaccurate estimations. This study explores nonlinear scenarios, assuming $p_i$ follows a flexible logistic single-index model, expressed as

where $\mathit{\boldsymbol{\beta}} = (\beta_1, \beta_2, \ldots, \beta_{q_n})^{\top} \in \mathbb{R}^{q_n}$ represents the unknown parameters, and $g(\cdot)$ is an unknown smooth function capturing the relationship between covariates and risk probabilities.

In semiparametric single-index models, the true parameters are generally considered non-identifiable without imposed constraints. To ensure the identifiability of $\mathit{\boldsymbol{\beta}}$, we impose a classical constraint: $\beta_1 = \sqrt{1 - \|\mathit{\boldsymbol{\beta}}_{-1}\|_2^2}$, where $\mathit{\boldsymbol{\beta}}_{-1} = (\beta_2, \beta_3, \ldots, \beta_{q_n})^{\top} \in \mathbb{R}^{q_{n}-1}$, and $\|\cdot\|_2$ denotes the $L_2$-norm. Note that both the function $g(\cdot)$ and the coefficient $\mathit{\boldsymbol{\beta}}$ in the single-index model are unknown. The $L_2$-norm constraint $\|\mathit{\boldsymbol{\beta}}\|_2 = 1$ is crucial for the identifiability of $\mathit{\boldsymbol{\beta}}$ as shown by Carroll et al. ^[28], Zhu et al. ^[29], Lin et al. ^[30], Cui et al. ^[31], and Guo et al. ^[32]. We assume that the true parameter $\mathit{\boldsymbol{\beta}}^{*}$ is sparse, defining the true model as $\mathcal{M}^{*} = \big\{j \in \{1, 2, \dots, q_n\} : \beta_{j}^{*} \neq 0 \big\}$.

For $i \in \{1, 2, \ldots, n\}$, $\tilde{Y}_i$ follows a Binomial distribution with parameter $p_i$, denoted as $\tilde{Y}_i \sim \operatorname{Binom}(1, p_i)$. In traditional single-index model studies, the true status $\tilde{\mathcal{Y}} = \{\tilde{Y}_i, \, i = 1, 2, \dots, n\}$, is directly observable. However, in group testing, $\tilde{\mathcal{Y}}$ is unobservable ^[33]. This paper investigates parameter estimation and statistical inference of single-index models based on group testing data. Moreover, if a group test result is positive, further testing is required to identify infected individuals. These results may depend on shared characteristics, leading to correlations within group test outcomes, complicating the modeling.

In group testing, we partition $n$ individuals into $J$ groups, denoted as $\mathcal{P}_{1, 1}, \mathcal{P}_{2, 1}, \dots, \mathcal{P}_{J, 1}$. Here, $\mathcal{P}_{j, 1}$ represents the initial index set of individuals for the $j$-th group, ensuring $\cup_{j = 1}^{J}\mathcal{P}_{j, 1} = \{1, 2, \dots, n\}$. For $j \in \{1, 2, \dots, J\}$, if any testing result for $\mathcal{P}_{j, 1}$ is positive, further testing may be warranted. Define $\mathcal{Z}_j = \{Z_{j, l}, \, l = 1, 2, \dots, L_j\}$ as the set of testing outcomes for the $j$-th group, where $L_j$ denotes the total number of tests conducted within $j$-th group. Each $Z_{j, l} \in \{0, 1\}$, where $Z_{j, l} = 0$ indicates a negative result and $Z_{j, l} = 1$ indicates a positive result. If $Z_{j, 1} = 0$, then $L_j = 1$; otherwise, $L_j \geq 1$. Let $\mathit{\boldsymbol{\mathcal{P}}}_j = \{\mathcal{P}_{j, l}, l = 1, 2, \dots, L_j\}$, where $\mathcal{P}_{j, l}$ corresponds to the individuals associated with $Z_{j, l}$. Define $\tilde{\mathcal{Z}}_j = \{\tilde{Z}_{j, l}, l = 1, 2, \dots, L_j\}$ as the true status corresponding to $\mathcal{Z}_j$. The true statuses of individuals determine the group's true status, defined as $\tilde{Z}_{j, l} = I\left(\sum_{i \in \mathcal{P}_{j, l}} \tilde{Y}_i\right)$, where $I(\cdot)$ denotes the indicator function.

In practical applications, measurement error of the test kits exists. We define $S_{e} = \operatorname{Pr}(Z_{j, l} = 1 \mid \tilde{Z}_{j, l} = 1)$ as sensitivity, representing the probability of correctly identifying positive samples, and $S_{p} = \operatorname{Pr}(Z_{j, l} = 0 \mid \tilde{Z}_{j, l} = 0)$ as specificity, denoting the probability of correctly identifying negative samples, where $l \in \{1, 2, \dots, L_j\}$ and $j \in \{1, 2, \dots, J\}$. According to the definitions of $S_e$ and $S_p$, given the true status $\tilde{Z}_{j, l}$, the group's testing results satisfy $Z_{j, l} |\tilde{Z}_{j, l} \sim \text{Binom}\left(1, \mathrm{S_e} ^{\tilde{Z}_{j, l}}\big(1- \mathrm{ S_p} \big)^{1-\tilde{Z}_{j, l}}\right)$.

Our approach is based on two widely accepted fundamental assumptions in group testing. The first assumption is that $S_e$ and $S_p$ are independent of group size, supported by various studies ^{[34,35,36,37]}. The second assumption posits that, given the true statuses of individuals in the $j$-th group $\{\tilde{Y}_i, i \in \mathcal{P}_{j, 1}\}$, the group's true statuses $\tilde{\mathcal{Z}}_j$ are mutually independent, as supported by previous research ^[23,34,35].

We apply our method to four group testing methods: master pool testing, Dorfman testing, halving testing, and array testing. Figure 1 illustrates the process of four testing methods: (a) Master pool testing, where a group of individuals (e.g., $\mathcal{P}_{j, 1}$ consisting of individuals 1, 2, 3, and 4) is tested as a whole to obtain the group testing result $Z_{j, 1}$; (b) Dorfman testing, where initially the same group testing as in master pool testing is conducted, and if the result of the master pool testing is positive ($Z_{j, 1} = 1$), each individual in the group is then tested separately to obtain individual testing results $Z_{j, 2}$, $Z_{j, 3}$, $Z_{j, 4}$, and $Z_{j, 5}$; (c) Halving testing, where the entire group (e.g., $\mathcal{P}_{j, 1}$) is tested as a whole, and if the result is positive ($Z_{j, 1} = 1$), the group is divided into two subgroups (e.g., $\mathcal{P}_{j, 2}$ and $\mathcal{P}_{j, 3}$) for subgroup testing, and if the result of subgroup testing is positive (e.g., $Z_{j, 2} = 1$), individuals in the positive subgroup are then tested individually; and (d) Array testing, where multiple individuals (e.g., 16 individuals) are arranged in an array for group testing to obtain multiple group testing results such as $Z_{j, 1}$, and if a specific group testing result is positive (e.g., $Z_{j, 1} = 1$), further subgroup testing is performed (e.g., obtaining results $Z_{j, 2}$, $Z_{j, 3}$), and if the group testing results for both the row and column where an individual is located are positive (e.g., $Z_{j, 3} = Z_{j, 4} = Z_{j, 7} = 1$), the individuals (e.g., $6$-th individual and $10$-th individual) are then tested.

Due to the nature of group testing, the true status of individuals, denoted as $\tilde{\mathcal{Y}}$, remains unknown. Our objective is to estimate $\mathit{\boldsymbol{M}}^*$, $\mathit{\boldsymbol{\beta}}^*$, and $g(\cdot)$ based on observed data $\mathcal{Z} = \{\mathcal{Z}_{j}, \, j = 1, 2, \ldots, J\}$ and covariate information $\mathbb{X} = (\mathit{\boldsymbol{X}}_{1}, \mathit{\boldsymbol{X}}_{2}, \ldots, \mathit{\boldsymbol{X}}_{n})^T \in \mathbb{R}^{n \times q_n }$ to ascertain individual risk probabilities. The likelihood function based on the observed data $\mathcal{Z}$ is defined as

where

and $\tilde{\mathcal{Y}}_{\mathcal{P}_{j, l}} = \{\tilde{Y}_i, \, i \in \mathcal{P}_{j, l}\}$ represents the set of true statuses for individuals in $\mathcal{P}_{j, l}$. Furthermore, the conditional probability $P\left(Z_{j, l} \mid \tilde{\mathcal{Y}}_{\mathcal{P}_{j, l}}\right)$ is expressed as

The likelihood function for the true disease status $\tilde{\mathcal{Y}}$ can be written as

Combining this with the logistic single-index model defined in (2.1), we obtain the log-likelihood function for $\tilde{\mathcal{Y}}$:

Since the smooth function $g(\cdot)$ is unknown, we approximate it using B-spline functions. Let the support interval of $g(\cdot)$ be $[a, b]$. We partition $[a, b]$ at points $a = d_{0} < d_{1} < \ldots < d_{N} < b = d_{N+1}$ into several segments, referred to as knots or internal nodes. This division generates subintervals $I_{k} = [d_{k}, d_{k+1})$ for $0 \leq k \leq N-1$ and $I_{N} = [d_{N}, d_{N+1}]$, ensuring that

where $M \in (0, \infty)$. The B-spline basis functions of order $q$ are denoted as $\mathit{\boldsymbol{\Phi}}(\cdot) = \Big(\phi_1(\cdot), \phi_2(\cdot), \ldots, \phi_S(\cdot)\Big)^{\top} \in \mathbb{R}^S$, with $S = N + q$. Thus, $g(\cdot)$ can be approximated as

where $\gamma_s$ are the spline coefficients to be estimated ^[38]. Denote $\mathit{\boldsymbol{\gamma}} = (\gamma_{1}, \gamma_{2}, \ldots, \gamma_{S})^{\top} \in \mathbb{R}^S$. We approximate $g(\mathit{\boldsymbol{X}}_i^T \mathit{\boldsymbol{\beta}})$ as

where $\mathit{\boldsymbol{\Phi}}(\mathit{\boldsymbol{X}}_{i}^{\top}\mathit{\boldsymbol{\beta}}) = \Big(\phi_1(\mathit{\boldsymbol{X}}_{i}^{\top}\mathit{\boldsymbol{\beta}}), \phi_2(\mathit{\boldsymbol{X}}_{i}^{\top}\mathit{\boldsymbol{\beta}}), \ldots, \phi_S(\mathit{\boldsymbol{X}}_{i}^{\top}\mathit{\boldsymbol{\beta}})\Big)^{\top}$. Therefore, we approximate $p_i$ by using a spline function, and denote the spline approximation of $p_i$ as $p_{iB}$, which is defined as follows:

In the following, we use the spline approximation $p_{iB}$ of $p_{i}$ to construct the log-likelihood function and the objective function in the subsequent EM algorithm. Thus, the log-likelihood function (2.3) for $\tilde{\mathcal{Y}}$ can be reformulated as

Furthermore, the target likelihood function (2.2) can be represented as

By employing spline approximation, we transform the estimation problem of $\mathit{\boldsymbol{\beta}}^{*}_{-1}$ and $g(\cdot)$ into estimating $\mathit{\boldsymbol{\beta}}^{*}_{-1}$ and $\mathit{\boldsymbol{\gamma}}$.

For high-dimensional group testing data, we aim to estimate $\mathit{\boldsymbol{\beta}}^{*}_{-1}$ using the penalized approach within a single-index model framework. The penalized log-likelihood function is defined as follows:

where $P_{\lambda}(\cdot)$ is the penalty function and $\lambda$ is a tuning parameter. We consider three common penalty functions: LASSO ^[39], SCAD ^[40], and MCP ^[41]. Specifically, for LASSO, $P_{\lambda}(x) = \lambda |x|$. For SCAD, it is defined as

where $\delta > 2$. In MCP, the penalty function is given by

with $\delta > 1$. The following section will detail the parameter estimation process.

2.2.   EM algorithm for regularized single-index model in group testing

The penalized log-likelihood function (2.5) lacks the individual latent status $\tilde{\mathcal{Y}}$. The complete data penalized log-likelihood function can be expressed as

Notably, $\ln P(\mathcal{Z}|\tilde{\mathcal{Y}})$ depends solely on known parameters $S_e$ and $S_p$, allowing us to disregard it in computations. The presence of the latent variable $\tilde{\mathcal{Y}}$ complicates direct maximization of the complete data penalized log-likelihood function (2.6). Therefore, we employ the EM algorithm, comprising two steps: the Expectation (E) step, and the Maximization (M) step.

In the E step, given the observed data $\mathcal{Z}$ and the parameters from the $t$-th iteration $(\mathit{\boldsymbol{\beta}}^{(t)}_{-1}, \mathit{\boldsymbol{\gamma}}^{(t)})$, calculate the following function:

where $w_i^{(t)} = \mathbb{E}[\tilde{Y}_{i}|\mathcal{Z}, \mathit{\boldsymbol{\gamma}}^{(t)}, \mathit{\boldsymbol{\beta}}^{(t)}_{-1}], \; i = 1, 2, \dots, n$. The calculation of the $w_i^{(t)}$ varies among the four grouping testing methods, which will be discussed in Section 3.

In the M step, we update $\mathit{\boldsymbol{\beta}}_{-1}^{(t+1)}$ and $\mathit{\boldsymbol{\gamma}}^{(t+1)}$, respectively. Initially, we update $\mathit{\boldsymbol{\gamma}}^{(t+1)}$ by maximizing:

Subsequently, we maximize $S^{(t)}(\mathit{\boldsymbol{\beta}}_{-1}, \mathit{\boldsymbol{\gamma}}^{(t+1)})$ to update the parameters $\mathit{\boldsymbol{\beta}}_{-1}^{(t+1)}$:

Given that $\mathit{\boldsymbol{\beta}}_{-1}$ appears in each B-spline basis function $\phi(\mathit{\boldsymbol{X}}_{i}^{\top}\mathit{\boldsymbol{\beta}})$, direct iteration presents challenges. Let $\tilde{g}^{(t)}(\mathit{\boldsymbol{X}}_i^{T}\mathit{\boldsymbol{\beta}}) = \mathit{\boldsymbol{\Phi}}^{\top} (\mathit{\boldsymbol{X}}_i^{T}\mathit{\boldsymbol{\beta}}) \mathit{\boldsymbol{\gamma}}^{(t+1)}$. We apply the approach by Guo et al. ^[42], approximating $\tilde{g}^{(t)}(\mathit{\boldsymbol{X}}_i^{T}\mathit{\boldsymbol{\beta}})$ via a first-order Taylor expansion

where $\mathit{\boldsymbol{J}}(\mathit{\boldsymbol{\beta}}) = \partial \mathit{\boldsymbol{\beta}} / \partial \mathit{\boldsymbol{\beta}}_{-1} = \left(-\mathit{\boldsymbol{\beta}}_{-1} / \sqrt{1 - \|\mathit{\boldsymbol{\beta}}_{-1}\|_2^2}, \mathit{\boldsymbol{I}}_{q_n-1} \right)^{\top}$ represents the Jacobian matrix of size $q_n \times (q_n - 1)$ and $\mathit{\boldsymbol{I}}_{q_n-1}$ denotes the $(q_{n}-1)$-dimensional identity matrix. This approximation is incorporated into $S^{(t)}(\mathit{\boldsymbol{\beta}}_{-1}, \mathit{\boldsymbol{\gamma}}^{(t+1)})$ to maximize the expression and update $\mathit{\boldsymbol{\beta}}_{-1}^{(t+1)}$. Therefore, we approximate $S^{(t)}(\mathit{\boldsymbol{\beta}}_{-1}, \mathit{\boldsymbol{\gamma}}^{(t+1)})$ by $\tilde{S}^{(t)}(\mathit{\boldsymbol{\beta}}_{-1}, \mathit{\boldsymbol{\gamma}}^{(t+1)})$ as follows:

We employ stochastic gradient descent ^[43] and coordinate descent ^[44] to update $\mathit{\boldsymbol{\gamma}}$ and $\mathit{\boldsymbol{\beta}}$, respectively. Let $\hat{\mathit{\boldsymbol{\gamma}}}$ and $\hat{\mathit{\boldsymbol{\beta}}}_{-1}$ denote the estimated parameters, and $\hat{\mathcal{M}} = \big\{j \in \{1, 2, \dots, q_n\} : \hat{\beta}_j \neq 0 \big\}$ represent the estimated model. Furthermore, $\hat{\mathit{\boldsymbol{\gamma}}}$ and $\hat{\mathit{\boldsymbol{\beta}}}_{-1}$ can be used to calculate individual risk probabilities and guide subsequent testing strategies. In summary, the EM algorithm offers a structured approach to handle the latent variable $\tilde{\mathcal{Y}}$ and estimate model parameters. The detailed steps of this method are summarized in Algorithm 1.

3.   Calculation of conditional expectations

Implementing Algorithm 1 requires deriving formulas to calculate the conditional expectations of individuals' true statuses. These expressions are essential for the effective application of the EM algorithm in various testing scenarios. Common group testing methods include master pool testing, Dorfman testing, halving testing, and array testing. We have derived the conditional expectation formula of these methods under our methodological framework, which will facilitate our other calculations.

For master pool testing, samples are divided into $J$ distinct groups, with each sample assigned to only one group, and each group undergoes a single test without subsequent testing. When the $i$-th individual is assigned to the $j$-th group, consider two cases for $w_i^{(t)}$:

While $Z_{j} = 0$,

Due to

let $\Delta_{j} = S_e+(1-S_e-S_p)\prod\limits_{i\in \mathcal{P}_j}(1-p^{(t)}_{iB})$, where $p_{iB}$ is an approximate result of $p_{i}$ in (2.4). Therefore,

Then,

While $Z_{j} = 1$,

In conclusion,

We apply our method to four group testing algorithms: master pool testing, Dorfman testing, halving testing, and array testing. For other algorithms, detailed expressions can be found in Appendix C. Using these expressions, we apply the EM algorithm to estimate the model parameters.

4.   Simulation study

In this section, we assess the performance of the proposed method using simulated datasets. The generation of covariates follows the approach described by Guo et al. ^[42]. Specifically, covariates $\mathbb{X} \in \mathbb{R}^{n \times q_n}$ are drawn from a truncated multivariate normal distribution. We first generate covariates from $N(0, \mathit{\boldsymbol{\Sigma}})$, where $\mathit{\boldsymbol{\Sigma}} \in \mathbb{R}^{q_n \times q_n}$ and $\Sigma_{ij} = 0.5^{|i-j|}$ for $1 \leq i, j \leq q_n$. These covariates are then truncated to the range $(-2, 2)$ to obtain $\mathbb{X}$. We consider logistic single-index models to describe $p_i = \operatorname{Pr}(\tilde{Y}_i = 1 \mid \mathit{\boldsymbol{X}}_i)$, with the function $g(\mathit{\boldsymbol{X}}_i^{\top} \mathit{\boldsymbol{\beta}})$ in the model (2.1) defined as follows,

Example 4.1. We set $n = 500$ and $\mathit{\boldsymbol{\beta}}^* = \left(\frac{3}{\sqrt{15.25}}, \frac{2.5}{\sqrt{15.25}}, 0, \ldots, 0\right)^{\top}$. We consider two scenarios: $q_n = 50$ and $q_n = 100$. The model is described as follows:

Under this setting, the disease prevalence is approximately 8.93%.

Example 4.2. We set $n = 1000$ and $\mathit{\boldsymbol{\beta}}^* = \left(\frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}}, \frac{1}{\sqrt{3}}, 0, \ldots, 0\right)^{\top}$. We consider two scenarios: $q_n = 100$ and $q_n = 500$. The model is described as follows:

In this example, the disease prevalence is approximately 11.41%.

Example 4.3. We set $q_n = 50$ and $\mathit{\boldsymbol{\beta}}^* = \left(\frac{9}{\sqrt{181}}, \frac{8}{\sqrt{181}}, \frac{6}{\sqrt{181}}, 0, \ldots, 0\right)^{\top}$. We consider two scenarios: $n = 500$ and $n = 1000$. The model is described as follows:

In this example, the disease prevalence is approximately 9.42%.

Example 4.4. We set $q_n = 100$ and $\mathit{\boldsymbol{\beta}}^* = (0.5, 0.5, 0.5, 0.5, 0, \ldots, 0)^{\top}$. Two scenarios are considered: $n = 750$ and $n = 1000$. The model is described as follows:

In this scenario, the disease prevalence is approximately 10.32%.

In our simulation study, we employed four group testing algorithms: master pool testing (MPT), Dorfman testing (DT), halving testing (HT), and array testing (AT) to evaluate the model. For MPT, DT, and HT, the group size was set to 4, while in AT, individuals were arranged in a $4 \times 4$ array. Both sensitivity and specificity were fixed at $S_e = S_p = 0.98$. Based on the methodologies of Fan and Li ^[40] and Zhang ^[41], we set $\delta$ values of 3.7 and 2 for SCAD and MCP, respectively. Each scenario was simulated $B = 100$ times, where $\hat{\mathit{\boldsymbol{\beta}}}^{[b]}$ denotes the estimated $\mathit{\boldsymbol{\beta}}^{*}$ in the $b$-th simulation, with $b \in \{1, 2, \dots, B\}$.

Following the approach of Guan et al. ^[45], we measured the estimation accuracy of $\hat{\beta}_j$ ($j = 1, 2, 3, 4$) using the mean squared error (MSE), defined as

We utilized average mean squared error (AMSE) to assess the accuracy of $\hat{\mathit{\boldsymbol{\beta}}}$, consistent with methods employed by Wang and Yang ^[46]:

Average mean absolute error (AMAE) was used to evaluate the estimation performance of $g(\cdot)$ and individual risk probabilities $p_i$ ^[42]. The AMAE for $g(\cdot)$ is defined as

while the AMAE for $\hat{p}^{[b]}_{i} = \frac{e^{g(\mathit{\boldsymbol{X}}_i^{\top} \hat{\mathit{\boldsymbol{\beta}}}^{[b]})}}{1 + e^{g(\mathit{\boldsymbol{X}}_i^{\top} \hat{\mathit{\boldsymbol{\beta}}}^{[b]})}}$ is defined as

where $p^{*}_{i} = \frac{g(\mathit{\boldsymbol{X}}_i^{\top} \mathit{\boldsymbol{\beta}}^{*})}{1 + e^{g(\mathit{\boldsymbol{X}}_i^{\top} \mathit{\boldsymbol{\beta}}^{*})}}$.

To evaluate variable selection performance, we employed true positive rate (TPR) and false positive rate (FPR). The FPR represents the proportion of false positives among identified predictors, while the TPR indicates the proportion of true positives among relevant predictors. Table 1 shows the results of variable selection. TPR and FPR are defined as follows:

The simulation results are summarized in Tables 2 to 5. As shown in the tables, the TPR was approximately 97%, with a very low FPR. The result shows that the probability that $\mathcal{M}^{*}$ is contained in $\hat{\mathcal{M}}$ is very close to $1$. This demonstrates the notable performance of our model in variable selection. The AMAE for $g(\cdot)$ and $p_i$ was approximately 0.5 and 0.01, respectively. This shows that we have accurately captured the form of the unknown smooth function $g(\cdot)$ and are able to precisely predict the individual risk probability. The AMSE for the model parameters $\mathit{\boldsymbol{\beta}}$ was around $10^{-4}$, while the AMSE for significant variables $\beta_j$ was approximately $10^{-3}$. This demonstrates the accuracy of our model in parameter estimation.

We set up two different sample sizes ($n$) or covariate scenarios ($q_n$) for each example. Results of Examples 4.1 and 4.2 suggest that our method maintains robust estimation performance as dimensionality increases in small sample scenarios. Furthermore, results of Examples 4.3 and 4.4 demonstrate that estimation accuracy improves with increased sample size. Figure 2 illustrates the estimation performance of $g(\cdot)$ and individual risk probabilities $p_i$, confirming our method's efficacy in estimating unknown functions and risk probabilities.

Moreover, we aim to evaluate our method's performance under different group sizes. Using Example 4.4, we investigated group sizes of 2, 4, 6, and 8 with the Dorfman algorithm and LASSO penalty function. Results are presented in Table 6, reporting the means of $\hat{\beta}_j$ for $j = 1, 2, 3, 4$. The simulation results indicate that our method consistently delivers strong estimation performance across various group sizes. At the same time, we set up comparative experiments with different $S_e$ and $S_p$, and the simulation results are shown in Tables 8 to 11 in Appendix A. As shown in these tables, our model maintains a certain level of stability, ensuring that $\mathcal{M}^{*}$ is still contained within $\hat{\mathcal{M}}$.

5.   Application to real data

In this section, we validate the effectiveness of our method using the diabetes dataset from the National Health and Nutrition Examination Survey (NHANES) conducted between 1999 and 2004. NHANES is a probability-based cross-sectional survey representing the U.S. population, collecting demographic, health history, and behavioral information through household interviews. Participants were also invited to equip mobile examination centers for detailed physical, psychological, and laboratory assessments. The dataset is accessible at https://wwwn.cdc.gov/Nchs/Nhanes/.

The dataset comprises $n = 5515$ records and 17 variables, categorizing individuals as diabetic or non-diabetic. Covariates include age ($\mathit{\boldsymbol{X}}_1$), waist circumference ($\mathit{\boldsymbol{X}}_2$), BMI ($\mathit{\boldsymbol{X}}_3$), height ($\mathit{\boldsymbol{X}}_4$), weight ($\mathit{\boldsymbol{X}}_5$), smoking age ($\mathit{\boldsymbol{X}}_6$), alcohol use ($\mathit{\boldsymbol{X}}_7$), leg length ($\mathit{\boldsymbol{X}}_8$), total cholesterol ($\mathit{\boldsymbol{X}}_9$), hypertension ($\mathit{\boldsymbol{X}}_{10}$), education level ($\mathit{\boldsymbol{X}}_{11}$), household income ($\mathit{\boldsymbol{X}}_{12}$), family history ($\mathit{\boldsymbol{X}}_{13}$), physical activity ($\mathit{\boldsymbol{X}}_{14}$), gender ($\mathit{\boldsymbol{X}}_{15}$), and race ($\mathit{\boldsymbol{X}}_{16}$). Notably, nominal variables from $\mathit{\boldsymbol{X}}_{10}$ to $\mathit{\boldsymbol{X}}_{16}$ are transformed using one-hot encoding, resulting in $q_n = 47$ covariates per individual. The first nine variables are continuous, while the remainder are binary. A detailed explanation of the variables as well as the content of the questionnaire can be found at https://wwwn.cdc.gov/Nchs/Nhanes/search/default.aspx. For convenience, the nominal variables are explained in Table 12. in Appendix B.

For $i \in \{1, 2, \ldots, n\}$, we define $\tilde{Y}_i = 1$ for diabetes and $\tilde{Y}_i = 0$ for non-diabetes. Individual covariate information is represented as $\mathit{\boldsymbol{X}}_i = (X_{i1}, X_{i2}, \ldots, X_{iq_n})^{\top}$. We construct the following single-index model for the probability of diabetes risk for the $i$-th individual:

where the smooth function $g(\cdot)$ is unknown, and our objective is to estimate the coefficients $\mathit{\boldsymbol{\beta}}$.

To verify the accuracy of our method, we compare the results with those obtained from two other methods. The first method is penalized logistic regression (PLR), which uses the true individual status, $\tilde{Y}_i$. This method is implemented using the R package "glmnet". The second method is the adaptive elastic net for group testing (aenetgt) data, as introduced by Gregory et al. ^[23]. This approach utilizes group testing data and employs a penalized Expectation-Maximization (EM) algorithm to fit an adaptive elastic net logistic regression model. The R package "aenetgt" is used for implementation. We generate Dorfman group testing data with a group size of 6, setting both sensitivity and specificity at $S_e = S_p = 0.98$.

To ensure comparability, we adhere to the standardization techniques referenced in Cui et al. ^[31]. First, we center the covariates to facilitate the comparison of relative effects across different explanatory variables. Second, we normalize the PLR and aenetgt coefficients by dividing them by their $L_2$-norm, as follows:

thereby obtaining coefficients with unit norm. This enables a comparison of regression coefficients from PLR, aenetgt, and the single-index group testing model.

The estimated coefficients from the three models are summarized in Table 7, and the parameter estimation of our method is denoted as $\hat{\beta}_{our}$. In this study, the estimated coefficients for age, $\hat{\beta}^{norm}_{PLR}$ and $\hat{\beta}_{our}$, are 0.280 and 0.307, respectively, indicating that the risk of diabetes increases with age, consistent with the findings of Turi et al. ^[47]. However, the coefficient $\hat{\beta}^{norm}_{aenet}$ is close to zero. For waist circumference, the coefficients $\hat{\beta}^{norm}_{PLR}$, $\hat{\beta}_{our}$, and $\hat{\beta}^{norm}_{aenet}$ are 0.178, 0.194, and 0.271, respectively, suggesting a positive association between waist circumference and diabetes risk, which is supported by Bai et al. ^[48] and Snijder et al. ^[49]. In addition, all three methods also identified leg length ^[50], hypertension ^[51], race ^[52], family history ^[53], and sex ^[54] as variables associated with diabetes. These covariates are widely recognized as being related to diabetes in the biomedical field ^[55].

We found that the covariable physical activity is associated with diabetes, but the aenetgt method failed to identify this association. The results of a study by Yu et al. ^[55], which used the same dataset as ours, are consistent with this finding. In addition, we found that education level was also a covariable associated with diabetes ($\hat{\beta}^{norm}_{PLR}$ and $\hat{\beta}_{our}$ are -0.523 and -0.335). Evidence for this association can also be found in the study by Aldossari et al. ^[56], and in this dataset, the probability that these participants will not develop diabetes is 100%. We also identified that household income is associated with diabetes, which is consistent with the study by Yen et al. ^[57]. In this dataset, the probability of developing diabetes for those who refused to answer about their household income is 60%. Furthermore, our model yields results similar to those obtained by the PLR method, which uses individual observations ($\tilde{\mathcal{Y}}$), suggesting that our method is able to extract information from group observations.

6.   Conclusions and discussion

This study presents a group testing framework based on a logistic regression single-index model for disease screening in low-prevalence environments. By employing B-splines to estimate unknown functions and incorporating penalty functions, our approach achieves high flexibility in capturing the relationships between covariates and individual risk probabilities while accurately identifying important variables. To address potential computational challenges in individual disease status estimation, we implemented an iterative EM algorithm for model estimation. Our simulation experiments demonstrate the proposed method's performance in high-dimensional covariate contexts with limited sample sizes, while application to real data confirms its efficacy. Our framework offers a unified approach for various group testing methods, showcasing its practical application value.

Despite these promising outcomes, our study acknowledges several limitations. First, our model assumes that sensitivity and specificity of testing are independent of group size, which may not always hold in practical applications. Second, data quality and variations in the testing population can impact the model's applicability. Therefore, exploring how to integrate prior information to enhance model accuracy and practical value remains a critical research direction. Furthermore, the potential high dimensionality of individual covariates poses significant challenges, necessitating the development of models capable of handling ultra-high-dimensional data.

Future research could explore the following directions. Firstly, examining model performance under varying group testing configurations, such as changes in testing errors and group sizes, could yield valuable insights. Secondly, investigating methods to incorporate additional prior knowledge to improve estimation accuracy is a worthwhile endeavor. Additionally, considering computational efficiency, developing faster algorithms for processing large-scale datasets will be a key focus for future work.

Author contributions

Changfu Yang: Methodolog, formal analysis, writing-original draft; Wenxin Zhou: Methodology, formal analysis; Wenjun Xiong: Conceptualization, methodology, writing-original draft, funding acquisition; Junjian Zhang: Conceptualization, methodology, writing-review and editing, funding acquisition; Juan Ding: Conceptualization, formal analysis, writing-review and editin, funding acquisition. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (Grant Nos. 12361055, 11801102), Guangxi Natural Science Foundation (2021GXNSFAA220054), and the Fundamental Research Funds for the Central Universities (B240201095).

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

A.   Appendix-A: Simulation results with different sensitivity and specificity settings

In this part, we tested the performance of four examples at different sensitivity and specificity, using the Dofman algorithm and the LASSO penalty function. The simulation results are shown in Tables 8 to 11.

B.   Appendix-B: Meaning of the nominal variable

C.   Appendix-C: Calculation of conditional expectations

In this part, we derive the conditional expectation formulas for Dorfman testing, halving testing, and array testing within the framework of our method. Before proceeding, it is necessary to clarify some notations. Let $\mathcal{P}_j \setminus \{i\}$ represent the set of individuals in $\mathcal{P}_j$ excluding the $i$-th individual, and $|\mathcal{P}_j|$ denotes the number of individuals in $\mathcal{P}_j$. Let $Y_i$ represent the test result of the $i$-th individual and $\mathcal{Y}_{\mathcal{P}_{j, l}} = \{Y_i, i \in \mathcal{P}_{j, l}\}$ represent the set of testing results of individuals in $\mathcal{P}_{j, l}$.

C.1.   Dorfman testing

If the initial group testing result is negative, no re-testing is performed. However, if $Z_{j, 1} = 1$, each individual in the group needs to undergo a separate re-testing.

1) When $Z_{j, 1} = 0$, the result is the same as the master poor testing:

2) When $Z_{j, 1} = 1$, each individual in the group must undergo a separate re-test. In total, the group has undergone $|\mathcal{P}_j| + 1$ tests.

The denominator is

Thus, the numerator is

Therefore, the final expression is

C.2.   Halving testing

Assume that the maximum number of partitions required during testing is two. Let the test result of the first testing be $Z_{j, 1}$. At this time, the set of all unpartitioned individuals is $\mathcal {P}_{j, 1}$, which contains $|\mathcal{P}_j|$ individuals. After the first partition, the partitioning method is to divide into two equal parts, with the two subsets of individuals being $\mathcal{P}_{j, 2}$ and $\mathcal{P}_{j, 3}$, respectively. The responses of the second testing are $Z_{j, 2}$ and $Z_{j, 3}$. There are five types of testing results in halving testing.

1) When $Z_{j, 1} = 0$:

Only one testing is performed, and the process is the same as master pool testing. Since the result of one testing is negative, no further partitioning and testing are performed. At this time,

2) When $Z_{j, 1} = 1, Z_{j, 2} = 0, Z_{j, 3} = 0$:

That is, the result of the first testing is $Z_{j, 1} = 1$. Subsequently, the first partition is performed, dividing into two equal parts $\mathcal{P}_{j, 2}$ and $\mathcal{P}_{j, 3}$. Then, testings are performed on the two sets respectively, with the testing results being $Z_{j, 2} = Z_{j, 3} = 0$. At this time,

The denominator is

Since the placement of the $i$-th individual in the sets $\mathcal{P}_{j, 2}$ and $\mathcal{P}_{j, 3}$ is symmetric, assume that $i$-th individual is placed in the set $\mathcal{P}_{j, 2}$. Then, the numerator is

3) When $Z_{j, 1} = 1, Z_{j, 2} = 0, Z_{j, 3} = 1$:

At this time, the second partitions are performed. The first partition divides all individuals into two sets, $\mathcal{P}_{j, 2}$ and $\mathcal{P}_{j, 3}$, with testing results $Z_{j, 2} = 0$ and $Z_{j, 3} = 1$, respectively. Individual testings are performed separately on the individuals in $\mathcal{P}_{j, 3}$, and the set of testing results is $\mathcal{Y}_{\mathcal{P}_{j, 3}}$. At this time,

The denominator is

Since an $i$-th individual may belong to either set $\mathcal{P}_{j, 2}$ or $\mathcal{P}_{j, 3}$, the numerator is discussed accordingly.

(a) Assume that $i$-th individual belongs to set $\mathcal{P}_{j, 2}$. Then, the numerator is

(b) Assume that $i$-th individual belongs to set $\mathcal{P}_{j, 3}$. Then, the numerator is

4) When $Z_{j, 1} = 1, Z_{j, 2} = 1$, and $Z_{j, 3} = 0$, the process is the same as when $Z_{j, 1} = 1, Z_{j, 2} = 0$, and $Z_{j, 3} = 1$, and the numerator needs to be discussed accordingly. At this time,

First, the denominator is

Next, the numerator is discussed.

(a) Assume that $i$-th individual belongs to set $\mathcal{P}_{j, 2}$. Then, the numerator is

(b) Assume that $i$-th individual belongs to set $\mathcal{P}_{j, 3}$. Then, the numerator is

5) When $Z_{j, 1} = 1, Z_{j, 2} = 1$, and $Z_{j, 3} = 1$, two similar partitions are performed as above, and individual retests are conducted separately for all individuals in $\mathcal{P}_j$. At this time, $\mathcal{Y}_{\mathcal{P}_{j}} = \mathcal{Y}_{\mathcal{P}_{j, 2}} \cup \mathcal{Y}_{\mathcal{P}_{j, 3}}$, and we have

The denominator is

The results of $i$-th individual belonging to either set $\mathcal{P}_{j, 2}$ or $\mathcal{P}_{j, 3}$ are symmetric. Assume that $i$-th individual belongs to set $\mathcal{P}_{j, 2}$. Then, the numerator is

C.3.   Array testing

For convenience, assume that the set of all individuals is $G$, and all individuals can be arranged into an $R \times C$ array, that is, $G = \left\{(r, c), r \in R, c \in C\right\}$. Define $\mathcal{R} = (R_1, R_2, \cdots, R_R)$ and $\mathcal{C} = (C_1, C_2, \cdots, C_C)$ as the collections of row and column testing results, respectively. Let $R = \max R_{r}$ and $C = \max C_{c}$. Furthermore, define $\tilde{R}_r = \max_c \tilde{Y}_{rc}$ and $\tilde{C}_c = \max_r \tilde{Y}_{rc}$ as the true result sets for rows and columns, respectively. Let $Y_{rc}$ denote the testing result of the individual in the $r$-th row and $c$-th column of the array, and $\tilde{Y}_{rc}$ represents the true disease status of the individual in the $r$-th row and $c$-th column of the array. Let

$\mathcal{Y}_Q$ represents the collection of responses from all potentially positive individuals, and $\tilde{\mathcal{Y}}_Q$ denotes the true disease statuses of all potentially positive individuals. Let $\mathcal{Z}_G = (R, C)$ denote the group testing responses. Since $(r, c) \in G$, define

Then,

1) When $\mathcal{Z}_G = (0, 0)$, there is no need to retest individuals within the group. At this time,

The denominator is

The numerator is

2) When $\mathcal{Z}_G \neq (0, 0)$, $\mathcal{Z}_G = (R, C)$ has multiple scenarios, specifically $\mathcal{Z}_G = (R, C) = (1, 0)$, $(R, C) = (0, 1)$, and $(R, C) = (1, 1)$. Therefore, when $\mathcal{Z}_G \neq (0, 0)$, the following classifications can be discussed:

(a) When $(R, C) = (1, 0)$,

The denominator is

At this point, the numerator requires further discussion:

$(i)$ If $(r, c) \in Q$ and $R_r = 1$ and $C_r = 0$, then

$(ii)$ If $(r, c) \notin Q$, then $\mathcal{R} = 0$, $\mathcal{C} \neq 0$, but $R_r = 0$ and $C_r = 0$:

(b) When $(R, C) = (0, 1)$, the denominator is

The numerator requires further discussion:

(ⅰ) If $(r, c) \in Q$ and $R_r = 0$ and $C_c = 1$, then

(ⅱ) If $(r, c) \notin Q$, then $\mathcal{R} = 0$, $\mathcal{C} \neq 0$, but $R_r = 0$ and $C_r = 0$:

(c) When $(R, C) = (1, 1)$, the denominator is

For the numerator, we provide the following derivations:

(ⅰ) If $(r, c) \in Q$ and $R_r = 1$ and $C_c = 1$, then

(ⅱ) If $(r, c) \notin Q$, then $\mathcal{R} \neq 0$, $\mathcal{C} \neq 0$, but $R_r = 0$ and $C_r = 0$: