Identifying Genetic Variants for Obesity Incorporating Prior Insights: Quantile Regression with Insight Fusion for Ultra-high Dimensional Data

Abstract

Obesity is widely recognized as a critical and pervasive health concern. We strive to identify important genetic risk factors from hundreds of thousands of single nucleotide polymorphisms (SNPs) for obesity. We propose and apply a novel Quantile Regression with Insight Fusion (QRIF) approach that can integrate insights from established studies or domain knowledge to simultaneously select variables and modeling for ultra-high dimensional genetic data, focusing on high conditional quantiles of body mass index (BMI) that are of most interest. We discover interesting new SNPs and shed new light on a comprehensive view of the underlying genetic risk factors for different levels of BMI. This may potentially pave the way for more precise and targeted treatment strategies. The QRIF approach intends to balance the trade-off between the prior insights and the observed data while being robust to potential false information. We further establish the desirable asymptotic properties under the challenging non-differentiable check loss functions via Huber loss approximation and nonconvex SCAD penalty via local linear approximation. Finally, we develop an efficient algorithm for the QRIF approach. Our simulation studies further demonstrate its effectiveness.

Jiantong Wang^a, Heng Lian^b, Yan Yu^a, and Heping Zhang^c

Keywords: BMI; Penalized regression; SCAD; Ultra-high dimensional data; Variable selection

1 Introduction

The escalating prevalence of obesity worldwide presents a daunting challenge for public health, endangering individual well-being and placing an overwhelming burden on healthcare systems. According to the World Health Organization (WHO), more than one billion people worldwide grappled with obesity in 2022, and this number is still increasing. It is estimated that over four million people die each year as a result of obesity. In a March 2023 report from the World Obesity Atlas, the associated healthcare costs are projected to reach a staggering $4 trillion annually by 2035. It is imperative to delve into the underlying risk factors associated with obesity so that effective strategies and interventions can be developed. Despite notable progress in research and treatment strategies, many unknown risk factors remain to be discovered, partly due to the complex nature of obesity (Loos and Yeo (2022)).

In this paper, we aim to identify important genetic risk factors associated with obesity over hundreds of thousands of single nucleotide polymorphisms (SNPs), incorporating prior insights from previous studies via a novel method, Quantile Regression with Insight Fusion (QRIF), for ultra-high dimensional data. Identifying a sparse set of important genetic risk factors for obesity over a large number of SNPs is important, yet challenging.

To analyze obesity, we advocate a simultaneous variable selection and estimation method through penalized quantile regression. Obesity, as defined by the World Health Organization (WHO), is a condition where an individual’s body mass index (BMI) is 30 or higher. BMI is calculated by dividing an individual’s weight in kilograms by the square of their height in meters. For obesity research, it is vitally important and necessary to investigate genetic risk factors for high or abnormal quantiles of BMI. Quantile regression, which has shown great success (Koenker and Bassett (1978),Yu et al. (2003)), is naturally more suitable for characterizing the conditional high quantiles of BMI that tend to be heterogeneous. In addition, conditional quantiles can provide a more comprehensive picture of BMI level specific risk factors, which could potentially lead to discovery for hidden risk factors and build a better view of genetic mechanism behind various levels of obesity. With the better understanding of conditional quantiles, healthcare providers could potentially develop targeted treatments, especially for individuals in high quantiles of BMI. Simultaneous variable selection and estimation method is crucial because it allows for a more comprehensive analysis of complex dataset, particularly in ultra-high dimensional genetic data. Traditional genome-wide association studies (GWAS) rely on univariate mean regression approach, which only focuses on the marginal mean effects. By employing a simultaneous penalized quantile method, we can better capture the potential joint effect for heterogeneous data, which may lead to more accurate and reliable variable selection and estimation results, and ultimately provide us with a more comprehensive understanding of the underlying genetic architecture of obesity.

We further present a novel approach to leveraging the rich insights from extensive previous studies or expert knowledge for obesity via quantile regression with insight fusion (QRIF). Traditional GWAS have been effective in pinpointing key SNPs linked to numerous phenotypes, such as those involved in BMI studies (Yengo et al. (2018)). Naturally, it will be beneficial to learn from prior insights from such GWAS studies. If we can fully trust the correctness of the set of important genetic factors, intuitively we can simply remove the penalty term with respect to those known important factors in the penalized quantile regression framework, which achieves the simultaneous estimation and variable selection via minimizing a quantile check-loss function and a penalty term on the variables. However, prior insights may contain information that is biased or not generalizable to different populations.

QRIF can incorporate useful insights and stay robust with potentially false information. It has three components in the objective function: a quantile check-loss function for observed data, a penalty term for variable selection, and a quantile check-loss function for prior insights. A tuning parameter is introduced to further balance the trade-off between the information from prior insights and the observed data. A desirable feature of our approach is the flexibility to accommodate a variety of insight formats, including important variables identified in earlier research, established models, or predicted model coefficients from previous studies.

In our obesity analysis, we use the Framingham Heart Study data that include over 500,000 SNPs. We integrate findings from prior insights, e.g., a meta-analysis of the Genetic Investigation of ANthropometric Traits (GIANT) consortium and UK Biobank, along with important genes identified as influencing the metabolic pathway to obesity. The marginal distribution of BMI of the data we use shows clear heterogeneity and a heavy tail. Its 80th percentile has a BMI of about 30, the obesity demarcation.

Interestingly, we find some new SNPs at the 80th conditional quantile of BMI that have not previously been linked with obesity traits, using our QRIF approach. Notably, rs12703441 in the MGAM2 gene, is predicted to be involved in carbohydrate metabolic process, through which it may possibly indirectly affect body shape by affecting nutrient absorption. DSCAML1 is verified by Cortes and Wevrick (2018) to be associated with autism spectrum disorder (ASD), a condition with obesity as a common comorbidity, which may indicate a plausible pathway to have influence on obesity status. Rs6125186 in LOC101927457 is found by Rashid (2020) to be associated with congenital hyperinsulinism, of which early on-set obesity is a known feature according to Zenker et al. (2023). Moreover, SNPs such as rs1957902 in the PRKCH gene, confirmed by Wheeler et al. (2013) to be associated with early-onset extreme obesity, are identified at a high quantile of 80% but not at the median, indicating a quantile-specific genetic influence on obesity. Many of the genes identified in our study can be validated by the medical literature as being connected to obesity-related characteristics, such as BMI, body height, body weight, waist or hip circumference, or waist-hip ratio. For example, rs2033236 in CDH9, rs964841 in TAFA2, and rs1899689 in CADPS2 have been linked to BMI (Kichaev et al. (2019)), while rs2186948 in KCTD1 is associated with body height (Yengo et al. (2022)). Furthermore, rs1873691 in KCNMA1 has been verified to be connected with human obesity (Jiao et al. (2011)), and rs2569034 in SGCD is related to waist-hip ratio and body height (Pulit et al. (2019)).

One appealing characteristic of our QRIF approach is that different important genetic risk factors can be selected for different conditional quantile levels. Genetic associations can vary between average levels and extreme values of BMI. For example, Cotsapas et al. (2009), Schlauch et al. (2020) and Helgeland et al. (2022) use GWAS one-variable-at-a-time through case-control studies and find unique genes associated with extreme obesity, differing from genes found associated with (mean) BMI levels in the traditional GWAS. Notably, gene RFTN1 (aka RAFTLIN) in Cotsapas et al. (2009) is consistent with our findings for $\tau=0.9$. This distinction in genetic effect is further supported by Sag et al. (2023), who report polygenic risk score performs well in discriminating obesity versus non-obesity, but has weak effect size association with mean BMI, which indicates potentially different genetic effect for extreme BMI and mean level of BMI. Goodarzi (2018) points out that despite numerous overlapping SNPs across various traits including extreme obesity and BMI, distinct genetic loci are frequently observed, emphasizing the need for targeted genetic studies across different BMI levels. These works are largely based on GWAS, which only considers marginal association and focuses on extreme obesity, where the continuous BMI level reduces to a binary response variable of extreme obesity or not, which may inevitably lose some useful information.

In our obesity analysis, we incorporate prior insights from previous meta-analyses of the GIANT consortium and UK Biobank. Traditional meta-analysis has made significant achievements, particularly in identifying new genetic loci associated with common phenotypes, and the majority of today’s recognized genetic risk variants are products of extensive meta-analyses of GWAS (Evangelou and Ioannidis (2013)). However, the inherent constraints of meta-analysis, especially its emphasis on homogeneity in study design and analytical methods (Yuan et al. (2016)), highlight the need for more versatile methodologies. Compared to traditional meta-analyses, our approach is distinguished by its adaptability and flexibility. This flexibility is particularly critical when there is a necessity to integrate domain knowledge from domain experts, which is hard to include in traditional meta-analytical frameworks. Beyond its inherent capability, by proper selection of tuning parameters, our model shows robustness to potentially misspecified or biased information from previous studies. Our approach, featuring its adaptability and robustness, provides an alternative method to integrate pre-established insights in the exploration of genetic risk factors for obesity.

To the best of our knowledge, in application, this is the first attempt to simultaneously discover and estimate important genetic risk factors over hundreds and thousands of SNPs for obesity, or high conditional quantiles of BMI, and incorporate insights from previous studies using QRIF model. Furthermore, our study offers methodological, theoretical, and computational contributions. In methodology, we have developed a novel QRIF approach designed for quantile analysis in the context of ultra-high dimensional datasets, incorporating insights from previous studies in flexible forms. Our QRIF model is partly inspired by Jiang et al. (2016)’s pLASSO model that very nicely integrates prior information to LASSO (Tibshirani (1996)) in the mean regression context. However, the set of genetic variables used in Jiang et al. (2016) real data analysis is at a much smaller scale of a total of 916 SNPs for bipolar disorder disease. More importantly, while quantile regression is desirable, major challenges arise both in theory and computation because the quantile check-loss functions for the current observed data and the prior insights are non-differentiable. The QRIF solution is no longer analytical and cannot use similar LASSO formulations as in Jiang et al. (2016). Instead, we employ Huber-loss approximation (HLA) for the quantile check-losses. We adopt SCAD penalty function (Fan and Li (2001)) for its desirable oracle properties. We use local linear approximation (LLA) to tackle the difficulties introduced by non-concave nature of the SCAD penalty. In theory, we establish challenging asymptotic properties including oracle properties with these approximations for QRIF for ultra-high dimensional data.

Computationally, we further develop an efficient algorithm and utilize parallel computing in order to accommodate the large-scale data and complex computation in our application. We optimize and rewrite our original R code in C++ that can boost the computational speed substantially, over 1000 times, in our study. Finally, even though this work primarily focuses on obesity analysis, the proposed QRIF method can not only help analyze obesity but is also expected to be applicable across diverse contexts such as blood pressure, glucose levels, LDL cholesterol, and HDL cholesterol levels, where abnormal conditional quantiles are of most interest.

The structure of this article is as follows: Section 2 introduces our QRIF method for obesity or high conditional quantiles of BMI and provides an efficient algorithm along with implementation specifics. Section 3 presents our application to obesity analysis and describes newly identified SNPs associated with high quantiles of BMI. In Section 4, we present simulation studies to further demonstrate the performance of QRIF. In Section 5, we establish theoretical properties in ultra-high dimensional settings. The proofs and additional numerical studies are relegated to supplementary materials.

2 Quantile Regression with Insight Fusion (QRIF) for Obesity

2.1 Model

In our obesity analysis, the response variable we use is BMI. We use $\mathbf{y}=(y_{1},y_{2},...,y_{n})$ to denote the BMI measurements of participant $i=1,...,n$. We denote the covariate matrix as $\mathbf{X}_{n\times d_{n}}$, which includes genetic and other risk factors, and we use $\mathbf{x}_{i}^{\top}=(x_{i1},x_{i2},\dots,x_{id_{n}})$ to denote the covariate vector of participant $i$. Note that the dimension of the covariates, $d_{n}$, is permitted to increase and can even reach ultra-high dimensionality at the exponential rate of $n$ (Fan and Lv, 2008). We utilize a subset of a popular ongoing NIH study, Framingham Heart Study (FHS) (Dawber et al. (1951)) with $n=1,964$ participants. The number of SNPs we use is over 500,000, which is much higher than the number of participants.

Fortunately, there are many well-established studies investigating the genetic risk factors of obesity. Notably, some of these investigations have been conducted using large biobank studies, such as UK Biobank. Furthermore, international collaborations, such as the Genetic Investigation of ANthropometric Traits (GIANT) consortium, have played a vital role in advancing our understanding of obesity. The GIANT consortium, which involves a collaboration of researchers from various institutions, countries, and studies, has primarily focused on meta-analysis of genome-wide association data and other extensive genetic datasets (Yengo et al. (2018)). Based on these foundational works, we propose to utilize the valuable insights from these previous studies, especially focusing on the confirmed significant SNPs in earlier meta-analyses.

To simultaneously select and estimate important risk factors for high quantiles of BMI with ultra-high dimensional data, while incorporating insights from established studies, we propose and employ our QRIF method. For a given $\tau\in(0,1)$, we employ the conditional quantile model below, to estimate and select variables from ultra-high dimensional covariates,

$$\theta_{\tau}(\mathbf{y}_{i}|\mathbf{x}_{i})=\mathbf{x}_{i}^{\top}\bm{\beta}_{\tau}.$$ (1)

In our settings, supports and values of $\bm{\beta}_{\tau}$ can be different for different conditional quantile levels $\tau$, which allows different genetic risk factors selected for different conditional quantiles of BMI. We center $\mathbf{X}$ and $\mathbf{y}$, therefore we omit the intercept from the coefficient vector and have $\bm{\beta}_{\tau}=(\beta_{\tau 1},\beta_{\tau 2},\dots,\beta_{\tau d_{n}})^{\top}$. We omit the subscript $\tau$ for notational simplicity and denote the coefficient vector as $\bm{\beta}=(\beta_{1},\beta_{2},\dots,\beta_{d_{n}})^{\top}$.

We employ the framework of penalized quantile regression proposed by Wu and Liu (2009), whose objective function is:

$$\mathbf{Q_{\lambda}}(\bm{\beta};\mathbf{X},\mathbf{y}) ={1\over n}\sum_{i=1}^{n}\rho_{\tau}(y_{i}-\mathbf{x}_{i}^{\top}\bm{\beta})+\sum_{j=1}^{d_{n}}p_{\lambda}(|{\beta_{j}}|),$$

where $\rho_{\tau}(u)=u\{\tau-\mathbf{I}(u<0)\}$ is the quantile check loss function. $p_{\lambda}(\cdot)$ is the penalty function with tuning parameter $\lambda$. Various popular penalty functions can be used. In our application, we mainly utilize the Smooth Clipped Absolute Deviation (SCAD) penalty function for its appealing oracle property. The SCAD penalty is a concave and non-decreasing function defined by $p_{\lambda}(0)=0$ and $p^{\prime}_{\lambda}(\beta)=\lambda\left[I(\beta\leq\lambda)+{(a\lambda-\beta)_{+}\over(a-1)\lambda}I(\beta>\lambda)\right]$ for $|\beta|>0$, and $\lambda$ is a regularization penalty parameter. We adopt $a=3.7$ as recommended by Fan and Li (2001).

Our approach equips the flexibility to accommodate a diverse variety of formats of insights from previous studies, such as important variables identified in earlier research, established models, or predicted effect sizes. In our application, we utilize the most common type of prior insights, the identified important variables. Let $\mathbf{S}_{p}$ be the set of indexes for important variables from prior insights. We firstly obtain the initial coefficient estimator $\hat{\bm{\beta}}^{p}$ with prior insights by fitting a traditional penalized quantile regression model with no penalty on $\mathbf{S}_{p}$:

$$\hat{\bm{\beta}}^{p} =\underset{\bm{\beta}}{\text{arg min}}\{\mathbf{Q}_{\lambda,\mathbf{S}_{p}}(\bm{\beta};\mathbf{X},\mathbf{y})\}={1\over n}\sum_{i=1}^{n}\rho_{\tau}(y_{i}-\mathbf{x}_{i}^{\top}\bm{\beta})+\sum_{j\not\in\mathbf{S}_{p}}p_{\lambda}(|{\beta_{j}}|).$$ (2)

In our application of obesity, we adopt identified important SNPs from the meta-analysis of large-scale GWAS including UK Biobank and GIANT (Yengo et al. (2018)). Incorporating predicted values from established work, our objective function of the QRIF model is:

$$\mathbf{Q}_{\lambda,\xi}(\bm{\beta};\mathbf{X},\mathbf{y},\hat{\bm{\beta}}^{p}) ={1-\xi\over n}\sum_{i=1}^{n}\rho_{\tau}(y_{i}-\mathbf{x}_{i}^{\top}\bm{\beta})+\sum_{j=1}^{d_{n}}p_{\lambda}(|\beta_{j}|)+{\xi\over n}\sum_{i=1}^{n}\rho_{\tau}(\mathbf{x}_{i}^{\top}\hat{\bm{\beta}}^{p}-\mathbf{x}_{i}^{\top}\bm{\beta}),$$ (3)

where $\xi\in[0,1]$ is a tuning parameter balancing the weight between the quantile check-loss of the observed data and that of the prior insights. When $\xi=1$, our QRIF method solely relies on the prior insights; when $\xi=0$, our method reduces to a traditional penalized quantile regression model. Note that here $\bm{\beta}$, $\hat{\bm{\beta}}^{p}$, $\rho(\cdot)$, and $\mathbf{Q}(\cdot)$ are different for different $\tau$. Again, we omit the subscript $\tau$ for notational simplicity.

In the QRIF method, the selection and estimation of important variables by minimizing the objective function (3) present significant challenges, including the non-differentiable nature of the quantile check loss function and the nonconvex nature of the SCAD penalty function. We tackle these challenges by adopting Huber loss approximation (HLA) as in Yi and Huang (2017) and Sherwood and Li (2022) for quantile check loss functions and local linear approximation (LLA) as in Zou and Li (2008) for the SCAD penalty function.

The check loss function $\rho_{\tau}(u)=u\{\tau-\mathbf{I}(u<0)\}$ is equivalent to $\rho_{\tau}(u)={1\over 2}[|u|+(2\tau-1)u]$, which is non-differentiable at $u=0$. We adopt Huber Loss Approximation (HLA, Huber (1973)) that first approximates the non-differentiable absolute value $|u|$ by:

$$g_{\gamma}(u)=\begin{cases}{u^{2}\over 2\gamma},&\text{ if }|u|\leq\gamma\\ |u|-{\gamma\over 2},&\text{ if }|u|>\gamma.\end{cases}$$ (4)

The Huber-approximated quantile loss is defined as:

$$h_{\gamma,\tau}(u)={1\over 2}[g_{\gamma}(u)+(2\tau-1)u].$$ (5)

When $\gamma$ is sufficiently small, $h_{\gamma,\tau}(u)\approx\rho_{\tau}(u)$. In our application, we adopt $\gamma=0.01$ as suggested in Sherwood and Li (2022).

To tackle the nonconvex nature of SCAD penalty function, we adopt local linear approximation (LLA). For any $j$, suppose $\beta_{j}$ is close to and initial estimate $\hat{\beta}^{(0)}_{j}$, we can use $\bm{\phi}^{*}\left(\beta_{j}\mid\hat{\beta}^{(0)}_{j}\right)$ to approximate $p_{\lambda}\left(|\beta_{j}|\right)$ by

$$\bm{\phi}^{*}\left(\beta_{j}\mid\hat{\beta}^{(0)}_{j}\right)=p_{\lambda}\left(|\hat{\beta}^{(0)}_{j}|\right)+p_{\lambda}^{\prime}\left(|\hat{\beta}^{(0)}_{j}|\right)\left(|\beta_{j}|-|\hat{\beta}^{(0)}_{j}|\right).$$ (6)

We have the approximated objective function with HLA and LLA:

$$\begin{split}\mathbf{H}_{\lambda,\xi,\gamma}(\bm{\beta};\mathbf{X},\mathbf{y},\hat{\bm{\beta}}^{p},{\hat{\bm{\beta}}}^{(0)})&\approx{1-\xi\over n}\sum_{i=1}^{n}h_{\gamma,\tau}(y_{i}-\mathbf{x}_{i}^{\top}\bm{\beta})+\sum_{j=1}^{d_{n}}\phi^{*}\left(\beta_{j}\mid\hat{\beta}_{j}^{(0)}\right)+{\xi\over n}\sum_{i=1}^{n}h_{\gamma,\tau}(\mathbf{x}_{i}^{\top}\bm{\beta^{p}}-\mathbf{x}_{i}^{\top}\bm{\beta}),\end{split}$$ (7)

for $\bm{\beta}$ near $\bm{\hat{\beta}}^{(0)}.$ Let $h_{\gamma}(\bm{\beta})={1\over n}\sum_{i=1}^{n}h_{\gamma,\tau}(y_{i}-\mathbf{x}_{i}^{\intercal}\bm{\beta})$ and $h^{p}_{\gamma}(\bm{\beta})={1\over n}\sum_{i=1}^{n}h_{\gamma,\tau}(\mathbf{x}_{i}^{\top}\hat{\bm{\beta}}^{p}-\mathbf{x}_{i}^{\top}\bm{\beta})$. Again we omit the subscript $\tau$ and with some abuse of notation of $h$ for simplicity. Therefore equivalently we have the approximated objective function of QRIF as:

$$\mathbf{H}_{\lambda,\xi,\gamma}(\bm{\beta};\mathbf{X},\mathbf{y},\hat{\bm{\beta}}^{p},{\hat{\bm{\beta}}}^{(0)}) \approx(1-\xi)h_{\gamma}(\bm{\beta})+\xi h_{\gamma}^{p}(\bm{\beta})+\sum_{j=1}^{d_{n}}\phi^{*}\left(\beta_{j}\mid\hat{\beta}_{j}^{(0)}\right).$$ (8)

By minimizing $\mathbf{H}_{\lambda,\xi,\gamma}$ in Equation (8), we obtain the estimator by QRIF as:

$$\widehat{}\mbox{\boldmath$\beta$}^{QRIF}=\operatorname*{arg\,min}_{\mbox{\boldmath$\beta$}}(1-\xi)h_{\gamma}(\mbox{\boldmath$\beta$})+\xi h_{\gamma}^{p}(\mbox{\boldmath$\beta$})+\sum_{j=1}^{d_{n}}p^{\prime}_{\lambda}(|\hat{\beta}^{(0)}_{j}|)|\beta_{j}|.$$ (9)

In our implementation, the LASSO estimator is used as an initial estimate.

2.2 Algorithm

We adopt the cyclic coordinate descent method to minimize the approximated objective function defined in (7). Firstly, we use the estimation from the traditional penalized quantile regression model as the initial value. We use ${\hat{\bm{\beta}}}^{(k)}$ to denote the estimated $\bm{\beta}$ at the k-th iteration. Let us use $\nabla\mathbf{H}_{\lambda,\xi,\gamma}(\hat{\beta}^{(k)}_{j})$ to denote the first derivative of $\mathbf{H}_{\lambda,\xi,\gamma}$ on $\hat{\beta}^{(k)}_{j}$. Therefore, conditional on the k-th iteration, for the j-th coefficient, we can obtain $\hat{\bm{\beta}}$ iteratively as follows:

$$\bm{\hat{\beta}}_{j}^{(k+1)} =\underset{\bm{\beta}_{j}}{\arg\min}\left\{\mathbf{H}_{\lambda,\xi,\gamma}\left(\bm{\beta};\mathbf{X},\mathbf{y},\hat{\bm{\beta}}^{p},\hat{\bm{\beta}}^{(k)}\right)\right\}$$

$$\mathbf{H}_{\lambda,\xi,\gamma}\left(\hat{\beta}^{(k+1)}_{j};\hat{\beta}^{(k)}_{j}\right) \approx\mathbf{H}_{\lambda,\xi,\gamma}(\hat{\beta}^{(k)}_{j})+\left(\hat{\beta}^{(k+1)}_{j}-\hat{\beta}^{(k)}_{j}\right)\nabla\mathbf{H}_{\lambda,\xi,\gamma}(\hat{\beta}^{(k)}_{j})+{\zeta_{j}\over 2}\left(\hat{\beta}^{(k+1)}_{j}-\hat{\beta}^{(k)}_{j}\right)^{2},$$ (10)

where $\zeta_{j}={2\over n\gamma}[\mathbf{X}^{\top}\mathbf{X}]_{(j,j)}$ for ${\hat{\bm{\beta}}}^{(k+1)}$ close to ${\hat{\bm{\beta}}}^{(k)}$. Let us use $\nabla h_{\gamma}(\hat{\beta}_{j}^{(k)})$ to denote the first derivative of $h_{\gamma}(\mbox{\boldmath$\beta$})$ on $\hat{\beta}^{(k)}_{j}$, and $\nabla h^{p}_{\gamma}(\hat{\beta}^{(k)}_{j})$ to denote the first derivative of $h_{\gamma}^{p}(\mbox{\boldmath$\beta$})$ on $\hat{\beta}^{(k)}_{j}$. Take the first derivative of (2.2) and set it to zero. Then we have the estimated $\hat{\beta}_{j}^{(k+1)}$ as:

$$\hat{\beta}_{j}^{(k+1)} =\hat{\beta}_{j}^{(k)}-{1\over\zeta_{j}}\left[(1-\xi)\nabla h_{\gamma}(\hat{\beta}_{j}^{(k)})+\xi\nabla h^{p}_{\gamma}(\hat{\beta}_{j}^{(k)})+p^{\prime}_{\lambda}(|\hat{\beta}_{j}^{(k)}|)\text{sign}(\hat{\beta}_{j}^{(k)})\right],$$ (11)

for each coefficient $\beta_{j}$ of $\bm{\beta}$ for $j=1,\dots,d_{n}$. Then we repeat the iteration of fitting until convergence of $\bm{\beta}$. Pseudo code of this algorithm is shown in Algorithm 1.

To select tuning parameters $\xi$ and $\lambda$, we apply cross-validation. First, based on $\lambda$ selected by penalized quantile regression, we select $\xi$ with minimum objective value of $\mathbf{H}_{\lambda,\xi,\gamma}$. Then, similarly, fix the $\xi$ selected in the previous step, we select the $\lambda$ with minimal check loss.

Algorithm 1 Algorithm of QRIF estimation and variable selection

1:Input: The original data set $\mathbf{X},\mathbf{y}$, prior insight important variable set $\mathbf{S_{p}}$.

2:Obtain the estimator fully trusting the prior insights via penalized quantile regression with a penalty on variables except those in $\mathbf{S_{p}}$:

$$\hat{\bm{\beta}}^{p}=\underset{\bm{\beta}}{\text{arg min}}\{\mathbf{Q}_{\lambda,\mathbf{S}_{p}}(\bm{\beta};\mathbf{X},\mathbf{y})\}={1\over n}\sum_{i=1}^{n}\rho_{\tau}(y_{i}-\mathbf{x}_{i}^{\top}\bm{\beta})+\sum_{j\not\in\mathbf{S}_{p}}p_{\lambda}(|{\beta_{j}}|).$$

3:Obtain an initial estimator via traditional penalized quantile regression:

$$\bm{\hat{\beta}}^{(0)}=\operatorname*{arg\,min}{1\over n}\sum_{i=1}^{n}\rho_{\tau}(y_{i}-\mathbf{x}_{i}^{\top}\bm{\beta})+\sum_{j=1}^{d_{n}}p_{\lambda}(|{\beta_{j}}|).$$

4:Apply cyclical coordinate descending method to calculate $\bm{\hat{\beta}}^{QRIF}$ as following:

5:repeat

6:     Based on $\bm{\hat{\beta}}^{(k)}$

7:     for $j=1,\dots,d_{n}$ do

8:         Update $\hat{\beta}_{j}$ with:

$$\hat{\beta}_{j}^{(k+1)}=\hat{\beta}_{j}^{(k)}-{1\over\zeta_{j}}\left[(1-\xi)\nabla h_{\gamma}(\hat{\beta}_{j}^{(k)})+\xi\nabla h^{p}_{\gamma}(\hat{\beta}_{j}^{(k)})+p^{\prime}_{\lambda}(|\hat{\beta}_{j}^{(k)}|)\text{sign}(\hat{\beta}_{j}^{(k)})\right]$$

9:     end for

10:until Convergence or reaching maximum iteration.

11:Output: $\bm{\hat{\beta}}^{QRIF}$.

3 Obesity Analysis

3.1 Real Data Analysis

The central aim of our study is to identify the important genetic risk factors contributing to obesity while utilizing prior information from previous studies. We use an ultra-high-dimensional genetic dataset from the Framingham Heart Study. We propose and employ a novel quantile regression framework, the QRIF, specifically focusing on the high quantiles of BMI.

The Framingham Heart Study data we use include 1,964 participants and 500,568 SNPs, along with covariates such as age and sex, and the phenotype variable of interest, body mass index (BMI). In our study, the mean BMI stands at 27.94 and the median at 27.25, indicating that average participants fall into the overweight range with BMI greater than 25. Of particular concern are the observations from the higher BMI quantiles: the marginal 80th percentile is 29.79, which is close to the threshold of obesity defined as a BMI over 30, and the marginal 90th percentile reaches 32.42, which is above the threshold of obesity. Obesity or abnormal BMI is of the most interest, and hence our study especially focuses on the high quantiles of BMI.

In the preprocessing phase of our genetic data, we follow the literature and exclude the SNPs with minor allele frequency (MAF) lower than 0.1 and Hardy-Weinberg Equilibrium (HWE) test p-value lower than 0.001. We first conduct a univariate GWAS on BMI. We include SNPs, sex, and age as covariates, following the recommendations from Hoffmann et al. (2018). To account for ancestral diversity, we include the first 5 genetic principal components as covariates (Price et al. (2006)). Based on the results of GWAS, we use the top 4,000 SNPs for the following analysis.

In our study, we incorporate information from prior established studies using our proposed QRIF method. Specifically, we utilize the results from the meta-analysis of the GIANT cohort and the UK Biobank, which focused on traits related to BMI. This meta-analysis provides a comprehensive set of important SNPs, established through rigorous research. We aligned these SNPs with the data from the FHS data, selecting 120 SNPs that overlap between the two sources as our foundation of prior insights. By doing so, we ensure that our QRIF model is not just informed by the current data, but also by insights that have been filtered through the lens of prior established research. This QRIF methodology allows us to focus on the most relevant variables, thereby may help to identify important genetic factors with greater precision and reliability.

We use 10-fold cross-validation to choose the tuning parameters. At 80% conditional quantile of BMI using QRIF, the choice of $\xi$, the tuning parameter balancing the trade-off between the information from the current data and prior insights, is 0.1. This may indicate that findings from different studies using different data may not agree with each other. It also suggests that integrating prior insights into the current study may offer valuable benefits to some degree.

We first focus our analysis on the 80th conditional quantile of BMI to identify the genetic risk factors for obesity through our QRIF approach, which enables us to target risk factors that are most relevant for individuals with higher obesity risk. In addition, we also report the results of the conditional 50% quantile via QRIF.

Table 1: Obesity Analysis. Top loci identified for conditional quantiles at both $\tau=0.8$ and $\tau=0.5$ (Matched), and top loci identified exclusively for $\tau=0.8$ or $\tau=0.5$ (Unmatched). This table lists SNPs, the associated genes, the traits they affect, and references confirming these traits.

$\tau=0.8$				$\tau=0.5$
SNP	Gene	Trait	Reference	SNP	Gene	Trait	Reference
Matched
rs12703441	MGAM2			rs4260813	MGAM2
rs2015768	SEMA3C *	Height Waist-hip ratio	Yengo et al. (2022) Kichaev et al. (2019)	rs2015768	SEMA3C *	Height Waist-hip ratio	Yengo et al. (2022) Kichaev et al. (2019)
rs538079	DSCAML1			rs538079	DSCAML1
rs2033236	CDH9	BMI	Pulit et al. (2019) Zhu et al. (2020)	rs2033236	CDH9	BMI	Pulit et al. (2019) Zhu et al. (2020)
rs1899689	CADPS2	BMI	Pulit et al. (2019) Kichaev et al. (2019)	rs1899689	CADPS2	BMI	Pulit et al. (2019) Kichaev et al. (2019)
rs1873691	KCNMA1	Obesity	Jiao et al. (2011)	rs1873691	KCNMA1	Obesity	Jiao et al. (2011)
rs3767392	PPP1R12B	Estradiol level	Comuzzie et al. (2012)	rs3767392	PPP1R12B	Estradiol level	Comuzzie et al. (2012)
Unmatched
rs871664	ARHGAP29*	Height	Yengo et al. (2022)	rs1491486	-
rs4761670	KRT19P2*	Height	Yengo et al. (2022)	rs10002317	-
rs6125186	LOC101927457			rs9313258	-
rs2237834	PDGFRL	Height	Yengo et al. (2022)	rs4875330	CSMD1	BMI	Pulit et al. (2019) Huang et al. (2022)
rs2567305	PLPPR1	BMI Height	Huang et al. (2022) Schoeler et al. (2023) Kichaev et al. (2019)	rs963158	CFAP20DC-DT	BMI	Locke et al. (2015)
rs13250654	-			rs4435016	ETS1	Height
rs2140450	PPP2R3A	BMI Waist-hip ratio	Schoeler et al. (2023) Lotta et al. (2018)	rs340863*	PROX1-AS1	Obesity	Kim et al. (2013)
rs4558646	LINC01429	Height	Yengo et al. (2022)	rs7602917	NRXN1	BMI	Pulit et al. (2019) Wang et al. (2022)
rs1460181	KIAA0825	Waist circumference Lean body mass	Zhu et al. (2020) Tachmazidou et al. (2017)	rs1926617	GPC5	Waist circumference Height	Tachmazidou et al. (2017) Yengo et al. (2022)
rs7624324	TP63	Height	Yengo et al. (2022)
rs1303470	NPAS3	BMI waist circumference Height	Huang et al. (2022) Tachmazidou et al. (2017) Zhu et al. (2020)	rs2247393	DDX31
rs6081077	ZNF133*	BMI	Turcot et al. (2018)	rs10495699	-
rs4245513	GJA1*	Height	Yengo et al. (2022)	rs7189918	DYNLRB2-AS1	Height	Kichaev et al. (2019) Yengo et al. (2022)
rs6080334	KIF16B*	BMI	Zhu et al. (2020)	rs10865208	PRKCE	BMI Height Weight	Tachmazidou et al. (2017) Yengo et al. (2022) Tachmazidou et al. (2017)

Table 1 presents the top loci identified at both conditional quantiles of $\tau=0.8$ and $\tau=0.5$, along with those unique to each quantile level of BMI, presenting the top SNPs with their corresponding genes, substantiating their association with obesity-related traits, and related scientific literature references. At the 80th quantile, genes such as PDGFRL and TP63 are reported by Yengo et al. (2022) to be linked to body height, while PLPPR1 and PPP2R3A are associated with both BMI (Huang et al. (2022)) and body height (Schoeler et al. (2023), Kichaev et al. (2019)), with PPP2R3A also impacting the waist-hip ratio (Schoeler et al. (2023), Lotta et al. (2018)). KIAA0825’s relationship with BMI-adjusted waist circumference (Zhu et al. (2020)) and lean body mass (Tachmazidou et al. (2017)), along with NPAS3’s broad associations to BMI, waist circumference, and body height, underscores the genetic complexity in obesity. According to Andreacchi et al. (2021), BMI, waist circumference, and waist-hip ratio are recognized as robust indicators of obesity. Moreover, according to Gadekar et al. (2020), there is a strong correlation between visceral fat and waist-hip ratio, which might indicate that waist-hip ratio could bring us more information about obesity.

Notably, rs6125186 in LOC101927457 has not been identified in the previous literature in association with obesity-related traits. According to Rashid (2020), LOC101927457 is found to over-express in congenital hyperinsulinism of infancy (CHI) tissue. Early onset obesity is a known feature of congenital hyperinsulinism according to Zenker et al. (2023), and many children with CHI become overweight in their early life (Banerjee et al. (2022)). These new findings may suggest a potential direction for future researchers to explore further the genetic risk factors associated with obesity.

Among the loci common to both quantiles, rs2015768 close to SEMA3C is correlated with height Yengo et al. (2022), with SEMA3C additionally linked to waist-hip ratio Kichaev et al. (2019), highlighting their potential influence on obesity-related traits. PPP1R12B is confirmed by Comuzzie et al. (2012) to be associated with estradiol level, which is an obesity-related trait. Notably, DSCAML1 and MGAM2 have not been identified in the previous literature in association with obesity-related traits. According to Cortes and Wevrick (2018), DSCAML1 is a verified gene associated with autism spectrum disorder (ASD). Moreover, based on lab-based experiments, Ogata et al. (2021) also verify that DSCAML1 will affect the ability to regulate the number of synapses, and therefore potentially cause neurodevelopmental disorders including ASD. Obesity is one of the common comorbidities of ASD, which might be a plausible pathway and needs further study in the future. According to Kim et al. (2014) and Xu et al. (2019), MGAM2 is homology to MGAM, which is a carbohydrate utilization-related gene and is highly expressed in the small and large intestines. MGAM2 plausibly influences the body shape status through the process of degradation of starch or glycogen, which also deserves further study. Another plausible pathway is that MGAM2 is confirmed by Wray et al. (2018) to be associated with major depressive disorder, while according to Rajan and Menon (2017) obesity and depression have a significant and bidirectional association, which deserves further investigation.

At the 80th conditional quantile, our findings largely align with the existing literature, confirming associations with traits such as height, BMI, waist-hip ratio, waist circumference, and lean body mass. Correspondingly, loci uniquely identified at the median quantile corroborate the relationship between several genes and obesity traits. For instance, CSMD1 and NRXN1 are associated with BMI Pulit et al. (2019); Huang et al. (2022), and ETS1 with body height Pulit et al. (2019); Huang et al. (2022); Zhu et al. (2020); Kichaev et al. (2019); Wang et al. (2022); Yengo et al. (2022); Kim et al. (2013). Notably, rs1957902 in PRKCH identified only for $\tau=0.8$ is implicated in early-onset extreme obesity Wheeler et al. (2013), despite not being among the top findings (not shown in the table), highlighting its potential link to the risk of obesity.

Furthermore, Table 2 consolidates the common loci findings across several quantiles at $\tau=0.5,0.8,0.9$, identifying loci such as CDH9 and TAFA2 with direct associations to BMI confirmed by Huang et al. (2022) and Kichaev et al. (2019). LINC00954 and KCTD1 are reported by Yengo et al. (2022) to be linked to body height, while FOXN3 is verified by Tachmazidou et al. (2017) and Yengo et al. (2022) to be associated with both height and BMI. SGCD is implicated in waist-hip ratio according to Pulit et al. (2019), and TMEM132D with body fat mass according to Tachmazidou et al. (2017), illustrating the diverse genetic influence on various obesity.

Table 2: Common loci identified for conditional quantiles at $\tau$ = 0.5, 0.8 and 0.9. This table lists SNPs, their corresponding chromosomes, the associated genes, the traits they affect, and references confirming these traits.

SNP	Chr	Gene	Trait	reference
rs10495699	$2$	LINC00954*	Height	Yengo et al. (2022), Kichaev et al. (2019)
rs2033236	$5$	CDH9	BMI	Pulit et al. (2019), Zhu et al. (2020)
rs2569034	$5$	SGCD	Waist-hip ratio	Pulit et al. (2019)
rs2015768	$7$
rs2582382	$8$
rs13281810	$8$	LOC10537563
rs964841	$11$	TAFA2	BMI	Comuzzie et al. (2012), Tachmazidou et al. (2017)
rs10507046	$12$
rs10773715	$12$	TMEM132D; LOC124903085	Body fat mass	Tachmazidou et al. (2017)
rs7147927	$14$	FOXN3	Height	Yengo et al. (2022)
			BMI	Tachmazidou et al. (2017)
rs11641334	$16$
rs7207687	$17$
rs2186948	$18$	KCTD1	Height	Yengo et al. (2022)

Our analysis across different BMI quantiles focuses on the high quantiles of most interest while enabling a comprehensive view of the complex genetic factors contributing to different levels of obesity. Particularly, our analysis using quantile regression with insight fusion (QRIF) may reveal genes that are crucial to understanding obesity risk, guiding future genetic research, and the pursuit of targeted intervention strategies.

3.2 Mimic Data Simulations

We further conduct a mimic genetic data simulation study. For this simulation, 400 subjects and 2,000 SNPs for high-dimensional setting and 100,000 SNPs for ultra-high-dimensional setting, are randomly selected from the real data analysis used in Section 3. Age and sex are included as covariates. The response variable is generated according to Model (1), where the error follows a normal distribution with mean of 0 and standard deviation of 1. In the model, the coefficients for age, sex, and 10 important SNPs are assigned a value of 0.5, with all other coefficients set to 0. We employ four different scenarios for our prior insights: S1 (Partial): {7 correct SNPs}, S2 (High quality mixed): {7 correct SNPs, 3 wrong SNPs}, S3 (Low quality mixed): {3 correct SNPs, 7 wrong SNPs}, S4 (Biased): {7 wrong SNPs}. Three different conditional quantile levels $\tau=0.1,0.5,0.9$ are reported.

To assess the performance of the variable selection accuracy, we employ the following measurements: (a) TP: True positive number that counts the correctly selected variables. (b) FP: False positive number that counts the incorrectly selected variables. (c) Bias: $L_{1}$ norm between the estimated coefficients and the true nonzero coefficients. (d) F1: F1 score, a model accuracy measure defined as ${TP\over TP+{1\over 2}(FP+FN)}$, where FN is the false negative number counting the incorrectly predicted zeros. We replicate our simulation for 200 runs and summarize the variable selection results in Tables 3 and 4.

Figure 1: F1 scores of GWAS, traditional penalized quantile regression (QR), penalized quantile regression solely rely on prior insights (Prior), and QRIF in the mimic data simulation with $d=2,000$. Rows 1, 2, and 3, refer to $\tau$ = 0.1, 0.5, and 0.9. Columns 1, 2, 3, 4 refer to scenarios S1, S2, S3, and S4.

Figure 2: F1 scores of GWAS, traditional penalized quantile regression (QR), penalized quantile regression solely rely on prior insights (Prior), and QRIF in the mimic data simulation with $d=100,000$. Rows 1, 2, and 3 refer to $\tau$ = 0.1, 0.5, and 0.9. Columns 1, 2, 3, 4 refer to scenarios S1, S2, S3, and S4.

From the results of our mimic data simulation with $d=2,000$, as shown in Figure 1 and Table 3, we observe that at $\tau=0.5$ and $0.9$, the simultaneous estimation and variable selection approaches of our QRIF approach incorporating prior information and penalized quantile regression with no prior information have higher F1 scores than the traditional GWAS methods on the mean and one-variable-at-a-time. For $\tau=0.1$, though traditional penalized quantile regression has lower F1 scores than GWAS, QRIF outperforms GWAS in scenarios S1, S2, and S3. Across all three levels of $\tau$, when prior insights are informative, as evidenced in scenarios S1, S2, and S3, QRIF consistently outperforms the traditional penalized quantile regression (QR). Even when the prior information is all biased as in scenario S4, QRIF demonstrates robustness against such misspecification, delivering results comparable to those from the penalized quantile regression model and better than the results solely relying on prior information.

In ultra-high-dimensional settings with $d=100,000$, for all three conditional quantile levels, again QRIF outperforms GWAS in scenarios S1, S2, and S3. When confronted with biased prior insights, as in S4, QRIF’s performance closely matches that of the penalized quantile regression, demonstrating its robustness to potential misspecification.

Table 3: Summary of TP, FP, and F1 scores in the mimic data simulation with $d=2,000$

		S1		S2		S3		S4
	QR	Prior	QRIF	Prior	QRIF	Prior	QRIF	Prior	QRIF
$\tau=0.1$
TP	9.83	10.00	9.50	10.00	9.50	9.76	9.37	9.87	9.83
FP	8.165	0.12	0.18	2.75	0.28	12.60	0.67	15.18	8.17
Bias	1.39	0.77	1.24	0.77	1.27	1.23	1.58	1.34	1.39
$\tau=0.5$
TP	10.00	10.00	10.00	10.00	10.00	10.00	10.00	10.00	10.00
FP	1.01	0.03	0.03	2.46	2.46	6.43	0.00	7.02	0.00
Bias	0.89	0.59	0.59	0.58	0.58	0.75	1.21	0.92	1.22
$\tau=0.9$
TP	10.00	10.00	10.00	10.00	10.00	10.00	10.00	10.00	10.00
FP	2.23	0.21	0.24	2.75	0.27	6.64	0.27	8.66	0.98
Bias	0.97	0.78	0.85	0.78	0.87	0.81	0.83	0.98	0.86

Table 4: Summary of TP, FP, and F1 scores for in the mimic data simulation with $d=100,000$

		S1		S2		S3		S4
	QR	Prior	QRIF	Prior	QRIF	Prior	QRIF	Prior	QRIF
$\tau$ = 0.1
TP	6.67	9.69	9.69	9.67	9.67	7.91	6.78	6.69	6.67
FP	12.82	4.30	4.30	7.08	7.08	14.52	1.73	18.79	12.82
Bias	4.35	1.17	1.17	1.16	1.16	2.91	4.17	4.26	4.30
$\tau$ = 0.5
TP	8.81	9.87	8.73	9.85	8.74	9.45	8.73	8.79	8.81
FP	3.75	0.27	1.30	3.29	1.30	8.15	1.33	10.73	3.75
Bias	2.66	0.74	2.42	0.76	2.43	1.39	2.48	2.59	2.67
$\tau$ = 0.9
TP	7.91	9.88	9.88	9.86	9.86	8.92	8.75	7.95	7.91
FP	9.01	0.34	0.34	3.31	3.31	9.28	3.08	16.09	9.01
Bias	3.25	0.84	0.84	0.85	0.85	1.81	2.44	3.11	3.25

4 Numerical Simulations

To further assess the performance of our Quantile Regression with Insight Fusion (QRIF) method on variable selection and estimation for different quantiles in ultra-high dimensional data, we conduct the following numerical simulations.

4.1 Example 1

In our simulations, we first employ a model with the following data generation formula for $i=1,\dots,n$, $y_{i}=\sum_{j=1}^{d_{n}}x_{ij}{\beta_{j}}+\epsilon_{i},$ where the design matrix $\mathbf{X}$ follows the multivariate normal distribution with mean $\mathbf{0}$ and covariance matrix $\mathbf{\Sigma}$. Here $\mathbf{\Sigma}=(\sigma_{ij})$ follows an AR(1) covariance structure and $\sigma_{ij}=\rho^{|i-j|}$, $i,j=1,\dots,d_{n}$. In our simulations, we use $\rho=0.5$. The coefficient vector $\bm{\beta}$ is set as $(1,\dots,1,0,\ldots,0)^{\top}$, where $\beta_{1},\dots,\beta_{20}=1$. The error terms $\epsilon_{i}$ are generated from a Cauchy distribution of $Cauchy(1,0)$. We set the number of covariates at a high dimension of $d_{n}=2,000$ and a sample size of $n=200$. A simpler case using normal error can be found in supplementary materials.

To investigate the model’s performance across different quantiles, we apply our QRIF method at distinct quantile levels and report the results at $\tau=0.5,0.8,$ and $0.9$.

To assess the QRIF model’s performance to incorporate valuable prior insights as well as its robustness against potential erroneous prior information, we examine the following four scenarios: S1: $\{x_{1},x_{2},\ldots,x_{15}\}$; S2: $\{x_{1},\ldots,x_{10}\}\cup\{x_{21}\ldots,x_{25}\}$; S3: $\{x_{1},\ldots,x_{5}\}\cup\{x_{21}\ldots,x_{30}\}$; S4: $\{x_{21},\ldots,x_{40}\}$. To assess the performance of variable selection, we measure the TP, FP, Bias and F1 score as defined in Section 3.2 based on 200 simulation runs.

We compare the performance of our QRIF method and traditional penalized quantile regression without prior insights (QR). We also benchmark with the penalized quantile regression estimator fully trusting the prior information (Prior) based on Equation (2). We use 10-fold cross-validation to select tuning parameters.

Figure 3: Bias in simulation Example 1. Rows 1, 2, 3 refer to $\tau$ = 0.5, 0.8 and 0.9. Columns 1, 2, 3, 4 refer to scenarios S1, S2, S3, and S4.

Figure 4: F1 Scores in simulation Example 1. Rows 1, 2, 3 refer to $\tau$ = 0.5, 0.8 and 0.9. Columns 1, 2, 3, 4 refer to scenarios S1, S2, S3, and S4.

In Example 1, as shown in Figures 3 and 4, we observe notable performance differences between our QRIF method and traditional penalized quantile regression (QR). With partially correct insights in scenario S1, QRIF demonstrates superior performance over QR at all 3 quantile levels, achieving lower Bias and higher F1 score. In scenarios with high-quality mixed insights S2, the QRIF still outperforms QR, especially for high quantiles of $\tau=0.8$ and $0.9$. In scenarios with low-quality mixed insights S3, the QRIF model works better than QR for $\tau=0.9$, confirming its robustness in less-than-ideal conditions. Importantly, even with erroneous insights, QRIF’s performance is comparable to QR, indicating its robustness against misleading prior information.

4.2 Example 2

For a more challenging setting, based on a model employed by Wang et al. (2012), we first generate $\tilde{\mathbf{X}}$ follows the multivariate normal distribution with mean $\mathbf{0}$ and covariance matrix $\mathbf{\Sigma}_{ij}=(\sigma_{ij})=\rho^{|i-j|}$ with $\rho=0.5$. Then we generate the design matrix $\mathbf{X}$ as $X_{1}=\Phi(\tilde{X}_{1}),$ where $\Phi$ is the cumulative distribution function of a standard normal random variable and $X_{j}=\tilde{X}_{j}$ for $j\neq 1$. Then for $i=1,\dots,n,$ we generate data from

$$y_{i}=\sum_{j=1}^{d_{n}}x_{ij}{\beta_{j}}+x_{i1}\epsilon_{i}.$$

We also set the number of covariates at a high dimension of $d_{n}=2,000$ and a sample size of $n=200$. The coefficient vector $\bm{\beta}$ is set as $\beta_{j}=1$ for $j\in\{10,20,\dots,50\}$, and $\beta_{j}=0$ for others. The error terms $\epsilon_{i}\sim N(0,\sigma^{2})$ with $\sigma=1$. Similar to Example 1, we examine the following four scenarios: S1: $\{x_{1},x_{10},\dots x_{40}\}$, S2: $\{x_{1},x_{10},\dots x_{40}\}\cup\{x_{60},x_{70}\}$, S3: $\{x_{1},x_{10},x_{20}\}\cup\{x_{60},\dots,x_{90}\}$, S4: $\{x_{70},\ldots,x_{100}\}$.

Figure 5: Bias in simulation Example 2. Rows 1, 2, 3 refer to $\tau$ = 0.5, 0.8 and 0.9. Columns 1, 2, 3, 4 refer to scenarios S1, S2, S3, and S4.

Figure 6: F1 scores in simulation Example 2. Rows 1, 2, 3 refer to $\tau$ = 0.5, 0.8 and 0.9. Columns 1, 2, 3, 4 refer to scenarios S1, S2, S3, and S4.

In Example 2 as presented in Figures 5 and 6, the overall estimation and variable selection performance is notably enhanced when the prior insights are of high quality or informative. Specifically, the QRIF method demonstrates a notable improvement in variable selection performance, particularly in cases involving partially correct insights and high-quality mixed insights. With informative prior insights, QRIF shows efficacy in identifying important variables more accurately in complex conditions, a common scenario we often encounter when analyzing genetic data. In scenarios of low-quality and erroneous insights, the results indicate that the performance of the QRIF method is comparable to QR, which validates QRIF’s robustness to misspecified prior information. In summary, the simulation study highlights that QRIF method is effective in variable selection in both examples when prior insights are informative and robust to potentially misspecified prior information.

5 Theoretical Properties

Theoretically, we establish the oracle properties for the estimator from our quantile regression with insight fusion (QRIF) method with Huber loss approximation for the quantile check loss and local linear approximation of the SCAD penalty. Again, we omit subscript $\tau$ for different conditional quantiles throughout.

We first establish an oracle inequality for the QRIF estimator with Lasso penalty or $L_{1}$ penalty

$$\widehat{}\mbox{\boldmath$\beta$}^{QRIF0}=\operatorname*{arg\,min}_{\mbox{\boldmath$\beta$}}(1-\xi)h_{\gamma}(\mbox{\boldmath$\beta$})+\xi h_{\gamma}^{p}(\mbox{\boldmath$\beta$})+\lambda\|\mbox{\boldmath$\beta$}\|_{1}.$$

Then, starting from the QRIF estimator with Lasso penalty, we show the QRIF estimator with SCAD penalty $\hat{\bm{\beta}}^{QRIF}$ defined in Equation (9) has the appealing oracle property. Note that the values of $\lambda$ in the above two equations may be different. In the following, $C$ denotes a generic positive constant which may assume different values at different places. The following assumptions are imposed.

• ${\bf x}_{i}$ is sub-Gaussian in the sense that $E\exp\{{\bf x}_{i}^{\top}{\bf a}\}\leq C\exp\{C\|{\bf a}\|^{2}\}$ for any vector ${\bf a}$. The matrix $\mbox{\boldmath$\Sigma$}=E[{\bf x}_{i}{\bf x}_{i}^{\top}]$ has eigenvalues bounded and bounded away from zero.

• The conditional density of $\epsilon_{i}=y_{i}-{\bf x}_{i}^{\top}\mbox{\boldmath$\beta$}_{0}$, $f_{\epsilon}(.|{\bf x})$, is uniformly bounded by a constant $\bar{f}$, and its derivative is uniformly bounded by another constant $\bar{f^{\prime}}$. On any interval $[-C,C]$, $f_{\epsilon}(.|{\bf x})$ is bounded below by a constant $\underline{f}>0$ ($\underline{f}$ depends on $C$ but this dependence is suppressed in the notation for simplicity).

• The true parameter $\mbox{\boldmath$\beta$}_{0}$ is sparse with support $S=\{j:\beta_{0j}\neq 0\}$ and $|S|=s_{n}$.

Theorem 1.

Under assumptions (A1)-(A3), also assuming that $\|\hat{\mbox{\boldmath$\beta$}}^{p}-\mbox{\boldmath$\beta$}_{0}\|\leq c\gamma$ for a sufficiently small constant $c$, with $\lambda\geq 2(1-\xi)\|\nabla h_{\gamma}(\mbox{\boldmath$\beta$}_{0})\|_{\infty}+2\xi\|\nabla h_{\gamma}^{p}(\mbox{\boldmath$\beta$}_{0})\|_{\infty}$, $\gamma>>\sqrt{s_{n}\hbox{log}(d_{n}\vee n)/n}$, $\gamma>>\lambda\sqrt{s_{n}}$, we have with probability at least $1-(d_{n}\vee n)^{-C}$ that the estimator of Equation (5) satisfies the following oracle inequality

$$\|\widehat{}\mbox{\boldmath$\beta$}^{QRIF0}-\mbox{\boldmath$\beta$}_{0}\|\leq C\lambda\sqrt{s_{n}}\;\mbox{ and }\|\widehat{}\mbox{\boldmath$\beta$}^{QRIF0}-\mbox{\boldmath$\beta$}_{0}\|_{1}\leq C\lambda s_{n}.$$

Remark 1.

By Theorem 1, we can take $\lambda=C(1-\xi)\sqrt{{\hbox{log}(d_{n}\vee n)\over n}}+2\xi\|\nabla h_{\gamma}^{p}(\mbox{\boldmath$\beta$}_{0})\|_{\infty}$. For the standard LASSO penalized quantile regression without the $h_{\gamma}^{p}$ term for prior insights, we would take $\lambda=C\sqrt{{\hbox{log}(d_{n}\vee n)\over n}}$ which would also lead to $\|\widehat{}\mbox{\boldmath$\beta$}^{QRIF0}-\mbox{\boldmath$\beta$}_{0}\|\leq C\lambda\sqrt{s_{n}}$ (here $\widehat{}\mbox{\boldmath$\beta$}^{QRIF0}=\widehat{}\mbox{\boldmath$\beta$}^{qLASSO}$ is the quantile LASSO estimator without using $h_{\gamma}^{p}$, equivalent to setting $\xi=0$). Thus if $\|\nabla h_{\gamma}^{p}(\mbox{\boldmath$\beta$}_{0})\|_{\infty}$ is of a smaller order than $\sqrt{{\hbox{log}(d_{n}\vee n)\over n}}$, the quantile regression with insights fusion (QRIF) estimator uses a smaller $\lambda$ and thus the upper bound is smaller. This observation is similar to that made in Jiang et al. (2016).

The proof of Theorem 2 utilizes the following lemmas.

Lemma 1.

Let $c>0$ be a sufficiently small constant. With the same assumptions used in Theorem 1, uniformly over the set $\Omega:=\{(\mbox{\boldmath$\beta$}_{1},\mbox{\boldmath$\beta$}_{2}):\|\mbox{\boldmath$\beta$}_{1}-\mbox{\boldmath$\beta$}_{2}\|\leq c\gamma,\|\mbox{\boldmath$\beta$}_{2}-\mbox{\boldmath$\beta$}_{0}\|\leq c,\|\mbox{\boldmath$\beta$}_{2}-\mbox{\boldmath$\beta$}_{0}\|_{1}\leq 3\sqrt{s_{n}}\|\mbox{\boldmath$\beta$}_{2}-\mbox{\boldmath$\beta$}_{0}\|\}$, we have with probability at least $1-e^{-t}$,

$$h_{\gamma}(\mbox{\boldmath$\beta$}_{1})-h_{\gamma}(\mbox{\boldmath$\beta$}_{2})-\langle\nabla h_{\gamma}(\mbox{\boldmath$\beta$}_{2}),\mbox{\boldmath$\beta$}_{1}-\mbox{\boldmath$\beta$}_{2}\rangle$$

$$\geq \left({\underline{f}\lambda_{\min}(\mbox{\boldmath$\Sigma$})\over 4}-{C\sqrt{s_{n}}\over\gamma}\left(\sqrt{{\gamma\hbox{log}(d_{n})\over n}}+{\hbox{log}(d_{n})\over n}\right)-C\sqrt{{t\over n\gamma}}-C{t\over n\gamma}\right)\|\mbox{\boldmath$\beta$}_{1}-\mbox{\boldmath$\beta$}_{2}\|^{2},$$

where $\lambda_{\min}(\mbox{\boldmath$\Sigma$})$ is the minimum eigenvalue of $\mbox{\boldmath$\Sigma$}=E[{\bf x}_{i}{\bf x}_{i}^{\top}]$.

Lemma 2.

Assume $c$ is a sufficiently small positive constant. With the same assumptions used in Theorem 1, we have that, with probability at least $1-e^{-t}$, uniformly over the set $\Omega^{\prime}=\{(\mbox{\boldmath$\beta$}_{1},\mbox{\boldmath$\beta$}_{2}):\|\mbox{\boldmath$\beta$}^{p}-\mbox{\boldmath$\beta$}_{2}\|\leq c\gamma,\|\mbox{\boldmath$\beta$}_{1}-\mbox{\boldmath$\beta$}_{2}\|\leq c\gamma,\|\mbox{\boldmath$\beta$}_{1}-\mbox{\boldmath$\beta$}_{2}\|_{1}\leq 3\sqrt{s_{n}}\|\mbox{\boldmath$\beta$}_{1}-\mbox{\boldmath$\beta$}_{2}\|\}$, we have

$$h_{\gamma}^{p}(\mbox{\boldmath$\beta$}_{1})-h_{\gamma}^{p}(\mbox{\boldmath$\beta$}_{2})-\langle\nabla h_{\gamma}^{p}(\mbox{\boldmath$\beta$}_{2}),\mbox{\boldmath$\beta$}_{1}-\mbox{\boldmath$\beta$}_{2}\rangle\geq\left({\lambda_{\min}(\mbox{\boldmath$\Sigma$})\over 2}-C\left({1\over\gamma}\sqrt{{s_{n}\hbox{log}(d_{n})\over n}}+\sqrt{{t\over n\gamma^{2}}}+{t\over n\gamma}\right)\right)\|\mbox{\boldmath$\beta$}_{1}-\mbox{\boldmath$\beta$}_{2}\|^{2}.$$

Lemma 3.

With probability at least $1-e^{-t}$,

$$\|\nabla h_{\gamma}(\mbox{\boldmath$\beta$}_{0})\|_{\infty}\leq C\left(\gamma^{2}+\sqrt{{t+\hbox{log}(d_{n})\over n}}+{t+\hbox{log}(d_{n})\over n}\right).$$

We now establish the oracle property for our quantile regression with insight fusion (QRIF) estimator with Huber loss approximation of the quantile check loss and local linear approximation of the SCAD penalty.

Define $\widehat{}\mbox{\boldmath$\beta$}^{ora}=(\widehat{}\mbox{\boldmath$\beta$}^{ora}_{S},{\bf 0})=\operatorname*{arg\,min}_{\mbox{\boldmath$\beta$}:\mbox{\boldmath$\beta$}_{S^{c}}={\bf 0}}(1-\xi)h_{\gamma}(\mbox{\boldmath$\beta$})+\xi h_{\gamma}^{p}(\mbox{\boldmath$\beta$})$. Theorem 2 shows the oracle property that, $\hat{\bm{\beta}}^{QRIF}$ defined in Equation (9), using $\widehat{}\mbox{\boldmath$\beta$}^{QRIF0}$ as the initial estimator, is equal to $\widehat{}\mbox{\boldmath$\beta$}^{ora}$ with probability approaching one. Here we require the smallest signal to be sufficiently large, as is always required for the oracle property to hold in the SCAD-penalized model.

Theorem 2.

Assume $\|\widehat{}\mbox{\boldmath$\beta$}^{ora}-\mbox{\boldmath$\beta$}_{0}\|\leq a_{n}$ and $\|\widehat{}\mbox{\boldmath$\beta$}^{QRIF0}-\mbox{\boldmath$\beta$}_{0}\|\leq b_{n}$ for some positive sequences $a_{n},b_{n}=o(1)$. If $b_{n}<<\lambda<<\min_{j\in S}|\beta_{0j}|$ and $\lambda>(1-\xi)\|\nabla_{\gamma}(\hat{}\mbox{\boldmath$\beta$}^{ora})\|_{\infty}+\xi\|\nabla h_{\gamma}^{p}(\widehat{}\mbox{\boldmath$\beta$}^{ora})\|_{\infty}$, for sufficiently large $n$, we have that the QRIF estimator of Equation (9) is the oracle estimator

$$\widehat{}\mbox{\boldmath$\beta$}^{QRIF}=\widehat{}\mbox{\boldmath$\beta$}^{ora}.$$

Remark 2.

In the statement of Theorem 2, we condition on $\|\widehat{}\mbox{\boldmath$\beta$}^{ora}-\mbox{\boldmath$\beta$}_{0}\|\leq a_{n}$, $\|\widehat{}\mbox{\boldmath$\beta$}^{QRIF0}-\mbox{\boldmath$\beta$}_{0}\|\leq b_{n}$ and the size of $(1-\xi)\|\nabla h_{\gamma}(\widehat{}\mbox{\boldmath$\beta$}^{ora})\|_{\infty}+\xi\|\nabla h_{\gamma}^{p}(\widehat{}\mbox{\boldmath$\beta$}^{ora})\|_{\infty}$, thus the statement is entirely deterministic. The rate of $\|\widehat{}\mbox{\boldmath$\beta$}^{ora}-\mbox{\boldmath$\beta$}_{0}\|$ of course depends on that of $\|\widehat{}\mbox{\boldmath$\beta$}^{p}-\mbox{\boldmath$\beta$}_{0}\|$. If $\widehat{}\mbox{\boldmath$\beta$}^{p}$ is far away from the true $\mbox{\boldmath$\beta$}_{0}$, we cannot expect $\widehat{}\mbox{\boldmath$\beta$}^{ora}$ to be a good estimator unless $\xi=0$. On the other hand, if $\|\widehat{}\mbox{\boldmath$\beta$}^{p}-\mbox{\boldmath$\beta$}_{0}\|=O_{p}(\sqrt{s_{n}/n})$, it can naturally be expected that the oracle estimator also has rate $a_{n}\asymp\sqrt{s_{n}/n}$, and as proved in Theorem 1 we can expect $b_{n}\asymp\sqrt{s_{n}\hbox{log}(d_{n}\vee n)/n}$. Furthermore, Lemma 4 gives a bound for $\|\nabla h_{\gamma}(\widehat{}\mbox{\boldmath$\beta$}^{ora})\|_{\infty}$. On the other hand, $\|\nabla h_{\gamma}^{p}(\widehat{}\mbox{\boldmath$\beta$}^{ora})\|_{\infty}$ depends on how good $\widehat{}\mbox{\boldmath$\beta$}^{p}$ is as an estimator of $\mbox{\boldmath$\beta$}_{0}$.

Lemma 4.

Assume assumptions (A1)-(A3) hold and that $\|\widehat{}\mbox{\boldmath$\beta$}^{ora}-\mbox{\boldmath$\beta$}_{0}\|=O_{p}(a_{n})$, then with probability at least $1-(d_{n}\vee n)^{-C}$,

$$\|\nabla h_{\gamma}(\widehat{}\mbox{\boldmath$\beta$}^{ora})\|_{\infty}=O_{p}\left(\gamma^{2}+a_{n}+\sqrt{{a_{n}s_{n}\hbox{log}(d_{n}\vee n)\over n}}+{s_{n}\hbox{log}^{3/2}(d_{n}\vee n)\over n}\right).$$

The detailed proofs of the main theorems and lemmas are provided in the supplementary materials.

6 Discussion

The escalating challenge of obesity has long been a severe threat to public health. For example, in 2001 the Surgeon General called for action to prevent and decrease overweight and obesity. Responding to the high healthcare cost by obesity, in 2023, the Congressional Budget Office has called for innovative research into obesity, especially new research on the use of anti-obesity medications. Members of Congress have introduced legislation, most recently the Treat and Reduce Obesity Act of 2023 (H.R. 4818 and S. 2407).

We aim to simultaneously identify and estimate important SNPs for high quantiles of BMI of the most interest. Leveraging ultra-high dimensional data available with over 500,000 genetic risk factors and various phenotypes, our study seeks to explore the genetic mechanism behind obesity, while incorporating previous research insights. Our analysis provides a comprehensive picture of important genetic risk factors for different quantiles of BMI.

We discovered some new interesting SNPs plausibly associated with obesity, providing possible direction for future research. Many of the identified SNPs are also validated by the literature as well as lab-based biological research. We further demonstrate the performance of our proposed quantile regression with insight fusion (QRIF) approach through numerical simulations. In our simulation mimicking real genetic data, the QRIF utilizing information from established studies outperforms the traditional penalized quantile regression and GWAS and shows robustness to potential misspecification.

Our analysis indicates that our proposed QRIF approach works successfully to identify and estimate important variables for different quantiles of BMI, while adopting flexible formats of prior information. By incorporating prior insights into obesity, our study contributes to possible personalized obesity treatment strategies and identifying high-risk individuals for preventative intervention. Furthermore, the QRIF approach, complemented by our computational algorithms, potentially serves as a useful tool for exploring other conditions where specific quantiles of phenotypes are of most interest, including but not limited to blood pressure and glucose levels.