Jocelyn Chauvet
email: jocelyn.chauvet@umontpellier.fr

Catherine Trottier
email: catherine.trottier@univ-montp3.fr

Xavier Bry
email: xavier.bry@umontpellier.fr

Component-based regularisation of multivariate generalised linear mixed models

Keywords: generalised linear regression, supervised components, random effects, structural relevance

We address the component-based regularisation of a multivariate Generalised Linear Mixed Model (GLMM) in the framework of grouped data. A set Y of random responses is modelled with a multivariate GLMM, based on a set X of explanatory variables, a set A of additional explanatory variables, and random effects to introduce the within-group dependence of observations. Variables in X are assumed many and redundant so that regression demands regularisation. This is not the case for A, which contains few and selected variables. Regularisation is performed building an appropriate number of orthogonal components that both contribute to model Y and capture relevant structural information in X. To estimate the model, we propose to maximise a criterion specific to the Supervised Component-based Generalised Linear Regression (SCGLR) within an adaptation of Schall's algorithm. This extension of SCGLR is tested on both simulated and real grouped data, and compared to ridge and LASSO regularisations. Supplementary material for this article is available online.

Introduction

In the framework of regression models on a large number of explanatory variables with redundancies and collinearities, the search for a reduced number of relevant dimensions to model responses has been an ongoing research over the last decades. In particular, the case where the explanatory variables outnumber the observations tends to be a new standard.

Generalised Linear Models (GLMs) are the most widely used regression models, because they are easy to interpret and address a very large scope of applications with a variety of response distributions. For instance, Epidemiology, Biology and Social Sciences need to model binary outcomes, count data and survival times. All these fields often have to deal with both grouped data and multivariate responses combining variables of different types (e.g. one binary and another Poisson). In this work, we particularly aim at modelling abundances of several tree genera on plots of land grouped into forest concessions, using multiple redundant explanatory variables.

As far as dimension-reduction is concerned, two main approaches have been developed.

The first one is variable-selection, whereas the second one builds components, i.e. linear combinations of the explanatory variables, which synthesise the useful part of their information. As far as variable-selection is concerned, the most popular method is currently the LASSO, introduced by [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF], which combines the likelihood with a penalty based on the L 1 -norm of the coefficient vector. LASSO is one of the penalty-based regularisation methods, as are also ridge [START_REF] Hoerl | Ridge Regression: Biased Estimation for Nonorthogonal Problems[END_REF] and elastic-net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. This LASSO selection approach has proved efficient to explain the phenomenon of interest when some of the explanatory variables are the "true" ones, surrounded by a high number of irrelevant others. Nevertheless, it may be very unstable and helpless when the true explanatory dimensions are latent and indirectly measured through highly correlated proxies. This is where the component-based approach turns out to be useful. [START_REF] Bry | Supervised component generalized linear regression using a PLS-extension of the Fisher scoring algorithm[END_REF] have developed a new methodology named Supervised Component-based Generalised Linear Regression (SCGLR), later extended and refined in [START_REF] Bry | Extending SCGLR to multiple regressor-groups: The Theme-SCGLR method[END_REF][START_REF] Bry | The Multiple Facets of Partial Least Squares and Related Methods, chapter Supervised Component Generalized Linear Regression with Multiple Explanatory Blocks: THEME-SCGLR[END_REF][START_REF] Bry | Component-based regularisation of a multivariate GLM with a thematic partitioning of the explanatory variables[END_REF]. As in any PLS-type method, the construction of components in SCGLR is guided both by the correlation-structure of variables in the explanatory space and by the prediction quality of the responses. Nevertheless, unlike PLS, SCGLR involves a general and flexible criterion allowing to specify the type of structure components are wanted to align with in the explanatory space (e.g. variable bundles, principal components, other subspaces). Moreover, SCGLR searches for explanatory directions common to multiple responses with probability distributions in the exponential family, each response being entitled to their own distribution. The current SCGLR method is implemented in the R package SCGLR [START_REF] Cornu | SCGLR: Supervised Component Generalized Linear Regression[END_REF] available at https://CRAN.R-project.org/package=SCGLR and https://github.com/SCnext/SCGLR.

In the present work, we aim at modelling responses with a repeated or grouped design.

For this purpose, the use of mixed models with random effects is widespread. Research on variance-component estimation in Generalised Linear Mixed Models (GLMMs) has been very active since the 1980s. For the most general distribution assumptions in such models, parameter estimation faces the intractability of the likelihood expressed as an integral with respect to the random effects. Several numerical approximations of the integral have been proposed: Gaussian quadrature [START_REF] Anderson | Variance Component Models with Binary Response: Interviewer Variability[END_REF] or adaptive versions of it [START_REF] Pinheiro | Approximations to the Log-Likelihood Function in the Nonlinear Mixed-Effects Model[END_REF], Laplace approximation leading to the definition of the penalised quasi-likelihood [START_REF] Breslow | Approximate Inference in Generalized Linear Mixed Models[END_REF] or modified versions of it [START_REF] Shun | Laplace Approximation of High Dimensional Integrals[END_REF]. An alternative to this type of analytic approximation is a stochastic approximation of the integral calculation via MCMC techniques. In this approach, [START_REF] Zeger | Generalized Linear Models With Random Effects; A Gibbs Sampling Approach[END_REF] described an approximate Gibbs sampling for GLMMs, which was extended by [START_REF] Clayton | Generalized linear mixed models[END_REF] to more general Metropolis-Hastings algorithms. In parallel, [START_REF] Mcculloch | Maximum Likelihood Algorithms for Generalized Linear Mixed Models[END_REF] developed the Monte Carlo EM algorithm where the expectation is computed numerically through a Monte Carlo approximation, after generating random effects with a Metropolis-Hastings sampler. Mention can also be made of the recent work by [START_REF] Knudson | Monte Carlo Likelihood Approximation for Generalized Linear Mixed Models[END_REF]: her strategy is to approximate the entire likelihood function using random effects simulated from a parametrised importance sampling distribution depending on the data.

Unfortunately, these different approaches are not necessarily suitable for the same types of random effect designs (one-dimensional random effect, embedded random effects, etc).

In the wake of the first type of approximations, we here adopt the "Joint-Maximisation" strategy [START_REF] Mcculloch | Maximum Likelihood Algorithms for Generalized Linear Mixed Models[END_REF], as introduced for instance by [START_REF] Schall | Estimation in Generalized Linear Models with Random Effects[END_REF]. The model is iteratively linearised conditional on the random effects and variance components are then estimated using adapted linear mixed models methods. This strategy can be used for any random effect design and is less computationally intensive than Monte Carlo methods.

Moreover, it provides us with a linear setting more suitable for the computation of components. Once the components calculated, model parameters can be estimated using any of the aforementioned strategies (see Section 7).

Modelling grouped responses through a GLMM with a large number of explanatory variables is the focus of this paper. The need for dimension-reduction and regularisation has to accommodate the presence of random effects in the model, but our main purpose still remains to investigate the explanatory structure and link it to interpretable dimensions.

For Gaussian responses, [START_REF] Eliot | Ridge Regression for Longitudinal Biomarker Data[END_REF] proposed to extend the ridge regression to Linear Mixed Models (LMMs). Based on a penalised complete log-likelihood, the adaptation of the Expectation-Maximisation algorithm they suggest includes a new step to find the best shrinkage parameter using a generalised cross-validation scheme at each iteration. More recently, [START_REF] Schelldorfer | GLMMLasso: An Algorithm for High-Dimensional Generalized Linear Mixed Lodels Using 1 -Penalization[END_REF] -and also [START_REF] Groll | Variable selection for generalized linear mixed models by L 1 -penalized estimation[END_REF] -proposed an L 1penalised algorithm for fitting a high-dimensional GLMM, using Laplace approximation and an efficient coordinate gradient descent. In this work, we combine Schall's iterative model linearisation with regularisation at each step. However, we do not use a penalty on the coefficient vector's norm -as proposed by [START_REF] Zhang | Regression Models for Multivariate Count Data[END_REF] within the framework of multivariate count data. We rather propose to combine dimension-reduction and predictor-regularisation using supervised components aligning on the most predictive and interpretable directions in the explanatory space.

The paper is organised as follows. In Section 2, we formalise the model and set the main notations used throughout the paper. In Section 3, we present the key features of SCGLR.

Section 4 designs an extension of this methodology to mixed models, and particularly to grouped data. In Section 5, our extended method "mixed-SCGLR" is evaluated on simulations and compared to ridge-and LASSO-based regularisations. Finally, in order to highlight the power of mixed-SCGLR in terms of model interpretation, Section 6 presents an application to real data in the Poisson case.

Model definition and notations

In the framework of a multivariate GLMM, we consider q response-vectors y 1 , . . . , y q forming matrix Y n×q , to be explained by two categories of explanatory variables. The first category consists of few weakly correlated variables A n×r = a 1 | . . . | a r . These variables are assumed to be interesting per se and their marginal effects need to be precisely quantified. The second category consists of abundant and highly correlated variables

X n×p = x 1 | . . . | x p considered
as proxies to latent dimensions which must be found and interpreted. Since explanatory variables in A are few, non-redundant and of interest, they are kept as such in the model. By contrast, X may contain several unknown structurally relevant dimensions K < p important to model and predict Y , how many we do not know. X is thus to be searched for an appropriate number of orthogonal components that both capture relevant structural information in X and contribute to model Y . This work addresses grouped data: the n observations form N groups. Within each group, observations are not assumed independent. For each response y k , a N -level random effect ξ k is used to model the dependence of observations within each group. Hence, each y k is modelled with a GLMM assuming a conditional distribution from the exponential family.

Notations and conventions

All variables (namely the a i 's, x j 's and y k 's) will be identified with n-vectors.

We will use bold lowercase letters for vectors (e.g. u) and bold capital letters for matrices (e.g. M ).

M being any matrix, M T denotes the transpose of M .

I n denotes the identity matrix of size n.

M (u, v) = u | v M u M v M
, where

u M = u | u M .
The space spanned by vectors u 1 , . . . , u h is denoted by span {u 1 , . . . , u h }. U being any matrix, span {U } refers to the space spanned by the column-vectors of U .

Let R n be endowed with metric W and let Z be a matrix of size n × p. Then Π W span{Z} refers to the W -orthogonal projector onto span {Z}. Let b be a vector in R n . The cosine of the angle between b and span {Z} with respect to W is defined by cos W (b, span {Z}) = cos W b, Π W span{Z} b .

SCGLR with additional explanatory variables

In this section we consider the situation where each y k is modelled with a GLM (without random effect). For the sake of simplicity, we focus on the single-component SCGLR (K = 1). Section 3.1 briefly recalls some standards for univariate GLMs. Section 3.2 defines the linear predictors considered in the SCGLR methodology, in a multivariate GLM framework with additional explanatory variables. Finally, Section 3.3 introduces the criterion SCGLR maximises to compute the component.

Notations and main features of univariate GLMs

We refer the reader to [START_REF] Mccullagh | Generalized Linear Models, Second Edition[END_REF] for a thorough overview of GLMs.

This section is only intended to recall the classical iterative scheme performing maximum likelihood (ML) estimation. Let X denote the n × p matrix of explanatory variables and β the p-dimensional parameter vector. At iteration t + 1, the Fisher Scoring Algorithm (FSA) for ML estimation calculates

β [t+1] = X T W [t] X -1 X T W [t] z [t] , (1)
where z [t] and W [t] respectively denote the classical working variable and the associated weight matrix at iteration t. As pointed out by [START_REF] Nelder | Generalized Linear Models[END_REF], update (1) may be interpreted as a weighted least squares step in the linearised model M [t] defined by [t] :

M
z [t] = Xβ + ζ [t]
with:

E ζ [t] = 0 and V ζ [t] = W [t] -1 .
(2)

Linear predictors for SCGLR with multiple responses

We are now considering a multivariate GLM [START_REF] Fahrmeir | Multivariate Statistical Modelling Based on Generalized Linear Models[END_REF]. In this context, SCGLR searches for a component common to all the y k 's. This component will be denoted f and its p-dimensional loading-vector will be denoted u, so that f = Xu. The linear predictor associated with response-vector y k then writes

η k = (Xu) γ k + Aδ k , (3)
where γ k and δ k are the regression parameters associated respectively with component f and additional explanatory variables A. f being common to all the y k 's, predictors are collinear in their X-part. For identification purposes, we impose u T M -1 u = 1, where M may so far be any p × p symmetric positive definite matrix. Let us note y k,i the i-th observation of the k-th response-vector and H = {η k,i | 1 k q, 1 i n} the predictor set. We assume that the q responses are independent conditional on f , and that the n observations are independent. The log-density then writes

(Y |H) = n i=1 q k=1 k (y k,i |η k,i ) ,
where k is the log-density of the k-th response, conditional on its linear predictor. As a result, z k being the working variable associated with y k and W -1 k its variance matrix, the corresponding linearised model derived from the FSA at iteration t is

M [t] k : z [t] k = (Xu) γ k + Aδ k + ζ [t] k with: E ζ [t] k = 0 and V ζ [t] k = W [t] k -1 . (4)
Although linearised models (2) and (4) seem very similar, (4) is no longer linear, owing to the product uγ k . An alternate version of the FSA must therefore be used:

(i) Given current values of all the γ k 's and δ k 's, a new loading-vector u is obtained by solving an SCGLR-specific program (see Section 3.3 for details).

(ii) Given a current value of u, each z k is regressed independently on Xu | A with respect to weight matrix W k , yielding new regression parameters γ k and δ k .

Calculating the component maximising an SCGLR-specific criterion

For an easier reading of this part, we omit the [t] index. For each k ∈ {1, . . . , q}, consider model M k endowed with weight matrix W k . As suggested in [START_REF] Bry | Supervised component generalized linear regression using a PLS-extension of the Fisher scoring algorithm[END_REF], the best loading-vector in the weighted least-squares sense would be the solution of min

u: u T M -1 u=1 q k=1 z k -Π W k span{Xu,A} z k 2 W k ⇐⇒ max u: u T M -1 u=1 q k=1 Π W k span{Xu,A} z k 2 W k .
The maximisation program also writes max

u: u T M -1 u=1
ψ A (u), where

ψ A (u) = q k=1 z k 2 W k cos 2 W k z k , span {Xu, A} = q k=1 z k 2 W k cos 2 W k z k , Π W k span{Xu,A} z k . (5)
Now, ψ A is a mere goodness-of-fit (GoF) measure that does not take into account the closeness of component f = Xu to interpretable directions in X. The GoF measure, ψ A , must therefore be combined with a measure φ of structural relevance (SR).

Assume matrix X consists of p standardised numeric variables. Consider a weight system ω = {ω 1 , . . . , ω p } -e.g. ω j = 1 p ∀j ∈ {1, . . . , p} -reflecting the a priori relative importance of variables. Also consider a weight matrix P -e.g. P = 1 n I n -reflecting the a priori relative importance of observations. We define the most structurally relevant loading-vector as the solution of max

u: u T M -1 u=1 φ (u) , where φ (u) = p j=1 ω j Xu | x j 2 P l 1 l = p j=1 ω j u T X T P x j x T j P X u l 1 l , l 1, (6)
for the scalar product is commutative. Formula ( 6) is in fact a particular case of the SR criterion proposed by [START_REF] Bry | THEME: THEmatic Model Exploration through multiple co-structure maximization[END_REF]; [START_REF] Bry | The Multiple Facets of Partial Least Squares and Related Methods, chapter Supervised Component Generalized Linear Regression with Multiple Explanatory Blocks: THEME-SCGLR[END_REF]. It can be viewed as a generalised average version of the usual dual PCA criterion: 6) is called "Variable-Powered Inertia" (VPI). It should be stressed that for X T P X to be invertible, X must be a column full rank matrix. In case of strict collinearities within X, as it always happens in high-dimensional settings, we replace X with the matrix C of its principal components associated with nonzero eigenvalues. The component is then sought as f = Cu. We have C = XV , where V is the matrix of corresponding unit-eigenvectors. Then, f = Cu = X u with u = V u. [START_REF] Bry | Component-based regularisation of a multivariate GLM with a thematic partitioning of the explanatory variables[END_REF] show that among all loading-vectors t such that Xt = f , u is that which has the minimum L 2 -norm.

p j=1 cos 2 P (Xu, x j ) = p j=1 Xu | x j 2 P . For M = (X T P X) -1 , (
Tuning parameter l allows to draw component towards more (greater l) or less (smaller l) local bundles of correlated variables, as depicted on Figure 1 in the particular instance of four coplanar variables. Informally, a bundle is a set of variables correlated "enough" to be viewed as proxies to the same latent dimension. The notion of bundle is flexible, and parameter l tunes the level of within-bundle correlation to be considered: the higher the correlation, the more local the bundle. Overall, taking l = 1 draws the components towards global structural directions (namely the principal components) while taking l higher leads to more local ones (ultimately, the variables themselves). The goal is to focus on the most interpretable directions. Finally, let s ∈ [0, 1] be a parameter tuning the importance of the SR relative to the GoF. SCGLR attempts a trade-off between (5) and ( 6) by solving max

u: u T M -1 u=1 [φ (u)] s [ψ A (u)] 1-s or equivalently max u: u T M -1 u=1 s log [φ (u)] + (1 -s) log [ψ A (u)] . (7)
More detail can be found in [START_REF] Bry | Component-based regularisation of a multivariate GLM with a thematic partitioning of the explanatory variables[END_REF].

Extension to mixed models

We now propose to extend SCGLR to mixed models. This extension will be called "mixed-SCGLR". A particular focus is placed on grouped data, for which the independence as-Figure 1: Polar representation of the VPI according to the value of l in the elementary case of four coplanar variables, x 1 , x 2 , x 3 , x 4 , with ω j = 1 4 ∀j ∈ {1, 2, 3, 4}. Loading-vector u is identified with complex number e iθ , where θ ∈ [0, 2π). Curves z l (θ) := φ e iθ l e iθ are graphed for l ∈ {1, 2, 4, 10, 50}. The intersection of curve z l with f = Xu has a radius equal to φ e iθ l . The red line is the direction of maximum for l = 1, which is in fact the first principal component. These four variables are then regarded as a unique bundle. By contrast, the blue lines represent the two directions of maximum for l = 4. The variables are then seen as two bundles containing two variables each. Finally, when l = 50, each variable is considered a bundle in itself.

sumption of observations is no longer valid. The within-group dependence of each response is modelled with a random group-effect. Consequently, each y k is modelled with a GLMM.

As in SCGLR, the responses are assumed to be independent conditional on the components. Section 4.1 presents the single-component mixed-SCGLR method. The underlying algorithm is given in Section 4.2. Considering only one component is generally not enough to explain the responses making it necessary to search for K explanatory components, with 1 K rank (X). The way in which we extract higher rank components is explained in Section 4.3.

First component

The random group-effect is assumed different across responses. This leads to q randomeffect vectors ξ 1 , . . . , ξ q , which are assumed independent and normally distributed:

∀k ∈ {1, . . . , q} , ξ k ind. ∼ N N (0, D k ) ,
where N denotes the number of groups. In this paper, variance components models will be considered. We assume D k = σ 2 k I N , where σ 2 k is the group variance component associated with response y k . Linear predictors involved in mixed-SCGLR are expressed as

∀k ∈ {1, . . . , q} , η ξ k = (Xu)γ k + Aδ k + U ξ k , ( 8
)
where U is the known random effects' design matrix. Predictor η ξ k epitomises the way we capture the dependence between outcomes. Indeed, as component f = Xu does not depend on k, it captures a structural dependence between the various y k 's. By contrast, the random effect ξ k models the within-group stochastic dependence of outcomes forming response-vector y k .

Recall that the distribution of the data conditional on the random effects is supposed to belong to the exponential family. The FSA was adapted by [START_REF] Schall | Estimation in Generalized Linear Models with Random Effects[END_REF] to the GLMM dependence structure. The key idea is to extend Schall's algorithm to the component-based predictors in (8).

Linearisation step

Let g k denote the link function for response y k , g k its first derivative and µ ξ k the conditional expectation (i.e.

µ ξ k := E (y k | ξ k ))
. The working variable associated with y k,i is calculated through

z ξ k,i = g k µ ξ k,i + y k,i -µ ξ k,i g k µ ξ k,i = η ξ k,i + e k,i , where e k,i = y k,i -µ ξ k,i g k µ ξ k,i .
In view of the conditional independence assumption, the conditional variance matrix for

z ξ k is Var z ξ k | ξ k = W ξ k -1 = Diag g k µ ξ k,i 2 a k,i (φ k ) v k µ ξ k,i i=1,...,n
, where a k,i and v k are known functions, and φ k is the dispersion parameter related to y k .

At iteration t, the conditional linearised model for working vector z ξ k is then defined by

M ξ k [t] : z ξ k [t] = (Xu) γ k + Aδ k + U ξ k + e [t] k with: E e [t] k | ξ k = 0 and V e [t] k | ξ k = W ξ k -1 [t] . (9)
Besides the variance component estimation, an alternated estimation step has to be developed (as aforementioned in Section 3.2) to deal with the non-linearity of (9).

Estimation step

Calculating the component: Given current values of all the γ k 's, δ k 's, ξ k 's and σ 2 k 's, a new component f = Xu is calculated by solving a (7)-type program. However, (5) has to be adapted to conditional linearised models M ξ k 's, involving weight matrices W ξ k 's. The appropriate goodness-of-fit measure is

ψ A (u) = q k=1 z ξ k 2 W ξ k cos 2 W ξ k z ξ k , span {Xu, A} . (10)
Estimating the regression parameters and variance-components: Given a current value of component f , we apply Schall's method with the linear predictors given in ( 8). New values of parameters γ k and δ k as well as new prediction ξ k are obtained by solving the following Henderson system [START_REF] Henderson | Best Linear Unbiased Estimation and Prediction under a Selection Model[END_REF]:

     f T W ξ k f f T W ξ k A f T W ξ k U A T W ξ k f A T W ξ k A A T W ξ k U U T W ξ k f U T W ξ k A U T W ξ k U + D -1 k           γ k δ k ξ k      =      f T W ξ k z ξ k A T W ξ k z ξ k U T W ξ k z ξ k      .
Finally, as mentioned by [START_REF] Schall | Estimation in Generalized Linear Models with Random Effects[END_REF], given prediction ξ k for ξ k , the update of the ML

estimation of variance component σ 2 k is σ 2 k ←- ξ k T ξ k N -1 σ 2 k Trace U T W ξ k U + D -1 k -1 .

The algorithm

The conditional linearised models considered at iteration t are given by ( 9). Algorithm 1 describes the (t + 1)-th iteration of the single-component mixed-SCGLR. It is repeated until stability of parameters is reached.

Step 1: Computing the component. Set

u [t+1] = arg max u: u T M -1 u=1 [φ (u)] s ψ [t] A (u) 1-s
, where ψ A (u) is given by ( 10) and φ (u) by ( 6)

f [t+1] = Xu [t+1]
Step 2: Henderson systems. For each k ∈ {1, . . . , q}, solve the system

      f [t+1] T W ξ k [t] f [t+1] f [t+1] T W ξ k [t] A f [t+1] T W ξ k [t] U A T W ξ k [t] f [t] A T W ξ k [t] A A T W ξ k [t] U U T W ξ k [t] f [t] U T W ξ k [t] A U T W ξ k [t] U + D [t] k -1             γ k δ k ξ k       =       f [t+1] T W ξ k [t] z ξ k [t] A T W ξ k [t] z ξ k [t] U T W ξ k [t] z ξ k [t]       Call γ [t+1] k , δ [t+1] k
and ξ

[t+1] k the solutions.

Step 3: Updating variance-component estimates. For each k ∈ {1, . . . , q}, compute

σ 2 k [t+1] = ξ [t+1] k T ξ [t+1] k N -1 σ 2 k [t] Trace U T W ξ k [t] U + D [t] k -1 -1 and D [t+1] k = σ 2 k [t+1] I N
Step 4: Updating working variables and weighting matrices.

For each k ∈ {1, . . . , q}, compute

η ξ k [t+1] = f [t+1] γ [t+1] k + Aδ [t+1] k + U ξ [t+1] k µ ξ k,i [t+1] = g -1 k η ξ k,i [t+1] , i = 1, . . . , n z ξ k,i [t+1] = η ξ k,i [t+1] + y k i -µ ξ k,i [t+1] g k µ ξ k,i [t+1] , i = 1, . . . , n W ξ k [t+1] = Diag   g k µ ξ k,i [t+1] 2 a k,i (φ k ) v k µ ξ k,i [t+1] -1   i=1,...,n Incrementing: t ←-t + 1
Algorithm 1: Iteration of the single-component mixed-SCGLR

Extracting higher rank components

Let F h = f 1 | . . . | f h be the matrix of the first h components, where h < K. An extra component f h+1 must best complement the existing ones plus A, i.e. A h = F h | A . So f h+1 must be calculated using A h as additional explanatory variables. Moreover, we must impose that f h+1 be orthogonal to F h , i.e.

F T h P f h+1 = 0. Component f h+1 = Xu h+1 is thus obtained by solving    max s log [φ (u)] + (1 -s) log [ψ A h (u)] subject to: u T M -1 u = 1 and D T h u = 0, (11)
where

ψ A h (u) = q k=1 z ξ k 2 W ξ k cos 2 W ξ k z ξ k , span {Xu, A h } and D h = X T P F h .
In the online Supplementary Material, we give a simple tool to maximise, at least locally, any criterion on the unit sphere: the Projected Iterated Normed Gradient (PING) algorithm. In particular, PING solves (11)-type programs, which give all components of rank h > 1. The rank-one component is computed using the same program with A 0 = A and D 0 = 0.

Comparative results on simulated data

Five simulation studies have been implemented to assess our method. The first one (discussed in Sections 5.1 -5.4) focuses on LMMs. It compares the performances of mixed-SCGLR, LMM-ridge [START_REF] Eliot | Ridge Regression for Longitudinal Biomarker Data[END_REF] and GLMM-LASSO [START_REF] Groll | Variable selection for generalized linear mixed models by L 1 -penalized estimation[END_REF][START_REF] Schelldorfer | GLMMLasso: An Algorithm for High-Dimensional Generalized Linear Mixed Lodels Using 1 -Penalization[END_REF]. The second simulation (Section 5.5) extends the first one to binary and Poisson outcomes. All simulation studies have been performed using R (R Core Team, 2017). To compute LASSO regressions, we have used the R package glmmLasso [START_REF] Groll | glmmLasso: Variable Selection for Generalized Linear Mixed Models by L1-Penalized Estimation[END_REF]. The extension of SCGLR to mixed models is available at https://github.com/SCnext/mixedSCGLR. Three additional simulations are presented in the online Supplementary Material. The first one reproduces the simulation scheme of Section 5.5 with binomial and Poisson outcomes. The second one assesses the performance of mixed-SCGLR on a different bundle structure and presents results concerning variance component estimates. The third one deals with high dimensional data.

Data generation

To generate grouped data, we consider N = 10 groups and R = 10 observations per group (i.e. a total of n = 100 observations). The random effects' design matrix is then

U = I N ⊗ 1 R .
Explanatory variables X consist of three independent bundles: X 0 (15 variables), X 1 (10 variables) and X 2 (5 variables). Each explanatory variable is normally simulated with mean 0 and variance 1. Parameter τ ∈ {0.1, 0.3, 0.5, 0.7, 0.9} tunes the level of redundancy within each bundle: the correlation matrix of bundle X j is cor (

X j ) = τ 1 p j 1 T p j + (1 -τ )I p j ,
where p j is the number of variables in X j . In order to enable comparison with LASSO and ridge and to focus on regularisation, our simulations do not involve additional explanatory variables (A = 0). Two random responses Y = y 1 | y 2 are generated as

   y 1 = Xβ 1 + U ξ 1 + ε 1 y 2 = Xβ 2 + U ξ 2 + ε 2 , (12)
such that y 1 is predicted only by bundle X 1 , y 2 only by bundle X 2 , and bundle X 0 plays no explanatory role. Our choice for the fixed-effect parameters is Finally, for each k ∈ {1, 2}, random effect and noise vectors are simulated respectively from

β 1 = (
ξ k ∼ N N 0, σ 2 k I N and ε k ∼ N n 0, ω 2 k I n , where σ 2 k = ω 2 k = 1.
For each value of τ , B = 100 samples are generated according to Model (12).

Parameter calibration

In order to compare mixed-SCGLR with the ridge and LASSO regressions, we recall the regularisation parameters required by each method. For both LMM-ridge and GLMM-LASSO methods, a unique shrinkage parameter has to be calibrated: λ ridge and λ LASSO respectively. For mixed-SCGLR, three tuning parameters need to be calibrated: the number of components K and the trade-off parameter s, which are both regularisation parameters, and the bundle-locality parameter l. For greater clarity, the simulation focuses on the behaviour of K and s. As recommended by [START_REF] Bry | Supervised component generalized linear regression using a PLS-extension of the Fisher scoring algorithm[END_REF], we set l = 4. In case-studies, l has to be tuned to maximise the interpretability of components.

For both mixed-SCGLR and GLMM-LASSO, optimal regularisation parameters are obtained through a 5-fold cross-validation, withdrawing 2 observations from each group every time. This could be termed "leave-two-observations-out per group." The data are thus divided into five parts P 1 , . . . , P 5 , each P j containing 20 observations, 2 for each of the 10 groups. Let y b) , is defined as

E (b) = 1 2 2 k=1 E (b) k , (13)
where

E (b) k = 1 5 5 j=1 1 20 i∈P j y (b) k,i -y (b) k,i(-j) 2 .
In the b-th sample, the optimal number of components K (b) , the trade-off parameter s (b) , and the shrinkage parameter λ (b)

LASSO are selected to minimise the cross-validation error (13). We then define b) , K = mode K (1) , . . . , K (B) and

s = 1 B B b=1 s (
λ LASSO = 1 B B b=1 λ (b) LASSO .
By contrast, [START_REF] Eliot | Ridge Regression for Longitudinal Biomarker Data[END_REF] suggest to calibrate the ridge parameter at each step of their EM implementation, using the generalised cross-validation. We thus define

λ ridge = 1 B B b=1 λ (b) ridge ,
where λ

(b)
ridge denotes the average of the ridge parameter values obtained over all the iterations of the EM algorithm in the b-th sample.

Table 1 summarises the optimal regularisation parameters selected through crossvalidation. In both ridge and LASSO, the shrinkage parameter value increases with the level of redundancy τ . Whereas for mixed-SCGLR, when τ increases, K decreases towards the true number of predictive variable-bundles: the greater the value of τ , the better mixed-SCGLR focuses on the structures in X that contribute to model Y . Moreover, when τ increases, the trade-off parameter s increases, meaning that regularisation requires a greater importance of the structural relevance relative to the goodness-of-fit.

Comparison of the estimate accuracies

Once tuning parameters are obtained, we focus on the fixed-effect estimates' accuracy.

Since the response-vectors y 1 and y 2 are normally distributed and have comparable orders of magnitude, the fixed-effect relative errors are on the same scale. Then we consider a risk-averse comparison criterion called "Mean Upper Relative Squared Error" (MURSE) defined as

MURSE (β 1 , β 2 ) = 1 B B b=1 max      β (b) 1 -β 1 2 β 1 2 , β (b) 2 -β 2 2 β 2 2      , where β (b)
k is the estimate of β k associated with sample b. The MURSE values for mixed-SCGLR, LMM-ridge and GLMM-LASSO are presented in Table 2. The LMM results obtained without regularisation are also presented. They were computed using the R package lme4 [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF]. In the latter case, relative errors increase dramatically with τ . Those of ridge and LASSO increase less drastically (but increase anyway) because these methods suffer from the high correlations among the explanatory variables. Except for τ = 0.1, mixed-SCGLR provides the most accurate fixed effect estimates. Indeed, if there are no real bundles in X (τ 0), searching for structures in X may lead mixed-SCGLR to be slightly less accurate. Conversely, mixed-SCGLR takes advantage of the high correlations among the explanatory variables: the stronger the structures (high τ ), the more efficient the method.

Model interpretation

This section aims at highlighting the power of mixed-SCGLR for model interpretation.

Additional simulations involving non-Gaussian outcomes

This section aims at assessing our method in the case of Bernoulli (B) and Poisson (P) distributions of responses. We still consider N = 10 groups and R = 10 observations per group. We keep design matrices X and U defined in Section 5.1, as well as the values of

β 1 , β 2 , σ 2 1 and σ 2 2 .
The group variance components are given by ς 2 1 = 0.1σ 2 1 and ς 2 2 = σ 2 2 so that for each k ∈ {1, 2}, ξ k ∼ N N (0, ς 2 k I N ). Then given ξ 1 and ξ 2 , we simulate

Y = y 1 | y 2 as      y 1 ∼ B p = logit -1 Xθ 1 + U ξ 1 y 2 ∼ P λ = exp Xθ 2 + U ξ 2 , (14)
where θ 1 = 0.1β 1 and θ 2 = β 2 . Again, for each value of τ , B = 100 samples are generated according to Model (14). As in Section 5.2, tuning parameters are calibrated so as to minimise the cross-validation error (13). However, since y 1 and y 2 do not have the same range of values, the prediction errors have to be standardised. The cross-validation error for response y k in the b-th sample is now given by

E (b) k = 1 5 5 j=1 1 20 i∈P j y (b) k,i -y (b) k,i(-j) 2 var y (b) k,i(-j)
.

Unlike in Section 5.3, the response-vectors do not come from the same distribution and have different orders of magnitude. The fixed-effect relative errors are thus not comparable. To compare mixed-SCGLR with GLMM-LASSO and classical GLMM (without regularisation), we thus use the Mean Relative Squared Error (MRSE) defined as

MRSE (θ k ) = 1 B B b=1 θ (b) k -θ k 2 θ k 2 , k ∈ {1, 2} , where θ (b)
k is the estimate of θ k from the b-th sample. MRSE values for the GLMM, mixed-SCGLR and GLMM-LASSO are presented in Table 3. For all methods, estimating a Bernoulli model is obviously a more challenging task than estimating a Poisson model.

Regardless of the level of redundancy τ , both mixed-SCGLR and GLMM-LASSO outperform classical GLMM estimation. Compared with the Gaussian case (Section 5.3), the results deteriorate but (overall) the same behaviours are observed.

For τ = 0.1, fixed-effect estimates provided by mixed-SCGLR are less accurate than those provided by GLMM-LASSO. In this case, GLMM-LASSO has indeed a double advantage. First, many θ k,j 's are true zeros. Unlike mixed-SCGLR, GLMM-LASSO often shrinks their estimates to exactly zero. Second, since the level of redundancy is low, GLMM-LASSO also provides accurate coefficient estimates of active variables.

By contrast, for τ 0.3, mixed-SCGLR takes advantage of redundancies within the explanatory variables. Thus, mixed-SCGLR outperforms GLMM-LASSO in this case, despite the sparse structure of the θ k 's.

Even if the response variables are not Gaussian, the power of mixed-SCGLR for model interpretation is preserved. Graphical diagnoses similar to those provided in Section 5.4

are available in the Supplementary Material. Physical factors are many and redundant: monthly rainfalls are highly correlated, and so are photosynthesis activity indicators. By contrast, geology and anthropogenic interference are weakly correlated and interesting per se. These variables are then considered as additional explanatory variables and included in matrix A.

Model and parameter calibration

Abundances of species given in Genus are count data. For each k ∈ {1, . . . , 8}, we consider a Poisson regression with log link

y k ∼ P λ = exp K j=1 (Xu j ) γ k,j + Aδ k + U ξ k ,
where ξ k is the 22-level random-effect vector used to model the dependence between the observations of y k within concessions. The first cross-validations we performed -with different fixed values of parameters s and l -indicated that four components were sufficient to capture most of the information in X needed to model and predict responses. We therefore keep K = 4. The optimal values of trade-off and locality parameter s and l are then determined through another cross-validation. Using the same procedure and notations as in Section 5.2, the data are divided into five parts P 1 , . . . , P 5 . Let n j be the size of P j .

E = 1 8 8 k=1 E k ,
where

E k = 1 5 5 j=1 1 n j i∈P j y k,i -y k,i(-j) 2 var y k,i(-j) . ( 15
)
On Figure 3, we plot the errors E for parameter pairs (s, l) ∈ E s × E l , where 10, 12, 14, . . . , 30, 35, 40, 45, 50} .

E s = {0.025, 0.1, 0.2, . . . , 1} E l = {1, 2, . . . ,
Parameter grid E s × E l therefore contains 264 pair values. Selecting the best parameter pair from E s × E l through a 5-fold cross-validation requires a computation time of about 65 minutes (parallel computing on 6 CPU cores, Intel Core i7-6700HQ, 2.6GHz). It should be noted that there is a risk of non-convergence when the trade-off parameter s is too close to 0. Indeed, if we consider no structural information (s exactly equal to 0) in X, mixed-SCGLR merely performs classical GLMM estimation and does not converge with this data. When s = 0.025, our algorithm converges but leads to fairly unstable estimates and high cross-validation errors because regularisation is then very weak. By contrast, the components calculated with s ∈ {0.5, 0.6, . . . , 1} are close to principal components.

The associated errors are therefore stable in most cases, but rather high. Finally, s ∈ {0.1, . . . , 0.4} leads to the lowest cross-validation errors, but only for l 10. Indeed, when s is not too high, mixed-SCGLR may focus on the most predictive structures of X. However, parameter l must not exceed a certain value, in order to avoid being drawn towards too local variable-bundles. As can be seen, choosing (s , l ) = (0.1, 10) minimises the cross-validation error.

Prediction quality and interpretation results

This part evaluates the benefits obtained by taking within-group dependence into account.

The predictions we get with mixed-SCGLR and with initial version of SCGLR are compared with respect to the cross-validation criterion given by ( 15). The decomposition of linear predictors on interpretable components allows to detect the species that tend to share common explanatory dimensions and those which are more idiosyncratic. We can then identify the variable-bundles these dimensions are related to. The underlying goal is a better understanding of the bio-and ecosystem diversity with a view to preserve them. Species 1, 2, 5 and 6 are sensitive to the same rain-wind regime, and Species 4 and 8 are explained by the same photosynthetic pattern. On the contrary, Species 3 and 7 are clearly separated. Species 7 grows at high altitudes where the atmosphere is rather dry while the abundance of Species 3 is favoured by regular rainfall and high humidity.

Discussion and Conclusions

Like Sufficient Dimension Reduction (SDR) methods, mixed-SCGLR is based on the construction of a reduction function of dimension less than p which tries to capture all the relevant information that X contains about Y . However, the two approaches do not exactly pursue the same objectives. Indeed, SDR methods look for the "central subspace" containing the predictive information irrespective of the structures within X (e.g. dimensions capturing a large part of X's variance, or bundles of correlated variables). Mixed-SCGLR rather aims at basing the explanatory subspace on such structural dimensions so as to both gain interpretability and stabilise prediction. We think that extracting a hierarchy of strong and interpretable dimensions, and decomposing the linear predictor on them, is an essential asset in model-building. The difference in goals entails a difference in means: SDR is based on the sufficiency principle, which is enough to identify a subspace but not to track strong predictive dimensions in it. By contrast, in the wake of PLS regression, mixed-SCGLR uses a criterion combining goodness-of-fit and structural relevance of components.

The supervised-component paradigm has proved effective in situations where regularisation is necessary but where variable selection is inappropriate -for instance when the true explanatory dimensions are latent and indirectly measured through highly correlated proxies.

When l = 1, trade-off parameter s allows to continuously tune the attraction of components towards the principal components of explanatory variables. This results in a continuum between classical GLMM estimation (s = 0 is associated with no regularisation) and principal component generalised linear mixed regression (with s = 1).

When l > 1, we take better advantage of local predictive structures in X. The components we build are then usually closer to local gatherings of variables, thus easier to interpret.

Mixed-SCGLR is able to identify more or less local predictive structures common to all the y k 's and performs well on grouped data with Gaussian, Bernoulli, binomial and Poisson outcomes. Compared to penalty-based approaches as ridge or LASSO, the orthogonal components built by mixed-SCGLR reveal the multidimensional explanatory and predictive dimensions, and greatly facilitate the interpretation of the model.

However, a natural question arises as to the accuracy of our methodology under significant deviations from normality. With binary data for instance, variance component estimates are prone to some bias towards zero [START_REF] Mcculloch | Maximum Likelihood Algorithms for Generalized Linear Mixed Models[END_REF]. That is why other estimation strategies might be considered, especially Monte Carlo integration methods which have the advantage of being based on direct approximations of the likelihood. Some examples are the MCMC methods developed by [START_REF] Hadfield | MCMC Methods for Multi-Response Generalized Linear Mixed Models: the MCMCglmm R Package[END_REF] in the GLMM framework, and the Monte Carlo Likelihood Approximation (MCLA) proposed by [START_REF] Knudson | Monte Carlo Likelihood Approximation for Generalized Linear Mixed Models[END_REF]. Indirect maximisations of the likelihood are also available such as Monte Carlo Expectation-Maximisation (MCEM) and Monte Carlo Newton-Raphson [START_REF] Mcculloch | Maximum Likelihood Algorithms for Generalized Linear Mixed Models[END_REF]. We think that these methodologies and Schall's could be combined sequentially. Indeed we could first take advantage of the linear approximation of the model in order to build the components, and then use MC-based methods to estimate both fixed-effect parameters and variance components. This would lead to replacing the current iteration of mixed-SCGLR given by Algorithm 1 with the following steps (to keep things simple, we take the canonical link):

1. Compute components F = f 1 | . . . | f K via the PING algorithm on Schall's linearised models.

2. For each k ∈ {1, . . . , q}, consider the hierarchy

f y k |ξ k ,γ k ,δ k (y k |ξ k , γ k , δ k ) = exp y T k η ξ k -1 T c η ξ k + 1 T d (y k ) ξ k |D k ∼ N (0, D k ) ,
where Even though such MC-based methods are computationally much more intensive than the "Joint-Maximisation" and have intrinsic disadvantages (particularly in the assessment of convergence and in the choice of prior distributions), they could give better results in case of binary data.

η ξ k = F γ k + Aδ k + U ξ k ,
1 IMAG, Univ Montpellier, CNRS, Montpellier, France. jocelyn.chauvet@umontpellier.fr ; xavier.bry@umontpellier.fr 2 Univ Paul-Valéry Montpellier 3, Montpellier, France.

catherine.trottier@univ-montp3.fr

The Projected Iterated Normed Gradient (PING) is a basic extension of the iterated power algorithm, for solving any program of the form

   max J h (u) , subject to: u T M -1 u = 1 and ∆ T h u = 0. (16)
Note that putting

v = M -1/2 u, G h (v) = J h M 1/2 v and B h = M 1/2 ∆ h , Program (16)
is strictly equivalent to Program (17):

   max G h (v) , subject to: v T v = 1 and B T h v = 0. (17)
Denoting

Π span{B h } ⊥ = I -B h (B T h B h ) -1 B T h and Γ h (v) = ∇ v G h (v) ,
a Lagrange multiplier-based reasoning gives the basic iteration of the PING algorithm:

v [t+1] = Π span{B h } ⊥ Γ h v [t] Π span{B h } ⊥ Γ h (v [t] ) . (18)
Although Iteration (18) follows a direction of ascent, it does not guarantee that G h actually increases on every step. We therefore propose a generic iteration of PING (Algorithm 2) and an alternative one (Algorithm 3), which both ensure that the criterion increases.

while convergence of v non reached do

κ [t] = Π span{B h } ⊥ Γ h v [t] Π span{B h } ⊥ Γ h (v [t] )
A unidimensional Newton-Raphson maximisation procedure is used to find the maximum of G h (v) on the arc v [t] , κ [t] and take it as v [t+1] .

t ←-t + 1 end Algorithm 2: Generic iteration of the PING algorithm while convergence of v non reached do m ←- Π span{B h } ⊥ Γ h v [t] Π span{B h } ⊥ Γ h (v [t] ) while G h (m) < G h v [t] do m ←- v [t] + m v [t] + m end v [t+1] = m t ←-t + 1 end Algorithm 3: Alternative generic iteration of the PING algorithm First rank component. Component f 1 = Xu 1 is obtained by solving    max s log [φ (u)] + (1 -s) log [ψ A (u)] subject to: u T M -1 u = 1.
This corresponds to Program ( 16) with h = 0, where

J 0 (u) = s log [φ (u)] + (1 -s) log [ψ A 0 (u)],
A 0 = A (the matrix of additional explanatory variables), and

∆ 0 = 0.
In this particular case, we have B 0 = M 1/2 ∆ 0 = 0, and so: This corresponds to Program ( 16), where

Π span{B 0 } ⊥ = I.

Higher rank components. Let

J h (u) = s log [φ (u)] + (1 -s) log [ψ A h (u)],
A h = F h | A , and

∆ h = X T P F h .
Initialisation. To quickly find f 1 , algorithm PING is initialised with the first PLS component of the responses on X. In like manner, for h 2, PING is initialised with the first PLS component of the responses on X deflated on components The power of mixed-SCGLR in terms of model interpretation remains the same for non-Gaussian outcomes. Figure 6 (respectively Figure 7) presents an example of the first component planes output by mixed-SCGLR in the binomial/Poisson (respectively Bernoulli/Poisson) case. As for Gaussian outcomes, the component planes reveal that y 1 is explained by bundle X 1 and y 2 by X 2 . In the binomial/Poisson case with τ = 0.3 (Figure 6), predictive bundles X 1 and X 2 are captured respectively by the first and the second components. The third component aligns on nuisance bundle X 0 , despite its high inertia. Figure 7 illustrates what may happen when the level of redundancy is very high (τ = 0.7 here). Since the explanatory variables are highly correlated, mixed-SCGLR regularisation requires that the structural relevance be given a heavy weight with respect to the goodness-of-fit, which leads to a trade-off parameter s close to 1 (s = 0.9 here). Having the greatest structural strength, the nuisance bundle is captured by the second component despite its lack of explanatory power. This is sometimes the price to be paid for the tradeoff. In our example, the second explanatory bundle is captured by the third component, so that the predictive dimensions are accurately represented in component plane (1, 3). 9 A new structure for the bundles -Gaussian, Pois-

F h-1 = f 1 | . . . | f h-1 .

son and binomial distributions

This simulation study tests mixed-SCGLR on a slightly more complex bundle structure.

Results concerning variance component estimates are also presented.

We consider a fixed-effect design matrix X n×p partitioned into 3 blocks G 1 , G 2 and G 3 . Block G 1 contains 10 predictive explanatory variables structured about a latent variable ϕ 1 ∼ N n (0, σ 2 LV I n ). Thus for each j ∈ {1, . . . , 10}, x j = ϕ 1 + ε j , where ε j ∼ N n (0, σ 2 noise I n ) such that σ 2 LV + σ 2 noise = 1. The correlation within G 1 is tuned by signal to noise (StN) ratio σ 2 LV /σ 2 noise (chosen in 1 3 , 1, 3 in practice). G 2 contains a predictive variable ϕ 2 = x 11 ∼ N n (0, I n ). G 3 contains 20 unstructured noise variables: for each j ∈ {12, . . . , 31}, x j ∼ N n (0, I n ). For each k ∈ {1, 2, 3}, random-effect vectors are simulated as ξ k ind.

∼ N

N (0, σ 2 k I N ). Given ξ 1 , ξ 2 , ξ 3 , we simulate 3 responses having different distributions, Y = y 1 | y 2 | y 3 , as            y 1 ∼ N n µ = α 1 ϕ 1 + U ξ 1 , Σ = I n y 2 ∼ P λ = exp α 2 ϕ 2 + U ξ 2 y 3 ∼ Bin trials = 25 1 n , p = logit -1 α 3 (ϕ 1 + ϕ 2 ) + U ξ 3 .
In our simulations, we set α 1 = σ 2 1 = 2, α 2 = σ 2 2 = 1 and α 3 = σ 2 3 = 0.5.

We consider in turn N = 10 and N = 50 groups, and R = 10 observations per group (n = 100 and n = 500 observations in total). B = 100 samples are generated for each pair of values (N, StN). The main goal of the study is to assess the ability of mixed-SCGLR to track down both latent variable ϕ 1 and predictive variable ϕ 2 . For j = 1 and 2, we then define

cor j = 1 B B b=1 cor ϕ j , f (b) j
, where f (b) j is the component most correlated with ϕ j issued from mixed-SCGLR in the b-th sample. Consistency of fixed-effect estimates is assessed through criteria err 1 , err 2 and err 3 defined by

err j = 1 B B b=1 α j ϕ j -X β (b) j 2 α j ϕ j 2 , j ∈ {1, 2} err 3 = 1 B B b=1 α 3 (ϕ 1 + ϕ 2 ) -X β (b) 3 2 α 3 (ϕ 1 + ϕ 2 ) 2 ,
where β (b) j is the fixed-effect estimate related to response y j associated with sample b.

Table 6 summarises the values of the afore-defined criteria and presents biases and standard errors of variance components estimates. For a given value of N , cor 1 increases towards 1 with ratio σ 2 LV /σ 2 noise : the tighter the block G 1 is structured about its latent variable, the better mixed-SCGLR can reconstruct it. The associated criterion err 1 then naturally decreases towards 0. On the other hand, cor 2 and err 2 are very stable, which proves that mixed-SCGLR is able to detect an isolated predictive variable among a large number of irrelevant others. As err 3 depends on how accurately mixed-SCGLR recovers To cope with high dimensional data, the key idea is to replace the fixed-effect design matrix, X, with the matrix C of its principal components associated with non-zero eigenvalues.

More precisely, λ j being the eigenvalue associated with the j-th eigenvector v j , the last eigenvector we consider, v r , is such that

λ r r j=1 λ j > 1 p ,
where p is the number of columns of matrix X. The matrix of the corresponding unit- Mixed-SCGLR then solves

eigenvectors is denoted V = v 1 | . . . | v
   max s log [φ (u)] + (1 -s) log [ψ A (u)] subject to u T C T P Cu = 1,
where the goodness-of-fit measure, ψ A , is given by

ψ A (u) = q k=1 z k 2 W ξ k cos 2 W ξ k z k , span {Cu, A} = q k=1 z k 2 W ξ k cos 2 W ξ k z k , Π W ξ k span{Cu,A} z k ,
and the structural relevance by

φ (u) = p j=1 ω j Cu x j 2 P l 1 l = p j=1 ω j u T C T P x j x T j P C u l 1 l .
This idea is tested on simulated data where the number of explanatory variables p exceeds the number of observations n.

Data generation

To generate grouped data, we consider N = 10 groups and R = 10 observations per group (i.e. a total of n = 100 observations). The random effects' design matrix is then U = 10.3 Results is set to (4, 0.5, 4). is set to (4, 0.9, 4).

1

  m denotes the all-ones vector of size m. Let u and v be non-zero vectors in R d and let M be a symmetric positive definite matrix of size d × d. Then u | v M = u T M v refers to the Euclidean scalar product of u and v with respect to metric M . The cosine of the angle between u and v with respect to M is given by cos

  0, . . . . . . . . . . . . . . . . . . , β 2 = ( 0, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . , 0 25 times , 0.3, 0.3, 0.4, 0.5, 0.5 ) T .

  be the i-th observation of the k-th response vector in the b-th sample. Let also y (b) k,i(-j) be the fit for y (b) k,i with part P j removed. The cross-validation error in the b-th sample, E (

Figure 2 Figure 2 :

 22 Figure2presents an example of the first component planes obtained for τ = 0.5, with associated optimal parameter values s = 0.7 and K = 3. We still impose l = 4. The first two components obtained are the ones which explain the responses. It clearly appears that y 1 is explained by bundle X 1 and y 2 by X 2 . Interestingly, although bundle X 0 is the one with maximum inertia (26.83%), it appears only along the third component, for having no explanatory part.

  29 physical factors linked to topography, rainfall or soil moisture, 25 photosynthesis activity indicators (the Enhanced Vegetation Indices, EVI, the Near-InfraRed indices, NIR, and the Mid-InfraRed indices, MIR), 2 indicators which describe the tree height.

Figure 3 :

 3 Figure 3: Behaviour of the cross-validation error E for trade-off parameter s ∈ {0.025, 0.1, 0.2, . . . , 1}, as a function of locality parameter l ∈ [1, 50].

Figure 4 :Figure 5 :

 45 Figure 4: Abundance maps issued from mixed-SCGLR. The plots respectively show real abundance (left) and associated conditional predictions (right) of the tree species number 8. Each point represents a land plot (2615 in total).

  and c, d are the functions associated with the natural parametrisation of the GLM. For example, for the Bernoulli-logistic regression, we have: c(x) = log(1 + e x ) and d(x) = 0. 3. Apply MC-based methods such as MCMC, MCLA, MCEM or MCNR to update γ k , δ k , ξ k and D k , k ∈ {1, . . . , q}. 4. Update working variables and weight matrices to define the new Schall's linearised models.

  F h = f 1 | . . . | f h be the matrix of the first h components and A h = F h | A . Let P denote the weight matrix reflecting the a priori relative importance of observations (P = 1 n I n if all observations are of equal importance). Component f h+1 = Xu h+1 is obtained by solving    max s log [φ (u)] + (1 -s) log [ψ A h (u)] subject to: u T M -1 u = 1 and F T h P Xu = 0.

Figure 6 :

 6 Figure 6: Example of component planes given by mixed-SCGLR in the binomial/Poisson case for τ = 0.3, with parameter triplet (K, s, l) = (3, 0.5, 2).

Figure 7 :

 7 Figure 7: Example of component planes given by mixed-SCGLR in the Bernoulli/Poisson case for τ = 0.7, with parameter triplet (K, s, l) = (3, 0.9, 4).

ϕ 1

 1 and ϕ 2 , it slightly decreases when the StN ratio increases. Both variance component biases and standard errors seem rather stable regardless of the value of StN. Finally, when N increases, all the cor j 's increase towards 1 and all the err j 's decrease towards 0. As far as variance component estimates are concerned, the biases are getting slightly closer to 0 and the standard errors decrease significantly. Model interpretation is revealed by Figure 8 in the case of N = 10 groups and R = 10 observations per group. The first component aligns with block G 1 which alone explains response y 1 . The second aligns with G 2 (containing single explanatory variable x 11 ) which alone explains y 2 . Finally, note that the projection of the X-part of the linear predictor related to y 3 is well represented on component plane (1, 2). This indicates that y 3 is explained jointly by G 1 and G 2 . Table 6: Summary of cor j and err j values, and presentation of biases and standard errors of estimated variance components. N = 10, R = 10 (n = 100) N = 50, R = 10 (n = 500)

Figure 8 :

 8 Figure 8: Examples of the first two-component planes given by mixed-SCGLR when σ 2 LV /σ 2 noise = 1/3 (top left), σ 2 LV /σ 2 noise = 1 (top right), and σ 2 LV /σ 2 noise = 3 (bottom). When StN ratio = 1/3 (resp. StN ratio ∈ {1, 3}), only the variables having cosine greater than 0.4 (resp. 0.5) with component plane (1, 2) are represented.

  r , and C = XV . The component f is then sought as a combination of the principal components: f = Cu = X u, where u = V u.

Figure 9 :

 9 Figure 9: Component planes (1, 2), (1, 3) and (1, 4) given by mixed-SCGLR for 100 observations, 150 explanatory variables and τ = 0.3. The tuning parameter triplet (K, s, l)

Figure 10 :

 10 Figure 10: Component planes (1, 2), (1, 3) and (1, 4) given by mixed-SCGLR for 100 observations, 200 explanatory variables, and τ = 0.9. The tuning parameter triplet (K, s, l)

Table 1 :

 1 Optimal regularisation parameter values obtained through cross-validation over 100 simulations.

		GLMM-LASSO	LMM-ridge	mixed-SCGLR
		shrinkage	shrinkage	number of	trade-off
		parameter λ LASSO	parameter λ ridge	components K	parameter s
	τ = 0.1	65	24	15	0.50
	τ = 0.3	92	54	5	0.58
	τ = 0.5	124	73	3	0.70
	τ = 0.7	163	78	3	0.73
	τ = 0.9	175	85	2	0.80

Table 2 :

 2 Mean Upper Relative Squared Error (MURSE) values associated with the optimal

	parameter values.
	LMM (no regularisation)	GLMM-LASSO LMM-ridge mixed-SCGLR
	τ = 0.1	0.12	0.05	0.08	0.12
	τ = 0.3	0.33	0.12	0.13	0.10
	τ = 0.5	0.61	0.20	0.16	0.07
	τ = 0.7	1.32	0.25	0.20	0.06
	τ = 0.9	4.62	0.26	0.31	0.05

Table 3 :

 3 Mean Relative Squared Error (MRSE) values obtained with Bernoulli and Poisson responses.

		GLMM (no regularisation)	GLMM-LASSO	mixed-SCGLR
		Bernoulli	Poisson	Bernoulli	Poisson	Bernoulli	Poisson
	τ = 0.1	316.48	0.54	8.61	0.30	14.71	0.46
	τ = 0.3	398.78	0.64	9.23	0.36	7.21	0.21
	τ = 0.5	576.68	0.87	14.48	0.44	2.01	0.09
	τ = 0.7	886.04	1.28	17.37	0.47	1.50	0.07
	τ = 0.9	2840.10	3.72	17.24	0.59	1.31	0.05
	6 An application to forest ecology data
	6.1 Data description

The present study is based on the Genus dataset of the CoForChange project (see http: //www.coforchange.eu). The subsample we consider gives the abundance of 8 common tree genera on 2615 Congo Basin land plots. These plots are grouped into 22 forest concessions. To predict abundances, we have 56 environmental variables, plus 2 explanatory variables which code geology and anthropogenic interference. X consists of all environmental variables which are:

Table 4 :

 4 (4, 0.1, 10) is selected for both methods. For each k ∈ {1, . . . , 8}, mixed-SCGLR gives a lower cross-validation error than SCGLR: taking into account the within-group dependence has clearly improved prediction performances. Cross-validation errors for each response variable.

			E 1 cv	E 2 cv	E 3 cv		E 4 cv	E 5 cv	E 6 cv	E 7 cv	E 8 cv
		SCGLR	1.32 2.46 3.27 1.43 2.56 1.28 1.54 3.44
		mixed-SCGLR 1.24 1.95 2.92 1.32 2.27 1.15 1.31 3.01
	Moreover, mixed-SCGLR enables to correctly reconstitute observed abundance maps, as
	illustrated on Figure 4.
		Real abundance						Prediction
		Central African							Central African
		Cameroon	Republic					Cameroon	Republic
		4					4
		3					3
	Latitude	2					2
						Democratic				Democratic
		Congo				Republic				Congo	Republic
						of Congo				of Congo
		1					1
		Gabon							Gabon
		0					0
		14	16		18				14	16	18
		Longitude							Longitude
			0	100	200	300	400	500	600	700	800

Table 4 summarises the E k 's for both SCGLR and mixed-SCGLR methods. Optimal parameter value triplet (K , s , l ) =

Table 5 :

 5 Mean Relative Squared Error (MRSE) values obtained with binomial and Poisson distributions. The R package glmmLasso does not handle binomial outcomes but only Bernoulli ones, which precludes comparison in this case.

		GLMM (no regularisation)	GLMM-LASSO	mixed-SCGLR
		Binomial	Poisson	Poisson	Binomial	Poisson
	τ = 0.1	2.31	0.50	0.31	0.51	0.45
	τ = 0.3	3.07	0.60	0.33	0.28	0.18
	τ = 0.5	3.93	0.75	0.39	0.15	0.09
	τ = 0.7	6.50	1.07	0.40	0.10	0.07
	τ = 0.9	19.29	2.71	0.42	0.07	0.05

Table 7 and

 7 Table 8 present the results for respectively 150 and 200 explanatory variables. They give the Mean Relative Squared Error (MRSE) values for β k , k ∈ {1, . . . , 4}, as well as biases and standard errors of estimated variance components, obtained on 20 samples for each value of τ .

	Table 7: Mean Relative Squared Error (MRSE) values for fixed-effect estimates, and biases
	and standard errors for estimated variance components, obtained with 100 observations and
	150 explanatory variables.
	β 1	β 2	β 3	β 4	σ 2 1	σ 2 2	σ 2 3	σ 2 4
	τ = 0.3 0.06 0.26 0.19 0.13	-0.01 -0.03 -0.02 0.09 0.09 0.03	0.02 0.06
	τ = 0.5 0.03 0.20 0.10 0.07	0.01 -0.03 0.11 0.08	0.00 0.07	0.01 0.07
	τ = 0.7 0.01 0.10 0.05 0.04	0.01 -0.05 0.07 0.09	0.01 0.10	0.02 0.07

Table 8 :

 8 Mean Relative Squared Error (MRSE) values for fixed-effect estimates, and biases and standard errors for estimated variance components, obtained with 100 observations and 200 explanatory variables.

	β 1	β 2	β 3	β 4	σ 2 1	σ 2 2	σ 2 3	σ 2 4
	τ = 0.3 0.06 0.15 0.18 0.10	-0.04 -0.05 0.04 0.09	0.01 -0.02 0.05 0.05
	τ = 0.5 0.03 0.17 0.09 0.05	-0.05 0.06	0.00 -0.02 -0.01 0.19 0.04 0.04
	τ = 0.7 0.01 0.15 0.04 0.03	0.03 0.08	0.00 -0.01 -0.02 0.14 0.05 0.05

Some component planes are given on Figure

9

(150 explanatory variables) and Figure

10

(200 explanatory variables).

ACKNOWLEDGMENTS

The extended data Genus required the arrangement and the inventory of 140.000 developed plots across four countries : Central African Republic, Gabon, Cameroon and Democratic Republic of Congo. The authors thank the members of the CoForTips project for allowing the use of this data. We are also grateful to the editor, the associate editor and to the referees for their thorough and constructive review of this work.

Code for running simulations: We also provide the R codes required to reproduce most of the simulation results (R and Rdata files).

SUPPLEMENTARY MATERIAL

In this section, we simply extend the simulation scheme presented in Section 5.5 to binomial and Poisson outcomes. We maintain design matrices X and U as defined in Section 5.1.

Fixed-effect parameters θ k 's and random-effect vectors ξ k 's are defined in Section 5.5.

Given ξ 1 and ξ 2 , we then simulate

Table 5 gives the Mean Relative Squared Error (MRSE) values for θ 1 and θ 2 obtained on 100 samples for each value of τ .

For the Poisson distribution, the results in Table 5 are essentially identical to those in the article: mixed-SCGLR outperforms GLMM-LASSO [START_REF] Groll | glmmLasso: Variable Selection for Generalized Linear Mixed Models by L1-Penalized Estimation[END_REF], R package glmmLasso)

except for τ = 0.1. As for the binomial distribution, the regularisation provided by mixed-SCGLR improves the results obtained without regularisation (Bates et al., 2015, R package lme4), regardless of the level of redundancy within the explanatory variables. Unsurprisingly, the errors are much smaller than in the binary case.

Each explanatory variable is normally simulated with mean 0 and variance 1. Parameter τ ∈ {0.3, 0.5, 0.7} tunes the level of redundancy within each bundle: the correlation matrix of bundle X j is cor (X

where p j is the number of variables in X j . For each k ∈ {1, 2, 3, 4}, random-effect vectors are simulated as ξ k ind.

Response y 1 is predicted only by bundle X 1 , y 2 only by bundle X 2 , y 3 only by bundle X 3 , y 4 by both bundles X 2 and X 3 , and bundle X 0 plays no explanatory role. Our choice for the fixed-effect parameters is