Bayesian hierarchical quantile regression with application to characterizing the immune architecture of lung cancer

The successful development and implementation of precision immuno‐oncology therapies requires a deeper understanding of the immune architecture at a patient level. T‐cell receptor (TCR) repertoire sequencing is a relatively new technology that enables monitoring of T‐cells, a subset of immune cells that play a central role in modulating immune response. These immunologic relationships are complex and are governed by various distributional aspects of an individual patient's tumor profile. We propose Bayesian QUANTIle regression for hierarchical COvariates (QUANTICO) that allows simultaneous modeling of hierarchical relationships between multilevel covariates, conducts explicit variable selection, estimates quantile and patient‐specific coefficient effects, to induce individualized inference. We show QUANTICO outperforms existing approaches in multiple simulation scenarios. We demonstrate the utility of QUANTICO to investigate the effect of TCR variables on immune response in a cohort of lung cancer patients. At population level, our analyses reveal the mechanistic role of T‐cell proportion on the immune cell abundance, with tumor mutation burden as an important factor modulating this relationship. At a patient level, we find several outlier patients based on their quantile‐specific coefficient functions, who have higher mutational rates and different smoking history.


INTRODUCTION
Immunotherapy is a class of cancer treatments that fosters the patient's own immune system to fight cancer (Waldman et al., 2020). Although immunotherapies represent a major breakthrough in cancer treatment, offering immense clinical benefit to some patients with a lower toxicity burden than chemotherapy, resistance to immunotherapy remains a major challenge (Walsh & Soo, 2020). Tumors use various strategies to protect themselves from antitumor immunity which might vary across patients. Antitumor immune responses might also be mediated by several different mechanisms, driven by patient-specific immune architecture. Cancer immunotherapy therefore needs to be personalized to recognize patient-specific rate-limiting steps and employ strategies for overcoming these hurdles (Kakimi et al., 2017).
Identifying patient-specific influences on immune response requires integrating measures of immune activity and characterizing their dependence on upstream factors. Here, we consider a setting where the response variable is a continuous measure of immune activation, such as the abundance of CD8+ T-cells, which play a key role in directly killing cancer cells. The abundance of these cells is driven by "Level 1" covariates 1 , … , , which represent measures of upstream immune activity. Here we take these Level 1 covariates to be features derived from T-cell receptor (TCR) sequencing, which capture activity of the adaptive immune system. The TCR repertoire depends, in turn on "Level 2" influences further upstream, such as antigens that the T-cells have been exposed to. This includes tumor-specific antigens, which are abnormal proteins produced by tumor cells due to DNA mutations in the cell (see Figure 1A). We consider mutational variables 1 , … , as (hierarchical) Level 2 covariates in our model. Our goal is to develop a hierarchical modeling approach providing insights into the mechanistic relationships among these variables wherein the measures of immune activity may depend on the upstream factors in a complex nonlinear fashion.
We now briefly review prior work addressing the challenge of flexible regression modeling, which lays the foundation for the proposed approach. In order to estimate the subject-specific effect of variables on the outcome variable, Hastie and Tibshirani (1993) proposed the varying coefficient model (VCM). Since then, several variations of the VCM have been proposed (Fan & Zhang, 1999;Park et al., 2013), including approaches that incorporate shrinkage in estimation (Wang & Xia, 2009). However, existing VCM methods do not enable explicit multilevel variable selection, which is crucial in settings such as ours where there are a large number of explanatory variables within a hierarchy. Although VCMs allow more flexible coefficients than traditional regression models, they are still focused on estimation of the mean of the response variable. If one is interested in obtaining a comprehensive picture of the effect of the predictors on the response variable, mean regression might be insufficient. For example, if the dependent variable is multimodal or skewed, estimating the mean effect might be misleading. In such scenarios, median or more generally, quantile regression may be more appropriate (Koenkar & Bassett, 1978). Specifically, if the interest lies specifically at higher or lower quantiles of the response variable, quantile regression is more suitable.
Over the last couple of decades, methodological developments in quantile regression have been proposed in both the classical and Bayesian frameworks (Das & Ghosal, 2018;Kottas & Gelfand, 2001;Reich, 2012;Yu & Moyeed, 2001). Several articles emerged in the literature integrating quantile regression with VCMs (Honda, 2004;Kim, 2007;Tang & Wang, 2005). In Bayesian settings, the literature on VCMs is relatively sparse (Biller & Fahrmeir, 2001), and these proposals do not address quantile regression. Recently, Ni et al. (2018) proposed a VCM incorporating variable selection in the Bayesian framework, which was applied to characterize the relationship of cancer patient outcomes to proteomics and genomics data. Although there are approaches for Bayesian variable selection in multilevel models that assume linear covariate effects (Koslovsky et al., 2020;Stingo et al., 2013), to the best of our knowledge, no prior work has considered variable selection in VCM for quantile regression.
In order to understand the subject-specific effect of hierarchically structured covariates on the outcome variable we propose Bayesian QUANTIle regression for hierarchical COvariates (QUANTICO), where the regression coefficients are allowed to differ across patients for any given quantile level of the outcome variable. Among the hierarchically structured (multilevel) covariates, we consider the Level 1 covariates have a direct effect on the response variable, modulated by Level 2 covariates. Since we expect the covariate effects to be heterogeneous across patients, we want to perform individualized inference as well as allowing variable selection for both Level 1 and Level 2 covariates. As we also expect the effects to be different at different parts of the distribution of the response variable, we model across the quantiles, rather than only considering fixed moments (e.g., mean). This provides a richer, broader, and more flexible exploration of the relationship structure. We show an illustration of the proposed model in Figure 1B with subjects, Level 1 and Level 2 covariates. As shown in the figure, the selection of Level 1 and 2 covariates is allowed to vary across quantiles. Due

F I G U R E 1 (A) Illustration of how mutation within tumor cells and T-cell receptor (TCR) repertoire impact immune cells. (B)
Illustration of the proposed model for a scenario with patients, mutation variables 1 , … , , and TCR variables 1 , … , . The response variable is denoted 1 , … , . Different colored lines describe the estimated dependency structure at different quantile levels denoted by ( 1 , 2 , … ). (C) Directed acyclic graph (DAG) of the QUANTICO model. Parameters are shown in circles and the observed data are shown in boxes to the presence of two levels of covariates and due to the fact that the effect of Level 2 covariates is induced on the output variable via its effect on Level 1 covariates, we call it a hierarchical model. The rest of the paper is organized as follows: Section 2 describes the QUANTICO modeling framework, with a discussion of the priors in Section 3. We then describe the computational algorithm for performing posterior inference (Section 4), and provide a simulation study comparing the performance of the proposed method with existing alternatives (Section 5). In Section 6, we apply QUANTICO to characterize the relationship of the CD8 immune marker with TCR and mutation variables for lung cancer patients. We conclude with a brief discussion and possible extensions of our methodology in Section 7.

QUANTICO MODEL
In Section 2.1, we introduce quantile regression in a VCM setting with two levels of covariates. In Section 2.2, we describe variable selection procedures on Level 1 and Level 2 covariates. Section 2.3 summarizes the likelihood constructions.

Varying sparsity quantile regression model
Suppose there are subjects whose response variables (immune marker values) are denoted by = ( 1 , … , ) and Level 1 covariates (TCR variables) = ( 1 , … , ). For the th subject, the set of Level 1 covariates are denoted as = ( 1 , … , ) for = 1, … , . We assume a linear relationship between the response variable and the Level 1 covariates . Now, for the th Level 1 covariate , consider a set of Level 2 covariates (mutation variables) = ( 1 , … , ) for = 1, … , . A Level 2 covariate induces its effect on the response variable through its effect on the th Level 1 covariate, that is, for = 1, … , . Note that the same Level 2 covariate may influence multiple Level 1 covariates, and that the number of Level 2 covariates can be different for each Level 1 covariate. However, in the particular case study considered in this paper, the same set of Level 2 covariates (mutation variables) are considered for all Level 1 covariates (TCR variables). For the th subject, the set of Level 2 covariates is denoted by = ( 1 , … , ) for = 1, … , , = 1, … , . The total number of possible effects for the variables which can be identified as a part of the coefficient of the variables is given by ∑

=1
. Instead of estimating the effect of and on the mean of the response variable, our interest lies in estimating their relationship for different quantile levels. To obtain the quantile-specific estimates, we assume the effect of the covariates on to be dependent on the quantile level (0 < < 1). Let denote the th conditional quantile (0 < < 1) of the response variable at = , = . The relation between the response and the covariates at the th quantile is given by where 0 denotes the values of all distinct Level 2 covariates for the th subject, and denotes the values of the Level 2 covariates corresponding to the th Level 1 covariate for th subject. This equation shows the dependence structure of the th quantile of the dependent variable for the th patient on the corresponding Level 1 ( ) and Level 2 ( ) covariates. Note that for each Level 1 covariate we have many distinct variables, but since two different Level 1 covariates may share the same variables, the ∑ =1 many variables may not be distinct. Also note that both the intercept and slope terms are patient-specific (indexed by ). In order to incorporate the effect of the variables on the dependent variable, at any given quantile level , all the slope terms (of ) are semiparametric functions of the variables and likewise the intercept is a semiparametric function of the distinct variables. The explicit structure of ( , ) is discussed in Section 2.2.

Varying-sparsity coefficient modeling and selection
For any given quantile level , we consider ( , ) as a smooth function of . The main motivation for taking ( , ) as a smooth function of is so that the coefficients of any Level 1 covariate ( ) for two different subjects are similar if the values of the corresponding Level 2 covariates ( ⋅ ) are similar as well. Since for any Level 1 covariate, "neighboring" patients (with respect to the variables) are expected to have a similar slope coefficient, this assumption enables borrowing of strength, increasing our power to estimate the subject-specific slope coefficients. For modeling the slope and intercept terms, we use spline functions due to their flexible construction, interpretation, and the ease of incorporating penalization.
At any given quantile (0 < < 1), the slope and the intercept terms are estimated as the sum of spline functions given by where 0 ( , 0 ) denotes the (global) intercept term, and ( , ) for ≥ 1 denotes the slope term. The spline components ( ) ( ) = where denotes the cubic B-spline bases for and denotes the corresponding spline coefficient. Note that we consider the intercept term to be a function of all Level 2 covariates. The number of knots for the B-spline bases is taken to be sufficiently large to capture the nonlinear features. We do not perform knot selection; rather equally spaced quantile knots on each of the Level 2 covariates are considered.
Smoothing is induced via regularization and overfitting is controlled through a roughness penalty on the spline coefficients, details of which are provided in Section A of the Supporting Information.
In Equation (3), ( , ⋅) is modeled as a sum of a set of smooth spline functions. As discussed in Section 3, we construct a prior on the spline coefficients, so that, for any given variable, the linear and nonlinear effects corresponding to all associated variables can be identified. Under these assumptions, the estimated coefficients of the variables would be zero if both the linear and nonlinear effects of all associated variables are zero. However, for a larger number of variables, it is unlikely that both the linear and nonlinear effects of all variables corresponding to any variable are zero. If the number of Level 1 covariates ( variables) is also large, it becomes crucial to perform variable selection on them to enforce sparsity and interpretability.
In order to incorporate sparsity on the Level 1 covariates, one naïve approach could be to use discrete mixture priors such as spike-and-slab (George & McCulloch, 1993). But, since we consider the selection of the covariates to be patient-specific, to apply spike-and-slab prior, we would need to assign a latent indicator for each coefficient for every patient. This approach would therefore substantially increase the number of parameters in the model. In addition, as discussed below, the spike-and-slab prior is not well-suited for functional regression coefficients. Instead, we rely on a Bayesian hard-thresholding approach where we truncate the coefficients with smaller absolute values to zeros so that only the important Level 1 variables are selected. For truncation, we take the Bayesian hardthresholding function ℎ( , ) = (| | > ) to threshold the slope coefficients of the variables. We modify the values of the slope coefficients given in Equation (3) as where is a "minimum effect size" for the effect of to be considered as nonzero. Our use of the hard-thresholding prior instead of the more commonly adopted Bayesian spike-and-slab prior is due to the functional nature of the regression coefficients (⋅). The hard-thresholding prior allows selection at any given input whereas spikeand-slab is not able to handle an infinite-dimensional object. Another advantage of using the Bayesian hardthresholding is that this "minimum effect size" can be set at any reasonable value by the user based on intuition or prior experience, or estimated by assigning a prior (as we discuss in the next section). Note that no hard-thresholding is performed (or needed) on the intercept term as our main interest is to perform variable selection among Level 1 covariates only.

PRIOR FORMULATIONS
Note that the intercept and the slope terms ( ( ) 0 ( ) and ( ) ( ), respectively) are functions of the spline coefficients . In this section, we describe the prior on the spline coefficients , the prior used to induce selection of Level 2 covariates, and the prior assigned to the thresholding parameter to select Level 1 covariates.

Prior on spline coefficients
We propose the use of a penalized spline in order to have a flexible but smooth fit. Specifically, we choose a large number of knots placed at equally spaced quantiles of the covariates so that the local features can be captured, and penalize the roughness of (⋅) through an improper Gaussian random walk prior on given by ∼ ( , − ). The penalty matrix is constructed from the second-order differences of adjacent spline coefficients, which essentially penalizes the second derivatives of (⋅). Larger value of leads to a smoother fit, while smaller value of leads to an irregular fit (Ni et al., 2015). Note that Equation (3) is not identifiable since adding some quantity to any term of the summation and deducting the same quantity from any other term would yield the same summation value. In addition, is singular and therefore no penalty is imposed on the linear and constant trend of (⋅) (i.e., the null space of ). In order to alleviate these two issues, we consider a similar approach as in Scheipl et al. (2012) and transform the spline bases into orthonormal bases. Let = ( 1 , … , ). For a given quantile , consider the spectral decomposition of the covariance of where is the orthonormal matrix of the eigenvectors corresponding to the positive eigenvalues in the diagonal matrix . Let * = 1 2 . Now if we assume an independent proper prior * ∼ ( , 2 ), then the nonlinear part of (⋅) can be parameterized by * * which has a proper distribution proportional to the distribution of the original improper prior . Suppose we denote the effect size of the th covariate for the th subject before applying the hard-thresholding by * ( , ). Thus the full reparameterization of * ( , ) is now given by * ( , where is the global constant term absorbing all constant terms from splines. This parameterization adds the flexibility to separately shrink and estimate the linear, nonlinear, and constant effects of the Level 2 ( ) variables on the coefficients of the Level 1 ( ) variables. In order to make the proposed method computationally more efficient, we only consider the first several eigenvectors which explain at least 99.5% of the variation, a similar idea as in principal component analysis. In order to induce sparsity on the linear, nonlinear, and the constant effects, we consider a parameter-expanded normal mixture of inverse gamma (peNMIG) prior on each * , 0 , and separately. We opt for the peNMIG prior since it is known to provide a more efficient Markov chain Monte Carlo (MCMC) algorithm compared to traditional spike-and-slab priors given the multivariate nature of the spline coefficients (Scheipl et al., 2012). peNMIG multiplicatively expands * as * = , where is a scalar parameter and is a vector of the same length as * . Each of 0 and is also expanded in a similar fashion. A brief discussion on the choice of such a prior is provided in Section A of the Supporting Information.

Priors for selection of Level 2 covariates
Since the Level 2 coefficients are at population level, we can induce explicit selection using a spike-and-slab prior on , where ∼ ( ) and 0 is a fixed very small quantity close to zero. The selection of as a nonzero effect is indicated by the binary variable , and thus indicates the selection of * vector due to the multiplicative construction. In terms of interpretation, the binary variable = 1 indicates that has nonzero nonlinear effect on . Similarly, in the expansion of 0 and , = 1 indicates the presence of linear and constant effects of on , respectively. Thus, based on the estimated values of in the expansion of * and 0 , we can identify the presence of nonlinear and linear effects of Level 2 covariates, respectively. In essence, this construction allows the flexibility of considering only linear or only nonlinear or joint effects simultaneously. We choose conjugate hyperpriors for and , ∼ ( , ) and ∼ ( , ). As the number of Level 2 covariates increases, this Beta-Bernoulli prior automatically corrects for multiplicity by making the posterior distribution of concentrated at small values near 0 (Scott & Berger, 2010). We assign a mixture normal prior on each element of the vector = ( ( ) ), (1, 1) + 1 2 (−1, 1). The structure of this assumed prior discourages small effects. In a similar fashion, peNMIG priors are assumed for 0 and as well.

Prior on hard-thresholding for selection of Level 1 covariates
As mentioned before, to select the Level 1 covariates, we adopt a hard-thresholding approach where we truncate the nonzero coefficient of the th Level 1 covariate with absolute value less than to 0. Since sparsity is induced while estimating the linear, nonlinear, and constant effects of the Level 2 covariates on the Level 1 covariates, it is possible that (⋅) = 0 even before hard-thresholding (when * , 0 , and are zero for all = 1, … , ). In that case, can take any value without affecting the resulting estimates.
To resolve this identifiability issue, we take = for = 1, … , . In the presence of at least one nonzero (⋅), is now well-defined. We put a gamma prior on , ∼ ( , ). The values of the shape and scale parameters ( , ) can be taken so that the mean of the gamma distribution (i.e., ) is equal to the desired cutoff (based on intuition or prior experience). A brief discussion on how to choose the parameters of the gamma prior is provided in Supporting Information (Section A). Note that a more conventional variable selection prior, the spike-and-slab, is often used for finite-dimensional parametric models. We choose to use the random hard-thresholding because it is more suitable for infinite-dimensional functional selection, as discussed in Section 2.2. A schematic illustration of the proposed model and its key parameters is given in Figure 1C.

MCMC AND POSTERIOR INFERENCE
The posterior distributions of the model parameters are not analytically tractable. Therefore, an MCMC sampling algorithm is required to generate samples from the posterior distribution. We use a Gibbs sampling scheme for updating parameters using their full conditional distributions ( , , ) and Metropolis sampling scheme for the parameters without closed-form conditional distributions ( , , ). Details on the full conditional distributions are provided in Section B of the Supporting Information.

Inferential summaries
When using this sampling algorithm, the quantile level of interest is fixed at the desired level. In order to estimate multiple quantiles simultaneously, the algorithm can be run in parallel for faster computation. The selection of Level 1 and Level 2 covariates can be performed in several ways based on the marginal posterior probability of inclusion. For the selection of Level 1 covariates, we consider the cutoff of the posterior probability of inclusion to be 0.5. A similar approach of thresholding the posterior inclusion probabilities can be taken for the selection of Level 2 covariates. In the presence of a higher number of Level 2 covariates, a false discovery rate (FDR) controlling approach can also be considered (Baladandayuthapani et al., 2010); see Section C of the Supporting Information.
The proposed model allows the identification of both linear and nonlinear effect components of within the effect size of for all subjects. From the posterior samples, the posterior probabilities of both linear and nonlinear effect of Level 2 covariates within the effect size of Level 1 covariate can be calculated explicitly. Due to the induction of sparsity on the effect-size components of Level 2 covariates, the effect size of each Level 2 covariate can be categorized as one of four possible cases: linear, nonlinear, both, or none.
Based on the posterior mean estimate of the effect sizes of the Level 1 covariates over a grid of quantiles, patientspecific posterior credible intervals over the quantiles can be obtained using QUANTICO. To calculate the posterior 95% credible interval of the quantile function for the th patient, we calculate the posterior mean of ( | ⋅ , ⋅ ) over the quantile grid = 0.1 for = 1, … , 9. Then, for quantile level , we calculate the 95th percentile of the values | ( ) ( | ⋅ , ⋅ ) −ˆ( | ⋅ , ⋅ )| from the posterior sample, where ( ) ( | ⋅ , ⋅ ) denotes the value of ( | ⋅ , ⋅ ) obtained in the th posterior sample. Thus, we derive the width of the posterior 95% credible interval at all 's. Furthermore, by taking a dense grid of quantiles over [0,1], patient-specific uniform posterior credible intervals can be calculated. A detailed description on the calculation of posterior pointwise and uniform credible intervals is provided in Supporting Information Section D.

SIMULATION STUDIES
In this section, we evaluate the variable selection performance of QUANTICO across both Level 1 and Level 2 covariates and illustrate the subject-specific estimation at different quantile levels using simulation studies. The performance of QUANTICO is compared with varying coefficient quantile regression model (VCQRM) and variable selection in quantile regression using the lasso penalty (LASSO-QR) (Wu & Liu, 2009) as well as its Bayesian alternatives. All of these methods are variants of quantile-based VCMs and we provide explicit details of each method in the Section C of the Supporting Information.

Variable selection performance
To compare the variable selection performance across both Level 1 ( ) and Level 2 ( ) covariates, we calculate the true positive rate (TPR), false positive rate (FPR), and area under receiver operating characteristic (ROC) curve (AUC) separately for the and variables. A detailed description of the simulation design and computation details of the metrics is provided in Supporting Information Section C.

Simulation structure
To compare the performance of QUANTICO, VCQRM, and LASSO-QR, we consider two scenarios with sample sizes of = 100 and = 200. In order to understand how the selection of the variables is carried out in the presence of a higher number of covariates, we consider the cases = 5, 10 for sample size scenario = 100 and the cases = 5, 10, 20 for the sample size scenario = 200. For the cases where > 5, the true model remains as given by Equation (1) of the Supporting Information. So, the additional covariates considered in > 5 scenarios have no true effect on . The simulation is repeated 25 times, and each time new data are generated from the quantile function given in Equation (1) of the Supporting Information. The mean and the standard deviation of the TPR, FPR, and AUC for the and variables are computed separately at = 0.1, 0.5, 0.9 for all the methods. In addition, for a few high-dimensional scenarios the performance of QUANTICO is evaluated in Section F of the Supporting Information. We observe that QUANTICO maintains a high TPR and AUC for larger values of and as well, along with low FPR rates. We also evaluate the performance of QUANTICO compared to two other existing quantile regression methods, implemented in the R packages rqPen (Sherwood & Maidman, 2019) and Brq (Alhamzawi & Ali, 2020) which is also provided in Section F of the Supporting Information. It is observed that QUANTICO outperforms rqPen and Brq, in general.

Comparative performance evaluation
The comparative performance of the three methods is reported in Table 1. In general, QUANTICO and VCQRM perform better than LASSO-QR. QUANTICO results in better selection of the (Level 1) variables, and a large improvement in terms of FPR and AUC over VCQRM. In VCQRM, we do not incorporate any thresholding of the slope terms. Although it is possible to have a zero slope or intercept term without thresholding if all estimated linear, TA B L E 1 The result for the method(s) with the best performance is marked in bold, comparative performance study of QUANTICO, varying coefficient quantile regression model (VCQRM) and LASSO quantile regression (LASSO-QR) methods at quantile levels = 0.1, 0.5, 0.9 for different sample size and number of variable scenarios. lvl1TPR, lvl1FPR, and lvl1AUC denote the true positive rate, false positive rate, and area under receiver operating characteristic curve (AUC) of Level 1 covariates; and lvl2TPR, lvl2FPR, and lvl2AUC denote the same for Level 2 covariates ( , ) Measures = .
= . nonlinear, and constant effects of all Level 2 covariates are zero, in the simulation results, the FPR of Level 1 covariates came out to be 1. We do not observe any strong pattern of differences in the comparative performance of QUANTICO and VCQRM in terms of TPR, FPR, and AUC for (Level 2) variables, since they follow same mechanism for selection of variables. In terms of performance, in Table 1, we report that QUANTICO yields very low TPR as well as FPR for Level 2 covariates for sample size scenario = 100 despite a high AUC, which implies that QUANTICO is very conservative in its estimation. In QUANTICO, the selection of Level 2 covariates is performed using an FDR-based approach with cutoff = 0.2. The performance of QUANTICO improves as increases (as noted in Table 1).

QUANTICO
We further report the selection and estimation performance of QUANTICO evaluated at quantile levels = 0.1, 0.2, … , 0.9. Figure 2A illustrates the selection of variables. QUANTICO selected the correct variables across all quantiles. Specifically, it selected variable 3 for the quantiles greater than 0.5 as it is in the true model. In Figure 2B,C, the corresponding true variables have the highest posterior probability of selection for 1 and 2 , respectively. In Figure 2D-F, the true and estimated quantile functions are plotted for three randomly selected subjects along with the corresponding estimated 95% credible bounds.
We also compute the uniform posterior 95% credible intervals of the quantile functions. To improve coverage of the uniform credible intervals, we increase the width using an inflation factor. Although for point estimates, the observed coverage might come close to the percentage of the computed credible interval, for functions, in practice, it is a common phenomenon that the observed coverage may be less than the actual percentage of the credible interval. This undercoverage of the uniform posterior credible band of smooth functions is a well-known property and has been addressed in several articles (Cox, 1993;Das & Ghosal, 2017;Knapik et al., 2011;Szabo et al., 2015). In order to improve coverage, either undersmoothing or inflation of the obtained credible interval is required (Yoo & Ghosal, 2016). To improve coverage of the posterior uniform credible intervals, we increase the width using an inflation factor. As mentioned in Remark 5.4 and Theorem 5.3 in Yoo and Ghosal (2016), the inflation factor of the uniform credible interval should be taken such that it slowly increases to infinity as a function of sample size. Through experimentation, we observe that the inflation factor ( ) = 1.5 √ log( ) (which has a similar form to that considered in Das and Ghosal (2017) in the context of uniform credible bound over quantiles) works well in simulation under various settings for the sample size and number of Level 1 covariates. A detailed discussion regarding the computation of the uniform posterior 95% credible intervals of the quantile functions along with an extensive study on the computational time is provided in Supporting Information Section D.

Scientific problem and data description
Lung cancer is the second most common cancer and is one of the leading causes of cancer death. Though early stage patients can be treated with surgery, late stage lung cancer requires the use of systemic therapies (Ettinger et al., 2017). In recent years, immunotherapies have emerged as a successful treatment in a subset of late stage lung cancers, largely through their ability to boost the activity of T-cells, the subset of immune cells that target infected or malignant cells based on the detection of specific antigens. However, a large proportion of lung cancer patients still do not respond to immunotherapy (Doroshow et al., 2019).
The nature of the antigens specifically recognized by Tcells has been the object of intense focus, and recent work has suggested that somatic mutations harbored by the tumor can be presented to T-cells as neoantigens (Lee et al., 2018). Analysis of specific somatic mutations and overall tumor mutational burden (TMB) has shown a positive association with response to immunotherapy in patients with NSCLC. This supports the role of these mutations in aiding T-cell responses by increasing tumor immunogenicity. Recent technologies have emerged to sequence the TCR to gain insight into T-cell responses, and studies have confirmed that TCR sequencing can be used to monitor immune response in various types of cancer (Page et al., 2016). In order to develop patient-specific immunotherapeutic treatment strategies for NSCLC, it is of critical importance to understand the interplay between somatic mutations, TCR variables, and the immune microenvironment, illustrated schematically in Figure 1A.
We consider a cohort of 215 NSCLC patients recruited at UT MD Anderson Cancer Center. The tumors of these subjects were analyzed to obtain immune profiling, TCR sequencing, and mutational status of immune-related genes. We focus here on understanding the impact of mutation and TCR variables on the CD8 marker, which is the outcome variable in our analysis and is discussed in more detail below. In order to assess the effect of the explanatory variables on the outcome variable across different patients, we would like to estimate the patientspecific effect at different quantiles of the response.
As our Level 1 covariates, we take standard summary measures of TCR sequencing data, including T-cell clonality, T-cell entropy, T-cell productive proportion of cells, and T-cell richness. T-cell clonality is a measure of heterogeneity among T-cells and has been linked to patient outcomes . T-cell entropy is Shannon's entropy, and highly diverse samples have comparatively higher entropy. The T-cell productive proportion of cells denotes the proportion of the tumor that consists of T-cells. TCR richness is another measure of TCR heterogeneity, defined as the number of different unique sequences in the sample.
As our Level 2 covariates, we consider mutation counts for the top six most frequently mutated genes across the cohort, namely, CSMD3, MUC16, RYR2, TP53, TTN, and USH2A, along with total TMB, which is the total number of mutations observed per sample. As the mutation counts for individual genes have sparse values (i.e., zeros for a large proportion of patients), we only consider the linear and constant effects for those six variables. For TMB, linear, nonlinear, and constant effects are allowed. As our response variable, we focus on CD8 abundance as a key measure of immune activity. CD8 is a protein found on the surface of cytotoxic T-cells. CD8+ T-cells have the ability to mount a response against pathogens and defend against tumors by killing transformed tumor cells (Berg & Forman, 2006), and are therefore a vital part of cancer immunity. In precision oncology, understanding the patient-specific effect of T-cell architecture on CD8+ T-cell abundance is therefore crucial.
Out of the cohort of 215 NSCLC patients, we focus on the 87 patients for which all three data types (immune, TCR, and mutational profiling) are available. Since the set of TCR variables considered have differing magnitudes, it is crucial to transform them to a similar range of values to make any comparison of their effect sizes meaningful. The same applies for the mutation variables. Therefore, before applying QUANTICO, we transform the TCR variables and the CD8 immune markers to unit intervals using log-normal cumulative distribution function (cdf) transformation (see Supporting Information Section G). The mutation variables are transformed into the unit interval using a linear transformation.

Population level findings
Using the proposed method, the coefficients of the TCR variables are estimated for each patient at nine quantile levels = 0.1, 0.2, … , 0.9. We use the same values of the hyperparameters as in the simulation study except the total number of iterations and burn-in, which are taken to be 50,000 and 10,000, respectively. The average posterior probability of selection of all TCR variables is plotted in Figure 3A. The T-cell productive proportion variable has an average posterior probability greater than 0.5 at all quantile levels except at = 0.9, and T-cell entropy has an average posterior probability marginally higher than 0.5 at = 0.2. In general, the average posterior probability of the TCR variables having a nonzero effect on the CD8+ immune cell abundance decreases at higher quantile levels. This implies that for patients with a lower abundance of CD8+ cells, the number of CD8+ cells has a stronger dependence on the TCR variable measures. However, in patients with a higher density of CD8+ cells, this dependence is less prominent. To summarize our Level 2 findings, we assessed the effect of the mutation variables on the coefficient of the Tcell productive proportion, which was identified to be the most important among the TCR variables ( Figure 3B). The nonzero effects of the mutation variables in the coefficient of the T-cell productive proportion are not strong for lower quantile levels, while TMB is shown to have a marginally higher posterior probability of having a nonzero effect at higher quantiles. Thus, at higher quantiles, a large proportion of the effect of the T-cell productive proportion on CD8+ immune cells is due to TMB. This is consistent with previous studies that have shown a positive correlation between TMB and CD8+ in melanoma , where immunotherapy using checkpoint inhibitors have shown to be successful. Hence, our findings suggest that TMB could be used for predicting the response to anti-PD-1/PD-L1 therapies in NSCLC.

Patient-level findings
Our results also provide insights into patient-specific immune profiles offering the potential to guide the development of precision immuno-oncology treatment strategies based on patient-specific information such as patient smoking history and mutational profiles. As an illustration of this, in Figure 3C, we show the estimated coefficients of the T-cell productive proportions for each subject across quantiles in the rows of a heatmap, along with patientlevel covariates. In this figure, we rely on hierarchical clustering over the estimated coefficients at quantile levels = 0.1, … , 0.8 to group patients with similar coefficients. Focusing on the last few rows of the heatmap, it is apparent that patients with overall lower values of the coefficient of the T-cell productive proportion have higher recurrence rates of NSCLC. Layering this analysis with clinical covariates reveals that the effect size of T-cell productive proportion is in general lower for nonsmokers compared to recent and former smokers ( Figure S3). Reuben et al. (2020) found higher T-cell clonality in current and former smokers compared to never smokers.
In Figure 4A, we show a boxplot of the coefficients of the T-cell productive proportion for all patients, and we detect outlier patients (patient numbers 4, 45, 51, 56, and 76) at quantile level 0.1 and an overlapping set of outliers (patient numbers 40, 45, 56, and 76) at quantile level 0.3. Since this coefficient is modeled as a function of the mutation variables, we further compare the number of mutations for these six distinct patients ( Figures 4B-D). In general, these outlier patients have a higher number of mutations compared to the average value in the patient cohort; specifically for TMB, TP53, and MUC16 genes.

CONCLUSIONS AND FUTURE WORK
In this paper, we propose QUANTICO with multilevel covariates where, at any specific quantile level, selection over the direct (Level 1) covariates is performed for each subject. A novel feature of the proposed model is the development of a quantile-specific varying sparsity coefficient estimation approach which allows us to explicitly delineate how at different quantiles, the response variable depends on different covariates for each subject. The proposed method also enables selection of the indirect (Level 2) covariates, and can be used for obtaining patientspecific posterior credible bands over the quantile levels.
The proposed method is used to analyze how the CD8 immune marker depends on TCR and mutation variables at different quantiles. We find that T-cell productive proportion is the most important TCR variable, and influences CD8 immune cells for most quantile levels. Out of all mutation variables considered, total TMB is found to be most important. Based on the structure of the relationship between the TCR and immune variable, we identify outlier patients, who turn out to have a higher number of mutations across several critical genes of known clinical relevance in cancer. This information is potentially useful to devise effective immunologic therapies for such patients(s) based on their unique immune architectures.
There are several potential refinements that could be made to our modeling framework. We can extend our approach to simultaneous quantile regression, where instead of estimating the model parameters at specific quantiles, the entire quantile function and associated parameters are estimated simultaneously (Das & Ghosal, 2018;Yang & Tokdar, 2017). In terms of theoretical excursions, one could investigate for such hierarchical VCQRMs, building on some of the theoretical results proposed for Bayesian quantile regression by ourselves and others (Das & Ghosal, 2017;Yang & Tokdar, 2017) regarding posterior consistency, rates of convergence and posterior contraction. Given the nontrivial nature of these explorations, we leave them as future work.

A C K N O W L E D G M E N T S
VB was supported by NIH grants R01-CA160736, R01CA244845-01A1, R21-CA220299, and P30 CA46592, NSF Grant 1463233, and start-up funds from the U-M Rogel Cancer Center and School of Public Health. CBP was supported by NIH/NCI CCSG P30CA016672 and CPRIT Grant RP150521. KAD was partially supported by a CCSG NCI Grant P30 CA016672, NIH Grants UL1TR003167, 5R01GM122775, the prostate cancer SPORE P50 CA140388, CPRIT Grant RP160693, and the Moon Shots funding at MD Anderson Cancer Center. YN was partially supported by the NSF DMS-1918851 and NSF DMS-2112943.

O P E N R E S E A R C H B A D G E S
This article has earned an Open Materials badge for making publicly available the components of the research methodology needed to reproduce the reported procedure and analysis. All materials are available at http://re3data. org/.

D ATA AVA I L A B I L I T Y S TAT E M E N T Code and data:
Code for the simulations and real data analysis, along with a synthetic data set designed to closely resemble our example T-cell receptor (TCR) data, are available online at https://github.com/bayesrx/QUANTICO.