High-dimensional sparse vine copula regression with application to genomic prediction

High-dimensional data sets are often available in genome-enabled predictions. Such data sets include nonlinear relationships with complex dependence structures. For such situations, vine copula based (quantile) regression is an important tool. However, the current vine copula based regression approaches do not scale up to high and ultra-high dimensions. To perform high-dimensional sparse vine copula based regression, we propose two methods. First, we show their superiority regarding computational complexity over the existing methods. Second, we define relevant, irrelevant, and redundant explanatory variables for quantile regression. Then we show our method's power in selecting relevant variables and prediction accuracy in high-dimensional sparse data sets via simulation studies. Next, we apply the proposed methods to the high-dimensional real data, aiming at the genomic prediction of maize traits. Some data-processing and feature extraction steps for the real data are further discussed. Finally, we show the advantage of our methods over linear models in the real data application.


Introduction
Genomic prediction (GP) aims at predicting a breeding value using genotypic measurements. Then an unobserved trait can be predicted using its genotype information like single-nucleotide polymorphism (SNP). With rapid developments in genomic technologies, researchers have high-dimensional SNP data sets. However, it poses some challenges in prediction modeling, such as a small number of observations and a large number of explanatory variables, skewness in variables, irrelevant and redundant variables, interactions among variables, and nonconstant error variance.
To solve the drawbacks regarding the data dimensionality in the GP, statistical or machine learning based approaches have been applied (Li et al., 2018). Recently, quantile regression approaches, which model the conditional distribution of the response, have been utilized to deal with the skewness and outliers in the data (Montesinos-López et al., 2019;Pérez-Rodríguez et al., 2020). Still, the question has been how to model conditional quantiles flexibly while handling data dimensionality in GP.
It is also important for researchers to identify the SNPs relevant for predicting breeding values to design future genotype studies. For instance, since the human population has been growing, the stability of food supplies has gained much more importance. Hence, recent plant breeding efforts have been given to the genetic improvement of crops. The key task is to predict a breeding value well so crops can be generated based on future requirements. Hölker et al. (2019) provided agronomic measurements and more than 500 thousand SNPs to make European flint maize landraces available for such a task.
Vine copula based (quantile) regression is a significant tool for modeling the nonlinear relationship between explanatory variables and response. It considers higher-order explanatory variables and their interaction. In addition, it can deal with unknown functional error forms. However, the current vine copula based regression methods' computational complexity makes them infeasible to be applied in high-dimensional data sets (Kraus and Czado, 2017;Chang and Joe, 2019;Zhu et al., 2021;Tepegjozova et al., 2022).
In this work, we propose two vine copula based regression methods that perform well in analyzing high-dimensional sparse data sets, where sparsity means that many explanatory variables do not predict the response. We show that their computational complexity is significantly less than the existing methods. We define relevant, redundant, and irrelevant explanatory variables for quantile regression and assess the methods' prediction power in high-dimensional sparse simulated data sets. Our analyses regarding the inclusion of relevant variables and exclusion of irrelevant variables show our methods' capability to provide sparse models. We apply the methods for genomic prediction of maize traits, proposing data preprocessing and feature extraction steps on the data given by Hölker et al. (2019). Such steps can be further applied and improved in future studies. Overall, our methods are novel and powerful and can resolve data dimensionality issues in vine copula based regression. To the best of our knowledge, there has not yet been any study performing the genomic prediction using vine copula models, as well as assessing vine copula regression methods' performance in the presence of redundant and irrelevant variables.
The paper is organized as follows: Section 2 introduces vine copulas and new methods; Section 3 provides simulation studies. We present the real data application in Section 4, discuss our findings in Section 5, and conclude in Section 6.

D-vine copulas and prediction
We recall essential concepts of vine copulas and refer to Aas et al. (2009) and Czado (2019) for more details. In the remainder, we assume all random variables to be absolutely continuous.
Copulas are distribution generating functions, allowing us to separate the univariate margins and dependence structure expressed in the copula. Let X = (X 1 , . . . , X p ) ∈ R p be a p-dimensional random vector with the joint cumulative distribution function (cdf) F and the univariate marginal distributions F 1 , . . . F p . By Sklar's theorem (Sklar, 1959), the copula C, corresponding to F , is a multivariate cdf with uniform margins such that F (x 1 , . . . , x p ) = C F 1 (x 1 ), . . . F p (x p ) . Moreover, when the univariate marginal distributions are continuous, C is unique. In addition, the p-dimensional joint density f can be written as f (x) = c F 1 (x 1 ), . . . , F p (x p ) · f 1 (x 1 ) · · · f p (x p ), x ∈ R p , where c is the copula density of the random vector F 1 (X 1 ), . . . F p (X p ) ∈ [0, 1] p .
Standard multivariate copulas, such as the multivariate exchangeable Archimedean or Gaussian, often do not accurately model the dependence among the variables. Aas et al. (2009) developed a pair-copula construction or vine copula approach using a cascade of bivariate copulas to extend their flexibilities. Such a copula construction can be represented by an undirected graphical structure involving a set of linked trees, i.e., a regular (R-) vine (Bedford and Cooke, 2002). For instance, a R-vine for p variables consists of p − 1 trees, where the edges in the first tree are the nodes of the second tree, the edges of the second tree are the nodes of the third tree and so on. If an edge connects two nodes in the (t + 1)th tree, their associated edges in the tth tree must have a shared node in the tth tree for t = 1, . . . , p − 2. The nodes and edges in the first tree represent p variables and unconditional dependence for p − 1 pairs of variables, respectively. In the higher trees, the conditional dependence of a pair of variables conditioning on other variables is modeled. To get a vine copula or pair-copula construction, there is a bivariate (pair) copula associated with each edge in the vine.
One special class of R-vines are D-vines, whose graph structure is a path, i.e., all nodes' degree in the graph is smaller than three. A node in a path represents a variable, and an edge between a pair of nodes corresponds to dependence among the variables of the respective nodes expressed by a pair copula. The node having a degree of one is called a leaf node. Once the order of the nodes in the first tree is determined, the associated D-vine copula decomposition is unique. Moreover, if the pair copulas in the higher tree levels than t are independence, where t < p and t ≥ 1, representing conditional independence, a t-truncated vine copula is obtained.
Example 2.1 Figure 1 shows the graphical specification for a 4-dimensional D-vine in the form of a set of linked trees. The D-vine with four variables has three trees, and the jth tree has 5 − j nodes and 4 − j edges, j = 1, 2, 3. The first and fourth variables are leaf nodes in the first tree. To model the dependence between a pair of variables, each edge is associated with a pair-copula. Therefore, in the first tree, the pair-copula densities, c 1,2 , c 2,3 , and c 3,4 model the dependence of the three pairs of variables, (1,2), (2,3), and (3,4), respectively. In the second tree, two conditional dependencies are modeled: between the first and third variables given the second variable, the pair (1,3;2), using associated pair-copula density c 1,3;2 , and between the second and fourth variables given the third variable, the pair (2,4;3), using associated pair-copula density c 2,4;3 . Likewise, the third tree models the conditional dependence between the first and fourth variables given the second and third variables using the associated pair-copula density c 1,4;2,3 . In addition, if c 2,4;3 is 1 everywhere in the domain, i.e., independence copula, the resulting D-vine is a 2-truncated vine.  D-vine copulas also allow to express the conditional density of a leaf node in the first tree in a closed form. Here we choose the leaf node in the first tree of a D-vine copula as a response variable. In the rest of the paper, assume that (y i , x i,1 , . . . , x i,p ), i = 1, . . . , n, are realizations of the random vector (Y, X 1 , . . . , X p ), and Y denotes a response variable with its marginal distributions F Y and the others correspond to explanatory variables with their marginal distributions F 1 , . . . , F p . For the D-vine copula with the node order 0 − 1 − . . . − p corresponding to the variables Y − X 1 − . . . − X p , p ≥ 2, the conditional density of the response given the explanatory variables f 0|1,...,p can be obtained as: (1) where F 0|1,...,j−1 is the conditional distribution function of the random variable Y |X 1 = x 1 , . . . , X j−1 = x j−1 , which can be calculated recursively using h-functions (Joe, 1996). We assume that the pair copulas in the higher tree levels than one do not depend on the conditioning values, e.g., it holds that c 0,j;1,...,j−1 (., .) = c 0,j;1,...,j−1 (., .; x 1 , . . . , x j−1 ). From now on, we make this simplifying assumption. Still, the pair copula density of the high tree levels depends on the conditioning value through its arguments (Stoeber et al., 2013).
Further, we can calculate the conditional log-likelihood cll 0|1,...,p (F Y (y), F 1 (x 1 ) . . . , F p (x p )) based on the D-vine copula as follows: Moreover, as stated in Kraus and Czado (2017), the conditional quantile function F −1 0|1,...,p at quantile α can be expressed in terms of the inverse marginal distribution function F −1 Y of the response and the conditional D-vine copula quantile function C −1 0|1,...,p at quantile α: Finally, since p-dimensional D-vine copula's input is the marginally uniform data on [0, 1] p , the estimation of the D-vine copula follows a two-step approach called the inference for margins (Joe and Xu, 1996). First, each marginal distribution is estimated. Then the data is converted into the copula data by applying probability integral transformation (PIT). We estimate the marginal distributions using a univariate non-parametric kernel density estimator with the R package kde1d (Nagler and Vatter, 2022), i.e.,F Y andF X d , d = 1, . . . , p. Next, we have the pseudo copula data: . . , n, being realizations of the random vector (V, U 1 , . . . , U p ).
2.2 Proposed methods: vineregRes and vineregP arCor Kraus and Czado (2017) proposed a D-vine copula based quantile regression with a forward variable selection, and the R package vinereg (Nagler, 2022) provides its implementation. We summarize the main steps of vinereg in Appendix A. However, its computational complexity makes it infeasible to run in high dimensions as will be discussed in Section 2.5. Thus, we propose two methods to perform a D-vine copula regression on high-dimensional sparse data sets: vineregRes and vineregP arCor.
The method vineregRes performs the variable selection at a given iteration based on the residuals of the previous iteration, i.e., the pseudo-response. It finds the variable among the candidates, which provides the best bivariate copula conditional log-likelihood conditioned on the variable and conditioning the pseudo-response of the previous iteration. Assumeỹ i i = 1, . . . , n, denote the pseudo-response and its pseudo copula data in the sth iteration, respectively, which are realizations of the random variable Y (s) and V (s) , respectively. V (0) and V (s) have the indices 0 and 0 (s) , respectively, and are always a leaf node in the first tree. vineregRes Step 1 (initialization): For given data (y i , x i,1 , . . . , x i,p ) and (v i , u i,1 , . . . , u i,p ), define the initial pseudo-responseỹ Step 2 (variable selection): Fit a parametric bivariate copula to data {(ṽ (s) i , u i,d ), i = 1, . . . , n} for d ∈ p (s) cand and denote the copula, copula density, and its estimated parameters by CR d (s) , cr d (s) andθ d (s) , respectively. Then, find the variable for which the conditional log-likelihood of the copula CR Step 3 (D-vine extension): Extend the D-vine order by adding the variable with index d * (s+1) to get a D-vine order D (s+1) = (D (s) , d * (s+1) ). Select the parametric pair copula families and estimate the parameters in the extended D-vine structure, where the associated D-vine copula and its estimated parameters are denoted byĈ (s+1) andθ (s+1) , respectively.
Step 4 (Chosen variable indices and hyperparameter updates): Extend the chosen variable set, adding the new variable, I Step 5 (Pseudo-response update or stop): If a stopping condition (see Section 2.4) does not hold, estimate the median of the response variable based on the D-vine copulâ C (s+1) and update the pseudo-response, i.e., where i = 1, . . . , n and {p 1 , . . . , p d (s+1) } ⊆ I (s+1) var . Another method to perform an efficient D-vine copula regression for high-dimensional data is to use the partial correlation between the response and a candidate explanatory variable given the chosen variables at each iteration based on their empirical normal scores. Such scores are calculated as detailed in Section 1 of Joe (2014).
Step 4 (Chosen variable indices and hyperparameter updates): As given in vineregRes.
In Section 2.7, we give an illustrative example of both methods.

Bivariate copula selection
Step 3 of vineregRes and vineregP arCor selects parametric pair copulas and estimates their parameters associated with the extension of the D-vine structure. Also, Step 2 of vineregRes fits a parametric bivariate copula to the pseudo-response and a candidate explanatory variable. The selection of pair copulas can be decided by analyzing the expected behavior of conditional quantiles in the extreme values of the variable's space. Alternatively, a model selection criterion can select a good fit among the candidate bivariate copula families. First, we estimate the parameters that maximize the log-likelihood of a candidate bivariate copula family (Appendix B). Later we can select the one with the lowest AIC or the Bayesian information criterion. While extending the D-vine structure and adding new trees at Step 3 of the methods, the fit of parametric pair copulas can be performed sequentially from the lowest to the highest trees (Brechmann, 2010).

Stopping criteria
To decide if a chosen candidate explanatory variable in a given iteration should be in the model, we will consider the conditional AIC, which penalizes the conditional loglikelihood of the given D-vine copula by the effective degrees of freedom in the model. Thus, in the sth iteration of vineregRes and vineregP arCor, the conditional AIC of the D-vine copulaĈ (s) with its effective degrees of freedomΘ (s) is var . Since we estimate the marginal distributions of the given data to have its pseudocopula data non-parametrically, we consider the effective degrees of freedom as a penalization term, for which more details are provided in Section 7.6 of Hastie et al. (2009). We stop adding variables in D-vine copula regressions when the current iteration's conditional AIC is equal to or larger than the previous iteration's conditional AIC. If the conditional AIC always improves in each iteration, we stop after all explanatory variables are included in the model.

Complexity
Assuming that the data consists of p explanatory variables, the complexity of the existing method vinereg is O(p 3 ) in terms of the total number of bivariate copulas to be selected during the algorithm (Tepegjozova, 2019). Thus, we evaluate the complexity of vineregRes and vineregP arCor using the same criterion. We will consider the worst case scenario that the algorithms run until all explanatory variables are included in the model. Further, the total number of estimated parameters is linear in terms of the number of bivariate copulas.
In vineregRes, the selection of bivariate copulas exists at Steps 2 and 3.
Step 2 fits a bivariate copula to the data of the iteration's pseudo-response and a candidate explanatory variable. Therefore, the number of bivariate copulas to be chosen at Step 2 is the total number of candidate explanatory variables at the given iteration. At the sth iteration, there are p − (s − 1) candidate explanatory variables. Thus, the number of bivariate copulas to be selected at this step is given by In vineregP arCor, the variable selection is based on partial correlations; thus, its Step 2 does not perform any selection of bivariate copulas.
Step 3 of vineregRes and vineregP arCor extends the D-vine structure, adding the selected explanatory variable at Step 2. Thus, at the sth iteration, it results in the selection of s pair copulas. Hence, Step 3 selects the following number of bivariate copulas: Considering Equations (3) and (4), we calculate the total complexity of vineregRes: Considering Equation (4), we calculate the total complexity of vineregP arCor: Equations (5) and (6) show that the complexity of vineregP arCor and vineregRes in terms of the total number of selected bivariate copulas is O(p 2 ). As a result, our methods significantly reduce the computational complexity of vinereg.

Relevant, redundant, and irrelevant variables
Now we define relevant, redundant, and irrelevant variables for predicting the conditional quantile of a response variable Y given the set of explanatory variables X . We will denote the cdf of X by F X in the following.
where the vectors x M and x R include the variables in the sets M and R, respectively, R ⊆ X , M ⊆ X , R ∩ M = ∅.
M is relevant for Y . Since R is not independent of M, R can be relevant for Y . However, when M and R are considered together, R is redundant.
where the vectors x M and x I include the variables in the sets M and I, respectively,

Illustrative example
Now we illustrate how vinereg, vineregRes, and vineregP arCor select variables. We simulate the data with 450 observations from the model (Y, X 1 , X 2 , X 3 , X 4 ) ∼ N 5 (0, Σ) with (0.40, 0.70, 0.32, 0, 0, 0, 0, 0, 0, 0) vectorizing the upper triangular part of symmetric matrix Σ. Thus, there are a response variable and four explanatory variables in the data, and the third and fourth variables are assumed to be irrelevant for predicting the response. First, we convert the observations on the copula scale (v i , u i,1 , u i,2 , u i,3 , u i,4 ), i = 1, . . . , 450, using the non-parametric kernel density estimator. Then we define the initial pseudo-response on the copula scaleṽ For an illustration purpose, we will show the methods' variable selection step in the third iteration, where it holds that the D-vine order D (3) = (0, 2, 1) for the D-vine copulâ C (3) with the parametersθ (3) , the chosen variable index set I  (3) is 975.30. The next iteration for all methods is to decide which and if the third or fourth variable should be added to the model.
To make the variable selection, vinereg extends the D-vine order as D (4) = (0, 2, 1, 3) for the third variable and D (4) = (0, 2, 1, 4) for the fourth variable. Accordingly, three new pair-copula families are selected for each vine order, thereby estimating their parameters. Thus, vinereg selects six bivariate copulas in this iteration (one for the first/second/third tree level for each D-vine order) to choose one variable. For more details about vinereg, see Example 3.2 in Kraus and Czado (2017).
To perform the variable selection, first, vineregRes estimates the median of the response based on the D-vine copulaĈ (3) . Then it calculates the pseudo-response and estimates the pseudo-response's marginal distribution nonparametrically. Next, the associated pseudo-response on the copula scale is obtained by the PIT, i.e., y , vineregRes identifies the third variable as the candidate explanatory variable to be added to the D-vine copula. Thus, it extends the D-vine order as D (4) = (0, 2, 1, 3), selecting the new pair-copula families and estimating new parameters. The resulting D-vine copulaĈ (4) has the conditional AIC of 975.30. SinceĈ (4) does not have a better conditional AIC thanĈ (3) , vineregRes stops the iterations and returns the final D-vine copula asĈ (3) . Thus, vineregRes chooses five bivariate copulas here. vineregP arCor chooses the next variable based on its partial correlation with the response given the chosen explanatory variables using their normal scores. It estimateŝ ρ 0,3;2,1 = 0.05 for the third variable andρ 0,4;2,1 = 0.03 for the fourth variable. Since the third variable has a higher (absolute) estimated partial correlation than the other, vineregP arCor extends the D-vine order as D (4) = (0, 2, 1, 3). Next, it selects the new pair-copula families and performs new parameter estimation associated with D (4) . Since the D-vine copula with D (4) = (0, 2, 1, 3) does not have a better conditional AIC than D (3) = (0, 2, 1), vineregP arCor gives the final model fitĈ (3) . Thus, vineregP arCor selects three bivariate copulas in this iteration.

Simulation study
We show the flexibility and effectiveness of the proposed methods on simulated datasets being nonlinear and having different sparsity. We explore the following questions: (Q1) How do vineregRes and vineregP arCor work in situations with nonlinear explanatory variable effects on the response's quantiles in the presence of redundant and irrelevant variables for prediction accuracy and computational complexity? (Q2) How well do vineregRes and vineregP arCor identify relevant and irrelevant variables for predicting the response's quantiles? (Q3) How do vineregRes and vineregP arCor perform compared to vinereg?

DGP2: redundant variables
where the samples of explanatory variables are independently generated from a multivariate normal distribution with a Toeplitz correlation structure, i.e., (X i,1 , . . . , X i,p d ) ∼ N p d (0, Σ), i = 1, . . . , n, j = 1, 2, 3, with the (a, b)th element of the covariance matrix Σ a,b = ρ |a−b| . To represent a challenging but realistic scenario, we set ρ = 0.75. We sample i ∼ N (0, 1), i = 1, . . . , n (iid) independently from the explanatory variables and set σ i ∈ {0.5, 1}, i = 1, . . . , n. Here, all variables are predicting the response's quantiles. However, given the set of the first ten variables, the others are redundant. To investigate the impact of the redundant variables on the methods' predictive power and sparsity of the model in terms of the total number of variables, we change the number of redundant variables and set it to (p d − 10) in each case d = 1, 2, 3, 4: Case 1 with p 1 = 20 (50% of variables are redundant given the ten relevant ones), Case 2 with p 2 = 40 (75% of variables are redundant given the ten relevant ones), Case 3 with p 3 = 100 (90% of variables are redundant given the ten relevant ones), Case 4 with p 4 = 1000 (99% of variables are redundant given the ten relevant ones).
Based on the DGPs in Equations (7) and (8), we simulate samples with size 450 (n=450) each time with a random split of 300/150 observations for a training set/a test set. We replicate our procedure 100 times and average performance measures per sample. Appendix D shows nonlinear relationships between the response and others.

Performance measures
We consider the computation time, the number of chosen variables, true positive rate (TPR), and false discovery rate (FDR) as methods' performance measures on the training set. To evaluate the performance of a method on the test set, we apply the pinball loss (P L α ) at α = 0.05, 0.50, 0.95. True positive rate (TPR) is the ratio of the number of chosen relevant variables by a method M to the total number of relevant variables. False discovery rate (FDR) is the ratio of the number of chosen irrelevant variables to the total number of chosen variables by a method M . Formally: where X rel is the set of relevant variables, X irrel is the set of irrelevant variables, and X M chosen denotes the set of chosen variables by a method M . Higher TPR and smaller FDR are better. The pinball loss (P L α ) measures the accuracy of the quantile predictionsŷ α,M i at the level α by a method M compared to the given response y i , i = 1, . . . , n and has the form of The pinball loss is identical to check score, and smaller pinball loss values are better (Steinwart and Christmann, 2011).

Results
All computations are run on a single node CPU with Intel Xeon Platinum 8380H Processor with around 25 GB RAM, running R version 4.2.0. Variable selection and computational complexity results on the training set: Table 1 compares the average performance of the methods for cases 1-3 and 1-4 specified in Equations (7) and (8), respectively, concerning the variable selection and computational complexity. We analyze the TPR and FDR only for the DGP1 since all variables are relevant in the DGP2. Furthermore, vinereg was not complete within three days per sample for the fourth case of the DGP2, making it computationally infeasible.
In all cases of the DGP1, vinereg has a better performance than the others regarding the true positive rate. However, its false discovery rate is also higher than others, adding many irrelevant variables into a model. For vineregRes, we observe that it correctly identifies more than 75% of the relevant variables in all cases of the DGP1. Furthermore, the most promising is that its false discovery rate is less than 15% there, making it the best method with respect to the FDR. On the other hand, vineregP arCor's TPR is higher than 50% in all cases of the DGP1. Nonetheless, like others, its FDR increases as the number of irrelevant variables increases in the model, reaching more than 50% in the third case of the DGP1.
Concerning the number of chosen variables by a method in the DGP1 and DGP2: While vineregRes selects the lowest number of variables between five and six as reflected in its high TPR and low FDR in the DGP1, vinereg includes almost half of the total number of variables in the data in each case. This highlights the power of vineregRes regarding the exclusion of irrelevant variables in sparse data sets. vineregP arCor's number of chosen variables is between vinereg and vineregRes in all evaluated cases. In an ultra high-dimensional case with 1000 explanatory variables, i.e., the last case in the DGP2, the number of variables chosen by vineregP arCor is, on average, 33.05 with the empirical standard error of 1.07 while it is 7.52 with that of 0.58 for vineregRes. Hence, our methods perform high-dimensional sparse vine copula regression as proposed.
As the number of variables increases, the average running time for all methods increases, and vineregP arCor provides the fastest computation as expected from the results in Section 2.5. In the second case of the DGP2, where 75% of variables are redundant given the ten relevant ones, vineregRes is on average 15 times faster than vinereg. In the third case of the DGP2 with 100 explanatory variables, vinereg takes more than three hours, whereas our methods run less than 2 minutes on average. Likewise, our methods vineregRes and vineregP arCor run less than 15 minutes in the ultra high-dimensional case. Table 1: Comparison of the methods' performance on the training set over 100 replications under the cases 1-3 and 1-4 specified in Equations (7) and (8) Prediction accuracy results on the test set: Table 2 shows that vineregRes provides the lowest pinball loss at the medium quantile in all cases of both DGPs. Moreover, it gives the best accuracy at the levels 0.05 and 0.95 in all cases of the DGP2, except the upper quantile in the first case. Since vineregRes gives the most sparse models in the DGP2 in Table 1, we can infer that including many relevant but potentially redundant variables in vinereg and vineregP arCor is worsening the prediction accuracy in the DGP2. However, vinereg provides better fits concerning the upper quantile than others in the first two cases of the DGP1. A likely explanation can be that including irrelevant variables in addition to the most relevant ones in a model does not negatively impact the prediction accuracy as redundant variables do but increases its computational complexity unnecessarily as shown in Table 1. In addition, Table 1 shows that half of the chosen variables by vinereg are irrelevant in the second case of the DGP1, resulting in a trade-off between prediction accuracy and inclusion of irrelevant variables in the model. vineregP arCor's accuracy is better than vinereg in 11 cases out of 18 evaluated cases and vineregRes in three cases out of 21 evaluated cases. Table 2: Comparison of the average performance of the methods on the test set at different quantile levels over 100 replications under the cases 1-3 and 1-4 specified in Equations (7) and (8)

Application: the genomic prediction of maize traits
We describe a real-data application on the doubled-haploid (DH) lines from European flint maize landraces that motivates our methods' usage. Hölker et al. (2019) evaluated 899 DH lines whose data contains genotypic measurements with the SNP array technology and phenotypic measurements of agronomic traits across environments. We can regard a DH line as an observation, a SNP as an explanatory variable, and a trait as a response in regression.
In this study, we are interested in the relationship between a DH line's genotype encoded by its SNPs and its phenotypic outcome described by its traits, i.e., the genomic prediction of maize traits. More specifically, we would like to find relevant SNPs for a trait in a multivariate prediction model using our high-dimensional sparse vine copula regression methods, performing variable selection.

Data description and preprocessing
There are three landraces in the data, and we focus on the Kemater Landmais Gelb (KE) landrace, which has the largest number of observations (471 out of 899). There are 501,124 explanatory variables, SNPs, which contain only zero and two values, e.g., zero corresponds to the genotype TT, and two denotes the genotype CC.
In addition, we predict four responses of agronomic traits separately: early plant height at the fourth and sixth stages (PH V4/V6), female flowering time (FF), and male flowering time (MF). They are quantitative and continuous as shown in Figure 2. Plant breeders need to increase the early plant development and adopt maize genotype in future studies. Also, plant breeders may avoid decreasing or increasing female and male flowering times during adoption. Thus, the traits' prediction from the genotypic measurements is crucial.
To compare the performance of regression methods, we partition our data randomly into training (67%) and test (33%) sets. As a result, the training and test sets contain 314 and 157 observations, respectively.
Moreover, we remove the duplicate explanatory variables, retaining only one. Next, we remove the common explanatory variables. For instance, assume an explanatory variable in our training set contains 300 zero values and 14 two values. Then, such a variable does not differ among the observations and might not be expected to have predictive power on a response. As a result, the number of explanatory variables in the training and test sets decreases from 501124 to 44789, i.e., we retain around 9% of the initial explanatory variables.

Feature extraction
Since our explanatory variables are binary, and there can be associated latent variables with a prediction power on the response, we focus on estimating these latent variables and using them as extracted features (new explanatory variables) in a regression method. Our approach is to group the explanatory variables and estimate their weights in their groups so that such weights are used to estimate the latent variables representing each group.
Let y k and SN P p denote the response k and explanatory variable p, respectively.
5. Estimate the latent variables, i.e., create the features f eature G k,d k by grouping G explanatory variables in O k and using their coefficients from the first step: . . , n k , n k = |O k | G , k = 1, . . . , 4. As a result, we have 174 (87) features for FF, 92 (46) features for MF,198 (99) features for PH V4, and 183 (93) features for PH V6 by grouping G = 100 (G = 200). Figure 2 shows a scatter plot of a selected continuous feature and a trait. Figure 2: Scatter plots of an extracted feature, combining 100 SNPs in a feature, and the traits. Feature 1 corresponds to the combination of the SNPs whose p-value from the OLS of the associated trait is lower than others. Similar correspondence applies to other features.

Prediction
We have our data D G k = (y k , f eature G k,1 k , . . . , f eature G k,n k ) for each response k = 1, . . . , 4 and G ∈ {100, 200}. Now we apply our methods and linear quantile regression (Koenker and Hallock, 2001) penalized with a LASSO function (LQRLasso) for variable selection to eight different data sets' training set. For the latter, we use the R package rqPen (Sherwood and Maidman, 2020) and provide more details on its fit in Appendix E. Next, we apply the fitted models to the test sets to evaluate the models' predictive performance using the pinball loss defined in Section 3.2 at the quantile levels 0.05, 0.50, 0.95. Table 3 shows that vineregRes is the most sparse and accurate model for PH V6, independent of G and quantile levels. It chooses four features for G = 100 and two for G = 200, identifying less than 5% of the features as relevant. For FF and PH V4, whenever LQRLasso is more accurate than vineregRes, it includes more features, giving a trade-off between model sparsity and accuracy. Further, vineregRes performs worse than LQRLasso and vineregP arCor for MF. We observe that dependencies among MF and its selected features by vineregRes are more linear than those among other traits since it fits mostly the Gaussian copula in the first tree of the D-vine for MF (Appendix F). Moreover, we remark that linear quantile regression by LQRLasso may perform well if it can avoid crossing quantile curves, but there is no guarantee that the 95% quantile curve exceeds the 90% quantile curve everywhere as shown in Appendix G. We also provide selected feature indices by our methods in Appendix H.
Our SNP screening and feature extraction steps are similar to Qian et al. (2020). Yet, they fit a simple linear regression on the first feature, which includes linearly and marginally most important SNPs. Thus, we compare our models' performance for the pinball loss at the levels 0.05, 0.50, 0.95 by fitting a vine copula regression on the first feature for PH V4 and PH V6, i.e., bivariate analyses or D-vines with two nodes: response and the first feature. We call such a bivariate regression with the first feature bicopreg. We see in Appendix H that having the first and 88th features in the model as in vineregRes performs better than bicopreg using G = 200 for PH V6. However, including more features as in vineregP arCor, which selects 15 features with G = 200, including the first one, worsens the performance compared to bicopreg. Similar results using G = 100 for PH V6 hold.
Likewise, the mean of selected feature indices using G = 100 for PH V4 is 102.33 for vineregRes and 106.18 for vineregP arCor, and both methods do not select the first feature as relevant. Nevertheless, their prediction power is stronger than bicopreg. Hence, linearly and marginally most important SNPs' groups might not be considered relevant for prediction when considering nonlinear dependencies as in our methods.

Discussion
High-dimensional sparse vine copula regression is a significant tool for efficiently allowing nonlinear relationships between explanatory variables and response and selecting relevant variables. In genomic prediction, genotypic measurements like SNPs are often very high-dimensional, which might be reduced by considering some SNP groups and their interactions. Also, many groups may be irrelevant for prediction. Our methods can handle such situations and predict responses at different quantile levels. Their performance might be improved with bivariate copula families having more asymmetries, e.g., more than two parameters. For our application, consider the following question: How would selected SNPs impact the low and high quantiles of the trait PH V4? Our methods provide a good answer as shown in Appendices F and H. Thus, plant breeders can assess the chosen SNPs' impact on the trait's various quantile levels and identify the associated SNPs. Moreover, chosen SNP groups can be compared to other genome-wide association studies, helping plant breeders to decide on future genotype adoption. In our application, comparing the identified SNPs with those in Mayer et al. (2020) is high on the agenda.
Feature extraction is a vital step that may impact our methods' genomic prediction power. For instance, the choice of SNPs' weights estimating their latent variable is open to future research. Also, even though it offers a trade-off between a computational burden and prediction power, one can apply cross-validation for the choice of the SNP's group size G. In addition, some SNPs might affect the trait, not marginally only in the presence of certain other SNPs. Alternatively, one may remove the P -value screening of the SNPs at the 10% level described in Section 4.1 and consider all possible extracted features. Likewise, some SNPs might influence the trait marginally, but not when certain other SNPs are in the model. For such cases, some post-processing steps for feature extraction might be applied.
Finally, our variable selection steps can be adapted for the vine copula based regression using more flexible vine tree structures, such as R-vines. Thus, they extend the vine copula based regression usage to areas with high-dimensional data sets.

Conclusion
We propose two methods to perform high-dimensional sparse vine copula based regression and analyze their performance in the presence of redundant and irrelevant variables that we define through simulation studies. We apply our approaches to the genomic prediction of maize traits, providing guidelines on real-data preprocessing and feature extraction. Finally, we discuss future research directions regarding the further application of vine copulas in genomic prediction and the extension of our methods with high-dimensional data sets to more general vine copula classes.

Acknowledgments
The

A vinereg
Assume the response variable's index is denoted by 0, and the response variable is a leaf node in the first tree of a D-vine.
Step 1 (initialization): Given data (y i , x i,1 , . . . , x i,p ) and (v i , u i,1 , . . . , u i,p ), the initial D-vine structure order D (1) = (0), the initial chosen variable index set I Step 2 (variable selection): Extend the D-vine structure order by adding a variable whose index is d to have a D-vine structure order D (s+1) = (D (s) , d). Fit a parametric D-vine copula having the structure D (s+1) and denote it and its estimated parameters bŷ CD d (s+1) andθ d (s+1) , respectively. Then, find the variable with the index d * (s+1) for which the conditional log-likelihood of the D-vine copulaĈD Step 3 (D-vine extension): Extend the D-vine structure order by adding the variable whose index is d * (s+1) to have a D-vine structure order D (s+1) = (D (s) , d * (s+1) ). Select the new parametric pair copula families and estimate their parameters in the extended D-vine structure. Then, denote the associated D-vine copula and its estimated parameters bŷ C (s+1) andθ (s+1) , respectively.
C Redundant variables in a D-vine copula Proposition C.1 (Redundant variables in a D-vine copula) If a p-dimensional D-vine copula is a t-truncated vine, where the response is a leaf node represented by the first node, t > 1, p ≥ 3, t ≤ p − 1, there are p − t − 1 redundant or irrelevant variables.
Example C.1 Suppose that the first node and others represent the response and three explanatory variables in the D-vine in Figure 1, respectively (p = 4). Assume that the D-vine is a 2-truncated vine (t = 2) and the copula for the pair (3, 4) is not independence. Then, there is one redundant variable for the response, which is the variable represented by the fourth node, given the two relevant variables represented by the second and third nodes.

E Penalized quantile regression with LASSO function (LQRLasso)
It extends the linear quantile regression method, adding a LASSO penalty function in the quantile regression objective function to perform the variable selection. Thus, it solves the following optimization problem to estimate the coefficients at the quantile level α: where ρ α (r) = r(α − 1(r < 0)) denotes the check function, λ ∈ R corresponds to a tuning parameter for the penalization of the coefficients, (y i , x i,1 , . . . , x i,p ), i = 1, . . . , n, are realizations of the random vector (Y, X 1 , . . . , X p ), and Y denotes a response variable, and the others correspond to explanatory variables. We perform a five-fold cross-validation using the training set to optimize λ in each replication and set the associated coefficients below 10 −8 to zero. The variables with zero coefficients are not relevant for predicting the response's quantiles. Such procedures are implemented in the R package rqPen.
Since it is not straightforward to decide which transformation of variables or some interaction effects should be applied in a linear model in the high-dimensional data, we analyze the variables without any transformations. Since a LASSO penalty function is sensitive to the feature scales, we apply a standardization with zero mean and unit variance for each feature before its fit.

G Quantile crossing of LQRLasso
LQRLasso performs variable selection and linear quantile regression at a specified quantile level at the same time. However, its final model fit might result in the crossing of quantile curves. For an illustration, we choose the first and fourth features obtained using G = 100 as explanatory variables for predicting the trait MF at the levels 90% and 95% and let LQRLasso run for both levels. Then only the fourth feature is identified in both, thereby estimating its coefficient. However, as seen in the following plot, the coefficient estimated by LQRLasso leads to the quantile crossing.