Springer: An R package for bi-level variable selection of high-dimensional longitudinal data

In high-dimensional data analysis, the bi-level (or the sparse group) variable selection can simultaneously conduct penalization on the group level and within groups, which has been developed for continuous, binary, and survival responses in the literature. Zhou et al. (2022) (PMID: 35766061) has further extended it under the longitudinal response by proposing a quadratic inference function-based penalization method in gene–environment interaction studies. This study introduces “springer,” an R package implementing the bi-level variable selection within the QIF framework developed in Zhou et al. (2022). In addition, R package “springer” has also implemented the generalized estimating equation-based sparse group penalization method. Alternative methods focusing only on the group level or individual level have also been provided by the package. In this study, we have systematically introduced the longitudinal penalization methods implemented in the “springer” package. We demonstrate the usage of the core and supporting functions, which is followed by the numerical examples and discussions. R package “springer” is available at https://cran.r-project.org/package=springer.


A.1 Penalized Group QIF
The penalized group QIF implemented in package springer has the score equation defined below: U (β n ) = Q(β n ) + p k=1 ρ(||η nk || Σ k ; λ 1 , γ), where ρ denotes the minimax concave penalty with tuning parameters λ 1 and γ. Based on Section 2.2, the coefficient vector β n represents all the main and interaction effects. η n k denotes the (q+1) main genetic and G×E interaction effects with respect to the kth G factor. Regularized identification has been only conducted on the group level because the shrinkage is imposed on ||η nk || Σ k , the empirical norm of η k .
Again, the Newton-Raphson algorithm can facilitate model fitting to obtainβ n , the penalized QIF estimate. At the (g + 1)th iteration, we solve forβ g+1 n usingβ g n from the previous iteration as follows: , with P (β g n ) and V (β g n ) being defined as: In addition, H(β g n ) is a diagonal matrix with the derivatives of MCP being the diagonal elements as follows: where λ 1 is the tuning parameter controlling the sparsity of the group level effects, and γ is the regularization parameter determining the concavity of MCP. In Section 2.4, a factor of √ q + 1 is multiplied to λ 1 so the group size can be adjusted in bi-level selection. With only the group level penalty, it is no longer necessary because all the groups are of the same size.
We set the first (q + 1) entries on the main diagonal of H(β g n ) to zero, since the intercept and main environmental effects are not subject to selection. The first and second derivative functions of group MCP can be approximated by nH(β g n ) and nH(β g n )β g n , respectively. The iterative Newton algorithm proceeds andβ g n can be updated across iterations. The stopping criterion is the same as we have adopted for the bi-level QIF. That is, the difference between two penalized estimates from successive iterations is small in terms of L1 norm. In general, convergence can be achieved in a few iterations.

A.2 Penalized QIF
Penalized QIF conducts selection of main and interaction effects in the G×E studies without accounting for the group structure. The penalized score function can be expressed as: where regression coefficient η nhk denotes the kth element of η nk . The Newton-Raphson update ofβ n is derived as: where P (β g n ) and V (β g n ) are defined as follows: The main diagonal of H(β g n ) includes the first order derivative of MCP: where λ 1 and γ are the tuning and regularization parameters, respectively. The first q + 1 zeros on the main diagonal of H(β g n ) leave no shrinkage on the coefficients corresponding to the intercept and E factors. The Newton algorithm proceeds in a similar manner as fitting for the bi-level and group level models.

B Derivations of Alternative Methods based on GEE
In addition to penalized QIF methods, package springer also provides the individual-, groupand bi-level feature selection based on GEE. We briefly introduce these counterparts of the aforementioned penalized QIF methods below.

B.1 Penalized Bi-level GEẼ
where the first order derivative of MCP, ρ ′ (·), has been placed on both the empirical group norm ||η nk || Σ k , and η nkh that represents the individual level effects. Here,Q(β n ) is the score equation of GEE, defined asQ(β ). We refer the readers to Section 2.2 where the rationale of developing the above sparse group penalty in longitudinal G×E interaction studies has been described in details.
In general, the Newton-Raphson algorithm can be applied to the optimization problem outlined above. The iterative update ofβ n at the (g+1)th iteration can be derived as: whereṼ (β g n ) denotes the first order derivative of −Q(β n ). Besides, the diagonal matrix H(β g n ) corresponds to the effects from both the individual-and group-level penalty functions, that is,H whereH(β g n ) is the sum of two diagonal matrices including the derivatives of MCP and group MCP on the main diagonal, respectively. The first (1 + q) entries on the main diagonal are zeros, indicating the intercept and environmental main effects are not subjected to penalized selection.
The iterative update ofβ n terminates when the L 1 difference betweenβ g+1 n andβ g n is smaller than a predefined threshold. In R package springer, we have adopted the same stopping criterion for all the implemented methods in both the GEE and QIF framework.

B.2 Penalized Group GEE
The penalized group QIF implemented in package springer has the score equation defined below:Ũ where ρ ′ (·) is the first order derivative of the minimax concave penalty, and we define the score equation of GEE asQ(β . The coefficient vector β n corresponds to all the main and interaction effects, and for the kth genetic factor, η nk represents the (q+1) main genetic and G×E interaction effects. The group level shrinkage has been imposed on ||η nk || Σ k , the empirical norm of η k .
At the (g + 1)th iteration of the the Newton-Raphson algorithm, we solve forβ g+1 n usinĝ β g n from the previous iteration as follows: , withṼ (β g n ) being the first order derivative of −Q(β g n ). The diagonal matrixH(β g n ) is defined as:H where the first (q + 1) entries on the main diagonal are zeros corresponding the main effects not subject to selection, and the rest nonzero entries are the first order derivatives of group MCP, corresponding to the kth G factor, respectively. Therefore, nH(β g n )β g n denotes the second order derivative function of group MCP in the aforementioned Newton-Raphson update usingβ g n . InH(β g n ), λ 1 and γ are the tuning and regularization parameters determining the magnitude of group level shrinkage and the concavity of group MCP, correspondingly. Since the group size is a constant, it is not necessary to adjust the group size by multiplying a factor of √ q + 1 to the group norm as in the bi-level selection case.