Laurène Perthame
email: laurene.perthame@agroparistech.fr

Frédéric Rees

Xavier Cornilleau

Céline Richard-Molard

Christophe Pradal

Alexandra Jullien

SIMBAL: A structural-functional plant model to simulate C and N dynamics and shoot-root architecture of winter oilseed rape associated with legumes

Keywords: plant model, fababean, plant-plant interactions, FSPM for agroecology

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Species mixtures are of great interest to promote low-input agricultural practices while maintaining agricultural production. For example, winter oilseed rape sown with a frost sensitive legume reduces the use of nitrogen fertilizer and herbicide while maintaining yield at levels equivalent to pure cropping (Lorin et al., 2016). A disadvantage of such mixtures is that the legume might compete with the cash crop for resources and reduce yield. To optimize management of these mixtures, it is necessary to better understand the ecophysiological processes that drive the sharing of resources between plants (e.g. carbon and nitrogen). Mechanistic functional-structural modeling (FSM) is particularly suitable for such studies as it allows the analysis of the processes underlying plant-plant interactions and integrates both the architectural growth processes of the plant and its functioning (e.g. acquisition/allocation of resources) [START_REF] Louarn | [END_REF]. For example, the model Virtual Grassland [START_REF] Louarn | The future of European Grasslands[END_REF] simulates competition for C, N and water between several herbaceous and legume species. However, no FSM for mixtures of annual crops with frost sensitive legume that considers the feedback of the carbon/nitrogen functioning in both shoots and roots has been developed so far. Such a modeling approach is shown in this study, with the aim of filling knowledge gaps on the role of carbon-nitrogen interactions in plastic responses of plants in mixtures, by taking the case study of winter oilseed rape sown with frost sensitive fababean at vegetative stage.

Conceptual framework of C and N fluxes

The model SIMBAL (SIMulated Brassica Associated with Legumes) is based on the ARNICA FSM developed for Arabidopsis thaliana (Richard-Molard et al., 2009). The model simulates shoot and root architecture of winter oilseed rape in which C and N fluxes determine the growth of each organ (leaf, internode, root segment). C is preferentially given to shoot and N to roots. As showed in Fig. 1, the total internal N amount results from root uptake and remobilization from storage (1) and determine the increment of leaf and internode area (4). N in leaves is distributed along the light gradient within the plant, while N in internodes is distributed homogeneously (3). The total internal C quantity results from photosynthesis (6) and remobilization of C storage (9). Root growth results from the N available to roots (2) and the C available after satisfaction of shoot demand (7). C storage is filled if total internal C quantity is superior to C organ demands (calculated as the C quantity necessary to produce a unit of leaf, internode or root surface) (9). N storage is filled from leaf senescence (driven by a leaf lifespan) (5) and when N quantity for shoot is superior to N used for actual total photosynthetic area (8). When the crop is grown in association with fababean rather than in monoculture, the light availability to the plant is reduced and parameters value such as phyllochron are modified based on experiments. A first FSPM was built to study the N and C fluxes and organ growth within a whole plant of winter oilseed rape at the vegetative stage in OpenAlea (Pradal et al., 2008). This conceptual framework aims to be generic and adaptable to other species, in particular frost sensitive fababean. In the near future, we intend to develop the model upon the entire crop cycle of winter oilseed rape and to assess the effect of frost sensitive fababean on oilseed rape growth and CN functioning through direct (competition for light harvest and N uptake) and indirect (degradation of N-rich legume residues) effects.

  SIMBAL: A structural-functional plant model to simulate C and N dynamics and shoot-root architecture of winter oilseed rape associated with legumes Laurène Perthame, Frédéric Rees, Xavier Cornilleau, Céline Richard-Molard, Christophe Pradal, Alexandra Jullien To cite this version: Laurène Perthame, Frédéric Rees, Xavier Cornilleau, Céline Richard-Molard, Christophe Pradal, et al.. SIMBAL: A structural-functional plant model to simulate C and N dynamics and shoot-root architecture of winter oilseed rape associated with legumes. FSPM2023 -10th International Conference on Functional-Structural Plant Model, Mar 2023, Berlin, Germany. hal-04262992

Fig. 1 :

 1 Fig.1: Conceptual framework of C and N dynamics and organ growth underlying the SIMBAL model


Julien Gibaud
email: julien.gibaud@umontpellier.fr

Xavier Bry

Catherine Trottier

Generalized linear model based on latent factors and supervised components

Keywords: EM algorithm, Factor model, Latent variables, Multivariate Generalized Linear Model, Supervised components

In a context of component-based multivariate model we propose to model the residual dependence of the responses. Each response of a response matrix is assumed to depend, through a Generalized Linear Model, on a set of explanatory variables, as well as on a set of additional covariates. Explanatory variables are partitioned into conceptually homogeneous variable groups, viewed as explanatory themes. Variables in themes are supposed many and redundant. Thus, generalized linear regression demands dimension reduction and regularization with respect to each theme. By contrast, additional covariates contain few variables, selected so as not to be too redundant, thus demanding no regularization. Regularization is performed searching each theme for an appropriate number of orthogonal components that both contribute to predict the responses and capture relevant structural information in themes. A small set of latent factors completes the model so as to model the covariance matrix of the linear predictors of the responses conditional on the components. To estimate the multiple-theme model, we present an algorithm combining thematic component-based model estimation and factor model estimation. This methodology is tested on simulated data and then applied to an agricultural ecology dataset.

Introduction

In multivariate modeling, accounting for mutual dependencies between the responses is a rich ongoing research field in both theoretical and applied statistics (Meyer, 2009). In an ecological context for instance, the species co-occurrences that are not explained by the environmental variables demand to be modeled. With this aim in mind, we model the residual dependence of the responses assumed drawn from Generalized Linear Models (GLM, McCullagh and Nelder, 1989), by means of random numeric latent variables called factors, which we introduce in the linear predictors. Estimating the residual covariance matrix will allow us to identify groups of species. In a Bayesian context, [START_REF] Bystrova | Clustering species with residual covariance matrix in joint species distribution models[END_REF] profit from the clustering properties of the Dirichlet process to group the species. Indeed, the species in the same cluster have the same rows in the residual correlation matrix. In the same vein, we propose to group the species sharing strong (positive or negative) residual correlations. To do this, we estimate the residual variance-covariance matrix of the linear predictors, and then perform clustering on a dissimilarity matrix built from the estimated correlation matrix. Identifying a group of correlated species is equivalent to identifying a square block of high absolute values in the correlation matrix, once the rows and columns are suitably reordered.

An approach, named Generalized Linear Latent Variable Model (GLLVM), has been proposed by [START_REF] Skrondal | Generalized latent variable modeling: Multilevel, longitudinal, and structural equation models[END_REF] to combine GLM with random latent variables. Unfortunately, when factors are involved, the log-likelihood cannot be derived from the GLLVM analytically. Several works propose to maximize this log-likelihood but, some of them suffer from a high computation time. We may cite for instance the works using the adaptive quadrature [START_REF] Rabe-Hesketh | Reliable estimation of generalized linear mixed models using adaptive quadrature[END_REF], the Expectation Maximization algorithm (EM, [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] in conjunction with Monte Carlo integration [START_REF] Hui | Model-based approaches to unconstrained ordination[END_REF] or the works using Bayesian Markov Chain Monte Carlo (MCMC) [START_REF] Hui | boral-Bayesian ordination and regression analysis of multivariate abundance data in R[END_REF][START_REF] Tikhonov | Joint species distribution modelling with the R-package Hmsc[END_REF]. Few methods reduce the computation time by using a closed form approximation of the log-likelihood. The works employing a variational approximation [START_REF] Hui | Variational approximations for generalized linear latent variable models[END_REF], a Laplace approximation [START_REF] Niku | Generalized linear latent variable models for multivariate count and biomass data in ecology[END_REF](Niku et al, , 2019a) ) or an extended variational approximation [START_REF] Korhonen | Fast and universal estimation of latent variable models using extended variational approximations[END_REF] deserve mentioning. To better reduce the computation time, we propose to use a modeling approach based on [START_REF] Saidane | Generalized linear factor models: A new local EM estimation algorithm[END_REF], which assumes that maximization should be performed through the EM algorithm after linearizing the and assuming the linearized model is approximately Gaussian. As demonstrated by [START_REF] Wolfinger | Generalized linear mixed models a pseudolikelihood approach[END_REF], the linearization method proves to be an interesting alternative to estimate parameters in a context of intractable likelihood. However, modeling the responses also requires taking into account a large set of possibly highly correlated explanatory covariates, so that the GLLVM demands regularization. This can be carried out by means of component-based dimension reduction. This consists in assuming that there exists a small number of common latent explanatory dimensions, which we aim to capture through as many linear combinations of the explanatory variables, named components.

In the wake of [START_REF] Marx | Iteratively reweighted partial least squares estimation for generalized linear regression[END_REF], the Supervised Component-based Generalized Linear Regression (SCGLR) was proposed by [START_REF] Bry | Supervised component generalized linear regression using a PLS-extension of the Fisher scoring algorithm[END_REF] in order to bridge the multivariate GLM estimation with the component-based dimension reduction of the explanatory space. Moreover, SCGLR allows both to find strong interpretable explanatory directions, and to produce regularized predictors in the high-dimensional framework (when the covariates outnumber the statistical units). Indeed, SCGLR optimizes a general and flexible trade-off criterion between the Goodness-of-Fit (GoF) of the model and the Structural Relevance (SR, [START_REF] Bry | THEME: THEmatic Model Exploration through multiple co-structure maximization[END_REF] of relevant directions of the space spanned by the explanatory variables. To extend this method, Bry et al (2020b) developed THEME-SCGLR with the aim to search for the components in a thematic partitioning of the explanatory variables, hence referred to as "themes". Within each theme, the components are required to extract the information that is useful to predict the responses, when associated with the components of the other themes. An package SCGLR is freely available at https://github.com/SCnext/SCGLR. All the aforementioned extensions assume that the responses are independent conditional on the explanatory covariates. We now propose to overcome this limitation by allowing the responses to have some conditional dependence, which we model by introducing common factors into their linear predictors. Besides this dependence, we may also need to take into account a thematic partitioning of the explanatory variables. We thus have to extend THEME-SCGLR to the GLLVM.

The paper is organized as follows. In Section 2, we recall the principle of the original SCGLR and its thematic extension. Section 3 presents our extension of SCGLR to the factor models. Section 4 details several simulation studies that illustrate the interest and the good performances of the proposed algorithm. Section 5 presents the results it yields on an agricultural ecology dataset. Finally, Section 6 provides a conclusion and a discussion.

SCGLR with a partitioning of explanatory variables 2.1 Preliminary notations

The sequel contains mathematical developments which use notations listed hereafter: • Let w n be the weight of unit n, and W = diag(w n ) n=1,...,N . Let R N be endowed with metric W , and let A ∈ R N ×P be a matrix. The W -orthogonal projector onto span[A] is given by Π

• Let a, b ∈ R N be
W span[A] = A A T W A -1 A T W .
The cosine of the angle between a vector b ∈ R N and span[A] with respect to metric W is given by cos

W (b, span[A]) = cos W b, Π W span[A] b .

The original SCGLR context

In the framework of a multivariate GLM we consider K response-vectors measured on N statistical units, encoded in a response matrix Y = [y 1 , . . . , y K ] ∈ R N ×K , to be predicted through explanatory variables partitioned in two groups. The first one A = [a 1 , . . . , a Q ] ∈ R N ×Q is a group of covariates that are few and weakly or not redundant. These variables are assumed to be interesting per se, and their marginal effects have to be taken into account explicitly in the model. The second group X = [x 1 , . . . , x P ] ∈ R N ×P gathers numerous and possibly highly redundant covariates, considered as proxies to latent dimensions, which must be found and interpreted. Thus, the matrix X demands dimension reduction and regularization, contrary to A.

To achieve this, SCGLR searches for explanatory components in X jointly supervised by the response set. A component f ∈ R N writes f = Xu, where u ∈ R P is a loading vector. For a single component model, the linear predictor associated with response y k is given by

η k = (Xu) γ k + Aδ k ,
where γ k and δ k are regression parameters. Component f is common to all the responses and, for identification, we impose u T u = 1. It was assumed originally that the responses are independent conditional on the explanatory variables, and consequently on f .

The original SCGLR specific criterion

For parameter estimation, SCGLR takes advantage of the GLM background. Here, we make use of the Fisher Scoring Algorithm (FSA). Let h k denote the canonical link function associated with the response y k , h k its first derivative and µ nk the mean parameter for statistical unit n. The working variable w nk associated with y nk is then calculated as the first order development of h k at point µ nk

w nk = h k (µ nk ) + (y nk -µ nk ) h k (µ nk ) = η nk + ζ nk ,
where ζ nk = (y nk -µ nk ) h k (µ nk ). In the spirit of [START_REF] Nelder | Generalized Linear Models[END_REF], this development leads to the conditional linearized model expressed column-wise

w k = (Xu) γ k + Aδ k + ζ k , with E[ζ k ] = 0 and V[ζ k ] =: W -1 k .
Due to the product uγ k , this linearized model derived from the FSA is not linear and must be estimated through an alternated weighted least squares process, estimating in turn {γ k , δ k } and u.

Let Π Wk span[Xu,A] be the projection on span[Xu, A] with respect to W k . As suggested by [START_REF] Bry | Supervised component generalized linear regression using a PLS-extension of the Fisher scoring algorithm[END_REF], the loading vector u may be viewed as the solution of the optimization program max

u T u=1 ψ A (u), where ψ A (u) := K k=1 w k 2 W k cos 2 W k w k , Π Wk span[Xu,A] w k .
The sub-criterion ψ A is merely a Goodness-of-Fit (GoF) measure, and maximizing it does not lead to strong and interpretable components. The GoF measure must therefore be aptly combined with an other sub-criterion to achieve both meaningful and predictive dimension reduction, together with regularization. [START_REF] Bry | THEME: THEmatic Model Exploration through multiple co-structure maximization[END_REF] proposed such a sub-criterion, named Structural Relevance (SR) to measure the ability of a component to capture information in a set of variables containing structures such as variable-bundles. Informally, a bundle is a set of variables correlated "enough" to be viewed as produced by a common latent dimension. Assuming that X consists of P standardized numeric variables, the associated SR measure φ is defined as the following generalized average of quadratic forms of u φ(u) := 1

P P p=1 Xu, x p 2l W 1/l , ( 1
)
where W is the weight matrix (in practice, W = 1 N I N ). Components will line up with a more or less thin bundle depending on whether l ≥ 1 is greater or smaller respectively. The main objective is to focus on the most interpretable directions.

The SCGLR specific criterion, proposed by Bry et al (2020b), introduces a hyperparameter s ∈ [0, 1] to tune the importance of the SR relative to the GoF. SCGLR thus attempts a trade-off between φ and ψ A by solving max

u T u=1 s ln(φ(u)) + (1 -s) ln(ψ A (u)).
When s = 0, the criterion identifies with the GoF, while at the other end, taking s = 1 makes it identify with the SR. Increasing s intensifies both the focus of components on "strong" dimensions, and the regularization, at the cost of some GoF. The explicit expression of the criterion is given in Appendix A.

THEME-SCGLR

Bry et al (2020b) refer to the "thematic model" as the conceptual model stating that variables in Y are dependent on R themes X 1 , . . . , X R plus a set of covariates A, and that structurally relevant dimensions should be explicitly identified in the X r 's. For a single component in each theme, the linear predictor associated with response y k is then given by

η k = (X 1 u 1 ) γ k1 + • • • + (X R u R ) γ kR + Aδ k .
To achieve theme-specific regularization, the SCGLR criterion has to be adapted. Denoting f r = X r u r the (first) component of theme X r , we have

Π Wk span[f1,...,fR,A] = Π Wk span[fr ,Ar ] where A r = [f 1 , . . . , f r-1 , f r+1 , . . . , f R , A].
For each theme, the GoF measure thus becomes

ψ Ar (u r ) := K k=1 w k 2 W k cos 2 W k w k , Π Wk span[Xr ur ,Ar ] w k .
The SR measure remains the same φ (u r ) as given by Equation (1). Finally, the optimization program can be solved by iteratively maximizing in turn the trade-off criterion on every u r ∀r, max

u T r ur =1
s ln(φ(u r )) + (1 -s) ln(ψ Ar (u r )).

(2)

This combined criterion is quite general. Indeed, the GoF measure adapts any situation where a likelihood function is available for the model taking the components and A as covariates. Generally, this likelihood involves a vector of parameters Θ. The maximization is carried out alternating two steps: (i) Given Θ, maximize the criterion with respect to each u r using a dedicated algorithm: PING (for Projected Iterated Normed Gradient) recalled in Appendix B, designed to maximize, at least locally, any criterion on the unit sphere [START_REF] Chauvet | Component-Based Regularization of Multivariate Generalized Linear Mixed Models[END_REF]Bry et al, 2020a,b;[START_REF] Gibaud | Response mixture models based on supervised components: Clustering floristic taxa[END_REF]. (ii) Given all u r , maximize the criterion with respect to Θ. This step is performed using a classical likelihood maximization algorithm relevant to the situation.

Higher rank components

Consider step (i) of the combined criterion maximization, and suppose we want to extract a given number of components H r from each theme X r . Let f h r = X r u h r be the rank-h component of theme X r , and let

F h r = [f 1 r , . . . , f h r ]
, where h ≤ H r , be the matrix of the first h components of theme X r . The new component f h+1 r must best complement both the existing ones and A, that is

A h r := [F H1 1 , . . . , F H r-1 r-1 , F h r , F H r+1 r+1 , . . . , F HR R , A]. So f h+1 r
has to be calculated using A h r as the new set of additional covariates. Moreover, to avoid linear redundancy of components, we impose that f h+1 r be orthogonal to F h r , i.e. F h T r W f h+1 r = 0. We calculate every new component as the solution of program (2), with the additional constraint: ∆ h T r u h+1 r = 0, where ∆ h r = X T r W F h r , and loop on r until overall convergence of the component system. For all r = 1, . . . , R, the rank-1 component of theme X r is calculated using the same program with F 0 r = ∅ and ∆ 0 r = 0.

3 Extending SCGLR to a factor model

As mentioned above, step (ii) of the combined criterion maximization boils down to maximizing the likelihood of the component-based model. In this section, the components are thus taken as known. For the sake of simplicity, we shall consider the matrix F = [F H1 1 , . . . , F HR R ] as the new set of explanatory variables and γ k = (γ k1 , . . . , γ kR )

T its vector of regression parameters associated with the response y k .

SCGLR in a factor model context

Let Y = [y 1 , . . . , y K ] ∈ R N ×K be the response matrix. For unit n, each response is assumed to be linearly modeled using the components and additional covariates, plus J random latent factors g n = (g n1 , . . . , g nJ )

T η nk = f T n γ k + a T n δ k + g T n b k ,
where f n and a n are the vectors composed of the nth rows of matrices F and A respectively, and b k is the vector of regression parameters associated with g n . The factors are assumed drawn from a multivariate normal distribution g n ∼ N J (0, I J ) and independent across statistical units. This model is designed so that the J factors capture as much as possible of the covariance between the responses not accounted for by the components and additional covariates, hence their residual covariance. Denoting G ∈ R N ×J the matrix containing all the realizations of factors, the linear predictor associated with the response y k expressed column-wise becomes

η k = F γ k + Aδ k + Gb k . Let B = [b 1 , . . . , b K ] ∈ R
J×K be the loading matrix. [START_REF] Jöreskog | A general approach to confirmatory maximum likelihood factor analysis[END_REF] notices that the loading matrix B is defined up to an arbitrary orthogonal rotation. To guarantee the identification of the model, we choose to constrain the J × J sub-matrix of B to be an upper triangular matrix with positive diagonal elements [START_REF] Geweke | Measuring the pricing error of the arbitrage pricing theory[END_REF]. An advantage of the factor model is to yield the matrix Σ = B T B ∈ R K×K , storing the residual covariances of the responses, in a parsimonious manner. Indeed, the number of factors retained may remain small with respect to the size of the covariance matrix.

Estimating the parameters of a GLM with factors

Let Θ = {γ k , δ k , b k | k = 1, . . . ,
= N n=1 ln K k=1 L(y nk | g n ; Θ)L(g n ) dg n .
In a context of non-Gaussian responses, the maximization of this log-likelihood is not allowed. In the spirit of [START_REF] Saidane | Generalized linear factor models: A new local EM estimation algorithm[END_REF], the estimation of the parameters is performed in two steps: first, we linearize the model; then, we maximize the pseudolikelihood of the linearized model under a Gaussian assumption.

The linearization step

Temporarily considering the factors given, i.e. conditional on G, the above loglikelihood is that of a classic multivariate GLM. The working variable w k could be viewed as the response in the linearized model

w k = F γ k + Aδ k + Gb k + ζ k , where E[w k | G] = F γ k + Aδ k + Gb k and V[w k | G] = V[ζ k ] = W -1 k .

The estimation step

In this step, we assume that the distribution of the working variables given F , A and G is Gaussian, and view the factors as latent variables. The model pseudo-loglikelihood l(Θ; W), where W denotes the matrix of working variables, being difficult to maximize directly, we use the EM algorithm to estimate the model parameters.

Assuming in the wake of [START_REF] Wolfinger | Generalized linear mixed models a pseudolikelihood approach[END_REF] that the working variables have a Gaussian distribution, we calculate and then maximize the expectation of their complete log-likelihood l(Θ; W, G). Further details of the EM algorithm are given in Appendix C.

The overall algorithm

The overall algorithm, presented in Appendix D, consists in alternating the following steps: (i) Given the current set of parameters, calculate all the components of all the themes iteratively through the PING algorithm. (ii) Given the current components, calculate the working variables of the linearized model and their variance matrix. (iii) Given the working variables, estimate the factor model parameters through the EM algorithm. The method thus implemented is named F-SCGLR (for Factor-SCGLR).

The clustering steps

Recall that the final aim of this work is to group the responses according to their mutual dependencies. In other words, two responses having a high residual correlation (positive or negative) should be cast to the same group. To achieve this, we propose the following strategy 5. Perform a K-means algorithm (taking as a starting point the output of a hierarchical clustering procedure) on the coordinates obtained on the previous step. We use the factoextra package [START_REF] Kassambara | Package 'factoextra[END_REF] where the function hkmeans runs the K-means and the function fviz-nbclust optimizes the number of clusters using the silhouette criterion.

Simulation study

Several simulation studies have been implemented to assess the performance of F-SCGLR. The first one focuses on the identification of the right combination of components and factors. The combination was calibrated across the cross-product grid (H 1 , . . . , H R , J) ∈ {1, . . . , 4} R × {0, . . . , 5} by minimizing the Bayesian Information Criterion (BIC, [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. As shown by [START_REF] Chauvet | Component-Based Regularization of Multivariate Generalized Linear Mixed Models[END_REF], the hyper-parameters must be chosen to avoid the components to be too close to the principal components (s > 0.5) or to be drawn towards too local bundle (l > 10). Thus, the second simulation aims at studying the influence of the hyper-parameters s ∈ {0.1, 0.3, 0.5} and l ∈ {1, 2, 3, 4, 7, 10} in a situation of a more or less clearly separated cluster pattern. In this simulation, we use the Rand Index (RI, [START_REF] Rand | Objective criteria for the evaluation of clustering methods[END_REF]) and the Adjusted Rand Index (ARI, [START_REF] Hubert | Comparing partitions[END_REF] to assess the correctness of the classification steps detailed in Section 3.4. In addition, to measure the quality of the latent dimensions recovery, we calculate the maximum square correlation between each true dimension, represented by their direction vector ξ used for simulation, and the components

ρ 2 (ξ, .) = max r,h ρ ξ, f h r 2 ,
where f h r denotes the hth component of theme X r . Finally, as reference values for comparison, we also calculated the RI and ARI of the partitions output by a competing package in a context of binary data. For each simulation, one hundred samples have been generated. The package FactorSCGLR and the simulation codes are available at https://github.com/julien-gibaud/FactorSCGLR.

Simulation in a context of mixed distributions

Generation of the simulated data

The variables are simulated on N = 100 statistical units. Five latent dimensions ξ 1 , ξ 2 , ξ 3 , ξ 4 and ξ 5 are simulated independently. The X matrix consists in two themes:

X = [X 1 , X 2 ]. The first theme X 1 = [X 1 , X 2 , M 1 ] is made of three blocks: X 1 ∈ R N ×90 and X 2 ∈ R N ×60
are bundles of variables distributed about ξ 1 and ξ 2 respectively, and M 1 contains fifty unstructured noise variables drawn from a multivariate normal distribution. Likewise, the second theme X 2 = [X 3 , X 4 , X 5 , M 2 ] is made of four blocks: X 3 ∈ R N ×100 , X 4 ∈ R N ×80 and X 5 ∈ R N ×60 are bundles of variables distributed about ξ 3 , ξ 4 and ξ 5 respectively, and M 2 contains sixty unstructured noise variables drawn from a multivariate normal distribution. More formally, for all i = 1, . . . , 5, a variable x p within a bundle is simulated as x p = ξ i +ε p , where ε p ∼ N N (0, 0.1I N ). This generation leads to P = 500 explanatory variables.

The N realizations of the J = 3 factors, simulated through g n ∼ N J (0, I J ), are stored in matrix G ∈ R N ×J . The matrix B ∈ R J×K of factor loadings is generated so as to exhibit a three-cluster pattern

∀k = 1, . . . , 5, b k ∼ N J µ 1 , σ 2 B I J ∀k = 6, . . . , 10, b k ∼ N J -µ 1 , σ 2 B I J ∀k = 11, . . . , 20, b k ∼ N J µ 2 , σ 2 B I J ∀k = 21, . . . , 35, b k ∼ N J -µ 2 , σ 2 B I J ∀k = 36, . . . , 50, b k ∼ N J µ 3 , σ 2 B I J ,
where σ 2 B = 0.1, µ 1 = (2, 0, 0) T , µ 2 = (0, -1, 0) T and µ 3 = (0, 0, 1.5) T . Finally, the response matrix Y is simulated as a mix of Gaussian, Poisson and Bernoulli distributions, with

∀k = 1, . . . , 20, y k ∼ N N µ = γ 1k ξ 1 + γ 2k ξ 2 + Gb k , Σ = σ 2 k I N ∀k = 21, . . . , 40, y k ∼ P (λ = exp [0.5γ 1k ξ 4 + 0.5γ 2k ξ 3 + Gb k ]) ∀k = 41, . . . , 50, y k ∼ B p = logit -1 [γ 2k ξ 3 + γ 3k ξ 2 + Gb k ] ,
where for all k, σ 2 k , γ 1k , γ 2k and γ 3k are uniformly generated, with

σ 2 k ∈ [0.1, 0.2], γ 1k ∈ [-4, 4], γ 2k ∈ [-2
, 2] and γ 3k ∈ [-0.5, 0.5]. In the linear predictors, we range the latent dimensions following their regression parameters.

Identification of the true model

In this simulation, the hyper-parameters are first calibrated through the SCGLR package (e.g. without factors) and set to s = 0.3 and l = 4. Appendix E sums up the results on a cross-product grid. As expected, the combination which minimizes the BIC is given by the true combination (H 1 , H 2 , J) = (2, 2, 3). However, several points deserve mentioning. We observe, for all component combinations, that the values of BIC decrease dramatically when factors are involved in the model. This shows that, because mutual dependencies may generally exist, the residual covariance should be modeled. When the model involves too many factors (when J = 4 and J = 5), the number of useful components is underestimated. Indeed, the variability of the model captured by the factors then contains a part of the variability otherwise captured by the components. In the opposite situation, when J = 0 or J = 1, the BIC leads to overestimate the number of components.

Varying the hyper-parameters and the variance within the clusters

Henceforth, keeping the true combination found by the BIC, we focus on the influence of the hyper-parameters s and l on the classification decision and latent dimension recovery. In order to compare the results in a context of more or less distinct cluster pattern, we vary the variance within the cluster by taking σ 2 B ∈ {0.1, 0.2, 0.3}. Figure 1 shows the residual correlation matrices for the three values of σ 2 B .

-1 -0.5 0 0.5 1 Table 1 gives the results for σ 2 B = 0.1. For the s = 0.3 and s = 0.5 cases, the values of RI and ARI are slightly better than the s = 0.1 case. Moreover, for s = 0.3 and s = 0.5, the maximum value for RI and ARI is reached for l = 4. This is in accordance with the hyper-parameters calibrated through the SCGLR package. The Tables summing up the results for σ 2 B = 0.2 and σ 2 B = 0.3 are presented in Appendix F. As expected, the higher the variance within the cluster, the weaker the values of RI and ARI for all the combinations of s and l. We may also note that the difference between the values of RI and ARI across the hyper-parameters s and l tends to fade when σ 2 B increases. The main result about the square correlations is that the variance within the cluster does not have a relevant influence on the quality of the latent dimensions recovery. Indeed, the search for components is related to the deterministic part of the model, while σ 2 B is involved in the stochastic one. The square correlations, for s = 0.3 and s = 0.5 with l ≥ 2, are greater than for s = 0.1. This observation is consistent with [START_REF] Chauvet | Component-Based Regularization of Multivariate Generalized Linear Mixed Models[END_REF] who notice that the thinner the bundles, the greater the value of s has to be to recover the latent dimensions correctly. Here, indeed, the variance within the bundles is equal to 0.1 (thin bundles). However, the particular case of l = 1 deserves mentioning. The components calculated with l = 1 being close to the principal components, the two components of theme X 2 are drawn between the three bundles and so, produce low square correlations with the latent dimensions. The interest of tuning the locality is shown by the gap between the results obtained for l = 1 and l = 2: in the latter case, the square correlations are dramatically better. Furthermore, ξ 3 being the less explanatory latent dimension, ρ 2 (ξ 3 , .) is always lower than the other square correlations.

Figure 2 shows the correlation scatterplots in the component planes (1, 2) for the first two themes. The components are almost perfectly aligned with the explanatory bundles. However, as observed in Table 1, the bundle X 3 seems slightly less correlated with the component f 2 2 than the other bundles with their corresponding components.

Comparative study

To compare the different GLLVM implementations, we use the package gllvm (Niku et al, 2019b). This package offers three ways to perform GLLVM: using a variational approximation (VA,

Generation of the simulated data

The variables are simulated on N ∈ {100, 200, 300} statistical units. For the sake of simplicity, a bundle X of ten variables distributed about the latent dimension ξ is generated. One categorical variable with three levels is taken as only additional covariate A. In this simulation, J = 2 factors are simulated to model the residual covariance of the K ∈ {10, 30, 50} responses. The regression coefficients of the factors are generated in order to get a two-cluster design

∀k = 1, . . . , 0.4K, b k ∼ N J (-1) k µ 1 , 0.1I J ∀k = 0.4K + 1, . . . , K, b k ∼ N J (-1) k µ 2 , 0.1I J ,
where µ 1 = (0, 2) T and µ 2 = (1.5, 0) T . The comparable package not allowing to consider different distribution families for the responses, we restricted the comparison to binary outcomes

∀k = 1, . . . , K, y k ∼ B p = logit -1 [γ k ξ + Aδ k + Gb k ] ,
where for all k, γ k and δ k are uniformly generated, with γ k ∈ [-4, 4] and δ k ∈ [-1, 1]. Figure 3 shows the residual correlation matrices obtained for the three values of K.

Results of the comparison

In this simulation the package SCGLR calibrates the hyper-parameters to s = 0.5 and l = 1, while the BIC selects H 1 = 1 and J = 2. Table 2 sums up the RI, ARI, Procrustes error and computation time output by, on the one hand, our package FactorSCGLR performing the F-SCGLR method, and, on the other hand, the package gllvm implementing the VA, LA and EVA approaches. We observe that, for all combinations of N and K, F-SCGLR gives the best values of RI and ARI, followed by EVA, VA, the shrinkage approach and then LA. Indeed, the highest values obtained of the ARI are respectively: 0.973, 0.617, 0.556, 0.532 and 0.516. Unlike with our package, which performs better when the number of either statistical units or responses increase, a higher number of responses may cause a deterioration of the classification of the compared methods. Likewise for all cases, the lowest Procrustes errors on B and G are reached by F-SCGLR. The ordination performances are not ranked in the same order as those obtained for the classification. The F-SCGLR, VA, EVA, LA and shrinkage methods respectively have an error on B equal to 0.054, 0.614, 0.615, 0.632 and 0.676. However, we may note that all methods could well predict the latent variables. In their best combinations, the values of the Procrustes error on G are equal to 0.133 (F-SCGLR) and 0.168 (LA). When the number of responses increases, the Procrustes error on B for the F-SCGLR, shrinkage and LA methods may get bigger while the Procrustes error on G decreases for all methods. Across the simulations, F-SCGLR is the fastest of the compared methods. The longest computation time (less than 12 seconds) occurred for N = 300 and K = 50, while EVA, VA, shrinkage and LA ran for 89, 80, 112 and 386 seconds respectively in that case. However, contrary to [START_REF] Korhonen | Fast and universal estimation of latent variable models using extended variational approximations[END_REF], we did not observe that EVA ran faster than VA: in the K = 10 and (N, K) = (300, 50) cases, the computation time of VA was lower. The conclusions of our respective works agree nevertheless that LA is relatively slow.

5 Analysis of an agricultural ecology dataset

Data description

We apply F-SCGLR to the dataset available following the link https://doi.org/10.15454/AJZUQN. The sample we consider gives the observation of K = 12 agrobiodiversity variables over N = 54 winter cereal fields in the French Vallées et Coteaux de Gascogne. The agrobiodiversity is reported through three carabid beetle variables (two abundances and one Shannon index), three vascular plant variables (one richness, one relative cover and one Shannon index) and six axes of correspondence analyses (CA) performed on presence-absence data of carabid species and plant species respectively. The three abundance and richness responses are assumed to be samples of Poisson random variables while the other responses are considered normally distributed. To model the agrobiodiversity, we have P = 21 variables partitioned into four themes and Q = 1 additional covariable. The first theme X 1 characterizes the pest control through four variables. Six farming intensity variables make up the second theme X 2 . The third and fourth themes X 3 and X 4 gather six and five variables representing the landscape heterogeneity related to seminatural covers and to the crop mosaic respectively. The binary categorical variable coding the observation year (2016 or 2017) is considered as the additional covariate put into matrix A. For more information about this dataset, we refer the reader to [START_REF] Duflot | Farming intensity indirectly reduces crop yield through negative effects on agrobiodiversity and key ecological functions[END_REF].

Results and interpretation

As in Section 4, we need to calibrate the hyper-parameters. We first tune s and l through the SCGLR package, then we find the best combination of number of components and factors according to the BIC. However, due to the small number of explanatory variables in each theme, we only allow the number of components to reach H r = 3. We thus minimize the BIC on the cross-product grid (H 1 , H 2 , H 3 , H 4 , J) ∈ {0, . . . , 3} 4 × {0, . . . , 5} with the s and l values previously found.

The SCGLR package recommends tuning hyper-parameters to s = 0.5 and l = 1 for this agricultural ecology dataset. Henceforth, the component combination minimizing the BIC is (H 1 , H 2 , H 3 , H 4 ) = (0, 3, 0, 0) meaning that only the farming intensity theme was found relevant for the prediction of the agrobiodiversity. Duflot et al ( 2022) make the assumption that agrobiodiversity is predictable from the farming intensity (theme X 2 ) and the landscape heterogeneity (themes X 3 and X 4 ). The combination found by the BIC validates this hypothesis as to the effect of the farming intensity and the non-effect of the pest control in the prediction of the agrobiodiversity. However the landscape heterogeneity themes proved here useless to this prediction.

We henceforth try to interpret the components of the second theme. The first component f 1 2 is correlated (ρ = -0.924, ρ = -0.794 and ρ = -0.738) with a bundle of three variables, of which "TFI.total" and "TFI.h" represent a treatment frequency index of herbicides, and "nb.op" is the total number of operations conducted by the farmers. The second component f 2 2 is correlated (ρ = 0.779) with the variable "cum.till.depth" measuring the cumulative tillage depth. The quantity of nitrogen denoted "qtyN.kg" is the most correlated explanatory variable (ρ = -0.781) with the last component f 3 2 . Figure 4 represents the correlation plots of the second theme. In this agricultural ecology dataset, three factors are recommended, according to the BIC, to model the residual variance-covariance matrix. By applying the clustering steps given in Section 3.4, four groups of responses are identified. The first group is composed by the three measures of the carabids. The second group gathers the first axis of the CA of the carabids, the plant richness and the plant Shannon diversity index. The carabids' second CA axis, the plant cover and the first and third plants' CA axes make up the third group. Finally, the fourth group contains the carabids' third CA axis and the plants' second CA axis. Figure 5 shows the residual correlation values alongside the Euclidean representation of the responses output by the MDS. The carabids' measure group having very high residual correlations, the distances between the responses composing it are close to 0 leading to a very compact group (in red) on the first principal plane of the MDS. On the contrary on the graph, the weaker the residual correlations, the wider the groups are scattered.

Conclusion and discussion

The original SCGLR was designed to regularize GLM estimation and reduce the explanatory dimension through components, so as to decompose the linear predictor in an interpretable way. It allowed to find strong and interpretable supervised components common to response variables, by achieving a trade-off between Goodness-of-Fit and a Structural Relevance measure. THEME-SCGLR extends SCGLR to a thematic partition of the explanatory variables, allowing to make better use of the complementary between the explanatory themes, both statistically when fitting the model, and conceptually when interpreting the components. F-SCGLR refines THEME-SCGLR in a major way: through a factor model, it models the residual variance-covariance matrix of the responses using a small number of latent factors. This matrix can be used as a basis for clustering, enabling to identify groups of linked responses.

In our simulation study, F-SCGLR proved to behave as expected regarding response clusters. Whenever the clusters were reasonably distinct, the original partitions were recovered. Whatever the dispersion of the regression coefficients within the clusters, it provided components aligned with the simulated latent dimensions. Our package outperforms the package gllvm in three ways: allowed. (iii) The performances of our package are better in terms of computation time and of cluster detection. On the agricultural ecology dataset, we found four groups of responses. Due to very high residual correlations, the first one gathers the measures of the carabids. The other groups are composed by a mix between the plant variables and the axes of the correspondence analyses. However, even though a strong residual covariance between responses may hint at a biological interaction between species [START_REF] Pollock | Understanding co-occurrence by modelling species simultaneously with a Joint Species Distribution Model (JSDM)[END_REF], [START_REF] Poggiato | On the interpretations of joint modeling in community ecology[END_REF] recall that the residual correlations cannot distinguish the biotic from the abiotic effects. Moreover, performing F-SCGLR, we revealed that the treatment by herbicides, the operations conducted by the farmers, In this research, some limitations have been reached. The use of EM algorithm on each step of the overall algorithm involves a high number of iterations. Due to the absence of consensus about the maximization of the log-likelihood, we think that more researches in this topics need to be effected. SCGLR and its extensions suffer of a high number of hyper-parameters involving the use of heuristics to well calibrate the latter. Only Bernoulli, Binomial, Gaussian and Poisson distributions can currently be handled in the FactorSCGLR package. The package should be improved by adding different distributions as Negative Binomial, Zero Inflated Poisson, Tweedie, Gamma, Beta or Exponential, which are allowed in the gllvm package.

Since span[X

h k ] is orthogonal to span[A h ], Π Wk span[Xu,Ah] = Π Wk span[X h k u,Ah] = Π Wk span[X h k u] + Π Wk span[Ah] .
Consequently, by classical Euclidean statistical concepts, we have

cos 2 W k (w k , span[Xu, A h ]) = cos W k (w k , span[Xu, A h ]) cos W k (w k , span[Xu, A h ]) =    Π Wk span[Xu,Ah] w k W k w k W k       w k , Π Wk span[Xu,Ah] w k W k w k W k Π Wk span[Xu,Ah] w k W k    = w k , Π Wk span[X h k u] + Π Wk span[Ah] w k W k w k 2 W k = w k , Π Wk span[X h k u] w k W k w k 2 W k + w k , Π Wk span[Ah] w k W k w k 2 W k .
The goodness of fit measure ψ A h (u) then writes more explicitly

ψ A h (u) = K k=1 w k 2 W k cos 2 W k (w k , span[Xu, A h ]) = K k=1 w k , Π Wk span[X h k u] w k W k + w k , Π Wk span[Ah] w k W k . Now, w k , Π Wk span[X h k u] w k W k = w T k W k Π Wk span[X h k u] w k = w T k W k X h k u u T X h T k W k X h k u -1 u T X h T k W k w k = u T X h T k W k w k w T k W k X h k u u T X h T k W k X h k u . Let, a k := X h T k W k w k w T k W k X h k , b k := X h T k W k X h k and c k := w k , Π Wk span[Ah] w k W k
.

Finally, we have

ψ A h (u) = K k=1 u T a k u u T b k u + c k and ∇ u ψ A h (u) = 2 K k=1 u T b k u a k u -u T a k u b k u (u T b k u) 2 .

Appendix B The PING algorithm

The Projected Iterated Normed Gradient (PING) algorithm is an extension of the Power Iteration algorithm. To find the hth component, we use the PING algorithm which aims at solving any optimization program of the form

max u J h (u), s.t. u T M -1 u = 1 and ∆ T h u = 0, (B1)
where J h is a function of u to maximize and ∆ h an additional constraint matrix. In the SCGLR context, J h (u) is the specific criterion and ∆ h the orthogonal constraint matrix. We rewrite this optimization program by posing

v = M -1/2 u, G h (v) = J h (M 1/2 v) and E h = M 1/2 ∆ h . max v G h (v), s.t. v T v = 1 and E T h v = 0. (B2)
To solve (B2), we must equate to zero the gradient of the following Lagrangian

L(v, λ, η) = G h (v) -λ(v T v -1) -η T E T h v. Setting Γ h (v) = ∇ v G h (v), we have ∇ v L(v, λ, η) = 0 ⇔ Γ h (v) -2λv -E h η = 0 (B3) ⇔ v = 1 2λ (Γ h (v) -E h η) . (B4) Multiplying (B3) by E T h 2λ E T h v =0 = E T h Γ h (v) -E T h E h η ⇔ E T h Γ h (v) = E T h E h η ⇔ η = E T h E h -1 E T h Γ h (v). (B5) Substituting (B5) in (B4), we get v = 1 2λ Γ h (v) -E h E T h E h -1 E T h Γ h (v) = 1 2λ I -E h E T h E h -1 E T h Γ h (v) = 1 2λ Π span[Eh] ⊥ Γ h (v), where Π span[Eh] ⊥ = I -E h E T h E h -1 E T h . Finally, the constraint v 2 = 1 gives v = 1 2λ Π span[Eh] ⊥ Γ h (v) 1 2λ Π span[Eh] ⊥ Γ h (v) = Π span[Eh] ⊥ Γ h (v) Π span[Eh] ⊥ Γ h (v)
, which suggests the basic iteration of the PING algorithm

v (t+1) = Π span[Eh] ⊥ Γ h (v (t) ) Π span[Eh] ⊥ Γ h (v (t) ) . ( B6
)
Let us show that the basic iteration of the PING algorithm follows a direction of ascent. One way to do this is to show that the direction given by the arc (v (t) , v (t+1) ) is a direction of ascent. In other words, show that

v (t+1) -v (t) , Γ h v (t) ≥ 0.
By construction, we know that on every iteration t of the algorithm, t+1) .

v (t) is orthogonal to span[E h ]. Thus, since for all t, v (t) = Π span[Eh] ⊥ v (t) , we have v (t+1) -v (t) , Γ h v (t) = Π span[Eh] ⊥ v (t+1) -v (t) , Γ h v (t) = v (t+1) -v (t) , Π span[Eh] ⊥ Γ h v (t) . Now, Equation (B6) implies that Π span[Eh] ⊥ Γ h (v (t) ) = v (t+1) Π span[Eh] ⊥ Γ h (v (t) ) . So, sgn v (t+1) -v (t) , Γ h v (t) = sgn v (t+1) -v (t) , v (t+1) = sgn v (t+1) 2 -v (t) , v (t+1) = sgn 1 -cos v (t) , v (
Finally, v (t+1) -v (t) , Γ h v (t) ≥ 0.
Although iteration (B6) follows a direction of ascent, it does not guarantee that function G h actually increases on every step. Indeed, we may go too far in such a direction, and overshoot the maximum. However, let us consider

κ (t) = Π span[Eh] ⊥ Γ h (v (t) ) Π span[Eh] ⊥ Γ h (v (t) )
.

Staying close enough to the current starting point on the arc (v (t) , κ (t) ) ensures that function G h increases on every iteration. Indeed, let be the plane tangent to the unit sphere on v (t) and let w denote the unit-vector tangent to arc (v (t) , κ (t) ) on v (t) . Then, there exists τ > 0 such that, w = τ Π κ (t) , and

w, κ (t) = τ Π κ (t) , κ (t) = τ cos 2 (κ (t) , ) > 0.
However, staying too close to the current starting point can impact the convergence speed of the algorithm to reach the maximum. We avoid that by using a one dimensional maximization function (e.g. Gauss-Newton type) to find the maximum of G h on the arc (v (t) , κ (t) ), and take it as v (t+1) . Therefore, we propose two possible generic iterations for the PING algorithm, which deal with this problem. Algorithm 1 and Algorithm 2 present these alternatives. The first one should be preferred, but is less easy to program.

Algorithm 1 PING algorithm 1: while not convergence do 2:

κ (t) ← Π span[Eh] ⊥ Γ h (v (t) ) Π span[Eh] ⊥ Γ h (v (t) ) 3:
Use a Newton-Raphson unidimensional maximization procedure to find the maximum of G h (v) on the arc (v (t) , κ (t) ) and take it as v (t+1) 4: t ← t + 1 5: end while

Appendix C The EM algorithm

We are now dealing with the linearized model, where the factors are latent. So, we shall use the EM algorithm to estimate the parameters. The previous developments lead to the conditional linearized model

w k = F γ k + Aδ k + Gb k + ζ k , where E[w k | G] = F γ k + Aδ k + Gb k and V[w k | G] = V[ζ k ] = W -1 k = diag υ -1 nk n=1,...,N , Algorithm 2 Alternative PING algorithm 1: m ← Π span[Eh] ⊥ Γ h (v (t) ) Π span[Eh] ⊥ Γ h (v (t) ) 2: while G h (m) < G h (v (t) ) do 3: m ← v (t) + m v (t) + m 4: end while 5: v (t+1) ← m 6: t ← t + 1 with υ -1 nk := a nk (φ k )v k (µ nk ) h k (µ nk ) 2
, where a nk and v k are known functions and φ k is the dispersion parameter related to y k . The linearized model expressed row-wise thus writes

w n = Γ T f n + ∆ T a n + B T g n + ζ n , where Γ = [γ 1 , . . . , γ K ], ∆ = [δ 1 , . . . , δ K ], B = [b 1 , . . . , b K ],
and where w n , f n , a n and g n are the vectors composed of the nth rows of matrices W, F , A and G respectively. The expectation and the variance are given by E

[w n ] = Γ T f n + ∆ T a n and V[w n ] = B T B + Υ -1 n , where Υ -1 n = diag υ -1 nk k=1,...,K .
Denoting Θ = {Γ, ∆, B} the set of parameters, the complete log-likelihood writes

l(Θ; W, G) = ln (L(W, G; Θ)) = N n=1 ln (L(w n | g n ; Θ)) + ln (L(g n ; Θ)) = N n=1 -ln (2π) K/2 det Υ -1 n 1/2 - 1 2 w n -Γ T f n -∆ T a n -B T g n T Υ n w n -Γ T f n -∆ T a n -B T g n -ln (2π) J/2 - 1 2 g T n g n = - 1 2 N n=1 K k=1
ln υ -1 nk + g T n g n + (K + J) ln(2π)

+ K k=1 υ nk w nk -f T n γ k -a T n δ k -g T n b k 2 .

C.1 The expectation (E) step

We Thus, we need to first find the law of g n | w n . Since the random vector (w T n , g T n ) T is Gaussian, we have

w n g n ∼ N Γ T f n + ∆ T a n 0 , B T B + Υ -1 n B T B I J .
Thanks to the conditioning rule of the multivariate Gaussian , we get

g n | w n ∼ N α n w n -Γ T f n -∆ T
υ nk w nk -f T n γ k -a T n δ k -g T n b k 2 | w n ; Θ = - 1 2 N n=1 (K + J) ln(2π) + K k=1 ln υ -1 nk + E g T n g n | w n ; Θ + E K k=1 υ nk w nk -f T n γ k -a T n δ k 2 + b T k g n g T n b k - 2 w nk -f T n γ k -a T n δ k g T n b k | w n ; Θ = - 1 2 N n=1 (K + J) ln(2π) + K k=1 ln υ -1 nk + E g T n g n | w n ; Θ + K k=1 υ nk w nk -f T n γ k -a T n δ k 2 + b T k Rn b k - 2 w nk -f T n γ k -a T n δ k gT n b k = - 1 2 N (K + J) ln(2π) + N n=1 K k=1 ln υ -1 nk + N n=1 E g T n g n | w n ; Θ + K k=1 w k -F γ k -Aδ k 2 W k + b T k N n=1 υ nk Rn b k - 2 Gb k T W k (w k -F γ k -Aδ k ) ,
where the rows of the matrix G are composed by gT n 's.

C.2 The maximization (M) step

The maximization step maximizes the conditional expectation of the complete loglikelihood with respect to Θ, subject to the upper triangular constraint on matrix B. However, for all k, the parameters γ k and δ k are not concerned by the constraint. Denoting β T k = (γ T k , δ T k ) and X = [F , A], the first order conditions of the maximization yield

∇ β k E[l(Θ; W, G) | W, Θ ] = 0 ⇔ ∇ β k w k -Xβ k 2 W k -2 Gb k T W k w k -Xβ k = 0 ⇔ XT W k w k -Xβ k -XT W k Gb k = 0 ⇔ XT W k Xβ k = XT W k w k -Gb k ⇔ β k = XT W k X -1 XT W k w k -Gb k .
If a response is drawn from a Gaussian law y k ∼ N N Xβ k , σ 2 k I N , the residual variance σ 2 k must be estimated. Besides,

∇ σ 2 k E[l(Θ; W, G) | W, Θ ] = 0 ⇔ ∇ σ 2 k N ln σ 2 k + 1 σ 2 k w k -Xβ k 2 + b T k N n=1 Rn b k -2 Gb k T w k -Xβ k = 0 ⇔ N - 1 σ 2 k w k -Xβ k 2 + b T k N n=1 Rn b k -2 Gb k T w k -Xβ k = 0 ⇔ σ 2 k = 1 N w k -Xβ k 2 + b T k N n=1 Rn b k -2 Gb k T w k -Xβ k .

C.3 The algorithm

As a result of the aforementioned developments, we shall use Algorithm 3 to estimate the parameters of the factor model.

Algorithm 4 The F-SCGLR algorithm 1: while not convergence do 2:

Compute the components through the PING algorithm

  vectors and W ∈ R N ×N be a symmetric positive definite matrix. The Euclidean scalar product between a and b with respect to metric W is given by a, b W = a T W b. Likewise, cos W (a, b) = a, b W a W b W denotes the cosine of the angle between a and b with respect to metric W . • If a and b are centred and W = 1 N I N , the cosine defines the linear correlation coefficient, denoted ρ. • A = [a 1 , . . . , a P ] ∈ R N ×P and B = [b 1 , . . . , b Q ] ∈ R N ×Q being matrices. The space spanned by their column-vectors is denoted span[A, B].

  K} be the set of parameters. The marginal loglikelihood of the model is obtained by integrating over latent variables g n l(Θ; Y ) = N n=1 ln (L(y n ; Θ))

1.

  Estimate the residual covariance matrix Σ := B T B. 2. Calculate the corresponding residual correlation matrix C where C ij := Σ ij / Σ ii Σ jj . 3. Calculate the associated dissimilarity matrix D where D 2 ij := 2 1 -C 2 ij . The squared residual correlation is used in order to consider as close two responses highly correlated, be it positively or negatively. 4. Perform Multidimensional Scaling (MDS, Cox and Cox, 2008) on the matrix D to obtain a Euclidean representation of the responses (i.e. a set of coordinates in a Euclidean space) with respect to this distance structure. We use the function cmdscale of the stats package (R Core Team, 2021).

Fig. 1 :

 1 Fig. 1: Residual correlation matrices for different values of σ 2 B

Fig. 2 :

 2 Fig. 2: Correlation scatterplot of plane (1,2) for the two themes with s = 0.3 and l = 4 obtained by the F-SCGLR algorithm. The red arrows represent the bundles X 1 and X 2 which explain the first theme. The blue ones represent the bundles X 3 and X 4 which explain the second theme. The percentage of inertia captured by each component is given in parentheses.

Fig. 3 :

 3 Fig. 3: Residual correlation matrices for different values of K

Fig. 4 :

 4 Fig. 4: Correlation plots of F-SCGLR plane (1,2), (2,3) and (1,3) of the second theme (farming intensity). The black arrows represent the covariates while the red ones are the linear predictors of the responses. The plot displays only variables having a cosine greater than 0.75 with the plane. The percentage of inertia captured by each component is given in parentheses.

Fig. 5 :

 5 Fig. 5: The residual correlation matrix alongside the Euclidean representation of the responses output by the MDS.

  first calculate the expectation of the complete log-likelihood conditional on the data WE[l (Θ; W, G) | W; Θ ] = N n=1 ln (L (w n | g n ; Θ) L (g n ; Θ)) L (g n | w n ; Θ ) dg n .

  a n , I J -α n B T ,where α n = B(B T B + Υ -1 n ) -1 .The moments of the random variable g n | w n are given bygn := E [g n | w n ; Θ] = α n w n -Γ T f n -∆ T a n and Rn := E g n g T n | w n ; Θ = V [g n | w n ; Θ] + E [g n | w n ; Θ] E [g n | w n ; Θ] T = I J -α n B T + gn gT n .Finally, we have the explicit form of the expectation of the complete log-likelihoodE[l (Θ; W, G) | W, Θ ]

W

  Now, we need to estimate the vector b k under the upper triangular constraint. For each k = 1, . . . , J, let b T k = (b T 1:k,k , 0 T ) be the regression parameters, where b T 1:k,k = (b 1k , . . . , b kk ) is a vector of length k to be estimated and 0 is a null vector of length (J -k) a priori fixed. In this case, we define ( Rn ) 1:k,1:k as the sub-matrix of size k × k of Rn and G1:k as the matrix composed by the first k columns of G. The maximization yields∇ b 1:k,k E[l(Θ; W, G) | W, Θ ] = 0 ⇔ ∇ b 1:k,k b T k w k -Xβ k .Likewise, for k = J + 1, . . . , K, b k is given by W k w k -Xβ k .

Table 1 :

 1 Mean values of RI, ARI and square correlation over a hundred samples with σ 2 B = 0.1, s ∈ {0.1, 0.3, 0.5} and l ∈ {1, 2, 3, 4, 7, 10}. (ξ 1 , .) ρ 2 (ξ 2 , .) ρ 2 (ξ 3 , .) ρ 2 (ξ 4 , .)

	s ρ 2 0.1 l RI ARI 1 0.926 0.839 0.938 2 0.927 0.839 0.980 3 0.920 0.816 0.981 4 0.928 0.837 0.979	0.906 0.959 0.962 0.966	0.759 0.810 0.805 0.816	0.838 0.935 0.942 0.945
		7	0.926 0.830	0.966	0.954	0.798	0.946
		10 0.927 0.835	0.967	0.955	0.792	0.945
		1	0.944 0.876	0.973	0.936	0.735	0.753
		2	0.944 0.877	0.993	0.972	0.934	0.950
	0.3	3 4	0.946 0.881 0.947 0.882	0.987 0.985	0.974 0.974	0.938 0.927	0.965 0.962
		7	0.943 0.875	0.984	0.974	0.911	0.964
		10 0.945 0.878	0.984	0.974	0.911	0.964
		1	0.942 0.871	0.974	0.937	0.697	0.659
		2	0.944 0.875	0.994	0.972	0.943	0.946
	0.5	3 4	0.947 0.882 0.948 0.884	0.988 0.986	0.975 0.975	0.948 0.946	0.961 0.967
		7	0.945 0.879	0.985	0.975	0.917	0.975
		10 0.944 0.877	0.985	0.975	0.911	0.969

[START_REF] Hui | Variational approximations for generalized linear latent variable models[END_REF]

, a Laplace approximation (LA,

[START_REF] Niku | Generalized linear latent variable models for multivariate count and biomass data in ecology[END_REF]

Table 2 :

 2 Mean values of RI, ARI, Procrustes error of the latent variables G and their loadings B and computation time over a hundred samples with N ∈ {100, 200, 300} and K ∈ {10, 30, 50}.

	N		100			200			300
	K	10	30	50	10	30	50	10	30	50
					F-SCGLR
	RI	0.845 0.869 0.895 0.922 0.951 0.980 0.950 0.976 0.986
	ARI	0.679 0.740 0.792 0.839 0.903 0.960 0.898 0.952 0.973
	Proc B	0.156 0.165 0.173 0.098 0.082 0.084 0.067 0.055 0.054
	Proc G	0.429 0.216 0.156 0.429 0.202 0.138 0.422 0.196 0.133
	Time	0.532 1.721 4.812 0.834 3.187 7.982 1.221 3.960 11.11
					gllvm-EVA
	RI	0.687 0.775 0.785 0.721 0.792 0.762 0.726 0.809 0.763
	ARI	0.386 0.550 0.571 0.456 0.583 0.518 0.466 0.617 0.521
	Proc B	0.982 0.670 0.615 1.021 0.713 0.664 0.948 0.635 0.621
	Proc G	0.578 0.300 0.223 0.544 0.256 0.182 0.518 0.231 0.166
	Time	1.436 5.182 11.02 4.167 18.20 44.76 8.995 36.86 88.99
					gllvm-VA
	RI	0.657 0.702 0.704 0.727 0.739 0.716 0.725 0.779 0.770
	ARI	0.337 0.402 0.403 0.468 0.476 0.426 0.455 0.556 0.534
	Proc B	0.737 0.633 0.635 0.677 0.618 0.668 0.624 0.614 0.618
	Proc G	0.709 0.296 0.244 0.697 0.245 0.180 0.702 0.224 0.165
	Time	1.384	8.80	20.18 3.645 22.78 47.58 8.665 38.34 79.65
					gllvm-shrinkage
	RI	0.606 0.700 0.765 0.648 0.735 0.698	0.65	0.697 0.685
	ARI	0.238 0.405 0.532 0.331 0.475 0.393 0.327 0.398 0.369
	Proc B	1.028 0.826 0.682 1.055 0.724 0.694 1.051 0.710 0.676
	Proc G	0.622 0.311 0.209 0.638 0.249 0.170 0.615 0.241 0.145
	Time	1.237 3.584 7.493 3.640 21.97 49.61 8.464 52.66 111.9
					gllvm-LA
	RI	0.598 0.699 0.689 0.652 0.720 0.734 0.710 0.749 0.761
	ARI	0.231 0.398 0.375 0.323 0.440 0.461 0.433 0.496 0.516
	Proc B	1.068 1.291 1.326 1.224 0.851 0.755 1.165 0.688 0.632
	Proc G	0.650 0.352 0.292 0.557 0.257 0.186 0.501 0.230 0.168
	Time	32.88 69.09 127.0 46.70 124.8 262.8 53.95 182.5 386.2

Declarations

• There is no potential conflict of interest.

• This research was supported by the GAMBAS project funded by the French Agence Nationale de la Recherche (ANR-18-CE02-0025). • The dataset is available following the link https://doi.org/10.15454/AJZUQN.

• The package FactorSCGLR and the simulation codes are available at https://github.com/julien-gibaud/FactorSCGLR.

• All the authors are contributed equally to this work.

Appendix A Analytical expression of the SCGLR-specific criterion and its derivative

The specific criteria which SCGLR needs to maximize for computing the (h + 1)-th loading-vector write

To facilitate the computation of the loading-vector, we give below an analytical expression of each sub-criterion and its gradient.

A.1 The structural relevance measure

In practice, we take either the variance component or the variable power inertia (VPI). In the first case, the SR and its gradient are easily given by

The explicit expression of VPI is

To calculate the gradient we use the classical rules of derivation

A.2 The goodness of fit measure

We aim at expressing ψ A h (u) as a function of quadratic forms. To achieve that, we decompose the projection on the regression space as follows

Algorithm 3 The EM algorithm applied to factor models with GLM 1: while not convergence do 2:

Expectation step 3:

end for 8:

Maximization step

9:

for k = 1, . . . , K do 10:

end for 12:

if Gaussian then 13:

Appendix D The overall F-SCGLR algorithm

Algorithm 4 consists in alternating the following steps: (i) Given the current set of parameters, calculate all the components of all the themes iteratively through the PING algorithm. (ii) Given the current components, calculate the working variables of the linearized model and their variance matrix. (iii) Given the working variables, estimate the factor model parameters through the EM algorithm.

Appendix E Identification of the true model

Appendix F Additional simulation studies

The Tables summing up the results for σ 2 B = 0.2 and σ 2 B = 0.3 are presented in Table F2 and Table F3.

Compute the model parameter through the EM algorithm 10:

Increment 12: t ← t + 1 13: end while Table E1: Mean values of BIC over a hundred samples for (H 1 , H 2 ) ∈ {1, 2, 3, 4} 2 and J ranging from 0 to 5. The lowest values are in bold font.