domain: stat sdv

text:
Benjamin Heuclin

Marie Denis

Catherine Trottier

Sébastien Tisné

Frédéric Mortier
email: frederic.mortier@cirad.fr

Continuous shrinkage priors for fixed and random e← ects selection in

Keywords: Angular parametrization, Fixed and random e↵ects selection, Horseshoe prior, Linear mixed model, Quantitative genetics

Introduction

Linear mixed models (LMM) are flexible tools for modeling data from a wide range of data types in various applied fields including ecology and evolution [START_REF] Bolker | Generalized linear mixed models: a practical guide for ecology and evolution[END_REF][START_REF] Ives | Generalized linear mixed models for phylogenetic analyses of community structure[END_REF], quantitative genetics [START_REF] Lynch | Genetics and Analysis of Quantitative Traits[END_REF] or medical research [START_REF] Brown | Applied Mixed Models in Medicine[END_REF]. An important practical point when using linear mixed models is the choice of random e↵ect components. Choosing which random grouping factors to include is vital to model appropriate dependence structures within data. This issue compounds when the possible number of random e↵ects is large, which also leads to identifiability problems and estimation instability. With new technologies, high-resolution satellite images, high-throughput genotyping/phenotyping techniques, such contexts are now very common. For instance, in quantitative genetics, linear mixed models are commonly used for IBD-QTL mapping [START_REF] George | Mapping quantitative trait loci in complex pedigrees: a two-step variance component approach[END_REF][START_REF] Van Eeuwijk | Mixed model approaches for the identification of qtls within a maize hybrid breeding program[END_REF][START_REF] Tisné | Mixed model approach for ibdbased qtl mapping in a complex oil palm pedigree[END_REF] or gene-set analyses (GSA) [START_REF] Fridley | Gene set analysis of snp data: benefits, challenges, and future directions[END_REF], allowing to decompose the global genetic e↵ect as a sum of local/specific e↵ects related to each position/block considered on the genome. In such studies, the genetic e↵ects are associated to random e↵ects that can vary from tens to a few hundreds.

This work has been motivated by two original applications. The first one takes place in the IBD-QTL mapping context and aims at identifying QTLs related to the oil palm production (Elaeis guineensis, Jacq). In this application, 135 random e↵ects associated to 135 genetic positions are considered. In the functional mapping framework, the second application aims at modeling the e↵ects of 38 genetic markers assumed varying with time, and at selecting those involved in the dynamics of shoot growth of Arabidopsis thaliana (L. Heynh) [START_REF] Marchadier | The complex genetic architecture of shoot growth natural variation in Arabidopsis thaliana[END_REF]. While providing an alternative modeling approach to the usual functional mapping methods based on non-parametric strategies [START_REF] Ma | Functional Mapping of Quantitative Trait Loci Underlying the Character Process: A Theoretical Framework[END_REF], our work also allows to estimate potential dependencies between random e↵ects.

The usual solution for dealing with identifiability and/or inference instability problems is to reduce the number of variables, using model choice procedures based on information criteria [START_REF] Müller | Model selection in linear mixed models[END_REF]. An alternative strategy relies on regularization approaches (Bickel et al., This paper has been submitted for consideration for publication in Biometrics 2006). In linear models (LM), regularization procedures have been widely studied and a large set of penalty functions has been proposed (see for example [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]; [START_REF] Desboulets | A review on variable selection in regression analysis[END_REF]). Most regularization methods, initially developed in a frequentist context, have been proposed in the Bayesian framework [START_REF] Kyung | Penalized regression, standard errors, and Bayesian lassos[END_REF]. Prior distributions act as penalty terms in the frequentist approach. A set of priors has been extensively developed, among others: spike-andslab [START_REF] George | Variable selection via Gibbs sampling[END_REF], Bayesian Lasso [START_REF] Park | The bayesian lasso[END_REF], Elastic-Net [START_REF] Kyung | Penalized regression, standard errors, and Bayesian lassos[END_REF], normalgamma [START_REF] Gri N | Inference with normal-gamma prior distributions in regression problems[END_REF] or horseshoe [START_REF] Carvalho | Handling sparsity via the horseshoe[END_REF] priors.

In the LMM context, literature is less developed. The first approaches to simultaneously select fixed e↵ects and variance components used model choice criteria [START_REF] Rao | A strongly consistent procedure for model selection in a regression problem[END_REF][START_REF] Vaida | Conditional akaike information for mixed-e↵ects models[END_REF][START_REF] Müller | Model selection in linear mixed models[END_REF][START_REF] Delattre | An iterative algorithm for joint covariate and random e↵ect selection in mixed e↵ects models[END_REF]. As in the LM context, penalized likelihood approaches have also been developed as alternatives, especially when the number of predictors increases. For instance, [START_REF] Bondell | Joint variable selection for fixed and random e↵ects in linear mixed-e↵ects models[END_REF] and [START_REF] Ibrahim | Fixed and random e↵ects selection in mixed e↵ects models[END_REF] combined penalized likelihood techniques, using adaptive lasso or smoothly clipped absolute deviation (SCAD) penalties, with the Cholesky decomposition of the random e↵ects covariance matrix or its modified version. In the high dimensional context, [START_REF] Fan | Variable selection in linear mixed e↵ects models[END_REF] proposed a two steps approach, while [START_REF] Li | Doubly regularized estimation and selection in linear mixed-e↵ects models for high-dimensional longitudinal data[END_REF] proposed a doubly regularized estimation and selection of fixed and random e↵ects in longitudinal data. The reparametrization of the LMM by the use of a modified Cholesky decomposition of the random e↵ects covariance matrix has initially been put forward in a Bayesian framework by [START_REF] Chen | Random e↵ects selection in linear mixed models[END_REF]. Such a technique allows to consider the standard deviations of random e↵ects as regression parameters.

In their approach, a spike-and-slab prior is used for the fixed e↵ects and the standard deviations, with a truncated Gaussian distribution for the slab part associated to the standard deviations. [START_REF] Frühwirth-Schnatter | Bayesian parsimonious covariance estimation for hierarchical linear mixed models[END_REF] propose a related approach, modeling directly the Cholesky decomposition elements using a spike-and-slab prior where the slab distribution is an unconstrained Gaussian distribution. As for the fixed e↵ect selection, alternative approaches to spike and-slab priors have also been developed and compared for the selection of variance components. In particular, continuous shrinkage priors, or their mixture versions, such as Student and Cauchy distributions, Bayesian Lasso and normal-gamma priors have been studied in the context of random intercept models [START_REF] Frühwirth-Schnatter | Bayesian variable selection for random intercept modeling of gaussian and non-gaussian data[END_REF]. Finally, we note that the question of shrinking variance parameters towards zero does not raise only in the LMM context. Such objectives have been studied in di↵erent statistical contexts. In structured additive regression models [START_REF] Fahrmeir | Penalized structured additive regression for space-time data: a bayesian perspective[END_REF] for instance, groups of fixed e↵ects are selected using a spike-and-slab prior on specific group variance components based on a mixture of inverse gamma distributions [START_REF] Scheipl | Spike-andslab priors for function selection in structured additive regression models[END_REF]. In time varying parameters and state-space models, [START_REF] Bitto | Achieving shrinkage in a time-varying parameter model framework[END_REF] and [START_REF] Cadonna | Triple the gamma-a unifying shrinkage prior for variance and variable selection in sparse state space and tvp models[END_REF] propose the use of double or triple gamma priors extending Normal-Gamma [START_REF] Gri N | Inference with normal-gamma prior distributions in regression problems[END_REF] or more generally the scaled-Beta distribution class [START_REF] Pérez | The scaled beta2 distribution as a robust prior for scales[END_REF].

In this paper, we propose to combine the horseshoe prior with its folded version to simultaneously select fixed and random e↵ects. We study the performances of the proposed prior through the two applications and discuss results in comparison with the two commonly used alternative priors: the folded Cauchy prior and the folded spike-and-slab prior. In the second application, to model dependency structures between random slopes e↵ects, we apply polar parametrization [START_REF] Pinheiro | Unconstrained parametrizations for variance-covariance matrices[END_REF] using sinusoidal prior on angles [START_REF] Pourahmadi | Distribution of random correlation matrices: Hyperspherical parameterization of the cholesky factor[END_REF] to ensure symmetry and positive-definiteness of the unknown correlation matrix. The paper is organized as follows. Section 2 presents the general model, priors formulation, and the specific context associated with each of the two applications along with their dedicated models. In section 3, we present the computational aspects of the Bayesian inference in order to optimize the MCMC algorithms. In Section 4, we firstly discuss results obtained by the three priors from a statistical point of view, and we then interpret results from a biological point of view.

Model specification and priors formulation 2.1 General considerations

LMMs can be expressed in the following general form:

y = X + Z e u + " (1)
where y is a n-response vector, X a n ⇥ (p + 1)-matrix of p covariates with a first unitary column for the intercept, Z a n ⇥ s known sparse random e↵ects design matrix associated to ũ a s-vector of random e↵ects assumed to be distributed as a multivariate Gaussian distribution with null expectation and covariance matrix denoted by ⌦. In the following, ũ and " are assumed independent. Such a formulation encompasses a broad set of LMMs. Each model leads to consider di↵erent design matrices X and Z, and variance matrices ⌦.

For the variance component model, considering q independent random e↵ects ũl , with c l levels each, following a Gaussian distribution centered on zero with covariance matrix ⌦ l (l = 1, . . . , q), then Z = L q l=1 Z l with Z l the l th random e↵ect design matrix, ũ = (ũ 0 1 , . . . , ũ0 q ) 0 , ⌦ is a block diagonal matrix where each block is the c l ⇥ c l matrix ⌦ l (⌦ = bdiag(⌦1, . . . , ⌦q)) and s = P q l=1 c l . Here L denotes the column concatenation operator. Z ũ can be decomposed as a sum of random e↵ects Z ũ = Z1 ũ1 + • • • + Zq ũq. Two special cases of the variance component model may be considered. The first one is the usual variance component model where ⌦ l = 2 l Ic l with Ic l the identity matrix of size c l . In this case, levels of each random e↵ect are assumed independent. The second one is the animal model where c l = n, Z l = In and ⌦ l = 2 l A l , where A l is a known n ⇥ n-IBD matrix. In the random intercept and slope (RIS) context with one grouping factor with c levels and p covariates, Z = J • X where • is the face-splitting product (row-by-row Kronecker product, see web appendix B for more details), J corresponds to the n⇥c-0/1-design matrix associated to the random e↵ect and ⌦ is a block diagonal matrix such as ⌦ = Ic ⌦ ⌦, where ⌦ is a (p + 1) ⇥ (p + 1) unknown correlation matrix related to dependencies between random intercept and slopes. It is straightforward to extend to q random e↵ects with c l levels each.

When p and s are large, such models (see equation 1) must be regularized. The already proposed approaches are mainly based on the Cholesky decomposition of ⌦ [START_REF] Chen | Random e↵ects selection in linear mixed models[END_REF][START_REF] Bondell | Joint variable selection for fixed and random e↵ects in linear mixed-e↵ects models[END_REF]. In this paper, we propose to use this decomposition: ⌦ = ⇤R⇤, where ⇤ is a diagonal matrix and R the associated correlation matrix. A general LMM model (see equation 1) can then be reformulated as:

y = X + Z⇤u + " = ⇥ X, u 0 ⌦ Z P ⇤  + " (2)
where is the unique diagonal elements vector of ⇤ and P is the matrix that transforms to vec(⇤) [START_REF] Ibrahim | Fixed and random e↵ects selection in mixed e↵ects models[END_REF]. Calculation details are presented in web appendix B. Now, u is the vector of Gaussian random e↵ects Ns (0, R). Finally, when ⌦ is supposed to be unknown (RIS models), we propose to write the correlation matrix R using the polar parametrization [START_REF] Pinheiro | Unconstrained parametrizations for variance-covariance matrices[END_REF][START_REF] Pourahmadi | Distribution of random correlation matrices: Hyperspherical parameterization of the cholesky factor[END_REF]. Such an approach ensures that the correlation matrix R, sampled though the posterior distribution, is a valid symmetric and positive-definite matrix with 1's on the diagonal. More details are given in section 2.3.2 bellow. Finally, " is a multivariate Gaussian residual vector assumed to be independent of u.

Priors formulation

To achieve the selection of fixed e↵ects ( ) and scale parameters ( ), we consider in this work local-global priors [START_REF] Polson | Sparsity information and regularization in the horseshoe and other shrinkage priors[END_REF][START_REF] Piironen | Sparsity information and regularization in the horseshoe and other shrinkage priors[END_REF]. Such priors, initially used for the selection of fixed e↵ects only, consist in a scale mixture of Gaussian distribution on parameters aj subjects to selection:

aj|⌧ 2 , ! 2 j , 2 ⇠ N 0, 2 ⌧ 2 ! 2 j , j = 1, . . . , q, ( 3
)
where 2 is the residual variance, ⌧ 2 is the global parameter while ! 2 j are the local ones. ⌧ 2 allows to shrink all coe cients towards zero while local parameters ! 2 j highlight non-null parameters. Such a prior encompasses a large set of well known priors. Assuming ! 2 j ⌘ 1 leads to global priors such as the well known normal-inverse-gamma (NIG) prior [START_REF] Gelman | Bayesian Data Analysis[END_REF], while assuming ⌧ 2 ⌘ 1 leads to local priors such as the Laplace prior [START_REF] Park | The bayesian lasso[END_REF], the student prior [START_REF] Gelman | Prior distributions for variance parameters in hierarchical models (comment on article by browne and draper)[END_REF] or the normal-gamma prior (Gri n et al., 2010) (see table 1 in web appendix A). Among local-global priors, the horseshoe prior assumes that local, as well as global parameters, are distributed from folded-Cauchy distributions [START_REF] Carvalho | Handling sparsity via the horseshoe[END_REF]. The horseshoe prior has demonstrated high performances for the selection of fixed e↵ects, comparable to the spike-and-slab prior [START_REF] Van Erp | Shrinkage priors for bayesian penalized regression[END_REF]. In this article, we propose to investigate horseshoe priors to simultaneously select fixed e↵ects and standard deviations. We note that since j is positive, the Gaussian distribution is replaced by a folded-Gaussian distribution

N + 0, 2 ⌧ 2 ! 2 j .
The horseshoe prior is systematically applied to fixed e↵ect. For standard deviations, we evaluate the folded horseshoe prior (fHS) addressing two specific questions: (i) what are the performances of such priors to select random e↵ects, and (ii) which impacts on the fixed and random e↵ect estimations. These results are discussed relatively to the use of alternative priors: folded Cauchy (fC) and folded zero-inflated (spike-and-slab, fSS).

In the specific RIS context, random intercept and slopes are commonly assumed to be non-independent through an unknown correlation matrix R such that R = Ic ⌦ R. Different priors have been proposed [START_REF] Lewandowski | Generating random correlation matrices based on vines and extended onion method[END_REF][START_REF] Pourahmadi | Distribution of random correlation matrices: Hyperspherical parameterization of the cholesky factor[END_REF]. Here, we adopt the polar parametrization introduced by [START_REF] Pinheiro | Unconstrained parametrizations for variance-covariance matrices[END_REF]. It consists in the use of a hyperspherical parametrization of the Cholesky factors of the correlation matrix and an appropriate distribution on related angles. [START_REF] Pinheiro | Unconstrained parametrizations for variance-covariance matrices[END_REF] demonstrate that any correlation matrix R can be factorized as BB 0 , with B1,1 = 1, Bi,1 = cos(✓i,1), i = 2, . . . , q and Bj,i =

( Q j 1 m=1 sin(✓i,m) for i = j, cos(✓i,j) Q j 1 l=1 sin(✓ i,l ) for 2 6 j 6 i 1, (4)
Pourahmadi and Wang (2015) proposed the following ✓'s sinusoidal distribution :

✓i,j / sin(✓) 2k+(p+1) j 0<✓<⇡ , i = j + 1, . . . , p, (5)
where k is a non-negative constant. This distribution ensures that angles are centered on ⇡/2 or equivalently that the distribution of R is centered on the identity matrix. Moreover and interestingly, parameter k can be interpreted as a shrinkage parameter [START_REF] Ghosh | Bayesian estimation of correlation matrices of longitudinal data[END_REF]. For instance, if k = 0 then R is distributed as a uniform distribution on the set of all (p + 1) ⇥ (p + 1)positive-definite correlation matrices, while if k tends to infinity, the distribution of R tends to a point mass on the unit diagonal (p + 1) ⇥ (p + 1)-matrix (see Figure 1). In this work, k is chosen based on cross-validation procedure.

2.3 Specific applied contexts 2.3.1 The oil palm dataset. We analyse this data set within the animal model framework and give the results in section 4.1. Indeed, this first application aims at identifying the genomic positions involved in the variability of oil palm production traits. A total of 144 palm trees belonging to the breeding program of PalmElit, a Cirad subsidiary and leading oil palm breeding company (www.palmelit.com), were analyzed. Palm trees were genotyped with 226 molecular markers and 1, 007 IBD matrices were estimated on a grid of 3 centimorgan (cM) along the genome [START_REF] Tisné | Mixed model approach for ibdbased qtl mapping in a complex oil palm pedigree[END_REF]. Each genetic position l is associated to a random e↵ect u l with a variance equal to 2 u l and a correlation matrix A l equal to the identity-by-descent (IBD) matrix (see equation 6). Then, the identification of the genomic positions is equivalent to the variance components selection. Due to the genetic characteristics of the population, i.e. a moderate number of individuals and generations, a subset of 135 genetic positions, spaced 10 cM apart, was considered to avoid a redundant information between consecutive genetic positions. In the next section, we will present the results for the bunch number trait.

As previously explained in subsection 2.1, the animal model can be formulated as follows:

y = µ + u1 1 + • • • + uq q + ", (6)
where µ is an intercept and u l is now assumed to follow a Gaussian distribution Nn(0, A l ), l = 1, . . . , q and l is the standard deviation associated to u l . Finally, the matrix version of the animal model is given by:

y = µ + U + ", ( 7
)
where U is a n ⇥ q-matrix of the concatenation of the random e↵ects U = L q l=1 u l . In a fully Bayesian framework, the intercept µ is supposed to be proportional to one and the residual variance 2 is supposed to follow an inverse-gamma distribution IG(s 2 , r 2 ) (shape and rate parametrization). The Bayesian hierarchical model is presented in web appendix C.

2.3.2

The arabidopsis thaliana dataset. We analyse this data set within the RIS model framework and give the results in section 4.2. Indeed, in this second application, we are interested by disentangling the evolution over time of the complex genetic architecture of shoot growth of Arabidopsis thaliana (L. Heynh). Data consists of leaf compactness phenotypic trait measured over T = 21 time points on n = 358 individuals. We use genetic covariates X containing p = 38 markers [START_REF] Marchadier | The complex genetic architecture of shoot growth natural variation in Arabidopsis thaliana[END_REF][START_REF] Heuclin | Bayesian varying coe cient model with selection: An application to functional mapping[END_REF]. In the RIS model framework, we consider time as the grouping factor. This model is an alternative approach to the usual non-parametric functional mapping [START_REF] Ma | Functional Mapping of Quantitative Trait Loci Underlying the Character Process: A Theoretical Framework[END_REF]. It can be expressed as follows:

yi,t = xi + xi ũt + ↵i + "i,t (8)
where yi,t is the observation of individual i at time t (i = 1, . . . , n and t = t1, . . . , tT ). xi is a (p + 1)-row vector of p genetic markers (constant over time) associated to the i th individual. The first element is fixed to one and is related to the intercept. is a (p + 1)-vector of fixed e↵ects, ũt a (p + 1)-vector of random intercept and slopes e↵ects assumed to follow a Gaussian distribution Np+1(0, ⇤R⇤), where ⇤ is an unknown (p + 1) ⇥ (p + 1)-diagonal matrix of standard deviation and R is an unknown (p + 1) ⇥ (p + 1)-correlation matrix. In this application, R is assumed to be block diagonal where each block is related to one chromosome [START_REF] Ghosh | Bayesian estimation of correlation matrices of longitudinal data[END_REF]. The random intercept is also assumed independent from the random slopes. ↵i is a Gaussian individual random e↵ect (N 0, 2 ↵ ) not subject to selection. "i,t corresponds to the residual part such that "i = ("i,t 1 , . . . , "i,t T ) is distributed from a multivariate Gaussian distribution Nt T (0, 2 e ) where is a tT ⇥tT -correlation matrix of a first-order autoregressive structure with unknown parameter ⇢.

Let y = (y 0 t 1 , . . . , y 0 t T ) 0 be the concatenation of all measurements over time for all individuals where yt = (y1,t, . . . , yn,t) 0 . Since the genetic information varies between individuals but is constant over time, X can be simplified such that X = ( t T ⌦ X) where X is the n ⇥ (p + 1)-matrix containing the p genetic markers (and the intercept) of all individuals. Matrix J is here equal to It T ⌦ n. The random e↵ects design matrix Z can also be simplified as: 1 , . . . , " 0 t T ) 0 is the concatenation of all residuals over time and for all individuals, where "t = ("1,t, . . . , "n,t) 0 is a n-vector of residuals associated to all individuals at time t. " is supposed to follow a Gaussian distribution centered on zero with covariance 2 where = ⌦ In. While the introduction of time random e↵ects allows to capture dependencies between observations within the same time measurement and to model dynamics of genetic e↵ects through the dependence structure. Moreover, introducing a random individual e↵ect combined with a specific residual correlation structure allows to take into account dependencies between measurements over time. Finally, in a fully bayesian framework, the variance associated to the individual random e↵ect 2 ↵ is supposed to follows an inverse-gamma distribution IG(s↵, r↵), the residual variance 2 is supposed to follows an inverse-gamma distribution IG(s 2 , r 2 ) and the autoregressive parameter ⇢ is supposed to follows a uniform distribution U( 1, 1). The Bayesian hierarchical model is presented in web appendix C.

Z = It T ⌦ X. ( 9

Computational aspects of the Bayesian inference

Both applications raise computational challenges mainly due to the number of parameters, dependency structures but also the number of latent variables. In the animal model, the number of parameters is equal to 137 (the intercept, 135 standard deviations and the residual variance) and 19440 unobserved latent variables should be updated (number of elements of U ). In the RIS model, the number of parameters is equal to 214 (39 fixed e↵ects, 39 standard deviations, one individual and one residual variances, one autoregressive parameter and 132 angles) and 1176 unobserved latent variables should be updated (number of elements of U ). While the animal model looks simpler (with a simple additive form) compared to the RIS model (involving complex unknown dependency structures), both complexities are high and the di↵erence between them are not clear. MCMC algorithms have to be appropriately designed for optimization purposes. These optimizations are achieved by reparametrizing standard deviations and by proposing an e cient sampling scheme to avoid inversion of dense posterior covariance matrices.

The first di culty relies on sampling the standard deviations j according to their full conditional distributions. These distributions are proportional to a non-central multivariate folded-Gaussian distribution. Such a distribution does not have a closed form and cannot easily be sampled. To overcome this challenge, following Gelman's work, we propose to reparametrize l as sign(⇠ l )⇠ l where ⇠ l is a parameter which can be positive or negative. It follows that ⇠ l is distributed from a Gaussian distribution (and not from a folded-Gaussian):

⇠ l ⇠ N(0, 2 ⌧ 2 ! 2 l )
. Thus, to sample a standard deviation l from its full conditional distribution p( l |.), we can more simply (i) sample ⇠ l from its full conditional distribution p(⇠ l |.) which is a Gaussian distribution and then (ii) compute l = sign(⇠ l )⇠ l . Demonstrations are provided for both models in web appendix C. High dimensionality causes a second issue to arise. Indeed, at each iteration of the MCMC algorithms, the random e↵ects sampling step involves either the inversion of q n ⇥ n-dense covariance matrices for the animal model (complexity O(qn 3 )) or one tT (p + 1) ⇥ tT (p + 1)-matrix for the RIS model (complexity in O((tT (p + 1)) 3 )). However, these covariance matrices have the form ⌃u = (aA + bI) 1 (after a reparametrization under the RIS model), which is the inverse of the addition of a dense matrix A and a unit identity matrix (a and b are scalars, a, b and A depend on the specific context). This form is very convenient because SVD of the dense matrix A can be used to compute the Cholesky decomposition of ⌃u e ciently. Thus, to sample a random e↵ect u from its full conditional distribution of the form p(u|.) ⇠ N(⌃uh, ⌃u = (aA + bI) 1 ) (where h is a vector), we can (i) compute A = W DW , the SVD of A, where W is an orthogonal matrix of singular vectors and D is a diagonal matrix of singular values, (ii) compute L = W (aD + bI) 1/2 , the Cholesky decomposition of ⌃u, (iii) sample z from a standard Gaussian distribution and then (iv) compute u = L(z + L 0 h). For the animal model context, dense matrices A are known IBD matrices and SVD can be computed only once at the beginning of the algorithm. Thus, the complexity of the sample scheme is in O(qn 2 ). For the RIS model context, matrix A is unknown. However, using specific reformulations of matrices B, ⇤ and ⌦ as Kronecker products, matrix A can be reformulated as a Kronecker product of two matrices and SVD of A can be computed using SVDs of both matrices. Thus, the complexity of the sample scheme is O((p + 1) 3 ) if p + 1 > tT , O(t 3 T ) otherwise. Such algebraic simplifications considerably accelerate MCMC algorithms. The third challenge, specific to the RIS context, is related to the sampling of fixed e↵ects, , and of random individual e↵ects, from their full conditional Gaussian distributions. Again, algebraic simplifications based on reformulations of X, D and ⌦ matrices as Kronecker products allow the simplification of posterior covariance matrices and highly increase the speed of MCMC algorithms.

All these manipulations allow to deal with full conditional posterior distributions and to propose an e cient Gibbs sampler algorithm for the animal model (see web appendix C) or a faster Metropolis-within Gibbs algorithm in the RIS context. A Metropolis-Hasting step is proposed to update angle parameters associated to the correlation matrix between random intercept and random slopes (see web appendix C). All results presented in the next section, are based on 3 MCMC chains initialized at random starting values, each with 50,000 iterations, a burn-in of 10,000 iterations and a thinning of ten. All output statistics are based on the pooled 120,000 posterior samples. The Gelman and Rubin's Potential Scale Reduction Factors (PSRF) statistics [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF] is used to evaluate chains convergence. For standard deviation parameters, estimation is based on the posterior median.

Results

In the next subsections, we show that the fHS prior distribution is e cient to infer and select fixed e↵ects and variance component parameters. As expected, when the number of parameters is large compared to the number of observations (first application), the fC prior does not shrink enough parameters towards zero, leading to clear over-fitting. We highlight that fHS and fSS priors perform similarly to select variance components as it has been shown in the multivariate linear context. In the second application, where the number of parameters is low compared to the number of observations, we show that the three priors perform well and no criteria, based on cross-validation procedure, allows to favour one more than the other.

The oil palm dataset (animal model)

Statistical results. Considering the algorithm does not converge, we adopt a fC prior for standard derivations as an alternative to the commonly used inverse-Gamma prior for variance parameters. The fC prior is not dedicated to selection but should allow for better model regularization than the inverse-Gamma. However, results show that even this prior does not shrink enough towards zero leading to a systematic bias in the estimations (see figure 2), with posterior medians varying around 0.17. The fC prior leads to over-fitting, which is particularly noticeable when analyzing the residual variance: it is estimated around zero (see figure 3) and it has a notable impact on the converge of the Gibbs sampling algorithm by leading to a PSRF greater than 2 for a few continuous parameters. Comparatively, the fHS prior exhibits a very di↵erent behavior. It shrinks towards zero most standard deviations and let some of them be far from zero. Thus, it enables the selection of random e↵ects and improves the MCMC convergence (PSRFs are always close to one for all continuous parameters). In this application, we propose the selection of variance components representing at least 0.05 percent of the total phenotypic variance (0.0023 or equivalently a threshold of 0.048 on standard deviations). This threshold leads to select 10 random e↵ects (see figure 2 andtable 1). The use of the fSS prior, with marginal inclusion posterior probability threshold equal to 0.1, leads to the selection of 7 standard deviations (see figure 2 andtable 1). Six markers are commonly selected by fHS or fSS priors. The selection of variance components is comparable. Such similarities have already been observed in the selection of the fixed e↵ects. Interestingly, the four markers selected using the fHS prior that are not selected using the fSS prior, have also been reported to impact phenotypic variability in di↵erent studies.

Thus, the fHS prior seems to e ciently shrink towards zero the non-relevant random e↵ects while properly estimating relevant parameters. Moreover, it presents better computational performances than the fSS prior. Indeed, computational time for the fHS prior is twice faster than the fSS prior (40 and 80 minutes respectively for 50,000 iterations). Then, the fHS prior should clearly be promoted in a high dimensional quantitative genetic context.

Biological interpretation. We turn to biological interpretations by focusing on the results obtained by the fHS prior. Comparing with the [START_REF] Tisné | Mixed model approach for ibdbased qtl mapping in a complex oil palm pedigree[END_REF] study that analyzed the same data using maximum likelihood ratio tests combined with a forward approach, all but one position identified in the former study were found. Surprisingly, the common positions were all identified at the 0.1% threshold selection, but none for the 0.05% selection. This could be due to the genetic design of the population studied derived from a breeding pedigree with unequal contributions of contrasted genetic groups: among the 144 palm trees, 73% were from La Mé (LM) genetic background, 15% from Yangambi (YBI) and 3% from their combination. Several other studies analyzed both genetic backgrounds with di↵erent genetic designs and common genetic markers. [START_REF] Billotte | Qtl detection by multi-parent linkage mapping in oil palm (elaeis guineensis jacq.)[END_REF], with 25% LM and 25% YBI, found four common positions including two at the 0.05% threshold, [START_REF] Ukoskit | Oil palm (elaeis guineensis jacq.) linkage map, and quantitative trait locus analysis for sex ratio and related traits[END_REF], with 50% YBI, four common positions including two at the 0.05% threshold and [START_REF] Seng | Qtls for oil yield components in an elite oil palm (elaeis guineensis) cross[END_REF], three common positions including two at the 0.05% threshold. The ability of selecting positions corresponding to YBI QTL that were segregating in a minor fraction of the population indicates that the method evaluated in this study performs well even with unbalanced genetic designs and rare allele segregations. This result highlights that a multivariate approach increases the power of detection of subtle e↵ects.

The arabidopsis thaliana dataset (RIS model)

In this second application, markers are labelled by their chromosome numbers and their positions (within the whole dataset of 538 markers) separated by a dash, such that marker 1-2 corresponds to the second position on the first chromosome. This notation was used by [START_REF] Heuclin | Bayesian varying coe cient model with selection: An application to functional mapping[END_REF] and will be used for comparison purposes. We compare our results with those of [START_REF] Heuclin | Bayesian varying coe cient model with selection: An application to functional mapping[END_REF], which used a non parametric functional mapping method, but also with the approach of [START_REF] Marchadier | The complex genetic architecture of shoot growth natural variation in Arabidopsis thaliana[END_REF], which is based on a stepwise strategy. For the three approaches, PSRF statistics of all continuous parameters are lower than 1.1 indicating chains' convergence.

Selection of fixed e↵ects. Fixed e↵ects are considered selected if zero does not belong to their credible intervals, leading to three selected markers (2-32, 2-62 and 5-104). Whatever the prior for variance components, the selection of fixed e↵ects using HS priors performs well and provides the same results (see table 2).

Selection of variance components and impact of the correlation matrix between random e↵ects. Using fC or fHS priors, a threshold representing 0.1 percent of the total phenotypic variance (0.0068, or equivalently a threshold of 0.083 on standard deviations) is used. For the fSS prior, a threshold of 0.5 is considered. While the time random intercept is systematically included in the model, the number of selected random slopes varies according to priors but also to the correlation matrix prior we use between random e↵ects. In particular, selection appears sensitive to the assumption we make on the angle shrinkage global parameter k which plays an important role on the correlation prior distribution (see equation 5). For example, when k tends to infinity (identity case) the number of selected random slopes is equal to 13, 10 and 10 for fC, fHS and fSS priors respectively. When k is fixed to one, numbers increase to 24, 18 and 11 respectively. To choose the most appropriate k value for each prior, a 10 crossvalidation scheme is performed. The log pointwise predictive density (lppd, [START_REF] Gelman | Understanding predictive information criteria for bayesian models[END_REF]), related to k = 1, 3, 5, 7, 10 and for independence assumption, are reported on table 3. Small di↵erences can be observed. For parsimony reasons, random e↵ects are assumed independent. In this example, the fC prior leads to select only few more markers (13) than the fHS or fSS (10) (see figure 4 andtable 2). Di↵erences between priors are less pronounced than in the animal context where the use of the fC prior leads to an estimation of the residual variance close to zero and then to over-fitting problems. Here, residual variance is slightly lower using the fC prior (2), than using fHS or fSS priors (2.5). These di↵erences cannot be used to evidence one prior rather than another. To decide if a prior can be promoted, we compare lppd between models (see table 3). But results are very close and no clear conclusion can be drawn from these results.

Selection of variance components and impact of the residual correlation.

The selection of fixed e↵ects is not impacted by the residual correlation matrix, on the contrary this dependency structure impacts the selection of variance components. Such conclusions have already been observed in the functional mapping context [START_REF] Ma | Functional Mapping of Quantitative Trait Loci Underlying the Character Process: A Theoretical Framework[END_REF][START_REF] Li | Functional mapping of growth and development[END_REF][START_REF] Heuclin | Bayesian varying coe cient model with selection: An application to functional mapping[END_REF]. When we compare selection of random e↵ects taking into account an AR(1) residual correlation structure or assuming independence between residuals, the number of selected markers di↵ers. It considerably increases with the on the oil palm trees dataset. Vertical dotted lines correspond to the selected positions using the fSS prior with posterior marginal probability of inclusion upper than a threshold of 0.1. The horizontal red dashed line corresponds to a threshold of 0.048 which is the root of 0.05% of the response variance. The alternated white and grey areas delimit the 16 chromosomes.

• • ••• • • • • • • • • • • • • • • • • • • • • • • • • • • • •• •• • • • • • • • •• • •• • •• • • • • • •• • •• • •• • • •• • • • •• • • • • ••• • • •• • • • • •• • • • • • • • • • • • • • • • • • • • •• • • • ••• • • • •••• •••• • • • • • •• 0.0 0.2 0.4 0.6 1 -1 1 -4 1 -6 1 -9 1 -1 3 1 -2 0 2 -1 2 -8 2 -1 2 2 -1 7 3 -1 3 -6 3 -1 2 4 -2 4 -7 4 -1 2 4 -1 6 4 -2 0 4 -2 4 4 -2 8 4 -3 5 5 -2 5 -5 6 -1 6 -4 6 -7 6 -1 2 6 -1 4 7 -1 7 -7 7 -9 7 -1 5 8 -2 8 -5 8 -1 3 8 -1 8 8 -2 4 8 -2 9 9 -1 9 -7 9 -1 2 9 -1 7 1 0 -1 1 0 -5 1 0 -8 1 0 -1 1 1 0 -1 5 1 1 -1 1 1 -4 1 1 -7 1 1 -1 4 1 1 -1 6 1 2 -6 1 2 -1 0 1 2 -1 5 1 2 -2 0 1 3 -1 1 3 -7 1 3 -9 1 4 -1 1 4 -7 1 4 -9 1 5 -5 1 5 -1 1 1 5 -1 7 1 5 -2 3 1 6 -2 1
Table 1: Selected standard deviation parameters l for the oil palm trees dataset using folded horseshoe (fHS) and folded spike-and-slab (fSS) priors. independence assumption, leading to potential over-fitting problems. Indeed, the 10 cross-validation lppd, considering an AR(1) residual correlation structure and the fC prior, is equal to 924, while considering an independent residual structure, it is equal to 1148. Here, results are clear and the residual correlation structure has to be included in the model.

Comparison with previous studies. The initial study identified eight markers using the last time measurement and a forward likelihood ratio test approach [START_REF] Marchadier | The complex genetic architecture of shoot growth natural variation in Arabidopsis thaliana[END_REF]. In a recent work, [START_REF] Heuclin | Bayesian varying coe cient model with selection: An application to functional mapping[END_REF] reanalyses this data proposing a non-functional mapping technique combined with group spike-and-slab and taking into account the full phenotypic profile over time. They identified the same eight markers but also highlighted five more e↵ects. In our current analysis, all positions already identified by the previous approaches are selected except two on chromosome three (3-1 and 3-25), compared to Heuclin et al.'s approach. And we select two extra positions (2-47 and 5-46). Moreover, in our approach, decomposing e↵ects as fixed and random allows to more precisely dissociate the type of e↵ects (null, constant or varying e↵ects). For instance, the position 5-104 selected as random e↵ect and varying over time in Marchadier et al. or Heuclin et al. is correspond to the HS-fHS and HS-fC approaches. Vertical red dotted lines correspond to the selected positions using the HS-fSS approach with posterior marginal probability of inclusion upper than a threshold of 50%. The horizontal red dashed line corresponds to a threshold of 0.083 which is the root of 0.1% of the response variance. The alternated white and gray areas delimit the 5 chromosomes.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0
Table 3: Log pointwise predictive density, considering either an unknown RIS correlation matrix with di↵erent fixed shrinkage parameters k, or an identity matrix. using 10 cross-validation, favours the RIS model (-925 for the HS-fC approach with R = I39) to the VCM model (-931).

Conclusion

In this paper, we show that the folded horseshoe prior should be promoted as a prior distribution for regularization in linear mixed models. Based on two real applications, we demonstrate that the folded horseshoe prior seems insensitive to high dimensional problems and leads to unbiased estimation even in low dimension. In the first example, where the number of parameters is close to the number of observations, the folded horseshoe prior shows advantages compared to the folded Cauchy and to the folded spike-and-slab priors. In particular, where the folded Cauchy prior does not allow to shrink parameters towards zero inducing a clear overfitting, the folded horseshoe prior performs well. Compared to the folded spike-and-slab prior, the folded horseshoe prior presents similar e↵ectiveness in terms of selection but a much greater computational e ciency. In the second application, where the number of observations is much greater than the number of parameters, no prior seems to take advantage. Such results observed in multivariate linear regression [START_REF] Van Erp | Shrinkage priors for bayesian penalized regression[END_REF] can then be extended to the linear mixed model framework. However, the folded horseshoe prior does not lead to biased estimations or under or over-fitting of models compared to the two other priors. Be that as it may, we recommend to use local-global priors.

We also propose a polar reparametrization of the model random e↵ect correlation matrix. This approach has received little attention in the past few decades. While Pourahmadi and Wang were the first to develop a prior to generate highdimensional random correlation matrix, Ghosh et al. were the first to infer, in a Bayesian framework, a correlation matrix in a longitudinal context. In this article, we show how this approach can be used to infer RIS correlation matrix. We also show that assuming independence or not can impact variance components selection. However, the number of parameters (angles) is equal to the number of elements of the sub-diagonal correlation matrix. Appropriate priors for the selection of angles such as considered by [START_REF] Ghosh | Bayesian estimation of correlation matrices of longitudinal data[END_REF] should be studied in combination with standard deviations shrinkage priors.

From a biological point of view, in the palm oil context, the folded horseshoe prior allows to identify positions which were segregated in a minor fraction of the population due to the unbalanced genetic design, while the frequentist stepwise selection approach considered by [START_REF] Tisné | Mixed model approach for ibdbased qtl mapping in a complex oil palm pedigree[END_REF] does not. In the Arabidobsis context, as already noticed by [START_REF] Heuclin | Bayesian varying coe cient model with selection: An application to functional mapping[END_REF], we show that a longitudinal approach allows a better detection of relevant markers compared to an approach that analyzes a single time point as proposed by [START_REF] Marchadier | The complex genetic architecture of shoot growth natural variation in Arabidopsis thaliana[END_REF]. Both applications highlight that multivariate approaches increase the statistical power.

Figure 1 :

 1 Figure 1: Prior density distribution of the angle ✓i,j (see equation 5) according to di↵erent values of k: 0, 10 of 100.

)

  Calculation details are presented in web appendix B. Finally, ⇤ is decomposed as It T ⌦ ⇤ and R = It T ⌦ R. P is the matrix that transforms to Vec(⇤) (or equivalently ⇤ = diag( ) and Vec(⇤) = P ). Then, as proposed in section 2.1, this model can be expressed as:y = [ t T ⌦ X, (U 0 ⌦ X)P ]  + D↵ + ".(10) Calculation details are presented in web appendix B. U is a (p + 1) ⇥ tT -matrix of the collection of the tT reparametrized vectors of random intercept and slopes associated to each time U = L t T k u k (U follows a matrix Gaussian distribution MN (p+1)⇥t T (0, R, It T )). D = t T ⌦ In is the design matrix associated to the individual random e↵ects. " = (" 0

Figure 2 :

 2 Figure 2: Posterior median of standard deviation parameters l for folded horseshoe (fHS) and folded Cauchy (fC) priors

Figure 3 :

 3 Figure 3: Posterior density of the residual variance parameter for folded Cauchy (fC), folded horseshoe (fHS) and folded spike-and-slab (fSS) priors on the oil palm trees dataset.

Figure 4 :

 4 Figure 4: Posterior median of standard deviation parameters l on arabidopsis thaliana dataset. Bullets black and blue

Table 2 :

 2 mostly identified as fixed e↵ect (see table 2). Finally, comparing the lppd statistics Selection of fixed e↵ects and scale parameters on arabidopsis thaliana dataset using HS-fC, HS-fHS and HS-fSS approaches. Alternative methods proposed by[START_REF] Marchadier | The complex genetic architecture of shoot growth natural variation in Arabidopsis thaliana[END_REF] and[START_REF] Heuclin | Bayesian varying coe cient model with selection: An application to functional mapping[END_REF] are also indicated.

				Chromosome	1	2	3	4	5	nb
	Marchadier et al. (2019) 1-20	2-62	3-3, 3-91	4-45	5-76, 5-104	8
	Heuclin et al. (2020)	1-1, 1-20, 1-79, 1-97, 1-110	2-62	3-1, 3-25, 3-91	4-45	5-33, 5-76, 5-104 13
	Fixed	e↵ects	HS-fC HS-fHS HS-hSS		2-32, 2-62 2-32, 2-62 2-32, 2-62			5-104 5-104 5-104	3 3 3
	Scale	parameters	HS-fC HS-fHS HS-fSS	1-1, 1-20, 1-79, 1-97, 1-110 1-1, 1-20, 1-79, 1-97, 1-110 1-1, 1-20, 1-79, 1-97, 1-110	2-47, 2-62	3-91 3-91 3-91	4-45 4-45 4-45	5-33, 5-46, 5-76, 5-104 5-33, 5-46, 5-76 5-33, 5-46, 5-76	13 10 10

Received May 2021. Revised --. Accepted --.

Acknowledgments:

F. Mortier and C. Trottier were supported by the GAMBAS project funded by the French National Research Agency (ANR-18-CE02-0025). M. Denis was fully supported by the European Union's Horizon 2020 Research and Innovation programme under grant agreement No 840383. We thank all people from Cirad/PalmElit (France) who planned this trial. We acknowledge P.T. Socfin Indonesia (Indonesia) for planting, observing and collecting data, and authorizing use of the phenotypic data for this study.

Supporting Information

Web appendix A, B and C, referenced in Section 1, 2 and 3, are available with this paper at the Biometrics website on Wiley Online Library. Algorithms for animal and RIS models are available in the R language [START_REF] Team | R: A language and environment for statistical computing[END_REF] on GitHub https://github. com/Heuclin/variance_component_selection. The oil palm dataset is available on request. For the arabidopsis thaliana dataset, the complete phenotypic dataset is freely available at: https://data.inra.fr/dataset.xhtml?persistentId= doi:10.15454/OCOP9B [START_REF] Loudet | Raw phenotypic data obtained on the arabidopsis rils with the phenoscope robots (marchadier, hanemian[END_REF]. The genotypic dataset is freely available at: http://publiclines.versailles. inra.fr/page/8.