# Broken Adaptive Ridge Method for Variable Selection in Generalized Partly
Linear Models with Application to the Coronary Artery Disease Data

Christian Chan1 Xiaotian Dai2 Thierry Chekouo3 Quan Long4 Xuewen Lu1
1 Department of Mathematics and Statistics, University of Calgary
2 Department of Mathematics, Illinois State University
3 Division of Biostatistics, School of Public Health, University of Minnesota
4 Department of Biochemistry and Molecular Biology, University of Calgary

###### Abstract

Motivated by the CATHGEN data, we develop a new statistical learning method
for simultaneous variable selection and parameter estimation under the context
of generalized partly linear models for data with high-dimensional covariates.
The method is referred to as the broken adaptive ridge (BAR) estimator, which
is an approximation of the $L_{0}$-penalized regression by iteratively
performing reweighted squared $L_{2}$-penalized regression. The generalized
partly linear model extends the generalized linear model by including a non-
parametric component to construct a flexible model for modeling various types
of covariate effects. We employ the Bernstein polynomials as the sieve space
to approximate the non-parametric functions so that our method can be
implemented easily using the existing R packages. Extensive simulation studies
suggest that the proposed method performs better than other commonly used
penalty-based variable selection methods. We apply the method to the CATHGEN
data with a binary response from a coronary artery disease study, which
motivated our research, and obtained new findings in both high-dimensional
genetic and low-dimensional non-genetic covariates.

Keywords: Bernstein polynomials; BAR regression method; Generalized partly
linear models; High-dimensional data; Logistic partly linear model

##  1 Introduction

In the era of high technology and supercomputing power, the combination of
lower financial costs and greater accessibility to DNA sequencing technology
has contributed to the rapid rise in omics research (Mardis, 2011). The high-
throughput DNA sequencing equipment produces high-dimensional data, which
motivates researchers to identify the genetic variations in the genome that
are relevant to the phenotype. In our research, we develop new statistical
learning methods to further decode acquired data and help scientists to find
relevant genetic covariates. Variable selection is one of these statistical
learning methods and is an important task when building a statistical model.
Given a large number of explanatory variables in a particular study, we want
to select the variables that are relevant to the response variable. One way to
do this is by using the best subset selection, which is based on the
$L_{0}$-regularization. The best subset selection method directly penalizes
the cardinality of a model subject to an information criterion, like the AIC
(Akaike, 1974) and BIC (Schwarz, 1978). There are several disadvantages to the
$L_{0}$-regularization method, the most important of them is the computational
complexity at scales of $2^{p}$, where $p$ is the dimension of the covariates,
thus making it computationally expensive for even a moderately large number of
covariates. Additionally, Breiman (1996) also showed that the
$L_{0}$-regularization is unstable in terms of variable selection. Penalty-
based variable selection methods were introduced to solve the computation
inefficiency of $L_{0}$-regularization. The significance of penalty-based
variable selection is the reformulation of the sparse estimation problem into
a continuous and nonconvex or convex optimization problem with a fewer number
of candidate models. Such methods include the Least Absolute Shrinkage and
Selection Operator (LASSO) (Tibshirani, 1996), Smoothly Clipped Absolute
Deviation (SCAD) (Fan and Li, 2001), the Elastic-Net penalty (Zou and Hastie,
2005), the Adaptive LASSO (Zou, 2006) and the Minimax Concave Penalty (MCP)
(Zhang, 2010).

Recently, the Broken Adaptive Ridge (BAR) regression method has been
introduced as an approximation to the $L_{0}$-regularization for variable
selection. The BAR regression can be summarized as an iteratively reweighted
squared $L_{2}$-penalized regression, where the estimators of the BAR method
are taken at the limit of the algorithm. Liu and Li (2016) first considered
the implementation of the BAR method under generalized linear models (GLM).
Since then, many papers have investigated the BAR method for different models
and data types, including the Cox PH model with large-scale right-censored
survival data (Kawaguchi et al., 2020), the linear model with uncensored data
(Dai et al., 2018), the additive hazards model with recurrent event data (Zhao
et al., 2018), the Cox PH model with interval-censored data (Zhao et al.,
2019), the partly linear Cox PH model with right-censored data (Wu et al.,
2020), and the accelerated failure time model with right-censored data (Sun et
al., 2022), among others. Most recently, Mahmoudi and Lu (2022) incorporated
the BAR method for semi-competing risks data under the illness-death model.
Previous work (Dai et al., 2018; Zhao et al., 2018; Kawaguchi et al., 2020)
have also proved that the BAR method possesses two desired large-sample
properties: consistency for variable selection and asymptotic normality, which
are called the oracle properties in the literature.

Motivated by the CATHGEN data detailed below, our goal in this research is to
extend the BAR method to select important variables in generalized partly
linear models (GPLMs) with a large number of genetic covariates, in the
presence of some low-dimensional non-genetic covariates. Particularly, we
apply the proposed method to select important single nucleotide polymorphisms
(SNPs) in a logistic partly linear model in the presence of both categorical
and continuous low dimensional non-genetic covariates, which belongs to a
family of generalized partly linear models. In the CATHerization GENetics
(CATHGEN) study, the primary objective was to assess the association of
multiple genetic markers with cardiovascular disease phenotypes. The study,
conducted by Duke University Medical Centre, collected peripheral blood
samples from consenting patients between 2001 until 2012. The follow-up period
of the recruited patients was between 2004 and 2014. Aside from the high-
dimensional genetic data, low-dimensional baseline clinical and demographical
variables were also measured when patients were first recruited to the study.
The data can be downloaded from U.S. National Institute of Health dbGaP data
accession number phs000704.v1.p1. We will use the proposed method to analyze
the data to identify important SNPs and the associated genes relevant to
coronary artery disease (CAD).

The contributions of our work can be summarized from three main aspects.
First, we develop a new statistical learning method for simultaneous variable
selection and estimation under the context of generalized partly linear models
using the BAR method. GPLMs extend GLMs by adding a non-parametric component
to it to allow for flexible modeling of linear and non-linear covariate
effects. Our method extends the work by Li et al. (2021) which incorporates
the BAR method under GLM with sparse high-dimensional and massive sample size
data for variable selection. Since the low-dimensional covariates in our work
contain both categorical and continuous variables, our method treats them
separately, only the continuous variables are considered to possess potential
non-linear effects. Second, we focus specifically on the logistic partly
linear regression model, as motivated by the presence of the binary response
variable (CAD vs. no CAD) and various types of covariates in the CATHGEN data.
We apply our proposed method to the CATHGEN data to identify relevant genetic
markers (i.e., SNPs) in high dimensions that contribute to developing CAD. We
are also interested in the estimation of the low-dimensional relevant non-
genetic covariate effects, which can handle both linear and non-linear
covariate effects. Third, our method can be easily implemented using existing
R packages developed for GLM, since we are able to use Bernstein polynomials
to construct a linear sieve space for estimating the non-parametric functions
of the low-dimensional covariates so that the resultant model form mimics a
GLM and the existing GLM R packages can be used for estimation and variable
selection. We make our code available at https://github.com/chrischan94/GPLM-
BAR.

The rest of this article is organized as follows. In section 2, we give a
comprehensive introduction to GPLMs, and a detailed explanation of our
proposed method and its algorithm. In section 3, we present the results of our
extensive simulation studies, where we compare our method to a few common
variable selection methods. In section 4, we present the results of the real
data analysis of the CATHGEN study and the biological interpretation of the
results. Finally, in section 5, we conclude our findings in this article and
discuss possible future directions for research. The description of the
involved algorithm and more simulation results are relegated to the Appendix.

##  2 Models and Methods

###  2.1 Generalized Partly Linear Models

Working under the framework of GLM, consider a random sample
$\textbf{v}_{i}=(\textbf{x}_{i}^{\top},y_{i})^{\top},i=1,\ldots,n$, where
$\textbf{y}=\\{y_{1},\ldots,y_{n}\\}^{\top}$ makes a $n\times 1$ response
vector and matrix
$\textbf{X}=\\{\textbf{x}_{1},\ldots,\textbf{x}_{n}\\}^{\top}$ makes a
$n\times p$ design. The observations
$\textbf{v}_{i}=(\textbf{x}_{i}^{\top},y_{i})^{\top},i=1,\ldots,n$ are
mutually independent. The distribution of $y_{i}$ conditional on
$\textbf{x}_{i}$ is from the exponential family with the following density,

$f_{\textbf{y}}(y_{i};\theta_{i},\phi)=\exp\left\\{\frac{y_{i}\theta_{i}-b(\theta_{i})}{a(\phi)}+c(y_{i},\phi)\right\\},$
(2.1)

where $a(\cdot)$, $b(\cdot)$ and $c(\cdot,\cdot)$ are known specific
functions, $b(\cdot)$ is assumed to be twice differentiable, $\theta_{i}$ is
the canonical parameter and $\phi$ denotes the dispersion parameter. Model
(2.1) indicates that
$\mathbb{E}(y_{i}|\textbf{x}_{i})=\mu_{i}=b^{\prime}(\theta_{i})$ and
$\mathbb{V}(y_{i}|\textbf{x}_{i})=b^{\prime\prime}(\theta_{i})a(\phi)$.
Through a link function $g(\mu_{i})=\boldsymbol{\beta}^{\top}\textbf{x}_{i}$,
the canonical parameter $\theta_{i}$ is connected to $\textbf{x}_{i}$ by a
linear combination of the coefficient parameter vector
$\boldsymbol{\beta}=\\{\beta_{1},\ldots,\beta_{p}\\}^{\top}$. When
$g(\mu_{i})=\theta_{i}$, it is called the canonical link function. Commonly
used canonical link functions in GLMs include the identity function for linear
regression, logit link function for logistic regression and the log function
for Poisson regression. Given our observed data, the likelihood function of
$\boldsymbol{\beta}$ for GLMs is

$\mathcal{L}_{n}(\boldsymbol{\beta};\textbf{v}_{i})=\prod^{n}_{i=1}f_{\textbf{y}}(y_{i};\theta_{i},\phi)=\prod^{n}_{i=1}\exp\left\\{\frac{y_{i}\theta_{i}-b(\theta_{i})}{a(\phi)}+c(y_{i},\phi)\right\\},$

and the log-likelihood is

$\ell_{n}(\boldsymbol{\beta})=\log\mathcal{L}_{n}(\boldsymbol{\beta};\textbf{v}_{i})=\sum^{n}_{i=1}\log
f_{\textbf{y}}(y_{i};\theta_{i},\phi).$

GLM assumes a linear relationship between the independent variables
$\\{\textbf{x}_{1},\ldots,\textbf{x}_{n}\\}^{\top}$ and the canonical link
function. If this assumption is violated for a subset of variables that may
have a non-linear relationship with the response, an alternative model form is
desirable. Motivated by the CATHGEN data that contains both genetic and non-
genetic variables, the covariates can be broadly grouped into three distinct
sets: a $n\times p$ design matrix
$\textbf{X}=\\{\textbf{x}_{1},\ldots,\textbf{x}_{n}\\}^{\top}$, a $n\times
q_{w}$ design matrix
$\textbf{W}=\\{\textbf{w}_{1},\ldots,\textbf{w}_{n}\\}^{\top}$, and a $n\times
q_{z}$ design matrix
$\textbf{Z}=\\{\textbf{z}_{1},\ldots,\textbf{z}_{n}\\}^{\top}$. The design
matrix X contains the high-dimensional genetic covariates, W contains the low-
dimensional and non-genetic categorical covariates, and Z contains the low-
dimensional and non-genetic continuous covariates. Then we define the
generalized partly linear model as follows, which extends GLM by adding a non-
parametric component in the linear predictor,

$g\left\\{\mathbb{E}(y_{i}|\textbf{x}_{i},\textbf{w}_{i},\textbf{z}_{i})\right\\}=\boldsymbol{\beta}^{\top}\textbf{x}_{i}+\boldsymbol{\alpha}^{\top}\textbf{w}_{i}+\Psi(\textbf{z}_{i}),$

where $\boldsymbol{\beta}=\\{\beta_{1},\ldots,\beta_{p}\\}^{\top}$,
$\boldsymbol{\alpha}=\\{\alpha_{1},\ldots,\alpha_{q_{w}}\\}^{\top}$,
$\Psi(\textbf{z}_{i})=\sum^{q_{z}}_{j=1}\psi_{j}(z_{ij})$, $\psi_{j}(\cdot)$’s
are unknown smooth functions, which model possible non-linear effects as shown
in the analysis of the CATHGEN data. A special case of GPLM is the logistic
partly linear model. Let
$\pi(y_{i}|\textbf{x}_{i},\textbf{w}_{i},\textbf{z}_{i})=P(y_{i}=1|\textbf{x}_{i},\textbf{w}_{i},\textbf{z}_{i})$,
then the logistic partly linear model has the model equation given by

$\log{\frac{\pi(y_{i}|\textbf{x}_{i},\textbf{w}_{i},\textbf{z}_{i})}{1-\pi(y_{i}|\textbf{x}_{i},\textbf{w}_{i},\textbf{z}_{i})}}=\boldsymbol{\beta}^{\top}\textbf{x}_{i}+\boldsymbol{\alpha}^{\top}\textbf{w}_{i}+\Psi(\textbf{z}_{i}).$

For the observations
$\\{\textbf{u}_{i},i=1,\ldots,n\\}=\\{\\{y_{i},\textbf{x}_{i},\textbf{w}_{i},\textbf{z}_{i}\\},i=1,\ldots,n\\}$,
the likelihood function of the logistic partly linear model can be constructed
as

$\mathcal{L}_{n}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\psi})=\prod^{n}_{i=1}\pi(y_{i}|\textbf{x}_{i},\textbf{w}_{i},\textbf{z}_{i})^{y_{i}}\left\\{1-\pi(y_{i}|\textbf{x}_{i},\textbf{w}_{i},\textbf{z}_{i})\right\\}^{1-y_{i}},$

where
$\boldsymbol{\psi}=\\{\psi_{1}(\cdot),\ldots,\psi_{q_{z}}(\cdot)\\}^{\top}$.
From the likelihood function, the log-likelihood can be easily derived as
follows

$\begin{split}\ell_{n}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\psi})&=\sum^{n}_{i=1}\big{[}y_{i}\left(\boldsymbol{\beta}^{\top}\textbf{x}_{i}+\boldsymbol{\alpha}^{\top}\textbf{w}_{i}+\Psi(\textbf{z}_{i})\right)\\\
&-\log\left(1+\exp\left\\{\boldsymbol{\beta}^{\top}\textbf{x}_{i}+\boldsymbol{\alpha}^{\top}\textbf{w}_{i}+\Psi(\textbf{z}_{i})\right\\}\right)\big{]}.\end{split}$
(2.2)

Direct estimation of
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\psi})$
in (2.2) will not be possible because of the presence of the unknown functions
$\Psi(\textbf{z}_{i})$ which are infinitely dimensional. Hence, approximating
the non-parametric part is needed. As the unknown functions
$\Psi(\textbf{z}_{i})$ are infinitely dimensional, we propose to construct a
sieve space to linearize them. To apply the sieve method, we employ the
Bernstein polynomials to approximate $\Psi(\textbf{z}_{i})$, then, the
Bernstein polynomials approximation reduces the infinitely dimensional space
to a finitely dimensional space. Let $\boldsymbol{\Theta}$ denote the
parameter space of $\boldsymbol{\vartheta}$ where

$\boldsymbol{\Theta}=\left\\{\boldsymbol{\vartheta}=\left(\boldsymbol{\alpha},\boldsymbol{\beta},\psi_{1},\ldots,\psi_{q_{z}}\right)\in\mathcal{A}\otimes\mathcal{M}_{1}\otimes\cdots\otimes\mathcal{M}_{q_{z}}\right\\},$

where

$\mathcal{A}=\left\\{\left(\boldsymbol{\alpha},\boldsymbol{\beta}\right)\in\mathbb{R}^{q_{w}}\times\mathbb{R}^{p},\left\lVert\boldsymbol{\alpha}\right\rVert+\left\lVert\boldsymbol{\beta}\right\rVert\leq
M\right\\},$

$M$ is a positive constant, and $\mathcal{M}_{j}$ denotes the collection of
all bounded and continuous functions over the range of the observed
$\textbf{z}_{j}$ for $j=1,\ldots,q_{z}$. Subsequently, the sieve space is
defined as

$\boldsymbol{\Theta}_{n}=\left\\{\boldsymbol{\vartheta}_{n}=\left(\boldsymbol{\alpha},\boldsymbol{\beta},\psi_{1n},\ldots,\psi_{q_{z}n}\right)\in\mathcal{A}\otimes\mathcal{M}_{1n}\otimes\cdots\otimes\mathcal{M}_{q_{z}n}\right\\},$

where

$\mathcal{M}_{jn}=\left\\{\psi_{jn}(z_{ij})=\sum^{m_{j}}_{k=0}\gamma_{jk}B_{jk}(z_{ij},m_{j},c_{j},u_{j}):\sum_{0\leq
k\leq m_{j}}|\gamma_{jk}|\leq M_{jn}\right\\},$ (2.3)

for $j=1,\ldots,q_{z}$ and $z_{ij}\in[c_{j},u_{j}]$, $c_{j}<u_{j}$. The
Bernstein basis polynomial of $m_{j}$ degree, denoted by
$B_{jk}(z_{ij},m_{j},c_{j},u_{j})$ in (2.3), has the equation

$B_{jk}(z_{ij},m_{j},c_{j},u_{j})=\begin{pmatrix}m_{j}\\\
k\end{pmatrix}\Bigg{(}\frac{z_{ij}-c_{j}}{u_{j}-c_{j}}\Bigg{)}^{k}\Bigg{(}1-\frac{z_{ij}-c_{j}}{u_{j}-c_{j}}\Bigg{)}^{m_{j}-k},\quad
k=0,\ldots,m_{j}.$

Therefore, the sieve log-likelihood of the logistic partly linear model using
the Bernstein polynomial to approximate $\Psi(\textbf{z}_{i})$ is

$\begin{split}\ell_{n}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma})=&\sum^{n}_{i=1}\bigg{\\{}y_{i}\Big{(}\boldsymbol{\beta}^{\top}\textbf{x}_{i}+\boldsymbol{\alpha}^{\top}\textbf{w}_{i}+\sum^{q_{z}}_{j=1}\sum^{m_{j}}_{k=0}\gamma_{jk}B_{jk}(z_{ij},m_{j},c_{j},u_{j})\Big{)}\\\
&-\log\Big{(}1+\exp\big{\\{}\boldsymbol{\beta}^{\top}\textbf{x}_{i}+\boldsymbol{\alpha}^{\top}\textbf{w}_{i}+\sum^{q_{z}}_{j=1}\sum^{m_{j}}_{k=0}\gamma_{jk}B_{jk}(z_{ij},m_{j},c_{j},u_{j})\big{\\}}\Big{)}\bigg{\\}},\end{split}$
(2.4)

where
$\boldsymbol{\gamma}=\\{\gamma_{10},\ldots,\gamma_{1m_{1}},\ldots,\gamma_{q_{z}0},\ldots,\gamma_{q_{z}m_{q_{z}}}\\}^{\top}$.

###  2.2 Simultaneous Estimation and Variable Selection using GPLM-BAR method

To conduct simultaneous estimation and variable selection in GPLMs, we propose
the GPLM-BAR method, which is an iterative method. Following the BAR method by
Li et al. (2021) for GLM, starting from an initial value vector computed from
the following the ridge regression,

$\left(\widehat{\boldsymbol{\alpha}}^{(0)},\widehat{\boldsymbol{\beta}}^{(0)},\widehat{\boldsymbol{\gamma}}^{(0)}\right)=\operatorname*{arg\,min}_{\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma}}\left\\{-2\ell_{n}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma})+\xi_{n}\sum^{p}_{j=1}\beta^{2}_{j}\right\\}.$
(2.5)

For $s\geq 1$, the estimator is iteratively updated by a reweighted squared
$L_{2}$-penalized regression

$\left(\widehat{\boldsymbol{\alpha}}^{(s)},\widehat{\boldsymbol{\beta}}^{(s)},\widehat{\boldsymbol{\gamma}}^{(s)}\right)=\operatorname*{arg\,min}_{\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma}}\left\\{-2\ell_{n}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma})+\lambda_{n}\sum^{p}_{j=1}\frac{\beta^{2}_{j}}{(\widehat{\beta}^{(s-1)}_{j})^{2}}\right\\},$
(2.6)

where $\xi_{n}$ and $\lambda_{n}$ are non-negative tuning parameters. The
updated step in (2.6) is continued until a pre-specified convergence criterion
is reached, where the estimators are taken at the limit as

$\widehat{\boldsymbol{\beta}}=\lim_{s\rightarrow\infty}\widehat{\boldsymbol{\beta}}^{(s)},\quad\widehat{\boldsymbol{\alpha}}=\lim_{s\rightarrow\infty}\widehat{\boldsymbol{\alpha}}^{(s)},\quad\widehat{\boldsymbol{\gamma}}=\lim_{s\rightarrow\infty}\widehat{\boldsymbol{\gamma}}^{(s)}.$

The implementation of the proposed method indicates the variable selection is
done only on high-dimensional covariates $\textbf{x}_{i}$, $i=1,\ldots,n$,
since the penalty is imposed on $\boldsymbol{\beta}$ only. In addition, the
proposed variable selection method can be applied in a similar fashion to the
Poisson partly linear regression and the partly linear model for counts and
continuous responses, respectively, since they are in the family of GPLMs.

####  2.2.1 A note on choosing the tuning parameters

Choosing the optimal values of tuning parameters is crucial for penalty-based
variable selection methods, as it greatly affects the variable selection
accuracy. In the absence of an external validation set, common methods to find
the optimal values of the tuning parameters include the $k$-fold cross-
validation (CV) method. The optimal tuning parameter value is the one that
minimizes a criterion. Typically, this is the mean squared error for
continuous outcome, or deviance for binary outcome. However, doing this only
adds to the computational complexity, and it is not ideal for larger datasets.
In the GPLM-BAR algorithm, we have two tuning parameters: $\xi_{n}$ and
$\lambda_{n}$. Unless the value of $\xi_{n}$ chosen is large, it is
empirically shown that the value chosen is inconsequential on the estimation
of $\boldsymbol{\beta}$, as seen in Figure 1. Hence, $\xi_{n}$ is set to a
relatively small value.

Figure 1: Path plots for the logistic partly linear regression with varying
$\xi_{n}$ for a random sample of size $n=200$ and $p=10$, using both the BIC
(left panel) and AIC (right panel) penalties.

For $\lambda_{n}$ in the Cox-BAR regression, it has been argued by Kawaguchi
et al. (2020) that it can be fixed. One example is fixing
$\lambda_{n}=\log(n)$, which corresponds the BIC penalty. Another example is
to fix $\lambda_{n}=2$, which corresponds to the AIC penalty. In our method,
both the AIC and BIC penalties are considered.

####  2.2.2 Computational aspects for GPLM-BAR

Except under the linear model, numerical approximation methods such as the
Newton-Raphson algorithm are integrated into the implementation of the BAR
penalty for simultaneous variable selection and estimation. When both the
number of covariates and sample size are small, calculating the partial
gradient vector and Hessian matrix at each iteration of the BAR algorithm is
computationally feasible. However, when both the number of covariates and
sample size becomes moderately big, numerical approximation becomes not
scalable because of the high computational costs and the numerical
instability. Alternative optimization techniques for parameter estimation
under large-scale regularization and regression problems (Zhang and Oles,
2001; Azoury and Warmuth, 2001) have been developed. The algorithm by Zhang
and Oles (2001) called column relaxation of logistic loss (CLG) can be
classified as a cyclic coordinate descent algorithm.

The R package `BrokenAdaptiveRidge` (Kawaguchi et al., 2020) was created to
implement BAR regression for GLM and the Cox model, which are linear models.
Since we have reparameterized our GPLM into a form of GLM in (2.4), we are
able to directly use the package to conduct variable selection and estimation
under the context of GPLM. This package uses the R package `Cyclops` (Suchard
et al., 2013) for efficient implementation of the iterative method as
described in Kawaguchi et al. (2020). To do this, first we create a data frame
`B` for Bernstein polynomial basis functions based on the low-dimensional
covariates `Z`, each column in `B(Z)` represents one basis function. Let `X`
represent the design matrix for high-dimensional genetic covariates, and `W`
for low-dimensional non-genetic covariates, `Y` represents the $(0,1)$ binary
response vector. Finally, we make a combined data frame `D`. Then, for
example, the GPLM-BAR estimates selected by AIC can be computed from the
following R code

    D <- data.matrix(cbind(X,W,B(Z))
    penAIC <- createBarPrior(penalty = 2,
                             exclude = c(1, ((ncol(X)+2):(ncol(D)+1)),
                             initialRidgeVariance = 1)
           #penalty=2 indicates AIC, and log(n) indicates BIC;
    cyD <- createCyclopsData(Y ~ D, modelType = "lr")
           #lr indicates logistic regression
    BARfit <- fitCyclopsModel(cyD, prior = penAIC)
           #estimates of all of the coefficients

The computation in the package is done by the cyclic coordinate descent
algorithm. We describe this algorithm for the GPLM-BAR regression in the
Appendix.

##  3 Simulation Studies

In this section, we present the results of a comprehensive simulation study in
three scenarios to demonstrate the effectiveness of our proposed method. The
first and second scenarios assess the performance under strong signals and
weak signals, respectively, under the setting of the logistic partly linear
model. The third scenario uses the selected model in the real data analysis
section for the CATHGEN data as a basis to simulate data, then assesses the
performance of our method, and the fourth scenario shows the performance under
the setting of the Poisson partly linear model.

### Scenario 1: Strong signals in the logistic partly linear model

In this scenario, let $q_{z}=5$ and $q_{w}=4$, and the number of non-zero
elements $q=5$ in the true parameter $p$-vector $\boldsymbol{\beta}_{0}$, for
various values of $p$. We generate the $n\times p$ design matrix X from a
multivariate normal distribution with mean zero and variance-covariance matrix
$\Sigma_{\textbf{X}}$, where the $(i,j)^{th}$ entry of it is $\rho^{|i-j|}$.
We fix $\rho=0.25$. We first consider large effects, i.e., large values of
$\boldsymbol{\beta}$, where the true value of $\boldsymbol{\beta}$ is
$\boldsymbol{\beta}_{0}=\\{1,-1,0,\ldots,0,-1,0.75,0.75\\}^{\top}$. We also
generate a $n\times q_{w}$ design matrix W from independent Bernoulli
distributions, with the same probability of success $\pi=0.5$. And, the true
value of $\boldsymbol{\alpha}$ is
$\boldsymbol{\alpha}_{0}=\\{1,-0.5,-0.5,0.75,-1\\}^{\top}$. Independently from
X and W, we also generate a $n\times q_{z}$ design matrix Z, where we draw
$\textbf{z}_{1}$ from the uniform distribution over (1,5), $\textbf{z}_{2}$
and $\textbf{z}_{3}$ independently from the standard uniform distribution, and
$\textbf{z}_{4}$ from the uniform distribution over $(-3,1)$. By setting the
non-linear functions to be $\psi_{1}(z_{i1})=0.1(z_{i1}-3)^{2}$,
$\psi_{2}(z_{i2})=0.2(\cos(2\pi z_{i2})+1)$, $\psi_{3}(z_{i3})=0.2\sin(2\pi
z_{i3})$, and $\psi_{4}(z_{i4})=0.2(z_{i4}+1)^{3}$, respectively, we generate
$y_{i}$ from the Bernoulli distribution with probability $\pi_{i}$, where
$\pi_{i}=1/(1+\exp\\{-\boldsymbol{\beta}_{0}^{\top}\textbf{x}_{i}-\boldsymbol{\alpha}_{0}^{\top}\textbf{w}_{i}-\sum^{4}_{j=1}\psi_{j}(z_{ij})\\})$.
The chosen non-linear functions have two common properties: 1) they are
symmetric at the midpoint of the interval of their domains, 2) The values of
the functions are zero at the midpoint for the purpose of identifiability. We
consider two different sample sizes $n=600$ and $800$, and two different
numbers of high-dimensional covariates $p=300$ and $450$. Each combination is
replicated 200 times. The number of basis functions for all non-linear
functions is set at $m_{j}+1=3+1=4$, since more than four basis functions only
add to the computational complexity while only marginally improving the
approximation of $\psi_{j}(\cdot)$, conversely having fewer than four basis
functions will not approximate $\psi_{j}(\cdot)$ well.

In the simulation studies, we compare our method against the methods of LASSO
and Adaptive LASSO. We use the R package `splines2` (Wang and Yan, 2021) to
generate the Bernstein polynomials. The LASSO and Adaptive LASSO methods are
implemented using the R package `glmnet` (Friedman et al., 2010; Simon et al.,
2011). To evaluate the estimation accuracy, we compute the median mean squared
error (MMSE), where the mean squared error has the equation
$(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{0})^{\top}\Sigma_{\textbf{X}}(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{0})$.
For the GPLM-BAR method, we fix $\lambda_{n}$ to two values, $\lambda_{n}=2$
and $\lambda_{n}=\log(n)$, which corresponds to the AIC and BIC penalties
respectively. Since the value of $\xi_{n}$ was shown to have an
inconsequential effect on estimation, we set $\xi_{n}=1$. For the other
methods, we use $10$-fold CV method to select the optimal value. To evaluate
the selection accuracy, we compute the average number of true positives (TP),
average number of false positives (FP), total misclassification rate (MC),
frequency of true model (TM) selected, and the average estimated size of the
model (MS), where $\text{MS}=\text{TP}+\text{FP}$.

Table 1: Variable selection and estimation results over 200 replications for Scenario 1. Standard deviations of the MMSE are in parentheses. Method | MMSE | TP | FP | MS | MC | TM
---|---|---|---|---|---|---
$n=600,p=300$
BAR(AIC) | 0.17(0.10) | 5 | 1.31 | 6.31 | 1.31 | 24$\%$
BAR(BIC) | 0.28(0.27) | 4.74 | 0 | 4.74 | 0.26 | 74$\%$
LASSO | 1.03(0.22) | 5 | 3.48 | 8.48 | 3.48 | 16$\%$
ALASSO | 0.46(0.24) | 5 | 2.68 | 7.68 | 2.68 | 51$\%$
Oracle | 0.08(0.08) | 5 | 0 | 5 | 0 | 100$\%$
$n=800,p=300$
BAR(AIC) | 0.11(0.07) | 5 | 1.21 | 6.21 | 1.21 | 28$\%$
BAR(BIC) | 0.17(0.12) | 4.97 | 0 | 4.97 | 0.03 | 98$\%$
LASSO | 0.67(0.17) | 5 | 8.10 | 13.10 | 8.10 | 1$\%$
ALASSO | 0.21(0.11) | 5 | 4.10 | 9.10 | 4.10 | 32$\%$
Oracle | 0.05(0.06) | 5 | 0 | 5 | 0 | 100$\%$
$n=600,p=450$
BAR(AIC) | 0.20(0.14) | 5 | 1.83 | 6.83 | 1.83 | 17$\%$
BAR(BIC) | 0.32(0.31) | 4.70 | 0 | 4.70 | 0.30 | 74$\%$
LASSO | 0.89(0.22) | 5 | 10.47 | 15.47 | 10.47 | 2$\%$
ALASSO | 0.33(0.16) | 5 | 10.35 | 15.35 | 10.35 | 11$\%$
Oracle | 0.08(0.07) | 5 | 0 | 5 | 0 | 100$\%$
$n=800,p=450$
BAR(AIC) | 0.14(0.09) | 5 | 1.77 | 6.77 | 1.77 | 16$\%$
BAR(BIC) | 0.16(0.13) | 4.98 | 0 | 4.98 | 0.02 | 98$\%$
LASSO | 0.70(0.18) | 5 | 9.80 | 14.80 | 9.80 | 1$\%$
ALASSO | 0.22(0.15) | 5 | 7.62 | 12.62 | 7.62 | 19$\%$
Oracle | 0.06(0.06) | 5 | 0 | 5 | 0 | 100$\%$

From Table 1, one can observe that the GPLM-BAR method performs better than
LASSO and Adaptive LASSO by most measures of accuracy. Although the average
number of TP is not the best particularly with the BIC penalty, the average
number of FP is far lower as compared to the other methods, resulting in the
total misclassification rate being the lowest. The GPLM-BAR method also
produces the sparsest model. It is interesting to observe a trade-off between
the AIC and BIC penalties, where estimation accuracy is better with the AIC
penalty, contrasting the better variable selection results with the BIC
penalty. This is explained by the larger tuning parameter in the BIC penalty,
which shrinks the relatively smaller signals in $\boldsymbol{\beta}$ to zero,
thus causing a larger estimation bias. We also report the estimation results
of $\boldsymbol{\alpha}$ in Table 12 in the Appendix, where the estimates of
our method is the best.

We also are interested in the estimation of non-linear covariate effects
$\psi_{j}(\cdot)$ using the GPLM-BAR method. The estimated curves are shown in
Figure 2, which compares the averaged estimates of each of the four non-linear
functions to the true function. Two observations are made. First, the
Bernstein polynomial using three basis functions to approximate each
$\psi_{j}(\cdot)$ is satisfactory, where the general shape of each function is
captured well. Second, different $\lambda_{n}$ tuning methods give slightly
different estimates of $\psi_{j}(\cdot)$, as the BAR method with the BIC
penalty (blue curve) produces more biases than the AIC penalty (yellow curve).
The GPLM-BAR also performs well for the other three combinations of $n$ and
$p$ (Appendix Figures 10,11 and 12).

Figure 2: Estimated nonlinear covariate effects $\psi_{j}(\cdot),j=1,2,3,4$,
for Scenario 1, $p=300$ and $n=600$.

### Scenario 2: Strong and weak signals in the logistic partly linear model

We also perform another scenario where a few signals of the non-zero entries
in $\boldsymbol{\beta}$ are weaker. Here, we fix $q=5$, and the true value of
$\boldsymbol{\beta}$ is
$\boldsymbol{\beta}_{0}=\\{1,-0.5,0,\ldots,0,-1,0.4,0.75\\}^{\top}$ in this
case. The true values of $\boldsymbol{\alpha}$ and the non-linear functions
are the same as in Scenario 1.

Table 2: Variable selection and estimation results over 200 replications for Scenario 2. Standard deviations of the MMSE are in parentheses. Method | MMSE | TP | FP | MS | MC | TM
---|---|---|---|---|---|---
$n=600,p=300$
BAR(AIC) | 0.19(0.13) | 4.63 | 1.28 | 5.91 | 1.65 | 21$\%$
BAR(BIC) | 0.50(0.18) | 3.14 | 0 | 3.14 | 1.86 | 0$\%$
LASSO | 0.67(0.15) | 4.66 | 6.51 | 11.17 | 6.85 | 3$\%$
ALASSO | 0.33(0.18) | 4.59 | 6.52 | 11.11 | 6.93 | 8$\%$
Oracle | 0.07(0.05) | 5 | 0 | 5 | 0 | 100$\%$
$n=800,p=300$
BAR(AIC) | 0.13(0.09) | 4.82 | 1.32 | 6.14 | 1.50 | 24$\%$
BAR(BIC) | 0.41(0.14) | 3.46 | 0 | 3.46 | 1.54 | 5$\%$
LASSO | 0.52(0.14) | 4.85 | 8.06 | 12.91 | 8.21 | 2$\%$
ALASSO | 0.24(0.11) | 4.81 | 5.93 | 10.74 | 6.12 | 14$\%$
Oracle | 0.05(0.05) | 5 | 0 | 5 | 0 | 100$\%$
$n=600,p=450$
BAR(AIC) | 0.23(0.14) | 4.65 | 1.89 | 6.54 | 2.24 | 8$\%$
BAR(BIC) | 0.50(0.16) | 3.22 | 0 | 3.22 | 1.78 | 1$\%$
LASSO | 0.68(0.18) | 4.65 | 8.66 | 13.31 | 9.01 | 2$\%$
ALASSO | 0.34(0.23) | 4.54 | 10.20 | 14.74 | 10.66 | 2$\%$
Oracle | 0.07(0.07) | 5 | 0 | 5 | 0 | 100$\%$
$n=800,p=450$
BAR(AIC) | 0.15(0.09) | 4.87 | 1.90 | 6.77 | 2.03 | 10$\%$
BAR(BIC) | 0.41(0.13) | 3.48 | 0 | 3.48 | 1.52 | 5$\%$
LASSO | 0.52(0.14) | 4.87 | 9.15 | 14.02 | 9.28 | 1$\%$
ALASSO | 0.23(0.12) | 4.79 | 9.11 | 13.90 | 9.32 | 6$\%$
Oracle | 0.05(0.04) | 5 | 0 | 5 | 0 | 100$\%$

One is able to observe the GPLM-BAR method with the AIC penalty outperforms
the LASSO and Adaptive LASSO methods from the results in Table 2. However, in
comparison to the results in Scenario 1, the selection and estimation accuracy
become worse, because the weaker signals in $\boldsymbol{\beta}$ have a
greater tendency to be shrunk to zero. The estimation of the non-linear
covariate effects using the GPLM-BAR method are good (Figures 13,14,15 and
16), and the estimation results of $\boldsymbol{\alpha}$ are satisfactory
(Table 13).

### Scenario 3: CATHGEN-based simulations

By using the results in the real data analysis, we also perform a simulation
study to investigate the performance of our method under the correlation
structure of the CATHGEN data. The results are presented in the Appendix.

### Scenario 4: Poisson partly linear model

We also perform a simulation study for the Poisson partly linear model. The
simulated results are presented in the Appendix.

##  4 Real data analysis

Coronary artery disease (CAD) is a major disease that inflicts death, and is
one of the biggest causes of death globally (Abubakar et al., 2015).
Environmental factors that contribute to CAD are typically age, smoking
status, obesity and lifestyle choices. However, genetic factors play a role in
death due to CAD, especially in younger patients (Marenberg et al., 1994).

We apply our proposed method on the CATHGEN data, which was downloaded from
dbGaP, with accession number phs000704.v1.p1. The study collected peripheral
blood samples from consenting patients who were undergoing cardiac
catheterization at Duke University Medical Center from 2001 to 2011. A total
of 1327 patients were recruited and followed-up between 2004 until 2014. The
binary response variable is the affection status, where the stratification
criteria is defined in Shah et al. (2010). The high-dimensional design matrix
contains 13991 columns of SNPs belonging to 331 genes that have been
associated with CAD using Ingenuity Pathway Analysis (Krämer et al., 2014). In
addition to the SNPs, there are ten clinical and demographical variables in
the data. These variables include age (Mean = 57.0, SD = 11.6), BMI (Mean =
30.8, SD = 7.8), smoking status (671 cases out of 1327), race (897 Caucasian-
Americans, 274 African-Americans and 156 Asian-Americans), hypertension status
(900 cases out of 1327), diabetes status (379 cases out of 1327),
hypercholesterolemia status (745 cases out of 1327), sex (684 males and 643
females), number of diseased vessels and history of myocardial infarction
(HXMI) (277 cases out of 1327). All clinical and demographical variables of
each subject were measured when they are included to the study. We exclude the
number of diseased vessels from further analysis because of conversion issues
when fitting the univariate logistic regression model.

Figure 3: The histograms of the BMI on the original scale (top panel), BMI on the log-transformed scale (middle panel), and the histogram of age on the original scale (bottom panel).  Table 3: Univariate polynomial logistic regression models of age and $\log$(BMI) fitted separately, with their respective higher order terms. Variable | Estimate | Std. Error | z-value | p-value
---|---|---|---|---
Intercept | -3.575 | 1.168 | -3.061 | 2E-03
Age | 0.107 | 4E-02 | 2.645 | 8E-03
$\text{Age}^{2}$ | -8E-04 | 4E-04 | -2.220 | 2.6E-02
Intercept | -38.952 | 8.684 | -4.486 | 7.3E-06
$\log(\text{BMI})$ | 22.561 | 5.035 | 4.480 | 7.45E-06
$\log(\text{BMI})^{2}$ | -3.255 | 0.729 | -4.468 | 7.88E-06

In Figure 3, the distribution of age is symmetrical on the original scale.
However, the distribution of the BMI on the original scale is right skewed.
The natural logarithm transformation of it fixes the skewness. Thus, we decide
to use the BMI on the log scale for further analysis. From Table 3, when
fitting age and it’s second order polynomial term in the logistic regression
model, both terms are found to be statistically significant. Likewise, when
fitting the log-transformed BMI and its second order polynomial term in the
logistic regression model, both terms are also significant. The results in
Table 3 indicate that age and log-transformed BMI have a non-linear effect on
the odds ratio of developing CAD. However, the functional form of the effect
is unknown, and this motivates us to consider a logistic partly linear
regression model.

Before we apply our proposed method, it is clear that the dimension of the
design matrix needs to be reduced. To reduce the dimension, we first remove
SNPs with a minor allele frequency (MAF) of less than 0.1. We then further
reduce the number of SNPs through pre-screening the candidate SNPs, by
performing univariate logistic regression, only selecting SNPs with a
$p$-value less than 0.1. In total, 1242 SNPs with a $p$-value less than 0.1
are retained for further analysis.

To choose the tuning parameters in GPLM-BAR, we decide to use the AIC penalty
for $\lambda_{n}$, because the individual estimated effect sizes of the SNPs
are small as shown in Figure 4. The value of $\xi_{n}$ and number of basis
functions are kept the same as the simulation study. The tuning parameters
values for the LASSO and Adaptive LASSO methods are chosen by 5-fold cross
validation. In Table 4, the estimated effects of the categorical clinical
variables obtained from the GPLM-BAR method have a larger magnitude. The
results in Table 4 indicate a positive risk association for hypertension,
diabetes, hypercholesterolomia, smoking and HXMI, where HXMI is the strongest
clinical indicator on the risk of developing CAD. We use the bootstrap method
with 100 random bootstrap samples to obtain the estimated standard error in
parentheses in Table 4. The GPLM-BAR method identified the fewest number of
SNPs that contribute to CAD. Specifically, the GPLM-BAR identified 19
different SNPs that are associated to 17 unique genes, the LASSO and Adaptive
LASSO methods identified 199 SNPs and 228 SNPs, respectively.

Figure 4: Histogram of the estimated effect sizes of SNPs using univariate logistic regression. Table 4: Estimation results of the categorical clinical variables for the CATHGEN data. Variable | GPLM-BAR | LASSO | ALASSO
---|---|---|---
Hypertension | $0.368_{(0.315)}$ | $0.324_{(0.157)}$ | $0.259_{(0.093)}$
Diabetes | $1.180_{(0.269)}$ | $1.097_{(0.242)}$ | $1.105_{(0.261)}$
Hypercholesterolomia | $1.403_{(0.299)}$ | $1.459_{(0.275)}$ | $1.399_{(0.234)}$
Sex | $0.361_{(0.224)}$ | $0.343_{(0.193)}$ | $0.313_{(0.219)}$
Smoking | $0.860_{(0.335)}$ | $0.787_{(0.301)}$ | $0.797_{(0.221)}$
HXMI | $37.388_{(1.786)}$ | $10.620_{(0.745)}$ | $10.477_{(0.561)}$
Race (African) | $-0.030_{(0.489)}$ | $-0.163_{(0.341)}$ | $0.879_{(0.584)}$
Race (Caucasian) | $0.385_{(0.458)}$ | $0.544_{(0.528)}$ | $0.084_{(0.394)}$
Male and Asian-American are the reference categories for the variables Sex and
Race, respectively. The remaining variables have the non-cases as the
reference categories.

From the genes identified using GPLM-BAR, RBFOX1 is found to be associated
with blood pressure and heart failure through transcriptome profiling (Gao et
al., 2016). CDH13 has been shown to be associated with blood cholesterol and
CAD through a genome-wide association study undertaken in the British
population (Consortium, 2007). F10 is associated with the lowering levels of
coagulation factor X, which is protective against ischemic heart disease
(Paraboschi et al., 2020). GABRG3 has been shown to be associated with density
of dodecanedioic acid, which plays a role in regulating blood sugar level
(Wang et al., 2021). ABCA1 has been shown to be associated to altered
lipoprotein levels which results in a increased risk for CAD (Clee et al.,
2001). IL1B belongs to the wider family of IL1 genes which is associated to
coronary heart disease (Francis et al., 1999; Vohnout et al., 2003; Tsimikas
et al., 2014). Certain subtypes of the APOE gene are identified to lipid
levels and coronary risk (Bennet et al., 2007). We report the complete set of
selected SNPs and genes in Table 14.

In addition to the results in Table 4, one can observe our method using
Bernstein polynomials approximation has showed that the effects of age and BMI
are non-linear, as seen in Figure 5. The plot on the left in Figure 5 shows
the risk of developing CAD increases non-linearly with age, and the plot on
the right shows the risk of developing CAD increases with BMI on the natural
logarithm scale until 3.5. After this cutoff point, it then decreases. The
unusual trend seen for BMI can be partially explained by the lack of data when
$\text{BMI}>35$ or $\log(\text{BMI})>3.5$, as BMI of the majority of patients
recruited to this study falls between 15 and 35 on the raw scale.

Figure 5: Estimated covariate effects of age (left panel) and log(BMI) (right
panel).

##  5 Discussion and Conclusion

In this article, we have proposed a new approach for simultaneous variable
selection and estimation under the context of GPLM, with a focus on the
logistic partly linear regression model. Our proposed approach was motivated
by the CATHGEN study, where the data contains both high-dimensional genetic
covariates and low-dimensional clinical and demographical covariates. We
considered GPLM as it grants us the flexibility to model possible non-linear
covariate effects. We employed the Bernstein polynomials to approximate the
non-parametric component of the model, where it has several advantages over
other approximation methods like splines and piecewise functions. First,
unlike the piecewise functions the Bernstein polynomials are differentiable
and continuous everywhere. This is desirable as the first and second
derivatives are calculated in each iteration of our algorithm. Second, the
Bernstein polynomials possesses computational scalability and optimal shape-
preserving property for all approximating polynomials (Carnicer and Pena,
1993). Third, the Bernstein polynomials do not require specification of the
number of interior knots and their locations, unlike B-splines. From the
results of our comprehensive simulation studies, we observe that our proposed
method outperforms common variable selection methods under a few practical
scenarios. Our method incorporating the BAR penalty produced the lowest total
misclassification rate and the highest frequency of the true model selected,
which is consistent with other empirical studies conducted by authors who also
employed the BAR penalty(Dai et al., 2018; Zhao et al., 2019; Sun et al.,
2022). Our method was also able to accurately estimate the true non-linear
functions. As an application, we applied our proposed method to the CATHGEN
data, where certain SNPs and genes were found to have a relevant contribution
to CAD, which is consistent with other variable selection methods applied to
this data (Li and Chekouo, 2022; Dai et al., 2023).

There are several directions one can take from our research. In the simulation
study, we only examined scenarios where the number of high-dimensional
covariates $p$ are diverging with a rate less than the sample size $n$.
Suppose $p$ diverges with a rate greater than $n$, it would be of interest to
investigate the performance of our method under this scenario. The CATHGEN
data also contains right-censored survival information. Under the context of
survival models, it would be of interest to investigate which relevant genetic
markers affect the survival probability, and the possible homogeneity or
heterogeneity of the two sets of genetic markers selected based on the two
different response outcomes and their biological interpretations.
Additionally, choosing the optimal tuning parameter poses a significant
challenge to researchers. The mixture of weak signals with strong signals
poses a noteworthy problem to researches, and requires more thorough
investigation.

## Acknowledgement

The authors would like to acknowledge New Frontiers in Research Fund to Quan
Long (NFRFE-2018-00748) administered by the Canada Research Coordinating
Committee, the Canada Foundation for Innovation JELF grant (36605) awarded to
Quan Long, and the Discovery Grant to Xuewen Lu (RG/PIN06466-2018)
administered by the National Science and Engineering Research Council of
Canada.

## Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the
research, authorship, and/or publication of this article.

## References

  * Abubakar et al. (2015) II Abubakar, Taavi Tillmann, and Amitava Banerjee.  Global, regional, and national age-sex specific all-cause and cause-specific mortality for 240 causes of death, 1990-2013: a systematic analysis for the global burden of disease study 2013.  _Lancet_ , 385(9963):117–171, 2015.
  * Akaike (1974) H. Akaike.  A new look at the statistical model identification.  _IEEE Transactions on Automatic Control_ , 19(6):716–723, 1974.  doi: 10.1109/TAC.1974.1100705.
  * Azoury and Warmuth (2001) Katy S Azoury and Manfred K Warmuth.  Relative loss bounds for on-line density estimation with the exponential family of distributions.  _Machine Learning_ , 43:211–246, 2001.
  * Bennet et al. (2007) Anna M Bennet, Emanuele Di Angelantonio, Zheng Ye, Frances Wensley, Anette Dahlin, Anders Ahlbom, Bernard Keavney, Rory Collins, Björn Wiman, Ulf de Faire, et al.  Association of apolipoprotein e genotypes with lipid levels and coronary risk.  _Jama_ , 298(11):1300–1311, 2007.
  * Breiman (1996) Leo Breiman.  Heuristics of instability and stabilization in model selection.  _The Annals of Statistics_ , 24(6):2350–2383, 1996.
  * Carnicer and Pena (1993) Jesús M Carnicer and Juan Manuel Pena.  Shape preserving representations and optimality of the Bernstein basis.  _Advances in Computational Mathematics_ , 1:173–196, 1993\.
  * Clee et al. (2001) Susanne M Clee, Aeilko H Zwinderman, James C Engert, Karin Y Zwarts, Henri OF Molhuizen, Kirsten Roomp, J Wouter Jukema, Michel van Wijland, Marjel van Dam, Thomas J Hudson, et al.  Common genetic variation in abca1 is associated with altered lipoprotein levels and a modified risk for coronary artery disease.  _Circulation_ , 103(9):1198–1205, 2001.
  * Consortium (2007) W. T. C. C. Consortium.  Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls.  _Nature_ , 447(7145):661–678, 2007.
  * Dai et al. (2018) Linlin Dai, Kani Chen, Zhihua Sun, Zhenqiu Liu, and Gang Li.  Broken adaptive ridge regression and its asymptotic properties.  _Journal of Multivariate Analysis_ , 168:334–351, 2018\.
  * Dai et al. (2023) Xiaotian Dai, Xuewen Lu, and Thierry Chekouo.  A bayesian genomic selection approach incorporating prior feature ordering and population structures with application to coronary artery disease.  _Statistical Methods in Medical Research_ , page 09622802231181231, 2023.
  * Fan and Li (2001) Jianqing Fan and Runze Li.  Variable selection via nonconcave penalized likelihood and its oracle properties.  _Journal of the American Statistical Association_ , 96(456):1348–1360, 2001.
  * Francis et al. (1999) Sheila E Francis, Nicola J Camp, Rachael M Dewberry, Julian Gunn, Petros Syrris, Nicholas D Carter, Stephen Jeffery, Juan Carlos Kaski, David C Cumberland, Gordon W Duff, et al.  Interleukin-1 receptor antagonist gene polymorphism and coronary artery disease.  _Circulation_ , 99(7):861–866, 1999.
  * Friedman et al. (2010) Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _Journal of Statistical Software_ , 33(1):1–22, 2010.  doi: 10.18637/jss.v033.i01.  URL https://www.jstatsoft.org/v33/i01/.
  * Gao et al. (2016) Chen Gao, Shuxun Ren, Jae-Hyung Lee, Jinsong Qiu, Douglas J Chapski, Christoph D Rau, Yu Zhou, Maha Abdellatif, Astushi Nakano, Thomas M Vondriska, et al.  Rbfox1-mediated RNA splicing regulates cardiac hypertrophy and heart failure.  _The Journal of Clinical Investigation_ , 126(1):195–206, 2016.
  * Gorst-Rasmussen and Scheike (2012) Anders Gorst-Rasmussen and Thomas H Scheike.  Coordinate descent methods for the penalized semiparametric additive hazards model.  _Journal of Statistical Software_ , 47:1–17, 2012.
  * Kawaguchi et al. (2020) Eric S Kawaguchi, Marc A Suchard, Zhenqiu Liu, and Gang Li.  A surrogate $L_{0}$ sparse Cox’s regression with applications to sparse high-dimensional massive sample size time-to-event data.  _Statistics in Medicine_ , 39(6):675–686, 2020\.
  * Krämer et al. (2014) Andreas Krämer, Jeff Green, Jack Pollard Jr, and Stuart Tugendreich.  Causal analysis approaches in ingenuity pathway analysis.  _Bioinformatics_ , 30(4):523–530, 2014.
  * Li et al. (2021) Ning Li, Xiaoling Peng, Eric Kawaguchi, Marc A Suchard, and Gang Li.  A scalable surrogate $L_{0}$ sparse regression method for generalized linear models with applications to large scale data.  _Journal of Statistical Planning and Inference_ , 213:262–281, 2021.
  * Li and Chekouo (2022) Weibing Li and Thierry Chekouo.  Bayesian group selection with non-local priors.  _Computational Statistics_ , 37(1):287–302, 2022\.
  * Liu and Li (2016) Zhenqiu Liu and Gang Li.  Efficient regularized regression with penalty for variable selection and network construction.  _Computational and Mathematical Methods in Medicine_ , vol. 2016(Article ID 3456153):12 pages, 2016.
  * Mahmoudi and Lu (2022) Fatemeh Mahmoudi and Xuewen Lu.  Penalized variable selection with broken adaptive ridge regression for semi-competing risks data.  _arXiv preprint arXiv:2211.09895_ , 2022.
  * Mardis (2011) Elaine R Mardis.  A decade’s perspective on DNA sequencing technology.  _Nature_ , 470(7333):198–203, 2011.
  * Marenberg et al. (1994) Marjorie E Marenberg, Neil Risch, Lisa F Berkman, Birgitta Floderus, and Ulf de Faire.  Genetic susceptibility to death from coronary heart disease in a study of twins.  _New England Journal of Medicine_ , 330(15):1041–1046, 1994.
  * Mittal et al. (2014) Sushil Mittal, David Madigan, Randall S Burd, and Marc A Suchard.  High-dimensional, massive sample-size Cox proportional hazards regression for survival analysis.  _Biostatistics_ , 15(2):207–221, 2014.
  * Paraboschi et al. (2020) Elvezia Maria Paraboschi, Amit Vikram Khera, Piera Angelica Merlini, Laura Gigante, Flora Peyvandi, Mark Chaffin, Marzia Menegatti, Fabiana Busti, Domenico Girelli, Nicola Martinelli, et al.  Rare variants lowering the levels of coagulation factor x are protective against ischemic heart disease.  _Haematologica_ , 105(7):e365, 2020.
  * Schwarz (1978) Gideon Schwarz.  Estimating the dimension of a model.  _The Annals of Statistics_ , pages 461–464, 1978.
  * Shah et al. (2010) Svati H Shah, Christopher B Granger, Elizabeth R Hauser, William E Kraus, Jie-Lena Sun, Karen Pieper, Charlotte L Nelson, Elizabeth R Delong, Robert M Califf, L Kristin Newby, et al.  Reclassification of cardiovascular risk using integrated clinical and molecular biosignatures: Design of and rationale for the measurement to understand the reclassification of disease of cabarrus and kannapolis (murdock) horizon 1 cardiovascular disease study.  _American Heart Journal_ , 160(3):371–379, 2010\.
  * Simon et al. (2011) Noah Simon, Jerome Friedman, Trevor Hastie, and Rob Tibshirani.  Regularization paths for Cox’s proportional hazards model via coordinate descent.  _Journal of Statistical Software_ , 39(5):1–13, 2011.
  * Suchard et al. (2013) Marc A Suchard, Shawn E Simpson, Ivan Zorych, Patrick Ryan, and David Madigan.  Massive parallelization of serial inference algorithms for a complex generalized linear model.  _ACM Transactions on Modeling and Computer Simulation (TOMACS)_ , 23(1):1–17, 2013.
  * Sun et al. (2022) Zhihua Sun, Yi Liu, Kani Chen, and Gang Li.  Broken adaptive ridge regression for right-censored survival data.  _Annals of the Institute of Statistical Mathematics_ , 74(1):69–91, 2022.
  * Tibshirani (1996) Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.
  * Tsimikas et al. (2014) Sotirios Tsimikas, Gordon W Duff, Peter B Berger, John Rogus, Kenneth Huttner, Paul Clopton, Emmanuel Brilakis, Kenneth S Kornman, and Joseph L Witztum.  Pro-inflammatory interleukin-1 genotypes potentiate the risk of coronary artery disease and cardiovascular events mediated by oxidized phospholipids and lipoprotein (a).  _Journal of the American College of Cardiology_ , 63(17):1724–1734, 2014.
  * Vohnout et al. (2003) Branislav Vohnout, Augusto Di Castelnuovo, Roberto Trotta, Andria D’Orazi, Gaetano Panniteri, Anna Montali, Maria Benedetta Donati, Marcello Arca, and Licia Iacoviello.  Interleukin-1 gene cluster polymorphisms and risk of coronary artery disease.  _Haematologica_ , 88(1):54–60, 2003.
  * Wang and Yan (2021) Wenjie Wang and Jun Yan.  Shape-restricted regression splines with R package splines2.  _Journal of Data Science_ , 19(3):498–517, 2021\.  doi: 10.6339/21-JDS1020.
  * Wang et al. (2021) Zixian Wang, Qian Zhu, Yibin Liu, Shiyu Chen, Ying Zhang, Qilin Ma, Xiaoping Chen, Chen Liu, Heping Lei, Hui Chen, et al.  Genome-wide association study of metabolites in patients with coronary artery disease identified novel metabolite quantitative trait loci.  _Clinical and Translational Medicine_ , 11(2), 2021.
  * Wu et al. (2020) Qiwei Wu, Hui Zhao, Liang Zhu, and Jianguo Sun.  Variable selection for high-dimensional partly linear additive Cox model with application to Alzheimer’s disease.  _Statistics in Medicine_ , 39(23):3120–3134, 2020\.
  * Wu and Lange (2008) Tong Tong Wu and Kenneth Lange.  Coordinate descent algorithms for Lasso penalized regression.  _The Annals of Applied Statistics_ , 2(1):224–244, 2008.
  * Zhang (2010) Cun-Hui Zhang.  Nearly unbiased variable selection under minimax concave penalty.  _The Annals of Statistics_ , 38(2):894–942, 2010\.
  * Zhang and Oles (2001) Tong Zhang and Frank J Oles.  Text categorization based on regularized linear classification methods.  _Information Retrieval_ , 4(1):5–31, 2001.
  * Zhao et al. (2018) Hui Zhao, Dayu Sun, Gang Li, and Jianguo Sun.  Variable selection for recurrent event data with broken adaptive ridge regression.  _Canadian Journal of Statistics_ , 46(3):416–428, 2018.
  * Zhao et al. (2019) Hui Zhao, Qiwei Wu, Gang Li, and Jianguo Sun.  Simultaneous estimation and variable selection for interval-censored data with broken adaptive ridge regression.  _Journal of the American Statistical Association_ , 115(529):204–216, 2019.
  * Zou (2006) Hui Zou.  The adaptive Lasso and its oracle properties.  _Journal of the American Statistical Association_ , 101(476):1418–1429, 2006.
  * Zou and Hastie (2005) Hui Zou and Trevor Hastie.  Regularization and variable selection via the elastic net.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 67(2):301–320, 2005.

## Appendix

### Cyclic coordinate descent algorithm

The cyclic coordinate descent algorithm first sets all parameters to some
chosen initial value. It solves a one-dimensional optimization problem by
estimating the first parameter that minimizes the objective function, while
holding the other parameters constant. It then estimates subsequent parameters
by solving a one-dimensional problem and holding the other parameters
constant. When all parameters are estimated, the iteration is complete and
returns to the first parameter for the algorithm to be repeated. Multiple
iterations are done over the whole set of parameters until the pre-specified
convergence criteria is met. Computing the Hessian matrix or gradient vector,
inverting the Hessian matrix, and transposing the gradient vector are not
required in the CLG algorithm, as only one-dimensional updates are used. As a
result, the CLG algorithm easily scales to high-dimensional data (Wu and
Lange, 2008; Simon et al., 2011; Gorst-Rasmussen and Scheike, 2012). It has
then been implemented under GLM for massive sample size data (Suchard et al.,
2013; Li et al., 2021) and high-dimensional massive sample size Cox PH model
(Mittal et al., 2014). Suchard et al. (2013) and Mittal et al. (2014)
developed an R package `Cyclops` that incorporates ridge and LASSO
regularization using Gaussian and Laplace priors. More specifically, when
$s=0$, the following prior is needed to obtain the initial estimates
$\widehat{\boldsymbol{\beta}}^{(0)}:f(\beta_{j}|\xi_{n})\sim N(0,1/\xi_{n})$.
For each $s\geq 1$, the reweighted prior
$f(\widehat{\beta}^{(s)}_{j}|\lambda_{n})\sim
N(0,\widehat{\beta}^{(s-1)2}_{j}/\lambda_{n})$ is used to obtain
$\widehat{\boldsymbol{\beta}}^{(s)}$.

The CLG algorithm finds
$\\{\alpha^{(\text{new})}_{j},\beta^{(\text{new})}_{j},\gamma^{(\text{new})}_{jk}\\}$,
where $\\{\alpha^{(\text{new})}_{j},\beta^{(\text{new})}_{j}\\}$ are the
updated estimates of the $j^{th}$ entry of $\boldsymbol{\alpha}$ and
$\boldsymbol{\beta}$, and $\gamma^{(\text{new})}_{jk}$ is the updated
estimated of the $(j,k)^{th}$ entry of $\boldsymbol{\gamma}$, while keeping
the other values of $\alpha_{j}$’s,$\beta_{j}$’s and $\gamma_{jk}$’s constant.
Suppose we have a tuning parameter $1/\phi_{j}$. Therefore, when $s\geq 0$,
finding $\beta^{(\text{new})}_{j}$ for $j=1,\ldots,p$ is equivalent to finding
$u$ that minimizes the penalized negative log-likelihood

$\begin{split}f(u)=&-\sum^{n}_{i=1}\bigg{\\{}y_{i}\Big{(}ux_{ij}+\sum^{p}_{k=1,j\neq
k}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\Big{)}\\\
&-\log\Big{(}1+\exp\big{\\{}ux_{ij}+\sum^{p}_{k=1,j\neq
k}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\big{\\}}\Big{)}\bigg{\\}}\\\
&+\frac{u^{2}}{2\phi_{j}}.\end{split}$ (5.1)

At the current $\beta_{j}$, one could use the Taylor expansion to approximate
$f(u)$ by

$f(u)\approx
f(\beta_{j})+f^{\prime}(\beta_{j})(u-\beta_{j})+\frac{1}{2}f^{\prime\prime}(\beta_{j})(u-\beta_{j})^{2}.$

In (5.1), different penalty terms can used. For example, $\phi_{j}=1/\xi_{n}$
and $\phi_{j}=(\widehat{\beta}^{(s-1)}_{j})^{2}/\lambda_{n}$ in the GPLM-BAR
algorithm. Similarly, finding $\alpha^{(\text{new})}_{j}$ for
$j=1,\ldots,q_{w}$ is equivalent to minimizing

$\begin{split}g(u)=&-\sum^{n}_{i=1}\bigg{\\{}y_{i}\Big{(}\sum^{p}_{k=1}\beta_{k}x_{ik}+uw_{ij}+\sum^{q_{w}}_{k=1,j\neq
k}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\Big{)}\\\
&-\log\Big{(}1+\exp\big{\\{}\sum^{p}_{k=1}\beta_{k}x_{ik}+uw_{ij}+\sum^{q_{w}}_{k=1,j\neq
k}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\big{\\}}\Big{)}\bigg{\\}},\end{split}$

and finding $\gamma^{(\text{new})}_{jk}$ for $j=1,\ldots,q_{z}$ and
$k=1,\ldots,m_{j}$ is equivalent to minimizing

$\begin{split}h(u)&=-\sum^{n}_{i=1}\Bigg{\\{}y_{i}\Big{(}\sum^{p}_{k=1}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+uB_{jk}(z_{ij},m_{j},c_{j},u_{j})+\sum^{q_{z}}_{s=1,j\neq
s}\sum^{m_{s}}_{t=0,k\neq
t}\gamma_{st}B_{st}(z_{is},m_{s},c_{s},u_{s})\Big{)}\\\
&-\log\bigg{(}1+\exp\Big{\\{}\sum^{p}_{k=1}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+uB_{jk}(z_{ij},m_{j},c_{j},u_{j})+\sum^{q_{z}}_{s=1,j\neq
s}\sum^{m_{s}}_{t=0,k\neq
t}\gamma_{st}B_{st}(z_{is},m_{s},c_{s},u_{s})\Big{\\}}\bigg{)}\Bigg{\\}}.\end{split}$

Suppose the Taylor series is also used to approximate $g(u)$ and $h(u)$, then
values of
$\\{\alpha^{(\text{new})}_{j},\beta^{(\text{new})}_{j},\gamma^{(\text{new})}_{jk}\\}$
can be computed as

$\beta^{(\text{new})}_{j}=\beta_{j}-\frac{f^{\prime}(\beta_{j})}{f^{\prime\prime}(\beta_{j})},\;\;\alpha^{(\text{new})}_{j}=\alpha_{j}-\frac{g^{\prime}(\alpha_{j})}{g^{\prime\prime}(\alpha_{j})},\;\;\gamma^{(\text{new})}_{jk}=\gamma_{jk}-\frac{h^{\prime}(\gamma_{jk})}{h^{\prime\prime}(\gamma_{jk})}.$

The first-order and second-order derivations of $f(u)$, $g(u)$ and $h(u)$ are
computed as follows. When $s=0$, $\phi_{j}=1/\xi_{n}$, then the first-order
and second-order derivatives are

$\begin{split}f^{\prime}(u)&=\left.\frac{\partial}{\partial
u}f(u)\right|_{u=\beta_{j}}\\\
&=-\sum^{n}_{i=1}\left\\{x_{ij}y_{i}-\frac{x_{ij}\exp\\{ux_{ij}+\sum^{p}_{k=1,j\neq
k}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\}}{1+\exp\\{ux_{ij}+\sum^{p}_{k=1,j\neq
k}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\}}\right\\}\\\
&+2\xi_{n}u,\end{split}$

and

$\begin{split}f^{\prime\prime}(u)&=\left.\frac{\partial^{2}}{\partial
u^{2}}f(u)\right|_{u=\beta_{j}}\\\
&=-\sum^{n}_{i=1}\left\\{-\frac{x^{2}_{ij}\exp\\{ux_{ij}+\sum^{p}_{k=1,j\neq
k}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\}}{(1+\exp\\{ux_{ij}+\sum^{p}_{k=1,j\neq
k}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\})^{2}}\right\\}\\\
&+2\xi_{n},\end{split}$

respectively. When $s\geq 1$,
$\phi_{j}=(\widehat{\beta}^{(s-1)}_{j})^{2}/\lambda_{n}$, then

$\begin{split}f^{\prime}(u)&=\left.\frac{\partial}{\partial
u}f(u)\right|_{u=\beta_{j}}\\\
&=-\sum^{n}_{i=1}\bigg{\\{}x_{ij}y_{i}-\frac{x_{ij}\exp\\{ux_{ij}+\sum^{p}_{k=1,j\neq
k}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\}}{1+\exp\\{ux_{ij}+\sum^{p}_{k=1,j\neq
k}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\}}\bigg{\\}}\\\
&+\frac{2\lambda_{n}u}{(\widehat{\beta}^{(s-1)}_{j})^{2}},\end{split}$

and

$\begin{split}f^{\prime\prime}(u)&=\left.\frac{\partial^{2}}{\partial
u^{2}}f(u)\right|_{u=\beta_{j}}\\\ &=-\sum^{n}_{i=1}\left\\{\
-\frac{x^{2}_{ij}\exp\\{ux_{ij}+\sum^{p}_{k=1,j\neq
k}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\}}{(1+\exp\\{ux_{ij}+\sum^{p}_{k=1,j\neq
k}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\})^{2}}\right\\}\\\
&+\frac{2\lambda_{n}}{(\widehat{\beta}^{(s-1)}_{j})^{2}}.\end{split}$

Similarly, the first and second derivatives of $g(u)$ are

$\begin{split}g^{\prime}(u)&=\left.\frac{\partial}{\partial
u}g(u)\right|_{u=\alpha_{j}}\\\ &=-\sum^{n}_{i=1}\bigg{\\{}w_{ij}y_{i}\\\
&-\frac{w_{ij}\exp\\{\sum^{p}_{k=1}\beta_{k}x_{ik}++uw_{ij}+\sum^{q_{w}}_{k=1,j\neq
k}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\}}{1+\exp\\{\sum^{p}_{k=1}\beta_{k}x_{ik}+uw_{ij}+\sum^{q_{w}}_{k=1,j\neq
k}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\}}\bigg{\\}},\end{split}$

and

$\begin{split}g^{\prime\prime}(u)&=\left.\frac{\partial^{2}}{\partial
u^{2}}g(u)\right|_{u=\alpha_{j}}\\\
&=-\sum^{n}_{i=1}\bigg{\\{}-\frac{w^{2}_{ij}\exp\\{\sum^{p}_{k=1}\beta_{k}x_{ik}+uw_{ij}+\sum^{q_{w}}_{k=1,j\neq
k}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\}}{(1+\exp\\{\sum^{p}_{k=1}\beta_{k}x_{ik}+uw_{ij}+\sum^{q_{w}}_{k=1,j\neq
k}\alpha_{k}w_{ik}+\sum^{q_{z}}_{k=1}\sum^{m_{k}}_{l=0}\gamma_{kl}B_{kl}(z_{ik},m_{k},c_{k},u_{k})\\})^{2}}\bigg{\\}}.\end{split}$

Finally, the first and second derivatives of $h(u)$ are

$\begin{split}h^{\prime}(u)&=\left.\frac{\partial}{\partial
u}h(u)\right|_{u=\gamma_{jk}}\\\
&=-\sum^{n}_{i=1}\bigg{\\{}\ddot{B}_{jk}y_{i}-\frac{\ddot{B}_{jk}\exp\\{\sum^{p}_{k=1}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+u\ddot{B}_{jk}+\sum^{q_{z}}_{s=1,j\neq
s}\sum^{m_{s}}_{t=0,k\neq
t}\gamma_{st}\ddot{B}_{st}\\}}{1+\exp\\{\sum^{p}_{k=1}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+u\ddot{B}_{jk}+\sum^{q_{z}}_{s=1,j\neq
s}\sum^{m_{s}}_{t=0,k\neq
t}\gamma_{st}\ddot{B}_{st}\\}}\bigg{\\}},\end{split}$

and

$\begin{split}h^{\prime\prime}(u)&=\left.\frac{\partial^{2}}{\partial
u^{2}}h(u)\right|_{u=\gamma_{jk}}\\\
=&-\sum^{n}_{i=1}\bigg{\\{}-\frac{\ddot{B}^{2}_{jk}\exp\\{\sum^{p}_{k=1}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+u\ddot{B}_{jk}+\sum^{q_{z}}_{s=1,j\neq
s}\sum^{m_{s}}_{t=0,k\neq
t}\gamma_{st}\ddot{B}_{st}\\}}{(1+\exp\\{\sum^{p}_{k=1}\beta_{k}x_{ik}+\sum^{q_{w}}_{k=1}\alpha_{k}w_{ik}+u\ddot{B}_{jk}+\sum^{q_{z}}_{s=1,j\neq
s}\sum^{m_{s}}_{t=0,k\neq
t}\gamma_{st}\ddot{B}_{st}\\})^{2}}\bigg{\\}},\end{split}$

where $\ddot{B}_{jk}=B_{jk}(z_{ik},m_{k},c_{k},u_{k})$ and
$\ddot{B}_{st}=B_{st}(z_{is},m_{s},c_{s},u_{s})$.

## Additional results of the simulation studies

### Scenario 3: CATHGEN-based simulations

We perform a simulation study based on the CATHGEN data, where the motivation
is to examine the performance of our proposed method, in comparison with
existing methods, under the correlation structure of our data. Before
conducting the simulation study, we first have to screen the data to reduce
dimension. The final screened set is 1313 observations and 1242 SNPs. Then, we
use the GPLM-BAR to determine which SNPs are relevant to the binary response,
where 19 SNPs are selected. In the simulation set-up, we set the effect size
of the identified SNPs as four times the true estimated effect size. We do
this because the original effects are too small for the model to distinguish
them. We do not change the effect size of the categorical and continuous
variables. We also do not change the coefficients of the Bernstein
polynomials.

Table 5: Variable selection and estimation results over 200 replications for Scenario 3. Standard deviation of the MMSE are in parentheses. Method | MMSE | TP | FP | MS | MC | TM
---|---|---|---|---|---|---
BAR(AIC) | 0.57(0.25) | 18.84 | 2.22 | 21.06 | 2.38 | 9$\%$
BAR(BIC) | 6.75(0.98) | 12.41 | 0.07 | 12.48 | 6.66 | 0$\%$
LASSO | 5.24(0.65) | 18.97 | 52.73 | 71.70 | 52.76 | 0$\%$
ALASSO | 2.82(0.63) | 18.89 | 27.48 | 46.37 | 27.59 | 0$\%$
Oracle | 0.30(0.26) | 19 | 0 | 19 | 0 | 100$\%$

From Table 5, the GPLM-BAR method with AIC penalty has the lowest average
number of FP and MC. Using the GPLM-BAR with BIC penalty, the average number
of TP decreases due to the large size of the tuning parameter $\lambda_{n}$.
The Lasso and Adaptive Lasso methods produce very large false positives. In
Table 6, the GPLM-BAR with AIC penalty produces the least biased estimates of
the relevant SNPs. In Figure 6, we observe that the estimation of the non-
linear effects of age and log(BMI) is also better when the AIC penalty, as
compared to the BIC penalty, is used.

Table 6: Averaged estimates of the relevant SNPs identified in the real data
analysis, performed over 200 replication for Scenario 3, along with the
clinical and demographical variables. Standard deviations are in parentheses.

True BAR-AIC BAR-BIC LASSO ALASSO Oracle rs1407961 -1.024 -0.922(0.136)
-0.444(0.134) -0.432(0.075) -0.630(0.099) -1.088(0.144) rs8037353 1.017
0.898(0.121) 0.277(0.192) 0.376(0.072) 0.567(0.098) 1.099(0.157) rs7107322
1.015 0.909(0.129) 0.300(0.176) 0.387(0.078) 0.560(0.105) 1.097(0.139)
rs2387952 -1.336 -1.228(0.129) -0.709(0.077) -0.621(0.081) -0.843(0.102)
-1.410(0.153) rs6680365 -0.918 -0.830(0.118) -0.370(0.131) -0.401(0.069)
-0.555(0.096) 0.980(0.131) rs3136558 -0.906 -0.806(0.116) -0.363(0.154)
-0.372(0.067) -0.530(0.091) 0.961(0.135) rs4131888 0.512 0.419(0.137)
0.008(0.049) 0.180(0.070) 0.253(0.106) 0.543(0.118) rs3845439 1.106
0.936(0.276) 0.411(0.154) 0.366(0.123) 0.533(0.168) 1.182(0.136) rs769449
0.773 0.701(0.104) 0.200(0.174) 0.338(0.063) 0.479(0.083) 0.833(0.125)
rs9549675 -1.092 -0.968(0.127) -0.461(0.109) -0.438(0.075) -0.636(0.100)
-1.169(0.157) rs2805543 -1.203 -1.091(0.123) -0.538(0.093) -0.531(0.085)
-0.731(0.106) -1.275(0.144) rs821292 -0.781 -0.666(0.136) -0.033(0.096)
-0.247(0.076) -0.362(0.117) -0.840(0.119) rs12612481 0.770 0.698(0.116)
0.266(0.161) 0.329(0.075) 0.448(0.104) 0.824(0.108) rs7188981 0.905
0.823(0.120) 0.402(0.136) 0.398(0.074) 0.560(0.101) 0.975(0.128) rs9932172
0.997 0.903(0.120) 0.483(0.083) 0.435(0.075) 0.591(0.104) 1.044(0.138)
rs9282537 0.763 0.658(0.129) 0.045(0.117) 0.256(0.074) 0.358(0.119)
0.826(0.123) rs244072 -0.608 -0.509(0.113) 0 -0.162(0.066) -0.234(0.103)
-0.660(0.136) rs11859718 0.760 0.675(0.138) 0.279(0.176) 0.316(0.071)
0.440(0.010) 0.811(0.136) rs17585580 -0.603 -0.491(0.159) 0 -0.153(0.070)
-0.197(0.118) -0.648(0.136) Hypertension 0.368 0.324(0.211) 0.228(0.131)
0.224(0.133) 0.246(0.161) 0.409(0.254) Diabetes 1.180 1.136(0.216)
0.769(0.137) 0.712(0.137) 0.855(0.161) 1.260(0.299) Hypercholesterolomia 1.403
1.294(0.218) 0.867(0.128) 0.879(0.132) 0.992(0.158) 1.483(0.258) Sex 0.361
0.325(0.205) 0.149(0.123) 0.153(0.125) 0.205(0.151) 0.392(0.229) Smoking 0.860
0.798(0.240) 0.431(0.145) 0.431(0.143) 0.546(0.138) 0.938(0.235) HXMI 37.388
42.379(1.097) 38.919(1.716) 11.130(0.310) 12.368(0.435) 26.677(0.803)
Race(African) -0.030 0.031(0.436) 0.308(0.283) 0.114(0.263) 0.116(0.322)
-0.050(0.474) Race(Caucasian) 0.385 0.357(0.378) 0.016(0.225) 0.123(0.223)
0.236(0.274) 0.412(0.411)

Figure 6: Estimation of the non-linear effects of age and log(BMI) based on
the simulated CATHGEN data. The black curve represents the true effect, yellow
curve represents the GPLM-BAR with AIC penalty, and the blue curve represents
the GPLM-BAR with BIC penalty.

### Scenario 4: Poisson partly linear model

In this scenario, we present the details of a simulation study for the Poisson
partly linear model. Let $q_{w}=5$ and $q_{z}=4$. The true values of
$\boldsymbol{\beta}$ is
$\boldsymbol{\beta}_{0}=\\{1,-0.75,0,\ldots,0,-1,0.75,-0.75\\}^{\top}$, and
the true values of $\boldsymbol{\alpha}$ is
$\boldsymbol{\alpha}_{0}=\\{0.75,-0.5,-0.5,0.75,-1\\}^{\top}$. We generate
$\textbf{X},\textbf{W}$, and Z from the same distributions used in Scenarios 1
and 2. By setting $\psi_{1}(z_{i1})=0.1(z_{i1}-3)^{2}$,
$\psi_{2}(z_{i2})=0.2(\cos(2\pi z_{i2})+1)$, $\psi_{3}(z_{i3})=0.2\sin(2\pi
z_{i3})$, and $\psi_{4}(z_{i4})=0.1(z_{i4}+1)^{3}$, we generate $y_{i}$ from a
Poisson distribution with a rate
$\lambda_{i}=\exp\\{\boldsymbol{\beta}^{\top}\textbf{x}_{i}+\boldsymbol{\alpha}^{\top}\textbf{w}_{i}+\sum^{4}_{j=1}\psi_{j}(z_{ij})\\}$.
From Table 7, the selection accuracy of GPLM-BAR method is better than the
LASSO and Adaptive LASSO methods. We set the number of basis $m_{j}=3$, for
$j=1,2,3,4$. We also report the values of the individual estimates of
$\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ in Tables 9 and 8,
respectively, where the individual BAR estimates are also the least biased
among all methods.

Table 7: Variable selection results over 200 replications for Scenario 3. Standard deviations are in parentheses. Method | MMSE | TP | FP | MS | MC | TM
---|---|---|---|---|---|---
$n=600,p=300$
BAR(AIC) | 0.003(0.003) | 5 | 0.90 | 5.90 | 0.90 | 46$\%$
BAR(BIC) | 0.002(0.001) | 5 | 0 | 5 | 0 | 100$\%$
LASSO | 0.368(0.967) | 4 | 0.67 | 4.67 | 1.67 | 37$\%$
ALASSO | 0.106(0.418) | 4.87 | 0 | 4.87 | 0.13 | 96$\%$
Oracle | 0.002(0.001) | 5 | 0 | 5 | 0 | 100$\%$
$n=800,p=300$
BAR(AIC) | 0.003(0.002) | 5 | 1.09 | 6.09 | 1.09 | 38$\%$
BAR(BIC) | 0.001(0.001) | 5 | 0.01 | 5.01 | 0.01 | 99.5$\%$
LASSO | 0.365(0.572) | 4.7 | 0.34 | 5.04 | 0.64 | 66$\%$
ALASSO | 0.407(0.320) | 4.96 | 0 | 4.96 | 0.04 | 98$\%$
Oracle | 0.001(0.001) | 5 | 0 | 5 | 0 | 100$\%$
$n=600,p=450$
BAR(AIC) | 0.004(0.004) | 5 | 0.87 | 5.87 | 0.87 | 44$\%$
BAR(BIC) | 0.002(0.002) | 5 | 0 | 5 | 0 | 100$\%$
LASSO | 0.321(1.064) | 3.75 | 0.90 | 4.65 | 2.15 | 28$\%$
ALASSO | 0.100 (0.525) | 4.80 | 0 | 4.80 | 0.20 | 96$\%$
Oracle | 0.002(0.002) | 5 | 0 | 5 | 0 | 100$\%$
$n=800,p=450$
BAR(AIC) | 0.003(0.003) | 5 | 1.18 | 6.18 | 1.18 | 34$\%$
BAR(BIC) | 0.001(0.001) | 5 | 0 | 5 | 0 | 100$\%$
LASSO | 0.310(1.101) | 3.65 | 0.38 | 4.03 | 1.73 | 47$\%$
ALASSO | 0.102(0.305) | 4.95 | 0 | 4.95 | 0.05 | 98$\%$
Oracle | 0.001(0.001) | 5 | 0 | 5 | 0 | 100$\%$
Table 8: Estimation results of non-zero entries in $\boldsymbol{\beta}$ for Scenario 3. Standard deviations are in parentheses. Method | Bias($\widehat{\beta}_{1}$) | Bias($\widehat{\beta}_{2}$) | Bias($\widehat{\beta}_{p-2}$) | Bias($\widehat{\beta}_{p-1}$) | Bias($\widehat{\beta}_{p}$)
---|---|---|---|---|---
$n=600,p=300$
BAR(AIC) | -0.001(0.021) | 0.002(0.021) | 0.003(0.024) | -0.001(0.023) | 0.002(0.022)
BAR(BIC) | -0.005(0.020) | 0.006(0.021) | 0.007(0.023) | -0.006(0.022) | 0.006(0.022)
LASSO | -0.440(0.286) | 0.385(0.193) | 0.438(0.285) | -0.425(0.177) | 0.374(0.200)
ALASSO | -0.179(0.141) | 0.204(0.130) | 0.189(0.148) | -0.257(0.158) | 0.206(0.138)
Oracle | 0(0.021) | 0(0.020) | 0.002(0.023) | 0.001(0.022) | 0(0.022)
$n=800,p=300$
BAR(AIC) | 0.0006(0.018) | -0.002(0.019) | 0.002(0.02) | -0.001(0.018) | 0.002(0.017)
BAR(BIC) | -0.002(0.018) | 0.002(0.019) | 0.005(0.02) | -0.004(0.017) | 0.005(0.017)
LASSO | -0.355(0.173) | 0.327(0.125) | 0.349(0.171) | -0.378(0.118) | 0.331(0.124)
ALASSO | -0.244(0.127) | 0.278(0.148) | 0.245(0.125) | -0.349(0.179) | 0.291(0.151)
Oracle | 0.001(0.018) | -0.002(0.019) | 0.001(0.020) | 0.001(0.017) | 0.001(0.017)
$n=600,p=450$
BAR(AIC) | 0.003(0.022) | 0.001(0.020) | 0.003(0.025) | 0.001(0.025) | 0.003(0.021)
BAR(BIC) | -0.001(0.022) | 0.004(0.020) | 0.006(0.024) | -0.004(0.024) | 0.008(0.021)
LASSO | -0.459(0.318) | 0.399(0.213) | 0.454(0.320) | -0.43(0.197) | 0.396(0.213)
ALASSO | -0.177(0.175) | 0.195(0.128) | 0.179(0.175) | -0.238(0.136) | 0.198(0.132)
Oracle | 0.004(0.022) | -0.001(0.020) | 0.001(0.024) | 0.003(0.024) | 0.002(0.021)
$n=800,p=450$
BAR(AIC) | 0.001(0.020) | -0.001(0.017) | -0.001(0.019) | -0.001(0.019) | -0.001(0.019)
BAR(BIC) | -0.001(0.020) | 0.002(0.017) | 0.001(0.018) | -0.004(0.019) | 0.002(0.019)
LASSO | -0.469(0.327) | 0.393(0.224) | 0.462(0.330) | -0.430(0.203) | 0.391(0.224)
ALASSO | -0.186(0.117) | 0.215(0.134) | 0.188(0.117) | -0.284(0.166) | 0.220(0.133)
Oracle | 0.002(0.020) | -0.002(0.017) | -0.002(0.018) | 0.001(0.019) | -0.002(0.019)
Table 9: Estimation results of $\boldsymbol{\alpha}$ for Scenario 3. Standard deviations are in parentheses. Method | Bias($\widehat{\alpha}_{1}$) | Bias($\widehat{\alpha}_{2}$) | Bias($\widehat{\alpha}_{3}$) | Bias($\widehat{\alpha}_{4}$) | Bias($\widehat{\alpha}_{5}$)
---|---|---|---|---|---
$n=600,p=300$
BAR(AIC) | -0.001(0.021) | 0.002(0.021) | 0.003(0.024) | -0.001(0.023) | 0.002(0.022)
BAR(BIC) | -0.005(0.020) | 0.006(0.021) | 0.007(0.023) | -0.006(0.022) | 0.006(0.022)
LASSO | -0.440(0.286) | 0.385(0.193) | 0.438(0.285) | -0.425(0.177) | 0.374(0.200)
ALASSO | -0.179(0.141) | 0.204(0.130) | 0.189(0.148) | -0.257(0.158) | 0.206(0.138)
Oracle | 0(0.021) | 0(0.020) | 0.002(0.023) | 0.001(0.022) | 0(0.022)
$n=800,p=300$
BAR(AIC) | 0.0004(0.037) | 0.002(0.037) | -0.001(0.035) | 0.004(0.034) | -0.001(0.034)
BAR(BIC) | 0.0003(0.037) | 0.002(0.036) | 0(0.035) | 0.004(0.033) | -0.0002(0.035)
LASSO | -0.085(0.204) | 0.067(0.157) | 0.056(0.164) | -0.080(0.204) | 0.107(0.253)
ALASSO | -0.032(0.134) | 0.042(0.135) | 0.016(0.125) | -0.033(0.136) | 0.035(0.133)
Oracle | 0.001(0.037) | 0.002(0.036) | -0.0003(0.035) | 0.005(0.032) | -0.001(0.035)
$n=600,p=450$
BAR(AIC) | 0.003(0.044) | -0.001(0.040) | 0.005(0.043) | -0.004(0.045) | 0.002(0.044)
BAR(BIC) | 0.0003(0.044) | 0.001(0.039) | 0.006(0.042) | -0.003(0.045) | 0.003(0.043)
LASSO | -0.242(0.314) | 0.152(0.220) | 0.146(0.228) | -0.228(0.319) | 0.291(0.424)
ALASSO | -0.066(0.161) | 0.040(0.119) | 0.037(0.123) | -0.065(0.161) | 0.071(0.208)
Oracle | 0.002(0.044) | 0.0004(0.040) | 0.005(0.042) | -0.002(0.045) | 0.002(0.043)
$n=800,p=450$
BAR(AIC) | 0.002(0.035) | -0.005(0.034) | 0.0003(0.035) | 0.001(0.036) | -0.006(0.038)
BAR(BIC) | 0.002(0.034) | -0.005(0.034) | 0.001(0.034) | 0.001(0.036) | -0.005(0.038)
LASSO | -0.239(0.322) | 0.152(0.229) | 0.159(0.223) | -0.228(0.331) | 0.304(0.433)
ALASSO | -0.021(0.119) | 0.018(0.104) | 0.018(0.115) | -0.011(0.123) | 0.023(0.124)
Oracle | 0.003(0.034) | -0.005(0.034) | 0.001(0.034) | 0.002(0.036) | -0.006(0.038)
Figure 7: Estimated nonlinear covariate effects $\psi_{j},j=1,2,3,4$ for
$p=300$ and $n=600$ for Scenario 3. The dashed line represents the GPLM-BAR
with AIC penalty, and the dotted line represents GPLM-BAR method with BIC
penalty. Figure 8: Estimated nonlinear covariate effects $\psi_{j},j=1,2,3,4$
for $p=300$ and $n=800$ for Scenario 3. The dashed line represents the GPLM-
BAR with AIC penalty, and the dotted line represents GPLM-BAR method with BIC
penalty. Figure 9: Estimated nonlinear covariate effects $\psi_{j},j=1,2,3,4$
for $p=450$ and $n=800$ for Scenario 3. The dashed line represents the GPLM-
BAR with AIC penalty, and the dotted line represents GPLM-BAR method with BIC
penalty.

### More simulation results for the logistic partly linear model

In this section, we present more simulation results that could not be included
in the main article, by presenting the following tables, where Scenarios 1 and
2 are for the logistic partly linear model and Scenario 3 is for the Poisson
partly linear model.

Table 10: Estimation results of non-zero elements in $\boldsymbol{\beta}$ in Scenario 1. Standard deviations are in parentheses. Method | Bias($\widehat{\beta}_{1}$) | Bias($\widehat{\beta}_{2}$) | Bias($\widehat{\beta}_{p-2}$) | Bias($\widehat{\beta}_{p-1}$) | Bias($\widehat{\beta}_{p}$)
---|---|---|---|---|---
$n=600,p=300$
BAR(AIC) | -0.01(0.15) | 0.005(0.14) | 0.01(0.13) | -0.03(0.14) | -0.01(0.14)
BAR(BIC) | -0.24(0.15) | 0.23(0.15) | 0.26(0.13) | -0.27(0.24) | -0.19(0.22)
LASSO | -0.49(0.11) | 0.48(0.11) | 0.49(0.10) | -0.40(0.10) | -0.32(0.11)
ALASSO | -0.23(0.16) | 0.22(0.16) | 0.23(0.15) | -0.22(0.16) | -0.14(0.15)
Oracle | 0.05(0.15) | -0.06(0.14) | -0.06(0.13) | 0.04(0.13) | 0.04(0.14)
$n=800,p=300$
BAR(AIC) | -0.01(0.11) | 0.01(0.11) | 0.01(0.11) | -0.02(0.12) | -0.02(0.12)
BAR(BIC) | -0.18(0.11) | 0.17(0.11) | 0.18(0.11) | -0.17(0.15) | -0.15(0.14)
LASSO | -0.43(0.09) | 0.43(0.09) | 0.43(0.09) | -0.35(0.09) | -0.29(0.09)
ALASSO | -0.20(0.13) | 0.19(0.13) | 0.20(0.12) | -0.18(0.12) | -0.14(0.12)
Oracle | 0.04(0.12) | -0.04(0.12) | -0.04(0.11) | 0.03(0.12) | 0.02(0.11)
$n=600,p=450$
BAR(AIC) | -0.01(0.15) | 0.004(0.16) | 0.01(0.14) | -0.03(0.14) | -0.01(0.13)
BAR(BIC) | -0.23(0.15) | 0.23(0.17) | 0.27(0.17) | -0.29(0.25) | -0.19(0.20)
LASSO | -0.50(0.11) | 0.49(0.12) | 0.50(0.10) | -0.41(0.10) | -0.33(0.09)
ALASSO | -0.22(0.15) | 0.22(0.16) | 0.22(0.14) | -0.21(0.15) | -0.14(0.13)
Oracle | 0.06(0.15) | -0.06(0.16) | -0.06(0.14) | 0.03(0.13) | 0.04(0.13)
$n=800,p=450$
BAR(AIC) | -0.01(0.21) | 0.01(0.21) | -0.001(0.18) | -0.001(0.18) | -0.01(0.20)
BAR(BIC) | -0.17(0.19) | 0.18(0.19) | 0.17(0.16) | -0.16(0.17) | -0.14(0.18)
LASSO | -0.44(0.17) | 0.45(0.17) | 0.44(0.15) | -0.35(0.15) | -0.30(0.16)
ALASSO | -0.19(0.20) | 0.19(0.19) | 0.18(0.17) | -0.15(0.17) | -0.13(0.19)
Oracle | 0.04(0.12) | -0.04(0.13) | -0.05(0.13) | 0.04(0.12) | 0.02(0.11)
Table 11: Estimation results of non-zero elements in $\boldsymbol{\beta}$ for Scenario 2. Standard deviations are in parentheses. Method | Bias($\widehat{\beta}_{1}$) | Bias($\widehat{\beta}_{2}$) | Bias($\widehat{\beta}_{p-2}$) | Bias($\widehat{\beta}_{p-1}$) | Bias($\widehat{\beta}_{p}$)
---|---|---|---|---|---
$n=600,p=300$
BAR(AIC) | -0.08(0.16) | 0.08(0.15) | 0.09(0.15) | -0.11(0.18) | -0.07(0.18)
BAR(BIC) | -0.46(0.12) | 0.46(0.11) | 0.46(0.12) | -0.35(0.09) | -0.30(0.15)
LASSO | -0.33(0.11) | 0.34(0.10) | 0.33(0.10) | -0.28(0.09) | -0.21(0.10)
ALASSO | -0.22(0.15) | 0.22(0.14) | 0.21(0.14) | -0.21(0.13) | -0.12(0.13)
Oracle | 0.02(0.11) | -0.02(0.10) | -0.02(0.11) | 0.01(0.11) | 0.01(0.12)
$n=800,p=300$
BAR(AIC) | -0.06(0.11) | 0.05(0.11) | 0.06(0.11) | -0.07(0.15) | -0.04(0.13)
BAR(BIC) | -0.41(0.16) | 0.40(0.17) | 0.41(0.16) | -0.32(0.13) | -0.26(0.18)
LASSO | -0.30(0.09) | 0.29(0.09) | 0.29(0.08) | -0.24(0.09) | -0.18(0.08)
ALASSO | -0.19(0.12) | 0.18(0.12) | 0.18(0.12) | -0.17(0.12) | -0.10(0.11)
Oracle | 0.004(0.09) | -0.01(0.10) | -0.01(0.09) | 0.01(0.10) | 0.01(0.10)
$n=600,p=450$
BAR(AIC) | -0.06(0.15) | 0.06(0.14) | 0.10(0.15) | -0.10(0.19) | -0.06(0.17)
BAR(BIC) | -0.45(0.13) | 0.45(0.13) | 0.48(0.09) | -0.36(0.08) | -0.31(0.15)
LASSO | -0.34(0.10) | 0.34(0.10) | 0.35(0.09) | -0.28(0.08) | -0.22(0.10)
ALASSO | -0.20(0.15) | 0.20(0.14) | 0.22(0.13) | -0.19(0.14) | -0.12(0.14)
Oracle | 0.03(0.12) | -0.03(0.11) | -0.01(0.11) | 0.02(0.12) | 0.02(0.11)
$n=800,p=450$
BAR(AIC) | -0.04(0.11) | 0.04(0.11) | 0.05(0.12) | -0.05(0.14) | -0.05(0.15)
BAR(BIC) | -0.38(0.18) | 0.38(0.18) | 0.39(0.17) | -0.31(0.14) | -0.24(0.18)
LASSO | -0.30(0.09) | 0.30(0.09) | 0.30(0.09) | -0.24(0.09) | -0.19(0.09)
ALASSO | -0.16(0.12) | 0.16(0.12) | 0.17(0.12) | -0.15(0.13) | -0.10(0.13)
Oracle | 0.02(0.10) | -0.02(0.10) | -0.02(0.10) | 0.02(0.10) | 0.01(0.10)
Table 12: Estimation results of $\boldsymbol{\alpha}$ for Scenario 1. Standard deviations are in parentheses. Method | $\text{Bias}(\widehat{\alpha}_{1})$ | $\text{Bias}(\widehat{\alpha}_{2})$ | $\text{Bias}(\widehat{\alpha}_{3})$ | $\text{Bias}(\widehat{\alpha}_{4})$ | $\text{Bias}(\widehat{\alpha}_{5})$
---|---|---|---|---|---
$n=600,p=300$
BAR(AIC) | 0.02(0.26) | -0.04(0.23) | 0.01(0.22) | 0.01(0.22) | -0.001(0.24)
BAR(BIC) | -0.10(0.23) | -0.03(0.20) | 0.07(0.19) | -0.09(0.20) | 0.13(0.21)
LASSO | -0.21(0.20) | 0.09(0.18) | 0.12(0.17) | -0.16(0.18) | 0.23(0.19)
ALASSO | -0.09(0.22) | 0.02(0.20) | 0.07(0.19) | -0.08(0.20) | 0.12(0.22)
Oracle | 0.06(0.26) | -0.06(0.24) | 0(0.22) | 0.03(0.23) | -0.03(0.24)
$n=800,p=300$
BAR(AIC) | 0.01(0.20) | -0.02(0.19) | -0.005(0.19) | 0.03(0.19) | -0.02(0.21)
BAR(BIC) | -0.09(0.18) | 0.04(0.17) | 0.04(0.17) | -0.05(0.17) | 0.07(0.19)
LASSO | -0.20(0.16) | 0.09(0.15) | 0.10(0.15) | -0.14(0.15) | 0.18(0.17)
ALASSO | -0.09(0.18) | 0.04(0.17) | 0.05(0.17) | -0.05(0.17) | 0.08(0.20)
Oracle | 0.04(0.20) | -0.02(0.20) | -0.01(0.19) | 0.04(0.20) | 0.05(0.21)
$n=600,p=450$
BAR(AIC) | 0.06(0.23) | -0.02(0.22) | -0.01(0.25) | 0.02(0.24) | -0.07(0.24)
BAR(BIC) | -0.08(0.20) | 0.06(0.19) | 0.06(0.21) | -0.08(0.21) | 0.07(0.21)
LASSO | -0.19(0.18) | 0.11(0.17) | 0.11(0.18) | -0.16(0.19) | 0.18(0.19)
ALASSO | -0.06(0.21) | 0.04(0.19) | 0.05(0.21) | -0.06(0.22) | 0.05(0.22)
Oracle | 0.09(0.23) | -0.02(0.22) | -0.02(0.25) | 0.05(0.24) | 0.10(0.25)
$n=800,p=450$
BAR(AIC) | 0.02(0.21) | -0.01(0.21) | -0.005(0.18) | 0.03(0.18) | -0.03(0.20)
BAR(BIC) | -0.08(0.19) | 0.04(0.19) | 0.05(0.16) | -0.05(0.17) | 0.08(0.18)
LASSO | -0.20(0.17) | 0.10(0.17) | 0.11(0.15) | -0.14(0.15) | 0.19(0.16)
ALASSO | -0.08(0.20) | 0.04(0.19) | 0.05(0.17) | -0.05(0.17) | 0.07(0.19)
Oracle | 0.04(0.21) | -0.02(0.21) | -0.01(0.18) | 0.05(0.19) | -0.04(0.20)
Table 13: Estimation results of $\boldsymbol{\alpha}$ for Scenario 2. Standard deviations are in parentheses. Method | $\text{Bias}(\widehat{\alpha}_{1})$ | $\text{Bias}(\widehat{\alpha}_{2})$ | $\text{Bias}(\widehat{\alpha}_{3})$ | $\text{Bias}(\widehat{\alpha}_{4})$ | $\text{Bias}(\widehat{\alpha}_{5})$
---|---|---|---|---|---
$n=600,p=300$
BAR(AIC) | 0.04(0.21) | -0.02(0.22) | -0.03(0.20) | 0.005(0.22) | -0.004(0.20)
BAR(BIC) | -0.07(0.19) | 0.04(0.20) | 0.03(0.19) | -0.07(0.21) | 0.11(0.18)
LASSO | -0.06(0.19) | 0.03(0.20) | 0.02(0.19) | -0.06(0.20) | 0.09(0.18)
ALASSO | -0.02(0.19) | 0.01(0.21) | 0.005(0.19) | -0.04(0.21) | 0.05(0.18)
Oracle | 0.07(0.21) | -0.04(0.22) | -0.04(0.21) | 0.03(0.22) | 0.04(0.20)
$n=800,p=300$
BAR(AIC) | 0.005(0.18) | 0.01(0.17) | -0.02(0.18) | 0.02(0.17) | -0.03(0.17)
BAR(BIC) | -0.10(0.17) | 0.06(0.17) | 0.03(0.17) | -0.07(0.16) | 0.07(0.16)
LASSO | -0.09(0.17) | 0.05(0.16) | 0.02(0.17) | -0.06(0.15) | 0.06(0.15)
ALASSO | -0.05(0.17) | 0.03(0.16) | 0.004(0.17) | -0.03(0.15) | 0.03(0.15)
Oracle | 0.03(0.18) | -0.003(0.17) | -0.04(0.18) | 0.03(0.17) | -0.05(0.17)
$n=600,p=450$
BAR(AIC) | 0.02(0.21) | -0.01(0.21) | -0.01(0.21) | 0.03(0.19) | -0.05(0.21)
BAR(BIC) | -0.11(0.20) | 0.05(0.20) | 0.06(0.20) | -0.07(0.17) | 0.08(0.20)
LASSO | -0.10(0.19) | 0.05(0.20) | 0.06(0.19) | -0.07(0.16) | 0.08(0.19)
ALASSO | -0.05(0.20) | 0.02(0.20) | 0.03(0.19) | -0.03(0.18) | 0.03(0.20)
Oracle | 0.05(0.21) | -0.03(0.21) | -0.02(0.21) | 0.05(0.19) | 0.08(0.21)
$n=800,p=450$
BAR(AIC) | 0.02(0.21) | -0.01(0.21) | 0.005(0.18) | 0.03(0.18) | -0.03(0.20)
BAR(BIC) | -0.08(0.19) | 0.04(0.19) | 0.05(0.16) | -0.05(0.17) | 0.08(0.18)
LASSO | -0.22(0.16) | 0.11(0.16) | 0.12(0.14) | -0.16(0.14) | 0.22(0.16)
ALASSO | -0.12(0.20) | 0.06(0.18) | 0.07(0.16) | -0.08(0.16) | 0.12(0.18)
Oracle | 0.03(0.19) | -0.01(0.21) | 0.01(0.20) | 0.05(0.18) | 0.04(0.18)
Figure 10: Estimated nonlinear covariate effects $\psi_{j},j=1,2,3,4$ for
$p=300$ and $n=800$ for Scenario 1. The yellow curve represents the GPLM-BAR
with AIC penalty, and the blue curve represents GPLM-BAR method with BIC
penalty. Figure 11: Estimated nonlinear covariate effects $\psi_{j},j=1,2,3,4$
for $p=450$ and $n=600$ for Scenario 1. The yellow curve represents the GPLM-
BAR with AIC penalty, and the blue curve represents GPLM-BAR method with BIC
penalty. Figure 12: Estimated nonlinear covariate effects $\psi_{j},j=1,2,3,4$
for $p=450$ and $n=800$ for Scenario 1. The yellow curve represents the GPLM-
BAR with AIC penalty, and the blue curve represents GPLM-BAR method with BIC
penalty. Figure 13: Estimated nonlinear covariate effects $\psi_{j},j=1,2,3,4$
for $p=300$ and $n=600$ for Scenario 2. The yellow curve represents the GPLM-
BAR with AIC penalty, and the blue curve represents GPLM-BAR method with BIC
penalty. Figure 14: Estimated nonlinear covariate effects $\psi_{j},j=1,2,3,4$
for $p=300$ and $n=800$ for Scenario 2. The yellow curve represents the GPLM-
BAR with AIC penalty, and the blue curve represents GPLM-BAR method with BIC
penalty. Figure 15: Estimated nonlinear covariate effects $\psi_{j},j=1,2,3,4$
for $p=450$ and $n=600$ for Scenario 2. The yellow curve represents the GPLM-
BAR with AIC penalty, and the blue curve represents GPLM-BAR method with BIC
penalty. Figure 16: Estimated nonlinear covariate effects $\psi_{j},j=1,2,3,4$
for $p=450$ and $n=800$ for Scenario 2. The yellow curve represents the GPLM-
BAR with AIC penalty, and the blue curve represents GPLM-BAR method with BIC
penalty.

## More real data analysis results

Table 14: SNPs and its associated gene selected by the GPLM-BAR method. SNP | Gene
---|---
rs1407961 | ZMYM2
rs8037353 | GABRG3
rs7107322 | PRCP
rs2387952 | ADARB2
rs6680365 | CAMTA1
rs3136558 | IL1B
rs4131888 | SLC7A11
rs3845439 | CACNA1E
rs769449 | APOE
rs9549675 | F10
rs2805543 | ADARB2
rs821292 | GFOD1
rs12612481 | PDE11A
rs7188981 | RBFOX1
rs9932172 | RBFOX1
rs9282537 | ABCA1
rs244072 | ADA
rs11859718 | CDH13
rs17585580 | CNTN5