Julien Clavel
email: clavel@biologie.ens.frmorlon@biologie.ens.fr

Leandro Aristide

Hélène Morlon

Clavel Julien

Aristide Leandro

Morlon Hélène

A Penalized Likelihood Framework for High-Dimensional Phylogenetic Comparative Methods and an Application to New-World Monkeys Brain Evolution

Keywords: Evolutionary covariances, Phylogenetic signal, Phylogenetic PCA, Regularization, Ridge, LASSO, large p small n, Shrinkage, approximate LOOCV, Generalized Information Criterion, New World monkeys, brain shape

published or not. The documents may come

2017; [START_REF] Bastide | Inference of Adaptive Shifts for Multivariate Correlated Traits[END_REF]). These multivariate phylogenetic models have for example allowed us to analyze the evolutionary integration of traits -i.e. how the evolution of one trait is associated to the evolution of others [START_REF] Revell | Testing quantitative genetic hypotheses about the evolutionary rate matrix for continuous characters[END_REF][START_REF] Revell | Phylogenetic analysis of the evolutionary correlation using likelihood[END_REF][START_REF] Bartoszek | A phylogenetic comparative method for studying multivariate adaptation[END_REF][START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF]. Estimating the correlated evolution of traits is also a necessary step in several multivariate statistical methods accounting for the evolutionary history of species such as canonical correlation Analysis (CCA) [START_REF] Revell | PCCA: a program for phylogenetic canonical correlation analysis[END_REF], principal component analysis (PCA) [START_REF] Revell Liam | Size-correction and principal components for interspecific comparative studies[END_REF], two-block partial least-squares (PLS) [START_REF] Adams | Assessing trait covariation and morphological integration on phylogenies using evolutionary covariance matrices[END_REF], and tests associated with traditional versions of phylogenetic multivariate analysis of variance (MANOVA) (Garland 1992(Garland , 1993)), and multivariate generalized least squares regressions (GLS) [START_REF] Grafen | The phylogenetic regression[END_REF][START_REF] Martins | Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data[END_REF]. A first step common to these multivariate phylogenetic comparative methods involves computing the maximum likelihood estimate of traits evolutionary variance covariance matrices -often referred to as the rate matrix R (e.g., [START_REF] Revell | Testing quantitative genetic hypotheses about the evolutionary rate matrix for continuous characters[END_REF]) -under specific models of phenotypic evolution (often, but not always, the multivariate Brownian model).

In many respects working with multivariate methods is a challenging task as the computational time, memory requirement and the number of parameters to estimate for fitting these models grows up non-linearly with increasing dimensions [START_REF] Ho | A linear-time algorithm for Gaussian and non-Gaussian trait evolution models[END_REF][START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF][START_REF] Goolsby | Rphylopars: fast multivariate phylogenetic comparative methods for missing data and within-species variation[END_REF]). The challenge is even greater when working with highdimensional datasets -where the number of traits (p) is large compared to the number of species (n). In this case, complications arise to compute and/or invert the maximum likelihood estimate (MLE) of traits evolutionary variance covariance matrices. When p approaches n, the MLE of the covariance matrix is no longer a good approximation of the true (population) covariance matrix [START_REF] Stein | Estimation with Quadratic Loss[END_REF][START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Van Wieringen | On the mean squared error of the ridge estimator of the covariance and precision matrix[END_REF].

For instance, the leading eigenvalues are systematically overestimated while the last eigenvalues are underestimated [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF][START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Ledoit | Nonlinear shrinkage estimation of large-dimensional covariance matrices[END_REF]. As a consequence, multivariate phylogenetic comparative methods based on maximum likelihood are no longer reliable and lead to increased model misspecification or loss of statistical power [START_REF] Dunn | Phylogenetic Analysis of Gene Expression[END_REF]Adams 2014aAdams , 2014b;;[START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF][START_REF] Adams | Multivariate Phylogenetic comparative methods: evaluations, comparisons, and recommendations[END_REF]. Things get even worse when p equals or exceeds n. In this case, several eigenvalues equal zero and result in a non-invertible MLE of the covariance matrix, which precludes the use of most phylogenetic comparative methods (Adams 2014a(Adams , 2014c;;[START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF][START_REF] Adams | Multivariate Phylogenetic comparative methods: evaluations, comparisons, and recommendations[END_REF]. In such a situation, a common practice is to reduce the dimensionality of the dataset through PCA and then treat the first axes as independent traits in downstream comparative analysis. But doing so can lead to erroneous conclusions about the underlying evolutionary model [START_REF] Revell Liam | Size-correction and principal components for interspecific comparative studies[END_REF][START_REF] Uyeda | Comparative Analysis of Principal Components Can be Misleading[END_REF]; see also [START_REF] Bookstein | Random walk as a null model for high-dimensional morphometrics of fossil series: geometrical considerations[END_REF] for a similar problem with non independent observations in time-series ). Alternative dimension reduction methods have recently been proposed (Tolkoff et al. 2017), but they are currently limited to the Brownian motion model.

As large p small n datasets become more and more widespread, as for example in geomorphometric and gene expression studies [START_REF] Dunn | Phylogenetic Analysis of Gene Expression[END_REF][START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF][START_REF] Cross | The inside story of 20,000 vertebrates[END_REF], there have been some recent efforts to tackle the issue. First, Adams (2014a, b, c) developed an approach, implemented in the geomorph R package, to compare evolutionary rates across clades, measure phylogenetic signal and perform regression (PGLS) analyses for highdimensional data by transforming the data using analytical results known for the Brownian motion model. This approach is useful, yet restricted to Brownian motion processes (in fact to simple cases of BM that do not require computing a likelihood). More recently, [START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF] showed that Adams' method suffers from low statistical performances, and proposed an alternative approach based on pairwise composite likelihoods, implemented in the phylocurve R package. The approach considers pairs of traits at a time and computes the corresponding pair-level likelihoods under any evolutionary model (such as Brownian motion, Ornstein-Uhlenbeck and Early burst); a pseudo-likelihood is then computed by summing the pair-level likelihoods. This is again a helpful approach, yet using a pseudolikelihood prevents applying classical statistical techniques; for example the parametric evaluation of uncertainty around parameters estimates and model selection procedures based on Information Criterion cannot be performed, and we have to resort to simulations (although analogues to conventional Information criteria have been derived for composite-likelihoods; [START_REF] Varin | An overview of composite likelihood methods[END_REF]. In his paper, [START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF] focused on model comparison, not estimates of the covariance matrix; such estimates can be obtained, but as discussed earlier we expect the obtained covariance to be a poor estimate of the true covariance when the ratio p/n approaches 1 or to be singular when p>n. Moreover, Goolsby's approach has also recently been criticized for being sensitive to potential data transformation such as rigid rotations, which limits its use, for example in the case of geomorphometric data [START_REF] Adams | Multivariate Phylogenetic comparative methods: evaluations, comparisons, and recommendations[END_REF]. Hence, while the two developments above advanced the field by providing some new tools for high-dimensional datasets, they did not directly address the problem of estimating accurate covariance matrices, and thus do not allow extending the full arsenal of available phylogenetic comparative methods to small n large p datasets.

In this article we develop a general penalized likelihood framework for estimating and inverting trait covariance matrices that allows multivariate phylogenetic comparative methods to be extended to high-dimensional data. Penalized likelihood [START_REF] Huang | Covariance matrix selection and estimation via penalised normal likelihood[END_REF][START_REF] Levina | Sparse estimation of large covariance matrices via a nested lasso penalty[END_REF][START_REF] Warton | Penalized Normal Likelihood and Ridge Regularization of Correlation and Covariance Matrices[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF]) is a modern formulation of the regularization/shrinkage approaches commonly used in the statistical literature (Hoerl and Kennard 1970a;[START_REF] Dempster | Covariance Selection[END_REF][START_REF] Friedman | Regularized Discriminant Analysis[END_REF][START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF][START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF][START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF]. The goal of regularization approaches is to constrain the estimates of the parameters, such as the MLE, in order to make them more accurate and to solve ill-posed problems. For instance, the ridge regularization was one of the first methods used in large multiple regression problems to circumvent singularity issues due to over parameterization and to obtain estimates with reduced mean squared errors (Hoerl andKennard 1970a, 1970b). Regularization techniques have been successfully applied in phylogenetics for reconstructing trees [START_REF] Kim | Penalized likelihood phylogenetic inference: bridging the parsimony-likelihood gap[END_REF]), estimating molecular rates and divergence times [START_REF] Sanderson | Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach[END_REF][START_REF] Smith | treePL: divergence time estimation using penalized likelihood for large phylogenies[END_REF], and more recently for estimating lineage-specific evolutionary rates and ancestral states under univariate Brownian models of trait evolution [START_REF] Kratsch | RidgeRace: ridge regression for continuous ancestral character estimation on phylogenetic trees[END_REF] or detecting shifts in the optimum of Ornstein-Uhlenbeck models [START_REF] Khabbazian | Fast and accurate detection of evolutionary shifts in Ornstein-Uhlenbeck models[END_REF]. In all these studies, regularization has proven to be an efficient alternative to Bayesian techniques, generally used to integrate over parameter rich spaces. While the idea of using regularization approaches to address the large p small n problem has been suggested before in another context [START_REF] Dunn | Phylogenetic Analysis of Gene Expression[END_REF], it has not been developed in practice. In the case of covariance matrix estimation, the goal of the penalized likelihood approach is to constrain the MLE to a symmetric positive definite and thus invertible matrix with a simpler structure and improved statistical properties (e.g., correcting for the over dispersion of the eigenvalues or reducing the quadratic loss). This is achieved by adding a penalty term to the likelihood that favors solutions with the required structure and by finding the optimal solution that maximizes the penalized likelihood.

We consider four commonly used penalizations that fall in two broad types: the ridge penalty, also called ܮ ଶ regularization [START_REF] Huang | Covariance matrix selection and estimation via penalised normal likelihood[END_REF][START_REF] Warton | Penalized Normal Likelihood and Ridge Regularization of Correlation and Covariance Matrices[END_REF][START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF], and the LASSO penalty (Least Absolute Shrinkage and Selection Operator), also called ܮ ଵ regularization [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF][START_REF] Levina | Sparse estimation of large covariance matrices via a nested lasso penalty[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Bien | Sparse estimation of a covariance matrix[END_REF]. We propose a simple method for jointly estimating model parameters and the level of penalization. In addition, we show how to perform model selection using Information Criteria [START_REF] Akaike | A new look at the statistical model identification[END_REF][START_REF] Konishi | Generalised Information Criteria in Model Selection[END_REF].

We apply our penalized framework to the estimation of evolutionary trait covariance matrices, model fit and model comparison for the high-dimensional multivariate Brownian (BM), Early-burst (EB), Ornstein-Uhlenbeck (OU) and Pagel's lambda models. We perform a series of simulations to test the performance of our penalized likelihood approach in terms of estimation of the variance covariance matrix and of evolutionary model parameters, as well as model comparison. When alternative approaches were available (i.e. ML when p<n for both parameter estimation and model comparison based on information criteria, and Goolsby's pairwise likelihood when p>n for parameter estimation), we compared the performance of our approach to these alternatives. We also implemented the estimation of ancestral states [START_REF] Martins | Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data[END_REF][START_REF] Cunningham | Reconstructing ancestral character states: a critical reappraisal[END_REF]) and phylogenetic PCA [START_REF] Revell Liam | Size-correction and principal components for interspecific comparative studies[END_REF] under the above-mentioned processes.

We illustrate potential applications of our approach by analyzing an extremely highdimensional 3-D geometric morphometric dataset describing brain shape variation for 48 species of New World monkeys. Previous model-fitting analyses on a subset of PCs axes derived from this dataset, together with disparity-through-time plots, have indicated that brain shape has quickly diversified during the early history of the clade, according to the occupation of different ecological niches (i.e., an adaptive radiation; [START_REF] Aristide | Brain shape convergence in the adaptive radiation of New World monkeys[END_REF]). However, it has been suggested that, under certain conditions, comparative analyses conducted only over the first PC axes of multivariate datasets may tend to artificially favor "early burst" like processes and should instead be ideally performed using fully multivariate approaches that avoid dimensionality reduction steps [START_REF] Uyeda | Comparative Analysis of Principal Components Can be Misleading[END_REF]. Here we apply our penalized likelihood framework to this high-dimensional empirical dataset, estimating the fit of alternative multivariate evolutionary models (BM, EB, and OU) and comparing the results to those obtained previously. Additionally, we utilize the best fitting model parameters to explore the global patterns of covariation among shape variables (i.e. evolutionary integration) derived from a phylogenetic PCA, as well as to obtain an ancestral reconstruction of the New World monkey's brain shape and evolutionary trajectories of its change through time.

Finally, we discuss our results and the broad potential of penalized-likelihood approaches for extending the multivariate comparative toolbox to high-dimensional data.

MATERIALS AND METHODS

Below we detail how to fit a class of multivariate trait evolution models to high-dimensional comparative data using penalized likelihood. These models include the multivariate Brownian model, and all other multivariate models (e.g. EB, OU, and Pagel's λ) for which a single (between traits) evolutionary variance covariance matrix R applies across the phylogeny (R is often noted , the multivariate counterpart of ߪ ଶ , the diffusion parameter of the univariate Brownian model), and all traits share a common phylogenetic variance-covariance matrix C representing the expected scaled variance for each species (a function of the height ݏ above the root of the ith taxa and the parameters of the model) and between species scaled covariances (a function of the shared evolutionary history ݏ between taxa i and j and the parameters of the model ;[START_REF] Felsenstein | Phylogenies and the comparative method[END_REF]Rohlf 2006). For these "simple" models, we can efficiently compute the corresponding log-likelihood using the multivariate normal density in Equation (1). For more complex models with for instance distinct C matrices for each trait or with distinct R matrices for different clades (e.g., the "BMM" model in mvMORPH; [START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF], different algorithms are required for computing the corresponding log-likelihood efficiently (e.g., [START_REF] Freckleton | Fast likelihood calculations for comparative analyses[END_REF][START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF][START_REF] Caetano | ratematrix: an R package for studying evolutionary integration among several traits on phylogenetic trees[END_REF][START_REF] Goolsby | Rphylopars: fast multivariate phylogenetic comparative methods for missing data and within-species variation[END_REF]. The general penalized-likelihood framework presented here can be applied to such models, provided these algorithms are available.

The Likelihood

The log-likelihood of traits is given by the multivariate normal density:

ℒ = - ଵ ଶ ሼ݊ logሺ2ߨሻ + log|| + ݊log|| + trሾ ିଵ ሺ -ߚሻ ் ିଵ ሺ -ߚሻሿሽ (1a)
Where n is the number of species, p is the number of traits, Y is the n by p matrix of trait values, X the design matrix mapping ancestral states to species for each trait (a n dimensional column vector of one in what follows, but this matrix could have different column dimensions in general; e.g. for regression models such as PGLS), ߚ a vector of size p representing the ancestral states for each trait, C is the phylogenetic covariance matrix of size n by n (see above), and R the traits covariance matrix of size p by p (see details [START_REF] Felsenstein | Inferring Phylogenies[END_REF][START_REF] Revell | Testing quantitative genetic hypotheses about the evolutionary rate matrix for continuous characters[END_REF][START_REF] Freckleton | Fast likelihood calculations for comparative analyses[END_REF][START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF]. For a given matrix A, trሺሻ stands for the trace (the sum of the diagonal elements), || for the determinant, ିଵ for the inverse, ் for the transpose and ି் for the transpose of the inverse. The likelihood expression written in Equation (1a) can also be written as a function of a Kronecker product of R and C, and also as a function of the inverse of the traits covariance matrix R, which is commonly referred to as the precision matrix ષ. We report these expressions in Appendix 1, as they are commonly used in the multivariate statistics literature; the expression with ષ is also the one we use to obtain some of our results.

In addition to this classical log-likelihood, we consider the restricted (RE) loglikelihood, given by [START_REF] Felsenstein | Maximum-likelihood estimation of evolutionary trees from continuous characters[END_REF][START_REF] Felsenstein | Inferring Phylogenies[END_REF][START_REF] Freckleton | Fast likelihood calculations for comparative analyses[END_REF]:

ℒ = - ଵ ଶ ሼሺ݊ -1ሻ logሺ2ߨሻ + log|| + ሺ݊ -1ሻlog|| + trሾ ିଵ ሺ -ߚሻ ் ିଵ ሺ -ߚሻሿ + log| ் ିଵ |ሽ (1b)
Estimates based on the unrestricted likelihood are known to be biased [START_REF] Harville | Maximum Likelihood Approaches to Variance Component Estimation and to Related Problems[END_REF], and the restricted likelihood is thus more appropriate for estimating the residuals and variance-covariance parameters of the models [START_REF] Felsenstein | Maximum-likelihood estimation of evolutionary trees from continuous characters[END_REF][START_REF] Felsenstein | Inferring Phylogenies[END_REF][START_REF] Freckleton | Fast likelihood calculations for comparative analyses[END_REF]).

Here we consider both types of likelihoods because there are straightforward extensions of the approach where restricted likelihood won't be appropriate for model comparison (e.g. cases when the number of ancestral states -or fixed-effects -coded in the design matrix X differ across models Freckleton 2012; but see [START_REF] Gurka | Selecting the best linear mixed model under REML[END_REF].

The maximum likelihood estimate (MLE) of the evolutionary covariance matrix R is given by:

= ሺିఉሻ షభ ሺିఉሻ (2)
or by using n-1 in the denominator of Equation ( 2) for the restricted log-likelihood.

Unfortunately when the number of variables p approaches the number of species n, (the sample estimate) does not necessarily provide a good estimation of the true R, and when p exceeds n, is singular, meaning that we cannot compute its inverse nor the log of its determinant, and thus lose the ability to compute the associated log-likelihood.

The Penalized Likelihood

In order to solve the large p small n problem above, we use the penalized likelihood framework, which aims to provide estimators for the covariance matrix R that are well conditioned, positive definite, and thus invertible. We consider four penalties (Table 1): a linear ridge penalty hereafter referred to as the "archetypal ridge penalty", two types of quadratic ridge penalties, and the LASSO penalty. These penalties are motivated by the bias-variance tradeoff and aim to find a compromise between the high-variance and low-bias and a well-behaved low-variance (but higher bias) target matrix T. Typical T matrices used in the literature are the null matrix (a matrix full of zeros), the identity matrix (or a multiple of the identity matrix), and the diagonal matrix composed of the estimated variances for each trait, that is the diagonal elements of [START_REF] Hoffbeck | Covariance matrix estimation and classification with limited training data[END_REF][START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF]. The latter is often considered as the best choice as it has less bias than the other matrices [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Van Wieringen | On the mean squared error of the ridge estimator of the covariance and precision matrix[END_REF] The archetypal ridge estimator is given by (e.g., [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Warton | Penalized Normal Likelihood and Ridge Regularization of Correlation and Covariance Matrices[END_REF]; van Wieringen and Peeters 2016):

ሺߛሻ ௧௬ ௗ = ሺ1 -ߛሻ + ߛ (3)
where ߛ is a regularization parameter (often also called tuning parameter) bounded between 0 and 1 that controls the intensity of the penalty, and T is a target matrix. We took T to be the diagonal matrix composed of the estimated variances for each trait [START_REF] Hoffbeck | Covariance matrix estimation and classification with limited training data[END_REF][START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF]. This estimator reduces to the MLE when ߛ = 0 and to T when ߛ = 1 (a situation in which traits are estimated to be independent when T is diagonal); intermediate ߛ values provide a good tradeoff between bias and variance reduction. We explain later how the value of the regularization parameter is selected. The archetypal ridge estimator results from maximizing the following penalizedlikelihood [START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF], see Appendix 2):

ℒ = - ଵ ଶ ሼ݊ logሺ2ߨሻ + plog|| + ݊log|| + ሺ1 -ߛሻtrሾ ିଵ ሺ -ߚሻ ் ିଵ ሺ -ߚሻሿ + ݊ߛtrሾ ିଵ ሿሽ (4)
More generally, different types of regularizations are obtained by amending a penalty term to the log-likelihood (Eq. 1). The first type of quadratic ridge penalized-likelihood that we consider (hereafter referred to as ridge quad. (var)) is given by (van Wieringen and Peeters 2016):

ℒ = ℒ - ఊ ସ trሾሺ ିଵ -ିଵ ሻ ் ሺ ିଵ -ିଵ ሻሿ (5)
where again we take the target matrix T to be the diagonal elements of , and the regularization parameter ߛ can take any value in ሾ0, ∞ሾ.

Finally, we consider two penalized-likelihoods given by:

ℒ = ℒ - ఊ ସ ‖ ିଵ ‖ (6)
where the penalty ‖⋅‖ = ∑ |⋅| మ ୀଵ is the matrix norm and ߛ can take any value in ሾ0, ∞ሾ.

With q=2 [START_REF] Huang | Covariance matrix selection and estimation via penalised normal likelihood[END_REF][START_REF] Witten | Covariance-Regularized Regression and Classification for High Dimensional Problems[END_REF][START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF], this penalty writes ∑ ∑ ൫ ିଵ ൯ ଶ ୀଵ ୀଵ = trሾ ି் ିଵ ሿ and is equivalent to the quadratic penalty given by Equation ( 5), but with T chosen to be the null matrix. We refer to this penalty as ridge quad. (null). With q=1, we obtain the penalty well known as the LASSO [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Witten | Covariance-Regularized Regression and Classification for High Dimensional Problems[END_REF]. The LASSO also shrinks towards the null matrix, but does so by progressively assigning the least significant covariance terms to zero instead of homogeneously reducing all terms as in the ridge quad. (null).

These penalties have been successfully applied to the estimation of large covariance matrices in other contexts (e.g., [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF][START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Huang | Covariance matrix selection and estimation via penalised normal likelihood[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Witten | Covariance-Regularized Regression and Classification for High Dimensional Problems[END_REF][START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF].

The archetypal ridge estimator was historically considered because efficient methods can be used to compute and solve it [START_REF] Friedman | Regularized Discriminant Analysis[END_REF][START_REF] Hoffbeck | Covariance matrix estimation and classification with limited training data[END_REF][START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF][START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Warton | Penalized Normal Likelihood and Ridge Regularization of Correlation and Covariance Matrices[END_REF][START_REF] Theiler | The incredible shrinking covariance estimator[END_REF]. The quadratic ridge estimator resembles more to the original ridge penalty used in regression analysis and shows improved statistical properties (van Wieringen and Peeters 2016; van Wieringen 2017), but its analytical solution is more computationally demanding. The LASSO exhibits the same interesting shrinkage properties as the ridge penalizations, and in addition provides sparse matrix estimates, which provides an efficient way to perform covariance selection, but it requires even more intensive iterative algorithms [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Witten | New Insights and Faster Computations for the Graphical Lasso[END_REF][START_REF] Bien | Sparse estimation of a covariance matrix[END_REF]. An important difference between the various penalizations is that some of them (e.g. the ridge quad null) are rotation-invariant, meaning that they are robust to data orientation, while others (e.g. the ridge quad. var. and the LASSO) are not (Table 1, [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF][START_REF] Warton | Penalized Normal Likelihood and Ridge Regularization of Correlation and Covariance Matrices[END_REF]Ledoit andWolf 2012, 2015). When working on untransformed data, methods that are not rotation-invariant can be used. However in other cases, for example for variables generated from shape spaces in geomorphometrics which usually depend on the arbitrary choice of the orientation of a baseline shape [START_REF] Rohlf | Shape statistics: Procrustes superimpositions and tangent spaces[END_REF], only methods insensitive to these rotations should be used [START_REF] Rohlf | Shape statistics: Procrustes superimpositions and tangent spaces[END_REF][START_REF] Adams | Multivariate Phylogenetic comparative methods: evaluations, comparisons, and recommendations[END_REF]. The merits of the various types of penalizations are discussed at greater length in our Discussion and summarized in Table 1.

In what follows, we need to compute first and second order derivatives of the penalized-likelihood for the different types of penalizations. These derivations are detailed in Appendix 2.

<TABLE_1>

Solving the Penalized Likelihood

Unlike the archetypal ridge estimator (Eq. 3), which formulation preceded the corresponding penalized-likelihood [START_REF] Warton | Penalized Normal Likelihood and Ridge Regularization of Correlation and Covariance Matrices[END_REF][START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF] Appendix 2), the three other estimators are found by computing the argmax of the various penalized-likelihoods.

The two quadratic ridge estimators (with T diagonal or T null) can be obtained analytically by finding the R matrix for which the first derivative of Equation ( 5) (and Eq. 6 with q=2) equals zero [START_REF] Witten | Covariance-Regularized Regression and Classification for High Dimensional Problems[END_REF][START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF]. For a fixed value of the regularization parameter ߛ in ሾ0, ∞ሾ, the quadratic ridge estimator is given by:

ሺߛሻ ொ௨ௗ௧ ௗ = ቂߛ + ଵ ସ ൫ -ߛ ିଵ ൯ ଶ ቃ ଵ ଶ ⁄ + ଵ ଶ ൫ -ߛ ିଵ ൯ (7)
where is the identity matrix of dimension p.

The LASSO estimator, obtained by solving the penalized likelihood given by Eq. 6

with q=1, cannot be computed analytically. Several algorithms have been proposed to compute it iteratively [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF], such as the "graphical lasso" algorithm implementation in the R package "glasso" [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Witten | New Insights and Faster Computations for the Graphical Lasso[END_REF]. This algorithm provides a sparse estimate of ିଵ rather than . Other algorithms that directly introduce zeros in [START_REF] Bien | Sparse estimation of a covariance matrix[END_REF]) are more computationally demanding and have not been as well tested. As we encountered non-termination convergence issues with the "glasso" package, we instead implemented the GLASSOFAST algorithm proposed by [START_REF] Sustik | GLASSOFAST: An efficient GLASSO implementation[END_REF], which is also generally faster than the one proposed by [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. We made our implementation of this general-purpose algorithm publicly available on github (https://github.com/JClavel/glassoFast) and in an R package called "glassoFast" on the CRAN repository (R Development Core Team 2016).

Estimating the regularization and model parameters

We used cross-validation to obtain the regularization parameter ߛ [START_REF] Hastie | The elements of statistical learning[END_REF]). This parameter cannot be computed jointly with the model parameters by directly maximizing the penalized log-likelihood (Eqs. 4-6), because it would then always be estimated to 0. This is because the true covariance R in Equations (4-6) is generally not available and is instead replaced by a covariance estimated from the data by taking the MLE. Hence, the same data is used both to evaluate the likelihood of each observation and to estimate the covariance matrix. One way to solve this issue is to approach the likelihood by an expression that uses independent data to evaluate the likelihood and to estimate the covariance matrix. Here we use a leave-one-out cross-validation (LOOCV) of the penalized likelihood for which the general idea is to average the likelihoods over each observation (testing samples) with a covariance matrix estimated without the observation (training samples) [START_REF] Hoffbeck | Covariance matrix estimation and classification with limited training data[END_REF][START_REF] Hastie | The elements of statistical learning[END_REF][START_REF] Theiler | The incredible shrinking covariance estimator[END_REF]. We develop a leave-one-out crossvalidation scheme that consists first in directly removing the phylogenetic correlation between species in the data; this avoids defining arbitrary blocks of independent species based on phylogenetic distance (e.g., [START_REF] Roberts | Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure[END_REF]. We note

෩ = ି భ మ , ෩ = ି భ మ
the transformed, decorrelated data, where ି భ మ is the square root of ିଵ . ߚ, the generalized least squares (GLS) estimate of the ancestral states (ߚ = ሺ ் ିଵ ሻ ିଵ ் ିଵ , e.g. Rohlf 2001), is then directly given by ߚ = ൫ ෩ ் ෩ ൯ ିଵ ෩ ் ෩ . The leave-one-out cross-validated penalized log-likelihood is expressed as:

ℒ = - ଵ ଶ ቂ݊ logሺ2ߨሻ + log|| + ∑ ቀ݈݃ห ෩ ሺߛሻ ሺିሻ ห + ݎݐ ቂ ෩ ሺߛሻ ሺିሻ ିଵ ൫ ෩ -෩ ߚ൯൫ܻ ෨ - ୀଵ ܺ ෨ ߚ൯ ் ቃቁቃ, (8)
Where ෩ and ෩ are the p by 1 column vectors made of the i th row of the corresponding matrix (i.e. the row related to taxon ݅; ൫ ෩ -෩ ߚ൯൫ ෩ -෩ ߚ൯

் is thus a p by p matrix) and ෩ ሺߛሻ ሺିሻ is the p by p penalized likelihood estimator computed as described in Eq. 3 for the archetypal ridge, Eq. 7 for the quadratic ridge, and the solution of Eq. 6 (with q=1) for the LASSO, with replaced by:

ሺିሻ = ൫ ෩ ሺିሻ -෩ ሺିሻ ߚ ሺିሻ ൯ ் ൫ ෩ ሺିሻ -෩ ሺିሻ ߚ ሺିሻ ൯ ݊ -1
where the subscript ሺ-݅ሻ refers to matrices computed without the row corresponding to taxon ݅.

We can jointly estimate the regularization parameter and the model parameters by maximizing the LOOCV (Equation 8). Calculating the LOOCV can be very computationally intensive, and so we consider faster implementations (for the archetypal ridge penalty,

Appendix 3) and approximations (for the quadratic ridge and LASSO penalties, Appendix 3).

Also, for the simple models considered here with X a n dimensional column vector of one (corresponding to the case when only one ancestral trait value is reconstructed for each trait, vs. models with specific subclades that can have distinct ancestral trait values for a given trait, such as when there are sudden variations in trait values during the evolutionary history of clades -O' Meara et al. 2006), we can compute the LOOCV directly on the phylogenetic independent contrasts scores that are, in this case, closely related to the residuals ෩ -෩ ߚ in Equation (8) (Appendix 3; see also [START_REF] Stone | Why the phylogenetic regression appears robust to tree misspecification[END_REF][START_REF] Blomberg | Independent contrasts and PGLS regression estimators are equivalent[END_REF]. We estimate the model and regularization parameters by maximizing the LOOCV, considering both the likelihood and the restricted likelihood, using the L-BGFS-B routine implemented in the optim function in R.

We also show in Appendix 2 how to compute standard errors around parameter estimates and provide some examples in the Supplementary Material (at http://datadryad.org; doi:10.5061/dryad.rf7317t).

Model selection

Penalized likelihoods cannot be compared using traditional information criteria such as the Akaike Information Criterion (AIC, [START_REF] Akaike | A new look at the statistical model identification[END_REF][START_REF] Burnham | Model selection and multi-model inference: a practical information-theoric approach[END_REF], but they can be compared using the Generalized Information Criterion (GIC), which is an extension of the AIC to Maximum likelihood-type estimators (or M-estimators, Konishi andKitagawa 1996, 2008). Consider a series of independent data ݔ ଵ , ݔ ଶ , … , ݔ and a fitted model ݂ሺߠ|ݔሻ with ߠ = ߠ ଵ , ߠ ଶ , … , ߠ the parameters of the model; the M-estimator ߠ is given by the solution of an equation of the type ∑ ߰ሺݔ , ߠሻ = 0 ୀଵ . Then, the GIC is defined by:

ܥܫܩ = -2 ∑ log ݂൫ݔ หߠ ൯ ୀଵ + 2trሺ ିଵ ሻ (9) where = - ଵ ∑ డటሺ௫ ,ఏሻ డఏ ୀଵ ห ఏୀఏ and = ଵ ∑ ߰൫ݔ , ߠ ൯ డሺ௫ |ఏሻ డఏ ୀଵ ห ఏୀఏ (10)
Applied to our case, f is the likelihood function and the GIC is obtained by taking ݔ = ෩ -෩ ߚ (i.e., the phylogenetically transformed residuals). In the case of ML estimation, ߰ሺ,ݔ ߠሻ is simply the derivative of the log-likelihood with respect to ߠ. In the case of penalized likelihood estimation, it is the derivative of the penalized log-likelihood. The bias term trሺ ିଵ ሻ in Equation ( 9) can be seen as an effective number of degrees of freedom -a monotonic function of the regularization parameter (e.g., [START_REF] Rondeau | Maximum Penalized Likelihood Estimation in a Gamma-Frailty Model[END_REF]Abbruzzo et al. 2014;[START_REF] Vinciotti | Model selection for factorial Gaussian graphical models with an application to dynamic regulatory networks[END_REF]. It can be computed as ∑ trሺ ିଵ ሻ ఛୀଵ , with [START_REF] Konishi | Information Criteria and Statistical Modeling[END_REF]. In our case we take the parameters of the model to be ߠ ଵ = ષ (to simplify the computation -Appendix 1), ߠ ଶ = ߚ, and ߠ ଷ the parameter of the trait evolution model (e.g. ߙ, r or ߣ). In the case of ML estimation, we compute the derivatives with respect to ષ and ߚ using the well-known closed form derivatives of the log-likelihood of the multivariate normal density (e.g., [START_REF] Mcculloch | Symmetric Matrix Derivatives with Applications[END_REF][START_REF] Anderson | Maximum-likelihood estimation of the parameters of a multivariate normal distribution[END_REF]; Appendix 2). In the case of the penalized likelihood estimation, we derive the derivatives and provide efficient solutions to compute them in Appendix 2. Finally, in order to speed up the computation of the GIC, we approximate the third term in the sum (i.e. the term corresponding to k=3) by the number of parameters of the trait evolution model, which corresponds to the number of degrees of freedom. Model selection with information criteria under REML estimation is often achieved with different formulations for the bias term (see discussions in [START_REF] Gurka | Selecting the best linear mixed model under REML[END_REF]). Here we follow [START_REF] Gurka | Selecting the best linear mixed model under REML[END_REF] and use the restricted log-likelihood with the bias term estimated as above, that is, by accounting for the estimation of all the parameters in the model.

= - ଵ ∑ డటሺ௫ ,ఏሻ డఏ ഓ ୀଵ ห ఏୀఏ and = ଵ ∑ ߰൫ݔ , ߠ ൯ డሺ௫ |ఏሻ డఏ ഓ ୀଵ ห ఏୀఏ

Implementation

We implemented our penalized-likelihood fit of multivariate evolution for the BM, EB, OU and Pagel's lambda model. The implementation is publicly available on github (https://github.com/hmorlon/PANDA/tree/PenalizedMLPhylo) and in the R packages RPANDA [START_REF] Morlon | RPANDA: an R package for macroevolutionary analyses on phylogenetic trees[END_REF], function fit_t_pl) and mvMORPH (Clavel et al. 2015, function mvgls). We allowed (as an option in the function) for the possibility to account for intraspecific variation and measurement errors (they may otherwise bias the model selection, [START_REF] Silvestro | Measurement errors should always be incorporated in phylogenetic comparative analysis[END_REF] by directly estimating them from the data in the fitting process [START_REF] Housworth | The phylogenetic mixed model[END_REF][START_REF] Ives | Within-species variation and measurement error in phylogenetic comparative methods[END_REF][START_REF] Hansen | Interpreting the evolutionary regression : the interplay between observational and biological errors in phylogenetic comparative studies[END_REF], see also the Supplementary Material). Three additional functions compute the corresponding GIC scores (function GIC), perform the reconstruction of ancestral states [START_REF] Martins | Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data[END_REF]Cunningham et al. 1998, function ancestral), and output the phylogenetic PC axes based on the estimated variance-covariance matrix (Revell 2009, function phyl.pca_pl).

Testing the performance of the PL approach with simulations

We use simulations to assess the ability of our penalized likelihood (PL) approach to recover parameter estimates and to select the proper model, for various multivariate evolutionary models. We compare the performance of our approach to that of maximum Likelihood (ML) and pairwise composite likelihood (PCL) [START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF] approaches. We focus on four simple and widely used models of phenotypic evolution: the Brownian Motion (BM), the Early-Burst (EB), the Ornstein-Uhlenbeck (OU) and Pagel's ߣ models (a measure of phylogenetic signal). We follow [START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF] and [START_REF] Bastide | Inference of Adaptive Shifts for Multivariate Correlated Traits[END_REF] and consider the pull matrix of the multivariate OU process (sometimes called ߙ or A; [START_REF] Bartoszek | A phylogenetic comparative method for studying multivariate adaptation[END_REF][START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF] to be diagonal -that is, the coevolution between traits is entirely modeled in R (i.e. in the stochastic, not deterministic part of the process); in addition, the same pull parameter ߙ applies to all the traits. For the OU and EB models we use parameter values (ߙ for OU and r for EB) corresponding to 0 (BM) and 0.5, 1, 3, 5 and 10 half-lives elapsing over the tree height ݐ( ଵ ଶ ⁄ = log ሺ2ሻ ሺݎ ݎ ߙሻሻ ⁄ , held identical across traits. For Pagel's ߣ transformation we use values of: 0, 0.2, 0.4, 0.6, 0.8 and 1 (BM), held identical across traits.

For each model and each parameter set, we simulate 1,000 trait datasets. We fix the size of the phylogeny to n=32 tips and vary the number of traits p from 2 to 50 (p=2, 5, 10, 25, 31, 50). Each simulation consists in the following steps: i) simulate a pure-birth tree with n=32 tips scaled to unit height using the pbtree function in "phytools" [START_REF] Revell | phytools: An R package for phylogenetic comparative biology (and other things)[END_REF], ii) transform the branch-lengths of the tree according to the phenotypic evolution model considered, using the stretching approach [START_REF] Pagel | Inferring the historical patterns of biological evolution[END_REF]; O'Meara 2012; this step is skipped for BM that does not require any transformation), iii) simulate a random R matrix using the approach of [START_REF] Uyeda | Comparative Analysis of Principal Components Can be Misleading[END_REF], which intends to reflect the structure of typical evolutionary rate matrices (e.g., [START_REF] Mezey | The dimensionality of genetic variation for wing shape in Drosophila melanogaster[END_REF]; this matrix is obtained by sampling the eigenvalues from an exponential distribution with rate ߣ = 1/100, assuring a certain degree of correlation between traits, and finally iv) simulate the traits on the transformed tree under a multivariate BM with evolutionary rate matrix R using the mvSIM function in "mvMORPH" [START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF]. All simulations were performed on a Linux platform with R version 3.2.3 (R Development Core Team 2016).

Next, to each simulated dataset, we apply our penalized likelihood approach with the archetypal ridge, the two quadratic ridge, and the LASSO penalties, as well as the PCL approach implemented in the evo.model function from the R package "phylocurve" [START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF]. Under conditions when it is applicable (that is, when p<n), we also apply the simple ML using the evo.model function. In addition, we apply our approximation of the LOOCV to the LASSO and ridge quad. (var) penalties; the approximation of ridge quad. (null) is expected to perform similarly as that of ridge quad. (var). We found that, while the approximations were very good for large n compared to p, they were not when p approached or exceeded n. Thus, we apply the approximations only for p<n; they then provide substantial computational advantage compared to the regular LOOCV (Equation 8). Finally, for each of these analyses, we use both the likelihood and the restricted likelihood formulations. In total, we thus compare the performance of a series of 8 approaches with the likelihood (and the same 8 with the restricted likelihood) when p<n, and a series of 5 approaches with the likelihood (and the same 5 with the restricted likelihood) when p>n. For each simulated dataset and each approach we obtain estimates of the parameters of the phenotypic model. In our optimization, we constrained r from the EB model to be negative, α from the OU process to be positive, and Pagel's λ to be between 0 and 1. We also obtain estimates of the traits evolutionary covariance matrix . For the ML and PCL approaches, the evo.model function directly provides estimates of the parameters of the phenotypic models. These can be used to compute C, which is then plugged into Equation 2 to obtain , which is not returned by default by the evo.model function. We compare to the true (i.e. simulated) matrix using the quadratic loss function (e.g., van Wieringen and Peeters 2016):

Δ ொ ൫܀ , ܀൯ = ฮ܀ ିଵ ܀ -ܫ ฮ ଶ (11)
This loss measure is zero when the matrices are perfectly similar. When p is equal or larger than n, non-regularized estimates of the R matrix (e.g., Eq. 2) are likely to be singular so we compute the estimate of ିଵ in the quadratic loss function using the pseudoinverse routine implemented in the R package "corpcor" [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF]). We also used instead the nearest positive definite matrix using the nearPD function in the "Matrix" R package [START_REF] Bates | Matrix: sparse and dense matrix classes and methods[END_REF] following the approach undertaken by [START_REF] Denton | A new phylogenetic test for comparing multiple highdimensional evolutionary rates suggests interplay of evolutionary rates and modularity in lanternfishes (Myctophiformes; Myctophidae)[END_REF] as well as [START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF]; this resulted in even larger quadratic loss for the non-regularized estimates (results not shown). Finally, we measure our ability to recover the proper model by estimating the proportion of time the generating model was selected by the GIC criterion. This was done for all approaches except the PCL approach for which we cannot compute the GIC. For data simulated under BM, EB and OU, these 3 models are compared. For data simulated under

Pagel's λ model, we compared Pagel's λ and BM (the typical test performed to assess phylogenetic signal).

An empirical example: the evolution of brain shape in New World monkeys

We analyzed the fit of three different evolutionary models (BM, OU and EB) to a comparative dataset of brain shape for 48 New World monkey species in order to better understand the diversification process of brain morphology in the clade. We used brain shape data from [START_REF] Aristide | Brain shape convergence in the adaptive radiation of New World monkeys[END_REF] which consists of a total of 26 anatomical landmark and 373 semi-landmark 3D coordinates (totaling 1197 variables) on digitized endocasts obtained by X-ray computed tomography and micro-computed tomography. Because we are considering 3D geometric morphometric landmark coordinates obtained after Generalized Procrustes

RESULTS

As expected, for all the methods we considered, we found a better performance of the restricted likelihood than of the unrestricted one. All the results shown in the main text are based on the restricted likelihood, while results based on the (unrestricted) likelihood are presented in the Supplementary Material.

Parameter estimation

With the restricted likelihood, parameters of the EB, OU and Pagel's λ models are estimated accurately by all methods for low dimensional traits (p=2-5) (Fig. 1 and S1, S3 and S5).

Estimators based on the unrestricted likelihood are on the other hand systematically biased, regardless of the method used (Fig. S2, S4, and S6), as is expected from unrestricted estimators [START_REF] Harville | Maximum Likelihood Approaches to Variance Component Estimation and to Related Problems[END_REF]: the parameter of the OU process (α) is overestimated, while the parameters of the EB (r) and Pagel (λ) models are underestimated. This bias is higher for the pairwise composite likelihood, and even more so for the maximum likelihood, than for the penalized likelihood approaches.

With both the restricted and unrestricted likelihoods, as the trait dimension p increases, parameters are estimated with increased precision when using the various non-approximated penalized approaches or the PCL approach, while the classical maximum likelihood estimate becomes completely inaccurate (Fig. S1-S6), and it cannot be computed at all when p>n.

The various non-approximated LOOCV approaches perform similarly, with no significant difference between the chosen penalties (Fig. 1 and S1-S6). The approximated LOOCV approaches perform as well as the non-approximated ones for small p values (Fig. 1); they lose precision as p increases (in particular in the case of the LASSO and when p approaches the number of species, Fig. S1-S6), but are still as efficient as the ML or pairwise composite ML when n is still larger than p (i.e., the situations where the use of the LOOCV is generally a major computational bottleneck).

<FIGURE_1>

Estimation of the traits variance-covariance matrix R

The error in the estimate of R, measured with the Q-loss function, is always smaller with the penalized likelihood approaches than with other approaches (Fig. 2 and Fig. S7-14). The Qloss is also slightly smaller when using the restricted versus unrestricted likelihood (Fig. S7- 14). Even with relatively small dimensions (p=2 or p=5), the penalized estimates outperform the PCL and ML estimates (Fig. S7-14). The difference in performance between the PL estimate and the other estimates increases as the dimension p increases (Fig. S7-14).

Among the penalized estimates obtained with the non-approximated LOOCV, the quadratic ridge penalty with a null target is the most efficient; next, the quadratic ridge with the diagonal target matrix composed of the traits' variances and LASSO penalties perform similarly. Finally, the archetypal ridge penalty performs slightly less well, consistent with observations from previous studies [START_REF] Ledoit | Nonlinear shrinkage estimation of large-dimensional covariance matrices[END_REF][START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF].

The penalized estimates obtained with the approximated LOOCV are not as good as those obtained with the non-approximated ones, but they still always outperform the PCL and ML approaches (Fig. 2 and Fig. S7-14).

<FIGURE_2>

Model comparison

Consistent with the other results, model comparison performs better with the restricted than with the unrestricted likelihood (Fig. S15 and S16); in particular, there is a high falsediscovery rate for the Ornstein-Uhlenbeck model with both the ML and the PL unrestricted approaches, which could be linked to the systematic overestimation of ߙ obtained in this case (Fig. S4).

For low trait dimensions (p=2-10), the ability to select the generating model is roughly similar with the penalized and maximum likelihood approaches (Fig. 3 and Figs. S15-17).

However when p>10 the penalized likelihood approaches clearly outperform the maximum likelihood (Fig. 3

<FIGURE_3>

Evolution of brain morphology in New World monkeys

Applying our penalized likelihood framework to the New World monkey brain shape dataset (Fig. 4a,b), we found substantial support for an Early-burst model and these conclusions were robust to the uncertainty in tree topology and dating (Table 2). The BM model was not supported in any of the 100 trees; the OU model was best supported for only 3 trees, and in these cases the difference in GIC with the EB model was very small (Table 2). The average parameter estimate for the Early-burst model ("r", Table 2) describes an evolutionary rate decay with a half-life -the time it takes for the rate to decay half of its initial value -of ~19.4

Ma for an average root height of ~25.5 Ma, indicative of a mild early-burst pattern of brain shape evolution. Considering the first two principal components of the trait evolutionary variance-covariance matrix obtained under the best fitting model (EB), we were able to extract the major patterns of evolutionary integration in brain shape across the clade.

Variation represented by the first PC axis (PC1) mainly reveals a pattern of correlated evolution among the stem, the base, and the parietal regions of the brain (Fig. 4c). Noticeably, variation along PC2 broadly indicates that the temporal, occipital and frontal regions evolve in a concerted fashion but independently from those aspects represented by PC1 (Fig. 4c).

Finally, we generated a hypothesis about the ancestral New World monkey brain shape, and the potential shape evolutionary trajectories leading to two species at opposite extremes in the brain morphospace (Alouatta macconelli and Saimiri boliviensis; Fig. 4b,d).

The ancestral brain shape presents a relatively reduced occipital lobe compared to S.

boliviensis, but a relatively expanded neocortex region and a more flexed base compared to A. macconnelli (Fig. 4d). Moreover, in agreement with the recovered EB model, the reconstructions reveal that for both considered lineages most changes occur at the intergeneric level rather than within each genus (Fig. 4b,d). This suggests that the brain shape of extant species was attained relatively early during the diversification process of the clade.

<FIGURE_4> <TABLE_2>

DISCUSSION

Phenotypic evolution is multivariate by nature given that traits covary due to pleiotropic effects, genetic linkages, developmental and functional constraints or because of correlated selection on multiple traits [START_REF] Felsenstein | Phylogenies and quantitative characters[END_REF][START_REF] Armbruster | Causes of covariation of phenotypic traits among populations[END_REF][START_REF] Walsh | Escape from flatland[END_REF][START_REF] Walsh | Abundant genetic variation + strong selection = multivariate genetic constraints: a geometric view of adaptation[END_REF][START_REF] Armbruster | Integrated phenotypes: understanding trait covariation in plants and animals[END_REF]. In this article we developed a penalized likelihood approach for phylogenetic multivariate comparative methods. We demonstrated through simulations that this approach performs well in estimating parameters of trait evolution models, and performs better than maximum likelihood (when p<n) and current alternatives (when p>n) in estimating model parameters. Even a rough penalization using approximations of the cross-validated log-likelihood to select the regularization/tuning parameter outperforms standard estimates when p<n. The approach allows the use of generalized information criteria (GIC) for model selection, therefore avoiding computationally intensive simulation-based model comparison techniques.

Parameter estimation

Our penalized likelihood approach provides an efficient way to estimate parameters of various multivariate phenotypic evolution models as well as to estimate corresponding evolutionary covariance matrices in high-dimensional datasets (with p large compared to n).

Maximum likelihood estimates are highly biased when p approaches n, and they can no longer be computed when p exceeds n. Penalized estimates generally perform better than maximum likelihood for estimating covariance matrices or their inverses [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF], 2012, 2015; van Wieringen 2017), and we have illustrated here that it is also the case for the evolutionary covariance in phylogenetic comparative methods. To date, the only alternative comparative method is the pairwise composite likelihood (PCL) approach proposed by [START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF]. This method performs as well as the penalized likelihood for estimating the parameters of various classical phenotypic evolution models; however, the penalized likelihood is much more accurate for estimating the evolutionary covariance matrix.

Thus, the penalized likelihood approach developed here is particularly useful for studying the correlated evolution of traits in high-dimensional settings. It should also be useful for extending comparative methods, most of which rely on accurate estimates of the traits evolutionary variance-covariance matrix, to such high-dimensional situations.

Model selection

The penalized likelihood approach also provides an efficient way to compare the relative fit of various multivariate phenotypic evolution models to high-dimensional trait datasets using information criteria. Model selection based on maximum-likelihood is inaccurate when p approaches n, and it is non applicable when p exceeds n. Previous attempts to deal with such situations have used parametric bootstrap to compare pairs of models [START_REF] Goolsby | Likelihood-Based Parameter Estimation for High-Dimensional Phylogenetic Comparative Models: Overcoming the Limitations of "Distance-Based[END_REF]. Parametric bootstrap approaches have appropriate statistical performances (e.g., type I and type II errors; Good 2005) but are computationally intensive and may be impractical when comparing several alternative models. In contrast, within the penalized likelihood framework, multiple model comparison as well as model averaging can be easily and efficiently performed using the Generalized Information Criterion (GIC). Information theoretic approaches (e.g., AIC, BIC, GIC) are commonly used in comparative studies for model selection and/or model averaging and we showed in our simulations that the GIC combined with penalized likelihoods clearly outcompetes information criteria based on ML when p approaches n and still performs well when p exceeds n.

Performing model selection using GIC requires computing a bias correction term for the log-likelihood. We have shown here how to derive this term analytically. However, this might be hard to achieve for more complex models (e.g. models with distinct trait covariation structures and/or evolutionary regimes in distinct subclades). In addition, we have derived a first-order correction term, and higher-order corrections could still improve the statistical performance of the model selection scheme [START_REF] Konishi | Asymptotic theory for information criteria in model selectionfunctional approach[END_REF]. Future developments could use numerical approximations to compute higher level bias correction terms even under complex models (e.g., [START_REF] Ueki | Optimal tuning parameter estimation in maximum penalized likelihood method[END_REF]. Alternatively, parametric bootstrap approaches such as the Bootstrap Information Criterion (or EIC, Efron Information Criterion; [START_REF] Ishiguro | Bootstrapping log likelihood and EIC, an extension of AIC[END_REF][START_REF] Konishi | Information Criteria and Statistical Modeling[END_REF]) could be used in combination with our penalized likelihood framework to do multiple model comparison. This would be more computationally demanding but would automatically achieve higher-order bias correction, and efficient bootstrap strategies have already been developed [START_REF] Ishiguro | Bootstrapping log likelihood and EIC, an extension of AIC[END_REF][START_REF] Konishi | Information Criteria and Statistical Modeling[END_REF][START_REF] Kitagawa | Bias and variance reduction techniques for bootstrap information criteria[END_REF].

Properties of the various penalizations

There are several regularization/shrinkage estimators for covariance matrices (reviewed in [START_REF] Engel | An overview of large-dimensional covariance and precision matrix estimators with applications in chemometrics[END_REF], none of which performs best in all situations; the choice must be guided by practical use (Table 1). The first choice concerns the type of penalization (linear versus quadratic). The archetypal ridge estimator shrinks the eigenvalues of the covariance matrix linearly, and is thus expected to perform well when the eigenvalues are similar in magnitude [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF]. The quadratic ridge estimators shrink the eigenvalues non-linearly, and are thus expected to perform better when the difference in magnitude between the leading and other eigenvalues is high (i.e. when the eigenvalues are dispersed, [START_REF] Ledoit | Nonlinear shrinkage estimation of large-dimensional covariance matrices[END_REF]. As the dispersion of eigenvalues partly reflects the numbers of clusters of integrated traits, their respective numbers of traits, and the average correlation between the traits within clusters (diffuse correlation when eigenvalues are of same magnitude and strong correlation within few clusters when eigenvalues are dispersed), we expect linear estimators to perform well when traits tend to evolve independently and their covariance matrix is not tightly integrated; the LASSO is also expected to perform well when many pairs of traits are not correlated. Quadratic estimators, on the other hand, should perform better when traits are highly correlated within few clusters, as is often observed in empirical data [START_REF] Wagner | On the eigenvalue distribution of genetic and phenotypic dispersion matrices: Evidence for a nonrandom organization of quantitative character variation[END_REF]. Unsurprisingly, in our simulations where eigenvalues are sampled from an exponential distribution, the quadratic ridge estimators clearly outperform the other penalties. In practice, when studying the evolution of traits that are thought to coevolve, quadratic estimators are probably more reliable; however, they are more computationally demanding. There is thus a tradeoff between reliability and computational load (see Table 1). In very high dimensions where extreme values of the regularization parameter may be necessary, archetypal and quadratic ridge estimators are expected to perform equally as they both tend towards the target matrix (van Wieringen and Peeters

2016).

The second choice to make when choosing a penalized likelihood estimator concerns the target matrix. In our study we analyzed the performance of the Quadratic ridge estimator when using either the null or the diagonal matrix of traits' variances. Here we obtained better R estimates in terms of quadratic loss when using the null target matrix, consistently with previous studies (van Wieringen and Peeters 2016); however, these results could vary with the choice of the metric for measuring statistical loss. In addition, each target matrix has a different variance-bias tradeoff, and so other targets may be a better choice depending on the goal of the study [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Lancewicki | Multi-Target Shrinkage Estimation for Covariance Matrices[END_REF].

Maybe a more important thing to consider is whether a rotation-invariant approach is needed. Approaches affecting the eigenvectors of the sample covariance matrix (e.g. Equation 2) are not rotation-invariant, while those that leave the eigenvectors unchanged are [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF][START_REF] Warton | Penalized Normal Likelihood and Ridge Regularization of Correlation and Covariance Matrices[END_REF]Ledoit andWolf 2012, 2015). As discussed above, rotationinvariant approaches should be used when there is no natural orientation of the data, i.e. when an arbitrary choice of orientation is made such as in geomorphometrics [START_REF] Rohlf | Shape statistics: Procrustes superimpositions and tangent spaces[END_REF]. In most cases however, there is a natural orientation of the data. In this case it is not necessary to apply rotation-invariant methods, and this can even yield better estimates [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Warton | Penalized Normal Likelihood and Ridge Regularization of Correlation and Covariance Matrices[END_REF][START_REF] Lancewicki | Multi-Target Shrinkage Estimation for Covariance Matrices[END_REF][START_REF] Van Wieringen | On the mean squared error of the ridge estimator of the covariance and precision matrix[END_REF] Wieringen and Peeters 2016); indeed, shrinking the eigenvectors can improve the estimation of covariances [START_REF] Pourahmadi | Covariance estimation: the GLM and regularization perspectives[END_REF]), which can for example be particularly interesting for ordination analyses [START_REF] Engel | An overview of large-dimensional covariance and precision matrix estimators with applications in chemometrics[END_REF].

Although computationally intensive, the LASSO penalty has the advantage over alternative penalizations that it directly identifies independencies between traits. Indeed, in the estimates of or ିଵ , the least significant entries are forced to zero. Such zeros in the covariance or its inverse correspond, respectively, to marginal independency between traits (pairs of traits are not correlated even when considering potential indirect factors) or conditional independency (a weaker form of independence where the pairs of traits are not directly related only when considering one or several other putative traits) [START_REF] Dempster | Covariance Selection[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Bien | Sparse estimation of a covariance matrix[END_REF]. Identifying such independencies is of major interest, for instance, in the study of patterns of phenotypic integration and modularity [START_REF] Magwene | New tools for studying integration and modularity[END_REF][START_REF] Magwene | Using Correlation Proximity Graphs to Study Phenotypic Integration[END_REF][START_REF] Goswami | Methods for studying morphological integration and modularity[END_REF]. We considered here the graphical-LASSO penalization where the target of the optimization is the precision matrix ିଵ rather than the covariance [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]; an alternative algorithm will need to be used to target the covariance matrix [START_REF] Bien | Sparse estimation of a covariance matrix[END_REF] in order to analyze marginal independencies.

Computational considerations

A limit to the applicability of the penalized likelihood framework can be the computational cost associated with the estimation of the regularized covariance matrix and its inverse when p is large. The least to most computationally intensive of the estimators we have considered here are the archetypal ridge, the quadratic ridge with null target, the quadratic ridge with diagonal matrix of variances, and finally the LASSO estimator. In our simulations, the computational time needed to cross-validate the log-likelihood with the LASSO penalty was order of magnitude higher than with ridge penalties. These computational considerations can guide the choice of the penalization strategy if there is no specific reason to choose one estimator versus another; all the penalized estimates improve upon the standard estimates.

There are several strategies to gain in computational efficiency. We considered approximations of the LOOCV score and have shown that even this rough regularization generally outcompetes maximum-likelihood estimation, but only when p<n. Parallel computing and/or alternative cross validation strategies such as k-fold CV [START_REF] Hastie | The elements of statistical learning[END_REF] could potentially be used. Approaches avoiding cross-validation for estimating the regularization parameter could also be developed. For instance, the regularization parameter could be chosen such that the deviation in likelihood from the maximum is insignificant [START_REF] Meyer | Performance of penalized maximum likelihood in estimation of genetic covariances matrices[END_REF] or by using a priori criteria depending on n, p, and quantities related to the model such as the number of free parameters (e.g., [START_REF] Foygel | Extended Bayesian Information Criteria for Gaussian graphical models[END_REF].

New World monkeys brain shape evolution

We demonstrated the applicability of our framework by analyzing a highlydimensional empirical morphometric dataset describing brain shape variation in New World monkeys. Our analyses recovered as the best fit an Early-Burst model of evolution, a model which is generally considered to better represent the expectations of adaptive radiations than the BM or OU models [START_REF] Harmon | Early bursts of body size and shape evolution are rare in comparative data[END_REF]. [START_REF] Aristide | Brain shape convergence in the adaptive radiation of New World monkeys[END_REF] had already found evidence in favor of a model that supports an adaptive radiation scenario in the brain shape data, but using only the first PC axes. Given the potential issues arising from comparative analyses of PC variables [START_REF] Uyeda | Comparative Analysis of Principal Components Can be Misleading[END_REF], this result was questionable. The fully multivariate approach used here on the complete morphometric dataset thus provides a stronger support to this result. Brain shape evolution in the clade shows a pattern of early diversification consistent with the hypothesis of an adaptive radiation [START_REF] Aristide | Brain shape convergence in the adaptive radiation of New World monkeys[END_REF]). This result is noteworthy because an early-burst pattern is generally hard to detect even on univariate data [START_REF] Slater | Robust Regression and Posterior Predictive Simulation Increase Power to Detect Early Bursts of Trait Evolution[END_REF] and highlights both the singularity of the New World monkey radiation [START_REF] Aristide | Modeling lineage and phenotypic diversification in the New World monkey (Platyrrhini, Primates) radiation[END_REF] and the strength of the multivariate framework we developed.

Additionally, we were able to extract the main patterns of evolutionary integration in brain shape, which indicated that most correlated changes are concentrated along well defined regions of the brain (e.g. frontal, occipital, parietal, etc.) and that some of these regions vary independently from the others. The functional and evolutionary implications of this remains to be assessed but ours constitutes the first attempt to analyze evolutionary integration patterns by using models of evolution in an extremely high-dimensional geometric morphometrics dataset. As stated previously, further developments that allow extracting these patterns from evolutionary covariance matrices of landmark data will help to provide a more detailed picture of brain evolution in the clade.

Finally, using the model parameters estimated for the fully multivariate dataset, we generated ancestral reconstructions for brain shape that can serve as hypotheses regarding various aspects of the evolution of the clade. For example, as brain shape is associated to social group size and other ecological traits in New World monkeys (e.g. [START_REF] Aristide | Brain shape convergence in the adaptive radiation of New World monkeys[END_REF], a careful consideration of reconstructed ancestral shapes may help to understand the evolution of traits that are in general not directly preserved in the fossil record.

Future directions

We see several directions for future developments. The models we have considered here are simple models with a common R matrix shared across the entire clade and a common phylogenetic covariance matrix C shared across traits. More complex models with multiple regime/clade specific covariance matrices could be developed. For such models, joint optimization of multiple penalties (each associated to a given R) could be envisioned (e.g., [START_REF] Guo | Joint estimation of multiple graphical models[END_REF][START_REF] Danaher | The joint graphical lasso for inverse covariance estimation across multiple classes[END_REF]. In even more complex multivariate models described by parameter rich stochastic processes with deterministic parts that are themselves nonindependent (e.g., in multivariate OU processes, if the pull matrix A is not diagonal, [START_REF] Bartoszek | A phylogenetic comparative method for studying multivariate adaptation[END_REF][START_REF] Reitan | Phenotypic evolution studied by layered stochastic differential equations[END_REF][START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF], penalties associated to parameters other than R (e.g. A) could be incorporated to reduce the variance in parameter estimates. Penalized approaches for such models have been shown to outperform alternative approaches in terms of model selection and parameter estimation in a non-phylogenetic context [START_REF] Wang | CGMM LASSO-type estimator for the process of Ornstein-Uhlenbeck type[END_REF]. In a phylogenetic context, the models are already available in a maximum likelihood framework, but they are not applicable with high-dimensional data and are generally estimated with large uncertainty [START_REF] Bartoszek | A phylogenetic comparative method for studying multivariate adaptation[END_REF][START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF].

Second, the regularized estimates provided by the penalized likelihood framework in high-dimensional settings can be used to improve and extend the multivariate comparative toolbox to high-dimensional datasets, beyond fitting multivariate trait evolution models.

Regularized versions of traditional (non-phylogenetic) multivariate statistics (such as the Wilks lambda or Pillai trace used in multivariate regressions and MANOVA) or multivariate classification methods (e.g., discriminant analysis, CCA, etc.) have already been shown to perform well and with increased power when the sample size is small compared to the number of variables [START_REF] Vinod | Canonical ridge and econometrics of joint production[END_REF][START_REF] Friedman | Regularized Discriminant Analysis[END_REF][START_REF] Warton | Penalized Normal Likelihood and Ridge Regularization of Correlation and Covariance Matrices[END_REF][START_REF] Ullah | Regularised Manova for High-Dimensional Data[END_REF][START_REF] Engel | Regularized MANOVA (rMANOVA) in untargeted metabolomics[END_REF]. Regularized estimates of evolutionary covariance matrices should likewise allow developing adequate phylogenetic equivalents of these multivariate statistics in high dimension. Similarly, penalized likelihood approaches for high-dimensional (nonphylogenetic) data have proven to be extremely efficient in estimating missing data [START_REF] Allen | Transposable regularized covariance models with an application to missing data imputation[END_REF]; missing data is also a common feature of high-dimensional phenotypic datasets [START_REF] Clavel | Missing Data Estimation in Morphometrics: How Much is Too Much?[END_REF][START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF][START_REF] Goolsby | Rphylopars: fast multivariate phylogenetic comparative methods for missing data and within-species variation[END_REF], and phylogenetic equivalents of these regularized missing values imputation techniques should solve this challenging task.

Finally, the regularization/shrinkage approaches we describe in this paper can be efficiently extended to Bayesian inferences (e.g., [START_REF] Caetano | ratematrix: an R package for studying evolutionary integration among several traits on phylogenetic trees[END_REF] where shrinkage priors are routinely used as penalties analogues -penalized likelihood estimates are often interpreted as maximum a posteriori (MAP) estimates [START_REF] Green | On the use of the EM for penalized likelihood estimation[END_REF]) -to reduce statistical risks or to obtain well-conditioned symmetric positive definite matrices [START_REF] Daniels | Shrinkage estimators for covariance matrices[END_REF][START_REF] Lu | Modeling between-trial variance structure in mixed treatment comparisons[END_REF][START_REF] Wang | Bayesian Graphical Lasso models and efficient posterior computation[END_REF][START_REF] Khondker | The Bayesian Covariance Lasso[END_REF]. Some of the computational tricks we developed here for computing the LOOCV should be useful for developing efficient proposals and sampling strategies in Bayesian MCMC approaches in high-dimensional parameter space.

CONCLUSION

Phylogenetic comparative methods for high-dimensional datasets are challenging, yet sorely needed. Regularization techniques are a powerful avenue to address this challenge. By providing the tools for properly estimating the evolutionary covariance matrix using penalized likelihood, we open the door to the extension of current multivariate phylogenetic comparative methods to high-dimensional datasets. This should allow addressing longstanding questions related to the correlated evolution of traits.

SUPPLEMENTARY MATERIAL

Our Supplementary Material as well as example codes for running analyses in the paper are available from the Dryad Digital Repository at http://datadryad.org, doi:10.5061/dryad.rf7317t.

APPENDIX 1. THE LIKELIHOOD WRITTEN WITH THE KRONECKER PRODUCT AND AS A FUNCTION OF THE PRECISION MATRIX

Expression with the Kronecker product

The log-likelihood formulation often used in studies of multivariate trait evolution is (e.g., [START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF]; Goolsby 2016):

ℒ = - 1 2 ሼ݊ logሺ2ߨሻ + log|⨂| + ܿ݁ݒሺ -ߚሻ ் ሺ⨂ሻ ିଵ ܿ݁ݒሺ -ߚሻሽ
Where ሺ⊗ሻ is the Kronecker product, and the vec operator stacks the columns of a matrix into a column vector [START_REF] Henderson | Vec and Vech Operators for Matrices, with Some Uses in Jacobians and Multivariate Statistics[END_REF]. This formulation is intuitive, but it is computationally prohibitive. The formulation we use in the paper (Eq. 1a) is equivalent to this formulation, using the following Kronecker product identities: |⨂| = || || where A is a p by p matrix and B is a n by n matrix, ሺ⨂ሻ ିଵ = ିଵ ⨂ ି and ሺݎݐ ሻ = ܿ݁ݒሺሻ ் ሺ ⊗ ሻܿ݁ݒሺሻ (see Magnus and Neudecker 2007, p. 32-35).

Expression with the precision matrix

We can re-express the log-likelihood in Equation (1a) as a function of the precision matrix ષ, defined as the inverse of the covariance matrix [START_REF] Anderson | Maximum-likelihood estimation of the parameters of a multivariate normal distribution[END_REF] as follows:

ℒ = - ଵ ଶ ൛݊logሺ2ߨሻ + log|۱| -݊log|ષ| + ݊tr൫ ષ൯ൟ (A1)
where is the MLE estimate of R given by = ሺିఉሻ షభ ሺିఉሻ (Equation 2). When the inverse of is used as an estimate for the precision matrix ષ, the last term in A1 is equal to ݊ ×  (see also Eq. 23.31 in Felsenstein 2004). Following [START_REF] Anderson | Maximum-likelihood estimation of the parameters of a multivariate normal distribution[END_REF], for computational reasons we use A1 (instead of Equation 1) in Appendices 2 and 3. The penalized estimate of ષ, noted ષ ࢽ in what follows, is obtained as ષ ࢽ = ሺߛሻ ିଵ (e.g., van

Wieringen and Peeters 2016). Likewise, ષ ෩ ࢽ ሺିሻ = ൣ ෩ ሺߛሻ ሺିሻ ൧ ିଵ in the equation of the LOOCV (Eq. 8), hence the LOOCV scores obtained maximizing with respect to ષ are equal to LOOCV scores obtained maximizing with respect to .

APPENDIX 2. COMPUTATION OF THE FIRST AND SECOND ORDER DERIVATIVES OF THE PENALIZED LIKELIHOOD AND APPLICATION TO THE COMPUTATION OF THE GIC

Model comparison using the generalized information criterion (Eq. 9) needs the computation of the first and second order derivatives of the penalized log-likelihood with respect to the parameters estimates (Eq. 10). This is also the case for the LOOCV approximation we describe in Appendix 3. The first and second order derivatives are called the gradient and the Hessian, respectively.

We recall first some matrix notations and differentiation rules (e.g., [START_REF] Dwyer | Some applications of matrix derivatives in multivariate analysis[END_REF]Magnus andNeudecker 1985, 2007). If ߶ is a differentiable scalar function of an ݊ × 1 vector

x, then the 1 × ݊ vector made of the partial derivatives of ߶ with respect to each element of x,

noted ߶ܦሺݔሻ = డథሺ௫ሻ డ௫ is called the derivative of ߶. If ݂ is an ݉ × 1 differentiable vector function of x, the derivative of ݂ is the ݉ × ݊ matrix noted ݂ܦሺݔሻ = డሺ௫ሻ డ௫
where the i th row of the matrix is given by ݂ܦ ሺݔሻ. Finally, if F is a differentiable ݉ ×  function of an ݊ × ݍ matrix X, the derivative of F at X (noted

ୢሺሻ ୢ ) is the ݉ × ݍ݊ matrix ܦሺሻ = ப௩ሼሺሻሽ பሼ௩ሺሻሽ .
We also use the following matrix differentiation rules (e.g., [START_REF] Dwyer | Some applications of matrix derivatives in multivariate analysis[END_REF][START_REF] Magnus | Matrix differential calculus with applications in statistics and econometrics[END_REF]:

ୢ୪୭||܆ ୢ܆ = ܿ݁ݒሺ ିଵ ሻ ் ; ୢ୲୰ሺሻ ୢ܆ = ܿ݁ݒሺሻ ் ; ୢ୲୰൫ ൯ ୢ܆ = ୢ୲୰൫ ൯ ୢ܆ = ܿ݁ݒሺ2ሻ ் ; ୢ ୢ܆ = మ ; ୢ ୢ܆ = మ ; ୢ షభ ୢ܆ = -ሺ ் ሻ ିଵ ⊗ ିଵ
where X and A are matrices of size p by p with A a constant matrix, మ and మ are respectively an identity matrix and a matrix full of zero both of size  ଶ by  ଶ .

Estimation of the first and second order derivatives with respect to ષ for the quadratic penalty

The penalized log-likelihood with quadratic ridge penalization is given by Equation 5( [START_REF] Witten | Covariance-Regularized Regression and Classification for High Dimensional Problems[END_REF][START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF]. Combined with Equation A1we have:

ℒ ൫ષ; ൯ = - ଵ ଶ ቄ݊logሺ2ߨሻ + log|۱| -݊log|ષ| + ݊tr൫ ષ൯ + ఊ ଶ trሾሺષ -ሻ ் ሺષ - ሻሿቅ (B1)
where T is the diagonal matrix corresponding to the diagonal elements of ିଵ . The first derivative of the penalized likelihood w.r.t. the precision matrix is:

ୢ ୢષ ℒ ൫ષ; ൯ = ଶ ܿ݁ݒሼષ ିଵ ሽ ் -ܿ݁ݒ൛ ൟ ் - ఊ ଶ ୢ ୢષ trሾሺષ -ሻ ் ሺષ -ሻሿ൨ (B2)
We can expand the derivative of the penalty term:

d dષ trሾሺષ -ሻ ் ሺષ -ሻሿ = d dષ ሾtrሺષ ் ષሻ -trሺષ ் ሻ -trሺ ் ષሻ + trሺ ் ሻሿ = d dષ trሺષ ் ષሻ - d dષ trሺષ ் ሻ - d dષ trሺ ் ષሻ = ܿ݁ݒሼ2ષሽ ் -2ܿ݁ݒሼ܂ሽ ் = 2ሾܿ݁ݒሼષሽ ் -ܿ݁ݒሼ܂ሽ ் ሿ
Which leads to:

d dષ ℒ ൫ષ; ൯ = ݊ 2 ቂܿ݁ݒሼષ ିଵ ሽ ் -ܿ݁ݒ൛ ൟ ் -ߛሺܿ݁ݒሼષሽ ் -ܿ݁ݒሼ܂ሽ ் ሻቃ ୢ ୢષ ℒ ൫ષ; ൯ = ଶ ܿ݁ݒ൛ષ ିଵ --ߛષ + ߛൟ ் (B3)
The second derivative of the penalized likelihood is:

ୢ మ ୢષ మ ℒ ൫ષ; ൯ = - ଶ ൛ષ ିଵ ⊗ ષ ିଵ + ߛ ൟ (B4)
Estimation of the first and second order derivatives with respect to ષ for the archetypal ridge penalty

The penalized log-likelihood with archetypal ridge penalization is given by Equation 4( Warton 2008, van Wieringen and[START_REF] Van Wieringen | Ridge estimation of inverse covariance matrices from high-dimensional data[END_REF]; expressed as a function of ષ, we have:

ℒ ൫ષ; ൯ = - 1 2 ൛݊ logሺ2ߨሻ + plog|| -݊log|ષ| + ݊ሺ1 -ߛሻtr൫ ષ൯ + ݊ߛtrሾષሿൟ (B5)
The first and second order derivatives with respect to ષ are:

ୢ ୢષ ℒ ൫ષ; ൯ = ଶ ܿ݁ݒ൛ષ ିଵ -ሺ1 -ߛሻ -ߛൟ ் (B6) ୢ మ ୢષ మ ℒ ൫ષ; ൯ = ଶ ሼ-ષ ିଵ ⊗ ષ ିଵ ሽ (B7)
Note that setting B6 to zero, we find that the solution of the penalized-likelihood for the archetypal ridge estimator (Equation 4) is indeed ሺ1 -ߛሻ + ߛ, as given in Equation 3.

Efficient inversion of the Hessian matrix and computation of the GIC

Computing the GIC criterion (Eq. 9) and the approximate LOOCV scores (Appendix 3) involves computing the inverse of the second order matrix derivative of the likelihood (or penalized likelihood) with respect to the covariance or the precision matrix (the Hessian, Eq.

B4 and B7). The Hessian matrix is of dimensions  ଶ ×  ଶ and its direct inversion and storage is computationally infeasible even for a dataset of moderate dimension. For instance with the New-World Monkeys brain dataset studied here the Hessian is a matrix with 2 × 10 ଵଶ entries, which would require more than 15 TB of memory to store (for double precision numerical data). Instead of directly inverting the Hessian, we use identities that reduce this operation to inversions and multiplications of matrices of size p by p.

In the case of maximum-likelihood, the negative Hessian ( in equation 9 and 10) is given by (e.g., [START_REF] Anderson | Maximum-likelihood estimation of the parameters of a multivariate normal distribution[END_REF]. ିଵ is easy to compute using the kronecker identity ሺ ⊗ ሻ ିଵ = ିଵ ⊗ ିଵ . In addition, an efficient computation of the GIC can be done by using the computational trick of Abbruzzo et al. (2014;Eq. 20) -see also [START_REF] Lian | Shrinkage tuning parameter selection in precision matrices estimation[END_REF] and [START_REF] Vujačić | A computationally fast alternative to cross-validation in penalized Gaussian graphical models[END_REF]:

ୢ మ ୢષ మ ℒ൫ષ; ൯ = ଶ ሼષ ିଵ ⊗ ષ ିଵ ሽ
trሺ ିଵ ሻ = ଵ ଶ ∑ ܿ݁ݒ൫ ൯ ் ܿ݁ݒ൛ ષ ષ ൟ ୀଵ - ଵ ଶ ܿ݁ݒ൫ ൯ ் ܿ݁ݒ൛ ષ ષ ൟ (B8) where = ൫ ෩ -෩ ߚ൯൫ ෩ -෩ ߚ൯
் and ષ = ିଵ .

For the archetypal ridge, the Hessian is also of type ⊗ (Eq. B7), and its inverse is easy to compute using the kronecker identity ሺ ⊗ ሻ ିଵ = ିଵ ⊗ ିଵ . The GIC can be computed using:

trሾ ିଵ ሿ = trൣ ିଵ ሺ ் ݊ ⁄ ሻ൧ = ଵ ∑ trൣ ିଵ ܿ݁ݒሺ ሻܿ݁ݒሺ ሻ ் ൧ ୀଵ = ଵ ∑ ܿ݁ݒሺ ሻ ் ିଵ ܿ݁ݒሺ ሻ ୀଵ (B9)
where is the matrix of gradient scores, i.e. the n by  ଶ matrix which i,j element is given by ߰ሺݔ , ൣષ ఊ ൧ ሻ (Eq. 10). With J of the form ⊗ and using the identity ܿ݁ݒሺሻ = ሺ ் ⊗ ሻܿ݁ݒሺሻ [START_REF] Henderson | Vec and Vech Operators for Matrices, with Some Uses in Jacobians and Multivariate Statistics[END_REF][START_REF] Magnus | Matrix differential calculus with applications in statistics and econometrics[END_REF], each element in the sum above can be efficiently computed as:

trሺሺ ⊗ ሻ ିଵ ܿ݁ݒሺሻܿ݁ݒሺሻ ் ሻ = ܿ݁ݒሺሻ ் ሺ ିଵ ⊗ ିଵ ሻܿ݁ݒሺሻ = ܿ݁ݒሺሻ ் ܿ݁ݒሺ ିଵ ିଵ ሻ = ∑ሾ ⊙ ିଵ ିଵ ሿ (B10)
Where ⊙ is the Hadamard (element wise) product.

For the quadratic ridge, the Hessian is of type ൫ ⊗ + ߛ ൯ (Eq. B4) and its inversion can be done by using the eigen-decomposition of Kronecker products identity (e.g., Magnus and Neudecker 2007 p. 33; see also [START_REF] Stegle | Efficient inference in matrix-variate Gaussian models with iid observation noise[END_REF]):

൫ ⊗ + ߛ ൯ = ሺ ⊗ ሻ൫ ⊗ + ߛ ൯ሺ ் ⊗ ் ሻ
where = ் is the eigenvalue decomposition of A, and = ் the eigenvalue decomposition of B ( is the matrix of eigenvectors and is the diagonal matrix of eigenvalues d). It follows that:

൫ ⊗ + ߛ ൯ ିଵ = ሺ ⊗ ሻ݀݅ܽ݃ሺ1. ሺ ⊗ + ߛሻ ⁄ ሻሺ ் ⊗ ் ሻ (B11)
Where ݀݅ܽ݃ሺݒሻ is the diagonal matrix with diagonal elements given by vector v and ./ is the element wise inverse operator. In order to efficiently compute the GIC we make use of the various Kronecker product identities given above, and obtain each element in the last sum of Equation (B9) as:

tr ቀ൫ ⊗ + ߛ ൯ ିଵ ܿ݁ݒሺሻܿ݁ݒሺሻ ் ቁ = ܿ݁ݒሺሻ ் ൫ ⊗ + ߛ ൯ ିଵ ܿ݁ݒሺሻ = ∑ሾܿ݁ݒሺ ் ሻ ⊙ ܿ݁ݒሾ1./ሺ ⊗ + ߛሻሿ ⊙ ܿ݁ݒሺ ் ሻሿ (B12)
Note that in our case (Eqs. B4 and B7), the matrices A and B are identical and further simplify the computations.

Computation of the GIC for the penalized likelihood with LASSO

The LASSO penalty is not easily differentiable [START_REF] Lian | Shrinkage tuning parameter selection in precision matrices estimation[END_REF]Abbruzzo et al. 2014;[START_REF] Vujačić | A computationally fast alternative to cross-validation in penalized Gaussian graphical models[END_REF]. Instead of computing the derivatives of the LASSO penalized likelihood with respect to the covariance matrix, we follow Abbruzzo et al. (2014) and [START_REF] Vujačić | A computationally fast alternative to cross-validation in penalized Gaussian graphical models[END_REF] and use instead the derivatives of the log-likelihood with a sparse estimate assumption as an approximation. This allows us to estimate the bias term in Equation ( 9) efficiently using the following formula:

trሺ ିଵ ሻ ≈ ଵ ଶ ∑ ܿ݁ݒ൫ ⊙ ఊ ൯ ் ܿ݁ݒ൛ષ ఊ ൫ ⊙ ఊ ൯ષ ఊ ൟ ୀଵ - ଵ ଶ ܿ݁ݒ൫ ⊙ ఊ ൯ ் ܿ݁ݒ൛ ષ ఊ ൫ ⊙ ఊ ൯ષ ఊ ൟ (B13)
Where ષ ఊ is the regularized estimate that maximize the LASSO penalized likelihood for a given ߛ, ఊ is an indicator matrix where the entries are 1 when the corresponding entries are non-zero in ષ ఊ and zero otherwise (see the details for the derivation of this approximation in Abbruzzo et al. (2014)), and is as above = ൫ ෩ -෩ ߚ൯൫ ෩ -෩ ߚ൯

் . Note that when every element of the indicator matrix ఊ is equal to one we obtain the formula described above (Eq. B8) for the MLE (i.e., when ߛ = 0).

Estimation of the first and second order derivatives with respect to ߚ

The gradient and Hessian with respect to the ancestral states ߚ are the same whether we consider the log-likelihood or the penalized log-likelihood (except in the case of the archetypal ridge, where they are equal up to ሺ1 -ߛሻ, which does not affect model comparison using GIC), and they can be obtained analytically. From matrix calculus and differential rules (e.g., Magnus and Neudecker 2007, p. 201), we have, if A is a symmetric matrix:

d dߚ ݎݐሾሺ -ߚሻ ் ሺ -ߚሻሿ = d dߚ ሾtrሺ ் ሻ -trሺ ் ߚሻ -trሺߚ ் ் ሻ + trሺߚ ் ் ߚሻሿ = - d dߚ trሺ ் ߚሻ - d dߚ trሺߚ ் ் ሻ + d dߚ trሺߚ ் ் ߚሻ = -ܿ݁ݒሼ ் ሽ ் -ܿ݁ݒሼ ் ሽ ் + ܿ݁ݒሼ ் ߚ + ሾߚ ் ் ሿ ் ሻሽ ் = ܿ݁ݒሼ-2 ் + 2 ் ߚሽ ்
Therefore, we have:

ୢℒ ು ሺఉ;ષሻ ୢఉ = ୢℒሺఉ;ષሻ ୢఉ = - ଵ ଶ ܿ݁ݒ൛-2 ෩ ் ෩ ષ + 2 ෩ ் ෩ ߚષൟ ் = ܿ݁ݒ൛ ෩ ் ෩ ષ -෩ ் ෩ ߚષൟ ் (B14)
From equation (B14) and results from Magnus and Neudecker (2007, p. 198, 360), the second order derivatives are given by:

ୢ మ ℒ ು ሺఉ;ષሻ ୢఉ మ = ୢ మ ℒሺఉ;ષሻ ୢఉ మ = -൛ષ ⊗ ෩ ் ෩ ൟ (B15)
The inverse of the negative Hessian -ቂ ୢ మ ℒሺఉ;ષሻ ୢఉ మ ቃ ିଵ and the GIC criterion (Eqs. 9-10) can then be efficiently computed using the tricks based on kronecker product identities described above for the precision matrix.

Standard error of parameter estimates

Taking ℒ ሺߠሻ = ℒሺߠሻ -ܲሺߠሻ to be the penalized likelihood with respect to parameters ߠ with penalty term P, the standard error (SE) of the parameter estimates ߠ can be computed using either the observed Fisher information matrix -the negative Hessian -or the so-called "sandwich" formula [START_REF] Li | Variable Selection via Nonconcave Penalized Likelihood and Its Oracle Properties[END_REF][START_REF] Rondeau | Maximum Penalized Likelihood Estimation in a Gamma-Frailty Model[END_REF]. A first simple approximation of the asymptotic variance-covariance matrix for the parameter estimates ߠ is given by: ݒܿ

ෞ ൫ߠ ൯ ଵ = ܪ ିଵ ൫ߠ ൯ (B16)
where

ܪ ൫ߠ ൯ = - ୢ మ ℒ ൫ఏ ൯
ୢఏ మ is the negative Hessian of the penalized log-likelihood. Another robust estimator of the covariance matrix for the parameter estimates ߠ is given by the "sandwich" formula [START_REF] Li | Variable Selection via Nonconcave Penalized Likelihood and Its Oracle Properties[END_REF][START_REF] Rondeau | Maximum Penalized Likelihood Estimation in a Gamma-Frailty Model[END_REF]:

ݒܿ ෞ ൫ߠ ൯ ଶ = ܪ ିଵ ൫ߠ ൯ܫ መ ൫ߠ ൯ܪ ିଵ ൫ߠ ൯ (B17)
where ܪ ൫ߠ ൯ is as previously defined and ܫ መ ൫ߠ ൯ = -

ୢ మ ℒ൫ఏ ൯ ୢఏ మ or ܫ መ ൫ߠ ൯ = ݒܿ ቀ ୢℒ൫ఏ ൯
ୢఏ ቁ. The standard errors for the parameter estimates ߠ are then obtained by taking the square root of the diagonal elements of either ݒܿ ෞ ൫ߠ ൯ ଵ or ݒܿ ෞ ൫ߠ ൯ ଶ .

= ݎ + ߶ ݎ ଶ 1 -ݎ߶ = ݎ 1 -ݎ߶ With ݎ = ݒ ் ఊ
ିଵ ݒ (see also Theiler 2012, eq. 16-18).

In addition, we can compute ห ෩ ሺߛሻ ሺିሻ ห in Equation ( 8) using Sylvester's determinant theorem (see section 3.3 in Hoffbeck and Landgrebe 1996):

ห ෩ ሺߛሻ ሺିሻ ห = ห ఊ -ݒݒ߶ ் ห = ห ఊ หሺ1 -ݎ߶ ሻ
Using these expressions, the LOOCV (Eq. 8) of the log-likelihood for the archetypal ridge is given by:

ℒ = - ଵ ଶ ቂ݊ *  logሺ2ߨሻ + log|| + ݊ * logห ఊ ห + ∑ ቀ݈݃ሺ1 -ݎ߶ ሻ + ଵିథ ቁ ୀଵ ቃ (C1)
This expression is an exact equivalent to the LOOCV score in Equation ( 8), yet the computational cost is greatly reduced (from O(np 3 ) to O(p 3 )+O(np 2 ) operations) because the computation of the determinant and inverse of ఊ is done only once [START_REF] Hoffbeck | Covariance matrix estimation and classification with limited training data[END_REF][START_REF] Theiler | The incredible shrinking covariance estimator[END_REF]. Moreover, the determinant and the inverse can be cheaply obtained from the factorization of ఊ (see [START_REF] Clavel | mvmorph: an r package for fitting multivariate evolutionary models to morphometric data[END_REF]Appendix S1). We note that in their original paper, Hoffbeck and Landgrebe (1996; section 3.4) considered a slight variant of the LOOCV . However, in our case, we found Equation C1 to provide better estimates than those obtained with the formula used by [START_REF] Hoffbeck | Covariance matrix estimation and classification with limited training data[END_REF].

with tr ቂ ෩ ሺߛሻ ሺିሻ ିଵ ൫ ෩ -෩ ߚ ሺିሻ ൯൫ ෩ -෩ ߚ ሺିሻ ൯ ் ቃ instead of tr ቂ ෩ ሺߛሻ ሺିሻ ିଵ ൫ ෩ -෩ ߚ൯൫ ෩ - ෩ ߚ൯ ் ቃ where ߚ ሺିሻ = ൫ ෩ ሺିሻ ் ෩ ሺିሻ ൯ ିଵ ෩ ሺିሻ ் ෩ ሺିሻ ; in
To further speed up the computations, we use the phylogenetic independent contrasts scores [START_REF] Felsenstein | Maximum-likelihood estimation of evolutionary trees from continuous characters[END_REF][START_REF] Felsenstein | Phylogenies and the comparative method[END_REF][START_REF] Felsenstein | Inferring Phylogenies[END_REF]) instead of the ෩ and ෩ matrices in Equation ( 8). We note U the ሺ݊ -1ሻ by  matrix of phylogenetic independent contrasts scores. We further note ܸ the column vector of size ݊ containing the variance associated to the estimation of all contrasts and the root state [START_REF] Felsenstein | Maximum-likelihood estimation of evolutionary trees from continuous characters[END_REF][START_REF] Freckleton | Fast likelihood calculations for comparative analyses[END_REF]. Applying the same calculations as above, we can write the LOOCV as:

log ሺℒሻ = - ଵ ଶ ቂ݊ *  logሺ2ߨሻ +  ∑ logሺܸ ሻ ୀଵ + ݊ * logห ఊ ห + ∑ ቀlogሺ1 -ݎ߶ ሻ + ଵିథ ቁ ିଵ ୀଵ ቃ (C2)
where ݎ = ݑ ் ఊ ିଵ ݑ with ఊ = ݊ * ߶ + ߛ and ݑ the column vector made of the ݅ ௧ row of matrix U and ߶ = ଵିఊ * ିଵ

. Here is simply computed as * [START_REF] Felsenstein | Inferring Phylogenies[END_REF]). Note that because the contrasts scores have a zero mean, ఊ and ߶ are different from those used for Equation (C1) above [START_REF] Theiler | The incredible shrinking covariance estimator[END_REF]. For REML estimation, ݊ * = ݊ -1 and the variance term associated to the root state in V is omitted.

Approximation of the LOOCV with the quadratic ridge penalty

ℒ ൫ષ ࢽ ; ൯ ≈ - ଵ ℒ൫ષ ࢽ ; ൯ - ଵ ଶ ∑ ቄ ୢ൫ ;ષ ࢽ ൯ ୢષ ܿ݁ݒ ቀષ ෩ ࢽ ሺିሻ -ષ ࢽ ቁቅ ୀଵ (C3)
where ݃ is the (scaled) penalized log-likelihood, which, for the quadratic ridge penalty, is defined as:

݃൫ , ષ൯ = log|ષ| -tr൫ ષ൯ - ߛ 2 trሾሺષ -ሻ ் ሺષ -ሻሿ
We can approximate the term ܿ݁ݒ ቀષ ෩ ࢽ ሺିሻ -ષ ࢽ ቁ in Equation (C3) using the first order Taylor expansion of the function Vujačić et al. 2015, 4.1):

ୢ൫ ;ષ ࢽ ൯ ୢષ around ൫ , ષ ࢽ ൯, evaluated at ቀ ሺିሻ , ષ ෩ ࢽ ሺିሻ ቁ ( Lian 2011;
d݃ ቀ ሺିሻ , ષ ෩ ࢽ ሺିሻ ቁ dષ ≃ ቊ d݃൫ , ષ ࢽ ൯ dષ ቋ + d ଶ ݃൫ , ષ ࢽ ൯ dષ ଶ ܿ݁ݒ ቀષ ෩ ࢽ ሺିሻ -ષ ࢽ ቁ + d ଶ ݃൫ , ષ ࢽ ൯ dષd ܿ݁ݒ൫ ሺିሻ -൯ Since ୢ൫ ,ષ ࢽ ൯ ୢષ
= 0 (it is the derivative -Eq. B3 -of the penalized-likelihood ݃ evaluated at the maximum penalized-likelihood estimates and ષ ࢽ ) and equivalently

ୢቀ ሺషሻ ,ષ ෩ ࢽ ሺషሻ ቁ ୢષ = 0,
when ሺିሻ and ષ ෩ ࢽ ሺିሻ are sufficiently close from and ષሺߛሻ it follows that:

ܿ݁ݒ ቀષ ෩ ࢽ ሺିሻ -ષ ࢽ ቁ ≃ -ቀ ୢ మ ൫, ષ ࢽ ൯ ୢષ మ ቁ ିଵ ୢ మ ൫, ષ ࢽ ൯ ୢષୢ ܿ݁ݒ൫ ሺିሻ -൯ (C4)
This approximation is accurate when ሺିሻis sufficiently small.

Therefore, the term ୢ൫ ,ષ ࢽ ൯ ୢષ ܿ݁ݒ ቀષ ෩ ࢽ ሺିሻ -ષ ࢽ ቁ on the right hand side of Equation (C3) can be approximated by:

d݃൫ , ષ ࢽ ൯ dષ -ቆ d ଶ ݃൫ , ષ ࢽ ൯ dષ ଶ ቇ ିଵ d ଶ ݃൫ , ષ ࢽ ൯ dષd܀ ܿ݁ݒ൫ ሺିሻ -൯൩
From Appendix 2 we have Using the eigen-decomposition of ષ ࢽ ିଵ = ் and similar calculus as above (see Appendix

2), we find that this term is equal to:

∑ ቂܿ݁ݒ൫ ் ൫ષ ࢽ ିଵ -൫ -ߛ൯ -ߣષ ࢽ ൯൯ ⊙ ቀܿ݁ݒሾ1./ሺ݀ ⊗ ݀ + ߛሻሿ ⊙ ܿ݁ݒ൫ ் ൫ - మ ሺିሻ ൯൯ቁቃ (C5)
Finally, after simplification we obtain the first order Taylor approximation of the LOOCV score for the penalized log-likelihood with the quadratic ridge penalty (see also Eq. 6 in [START_REF] Vujačić | A computationally fast alternative to cross-validation in penalized Gaussian graphical models[END_REF]:

ℒ ≈ -ଵ ℒ൫ષ ࢽ ; ൯ + 1 2݊ሺ݊ -1ሻ ⁄ ∑ ቄ∑ ቂܿ݁ݒ൫ ் ൫ષ ࢽ ିଵ -൫ -ߛ൯ -ߛષ ࢽ ൯൯ ⊙ మ ୀଵ ቀܿ݁ݒሾ1./ሺ݀ ⊗ ݀ + ߛሻሿ ⊙ ܿ݁ݒ൫ ் ൫ -൯൯ቁቃቅ (C6) Further, we note that once ષ ࢽ ିଵ (= ሺߛሻ) has been computed (using Eq. 7), ષ ࢽ can be easily obtained from the eigen-decomposition of ષ ࢽ ିଵ or by noticing that ષ ࢽ = ଵ ఊ ൣષ ࢽ ିଵ -ሺ -ߛሻ൧

(van Wieringen and Peeters 2016; section 3.1).

Efficient approximation of the LOOCV with the LASSO penalty

For the LASSO penalty we use the LOOCV approximation of the multivariate normal loglikelihood with the assumption of a sparse estimate for ષ (i.e., ષ = ષ ࢽ ௌௌை ) given by Vujačić et al. (2015 -Eq. 6 and 9): LOOCV approaches are not reliable (and thus not shown), and the maximum likelihood approach is not applicable. For each trait model, we report results for only 4 of the 6 simulated parameters for presentation purposes.

ℒ ≈ - 1 ݊ ℒ൫ષ ࢽ ; ൯ + 1 2݊ሺ݊ -1ሻ ⁄ ܶ ୀଵ with ܶ = ∑൛൫ષ ࢽ ିଵ -൯ ⊙ ఊ ⊙ ષ ࢽ ൛൫ -൯ ⊙ ఊ ൟષ ࢽ ൟ (C7)
283x175mm (300 x 300 DPI) NO ++++ ++++ * Computational complexity refers to both the computation of the solution of the penalized likelihood and the LOOCV.** Flexibility refers to the robustness of the penalization procedure with respect to structure in the data ; more confidence in the estimates (and thus also in model selection) will be provided by more flexible penalizations, at the expense of computational efficiency. *** In some LASSO implementations, the diagonal elements of (or ିଵ ) are not penalized and (or ିଵ ) is shrunk towards ݀݅ܽ݃ሺሻ (or ݀݅ܽ݃ሺ ିଵ ሻ) which can be then seen as the target matrix (e.g., [START_REF] Bien | Sparse estimation of a covariance matrix[END_REF].

  and Figs.S15-17): the power to detect the true model increases with increasing trait dimensions for the penalized likelihoods while maximum likelihood fails to detect the generating model when p is close to n (i.e., p=31; Fig.3and Fig.S15-S17). The different types of penalizations perform similarly well (Fig.3and Fig.S15-17). When models are fitted using the approximated LOOCV, model misspecification increases when p is close to n, but not as badly as when models are fitted using maximum likelihood (Figs. S17).

Figure 1 .

 1 Figure 1. Performance of penalized likelihood for the estimation of model parameters and comparison with other approaches. Estimates of the parameters of the OU, EB and Pagel's λ models obtained with the various (restricted) penalized likelihood approaches, the (restricted) pairwise composite likelihood approach, and the (restricted) maximum likelihood. Boxplots represent the median, first and 3rd quartile and range of estimates obtained over 1000 simulated datasets. The red line represents simulated parameter values. Left panels (p=5) illustrate results with p<n (n=32), and right panels (p=50) illustrate results with p>n. Results with other p values and unrestricted likelihoods are presented in Figure S1-S6. For p>n, the approximatedLOOCV approaches are not reliable (and thus not shown), and the maximum likelihood approach is not applicable. For each trait model, we report results for only 4 of the 6 simulated parameters for presentation purposes.

Figure 2 .Figure 3 .

 23 Figure 2. Performance of penalized likelihood for the estimation of the variance covariance matrix and comparison with other approaches. Quadratic loss between the simulated and estimated traits evolutionary variance-covariance matrix R obtained with the various (restricted) penalized likelihood approaches, the (restricted) pairwise composite likelihood approach, and the (restricted) maximum likelihood. Boxplots represent the median, first and 3rd quartile and range of Q loss values obtained over 1000 simulated datasets. Left panels (p=5) illustrate results with p<n (n=32), and right panels (p=50) illustrate results with p>n. Results with other p values and unrestricted likelihoods are presented in Figures S7-S10. 289x177mm (300 x 300 DPI)

  Wieringen and Peeters 2016, but see Lancewicki and Aladjem 2014).

  Previous authors have shown that it is possible to efficiently approximate the LOOCV score for the multivariate normal log-likelihood and used these results to approximate the LOOCV score for the graphical LASSO[START_REF] Lian | Shrinkage tuning parameter selection in precision matrices estimation[END_REF][START_REF] Vujačić | A computationally fast alternative to cross-validation in penalized Gaussian graphical models[END_REF]; see below Efficient approximation of the LOOCV with the LASSO penalty). Building on these previous developments we give here the derivation of the approximate LOOCV score for the multivariate normal penalized likelihood with the quadratic ridge penalty.From[START_REF] Lian | Shrinkage tuning parameter selection in precision matrices estimation[END_REF] Appendix A), and Vujačić et al. (2015; section 4.1), the LOOCV score for the penalized likelihood can be approximated by:

Table 1 .

 1 Summary of various penalized likelihood approaches, associated target matrices, and their properties. The four penalties evaluated in the paper are highlighted in bold. I is the identity or multiple of the identity matrix. V is a diagonal matrix with distinct diagonal values, such as the estimated variances for each trait. N is the null matrix (a matrix full of zero).

	Penalty	Target matrix	Rotation-invariance Computational complexity* Flexibility**
	ridge archetypal	I	YES	+	+
		V	NO	+	++
	ridge quadratic	N	YES	++	++
		I	YES	+++	++
		V	NO	+++	+++
	LASSO	N ***

Table 2 .

 2 Support for BM, OU and EB models for the evolution of New World monkeys' brain shape over 100 trees from the Bayesian posterior distribution. The EB model is preferred(lowest GIC value) in 97% of the trees. The Brownian parameters are given in the multidimensional matrix (see also Figure4).

	Model	GIC (mean±2sd)	∆GIC (2.5%-97.5% range) % trees preferred Parameters* (mean±2sd)
	BM	-726505±410	[4.95 -258]	0	-
	OU	-726504±409	[1.03 -260.4]	3	ߙ = 9.34e-5±4.8e-4
	EB	-726618±379	[0 -0.30]	97	r = -0.91±0.69

*Downloaded from https://academic.oup.com/sysbio/advance-article-abstract/doi/10.1093/sysbio/syy045/5040209 by ENS Paris Biologie, Julien Clavel on 19 June 2018

Downloaded from https://academic.oup.com/sysbio/advance-article-abstract/doi/10.1093/sysbio/syy045/5040209 by ENS Paris Biologie, Julien Clavel on 19 June 2018

ACKNOWLEDGEMENTS

The authors wish to thank Ivan Vujačić and Wessel N. van Wieringen for precisions on the estimation of the regularization parameter in their papers. We thank Eric W. Goolsby for assistance with the R package "phylocurve" and Mátyás A. Sustik for supplying source codes for the graphical LASSO. Finally, we thank Renske Gudde, Eric Lewitus, Odile Maliet, Marc Manceau, Olivier Missa, Benoît Perez, Guilhem Sommeria-Klein, Thomas Near, Luke Harmon, Michael Landis and two anonymous reviewers for helpful comments on the manuscript.

FUNDING

This work was supported by European Research Council grant ERC 616419-PANDA (to HM).

Analysis (see details in [START_REF] Aristide | Brain shape convergence in the adaptive radiation of New World monkeys[END_REF], model fit must be performed with a rotationinvariant procedure [START_REF] Rohlf | Shape statistics: Procrustes superimpositions and tangent spaces[END_REF][START_REF] Adams | Multivariate Phylogenetic comparative methods: evaluations, comparisons, and recommendations[END_REF]. We applied the fit_t_pl function to the data, using the rotation-invariant "ridge quadratic null" penalty and accounting for intra-specific variation and measurement errors in the model fit by setting the option SE to TRUE. The relative fit of the three models was assessed using the GIC criterion with the GIC function. We considered phylogenetic uncertainty by performing the analyses over a sample of 100 fossil-calibrated trees from the Bayesian posterior distribution of the phylogenetic analysis described in [START_REF] Aristide | Modeling lineage and phenotypic diversification in the New World monkey (Platyrrhini, Primates) radiation[END_REF] (pruned to the 48 species considered here).

Additionally, we studied the patterns of evolutionary integration and modularity (i.e., the correlated/independent change of traits through evolution; Klingenberg and Marugán-Lobón 2013) in brain shape by analyzing the estimated evolutionary variance-covariance matrix derived from the model fitting analysis. We performed a principal component analysis of this matrix using the phyl.pca_pl function, which allowed us to extract evolutionary independent axes of correlated change among variables that broadly represent evolutionary integration patterns across the structure [START_REF] Klingenberg | Morphological integration and developmental modularity[END_REF]). Finally, with the ancestral function, we utilized the parameters derived from the evolutionary model that best explained our data to obtain brain shape reconstructions through time that illustrate the hypothetical evolutionary changes leading to the two most different species in the dataset. To generate 3D model visualizations of the reconstructed shapes, a 3D surface model of the sample's mean shape was warped to each target ancestral shape by aligning the reconstructed and mean shape landmark data and using a TPS interpolation [START_REF] Wiley | Evolutionary morphing[END_REF], as implemented in the "Morpho" package for R [START_REF] Schlager | Morpho and Rvcg --Shape Analysis in {R}[END_REF]. Color maps depicting shape changes were generated by computing the radial distance between each pair of corresponding vertices in two given 3D surfaces.

When computing standard errors around the precision matrix ષ (or the covariance R), there are ሺ + 1ሻ/2 parameter estimates (not  ଶ ), because ષ (or R) is symmetric. Therefore, instead of computing the derivatives with respect to ષ, we compute the derivatives with respect to vech(ષሻ, where vech is an operator that stacks the columns of a square matrix starting at its diagonal elements into a column vector [START_REF] Henderson | Vec and Vech Operators for Matrices, with Some Uses in Jacobians and Multivariate Statistics[END_REF]. We have:

where , known as the "duplication matrix", is a  ଶ by ሺ + 1ሻ/2 matrix with indicator entries 0 or 1 that satisfies ܿ݁ݒሺሻ = ܿ݁ݒℎሺሻ [START_REF] Henderson | Vec and Vech Operators for Matrices, with Some Uses in Jacobians and Multivariate Statistics[END_REF][START_REF] Moneta | Implicit construction of McCulloch's G matrix for the numerical evaluation of Fisher information matrixes[END_REF]. For the archetypal ridge for example, from B7 we deduce (e.g., McCulloch 1982):

Codes for this calculation and the associated computation of standard error for the precision matrix are provided as an illustrative example in our Supplementary Material.

APPENDIX 3. SPEEDING-UP THE COMPUTATIONS OF THE LEAVE-ONE-OUT CROSS-VALIDATION (LOOCV)

The ordinary approach for computing the LOOCV score (Eq. 8) requires the inversion of n matrices of size p by p for a given value of the regularization parameter, which is computationally intensive, especially for large n. Here we develop approaches for reducing computation time. In the case of the archetypal ridge penalty, we use an analytical trick and illustrate the use of independent contrasts scores instead of the ෩ and ෩ matrices in Equation ( 8). In the case of the quadratic ridge and LASSO penalties, we derive an approximation of the LOOCV score.

Exact and efficient computation of the LOOCV with the archetypal ridge penalty

The LOOCV of the penalized log-likelihood (Eq. 8) with the archetypal ridge estimator, target matrix T, and regularization parameter ߛ can be efficiently computed using the rank-one update trick proposed by [START_REF] Hoffbeck | Covariance matrix estimation and classification with limited training data[END_REF]section 3.3-3.4). These authors showed that we could obtain ሺିሻ required to compute ෩ ሺߛሻ ሺିሻ in Equation ( 8) using a rank-one update of the covariance matrix :

with ݊ * = ݊ for the maximum likelihood estimator and ݊ * = ݊ -1 for the restricted maximum likelihood (see also [START_REF] Theiler | The incredible shrinking covariance estimator[END_REF]eq. 50). Hence, we can re-express the archetypal ridge estimate as:

෩

We deduce: