Charlotte Baey

P.-H Cournède

Estelle Kuhn
email: estelle.kuhn@inra.fr

Paul-Henry Cournède

Asymptotic distribution of likelihood ratio test statistics for variance components in nonlinear mixed effects models

Keywords: Chi-bar-square distribution, inference under constraints, hypothesis testing, likelihood ratio test, nonlinear mixed effects models, variance components

Introduction

Mixed effects models have been extensively used in population models in order to account for heterogeneity in populations and to describe the intra and inter-individual variability (see [START_REF] Pinheiro | Mixed-Effects Models in S and S-PLUS[END_REF]; [START_REF] Davidian | Nonlinear models for repeated measurement data: An overview and update[END_REF]; [START_REF] Lavielle | Mixed effects models for the population approach: models, tasks, methods and tools[END_REF]). In mixed effects models, parameters are of two types: on one side, fixed effects that are common to all the individuals of the population; on the other side, random effects that vary from one individual to the other. The last ones are also called individual parameters.

From a modelling point of view, a key question when adjusting a population model to a dataset is to identify the fixed and random effects of the model. Therefore, being able to compare a simpler model with less random effects to a larger one will help to choose a parsimonious model. From a statistical point of view, it can be rephrased as a test on the nullity of the variances of a subset of all the random effects (Silvapulle and Sen (2011), p14).

Several authors have been interested in likelihood ratio tests (LRT) in this context. We recall here some of the main existing results. Let us denote by Θ the parameter space, by Θ 0 the subset of Θ corresponding to the null hypothesis and by Θ 1 the subset corresponding to the alternative hypothesis, with Θ 0 ⊂ Θ 1 ⊂ Θ. [START_REF] Chernoff | On the Distribution of the Likelihood Ratio[END_REF], assuming that Θ is open, treated the case where the true value of the parameter lies on the boundary of Θ 0 and Θ 1 , which is assumed to be a proper set of Θ, i.e. strictly contained in Θ. He provided a representation of the asymptotic distribution of the LRT and proved that it is asymptotically equivalent to testing the mean of a multivariate Gaussian distribution based on one single observation. A few years later, [START_REF] Chant | On Asymptotic Tests of Composite Hypotheses in Nonstandard Conditions[END_REF] generalized these results by considering the case where the true value lies in a subset of Θ which may not be a proper subset. [START_REF] Shapiro | Asymptotic distribution of test statistics in the analysis of moment structures under inequality constraints[END_REF] studied the asymptotic distribution of a larger class of tests when the true value is on the boundary of Θ 0 and an interior point of Θ.

He established that the asymptotic distribution is a mixture of chi-square distributions. Simultaneously, [START_REF] Self | Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions[END_REF] obtained similar results in the case where the true value is on the boundary of Θ. They established a general convergence result and they derived the expression of the asymptotic distribution only in specific cases, assuming in particular that the parameter space Θ is equal to the product of a finite number of either closed, half-open or open intervals of R. They considered several specific cases for applications.

They identified in particular the limiting distribution of the LRT for testing that the variance of one single random effect is equal to zero as a mixture 1 2 δ 0 + 1 2 χ 2 1 , where δ 0 is the Dirac distribution at zero and χ 2 1 is the chi-square 3 distribution with one degree of freedom (Self and Liang, 1987, case 5). They also identified the limiting distribution of the LRT for testing simultaneously that the variance of one single random effect is equal to zero and that its mean is equal to a constant as a mixture 1 2 χ 2 1 + 1 2 χ 2 2 (Self and Liang, 1987, case 6). Nevertheless, their results do not allow to identify the asymptotic distribution of the LRT statistics in all cases, for example testing that one variance is equal to zero when considering two correlated random effects. Building upon those works, several authors have addressed the issue of variance components testing in mixed effects models. In the context of linear mixed effects models, Stram andLee (1994, 1995) considered the likelihood ratio test procedure and identified the limiting distribution of the LRT statistics in some particular cases. They suggested also that the limiting distribution might in fact be influenced by the presence of correlations between random effects but they do not investigate this issue. Some authors have also proposed finite sample test procedures for variance components testing in linear mixed models, either by deriving the finite sample distribution of test statistics, or using bootstrap and permutation tests. For example, the finite sample size distribution of the likelihood and restricted likelihood ratio test statistics was studied by [START_REF] Crainiceanu | Likelihood ratio tests in linear mixed models with one variance component[END_REF] for linear mixed models with one single random effect, and [START_REF] Greven | Restricted Likelihood Ratio Testing for Zero Variance Components in Linear Mixed Models[END_REF] extended these results to linear mixed models with more than one random effect.

Other approaches were also inquired. Several years later, [START_REF] Qu | Linear score tests for variance components in linear mixed models and applications to genetic association studies[END_REF] proposed a procedure based on the score test for testing several variance components in linear mixed models, and [START_REF] Wood | A simple test for random effects in regression models[END_REF] studied the finite sample distribution of a test based on the restricted likelihood for testing that one variance is null in generalized linear mixed models. Also in the context of linear mixed models, [START_REF] Sinha | Bootstrap tests for variance components in generalized linear mixed models[END_REF] studied a bootstrap test based on the score test for testing several variance components in a generalized linear mixed model, while [START_REF] Fitzmaurice | A note on permutation tests for variance components in multilevel generalized linear mixed models[END_REF], [START_REF] Samuh | The Use of Permutation Tests for Variance Components in Linear Mixed Models[END_REF] and [START_REF] Drikvandi | Testing multiple variance components in linear mixed-effects models[END_REF] considered permutation tests for testing several variance components in the context of linear and generalized mixed effects models. [START_REF] Molenberghs | Likelihood Ratio, Score, and Wald Tests in a Constrained Parameter Space[END_REF] proposed a review of the existing results for testing variance components in mixed effects models, and studied in particular the equivalence between the LRT, the Score test and the Wald test, based on results by [START_REF] Silvapulle | A score test against one-sided alternatives[END_REF] or [START_REF] Stram | Variance components testing in the longitudinal mixed effects model[END_REF].

They also exhibited the common limiting distribution in some specific cases.

However, to the best of our knowledge, there exist no results identifying the limiting distribution of the LRT for general tests on variance components in mixed effects models.

In this paper, we consider the LRT in nonlinear mixed effects models to test that the variances of any subset of the random effects are equal to zero and identify its asymptotic distribution as a mixture of chi-square distributions. In Section 2, we present the framework of nonlinear mixed effects models. Section 3 is devoted to the description of the proposed test and its theoretical properties. Practical implementation guidelines are presented in Section 4. Experimental results illustrate the performances of the procedure through simulation studies and real datasets analysis in Section 5. The paper ends with some discussion in Section 6. The technical proofs are given in Appendix.

Nonlinear mixed effects model

Description

We consider the following nonlinear mixed effects model [START_REF] Davidian | Nonlinear models for repeated measurement data[END_REF] p98, [START_REF] Pinheiro | Mixed-Effects Models in S and S-PLUS[END_REF] p306, Lavielle (2014) p24):

y i = g(ϕ i , x i ) + ε i ,
where y i denotes the vector of observations of the i-th individual of size J, 1 ≤ i ≤ N , g a nonlinear function, ϕ i the vector of individual parameters of individual i, x i a vector of covariates, and ε i the random error term. The vectors of individual parameters (ϕ i ) 1≤i≤N are modeled as:

ϕ i = U i β + V i b i , 1 ≤ i ≤ N,
where β is the vector of fixed effects taking values in R b , b i is the centered vector of Gaussian random effects with covariance matrix Γ, U i and V i are covariates matrices of individual i of sizes p × b and p × p respectively, and Γ is a covariance matrix of size p × p. The random vectors (b i ) are assumed to be independent. The vectors (ε i ) 1≤i≤N are assumed to be independent and identically distributed centered Gaussian vectors with covariance matrix Σ, and the sequences (ε i ) and (b i ) are assumed mutually independent.

Let us denote by θ = (β, Γ, Σ) the vector of all the model parameters and by q its dimension. Thus, the parameter space is defined as Θ = R p ×S p + ×S J + , where S p + is the set of symmetric, positive semi-definite p × p matrices.

Examples

A special but very common case is the one where the function g is linear.

The model can be rewritten in the following usual form [START_REF] Pinheiro | Mixed-Effects Models in S and S-PLUS[END_REF], p58):

y i = X i β + Z i b i + ε i ,
where y i is the observation vector for individual i, X i and Z i are matrices of known covariates, β is the vector of fixed effects, b i is the vector of centered random effects for individual i, with b i ∼ N (0, Γ), and ε i is a random error term, with ε i ∼ N (0, Σ).

Another famous example of a nonlinear mixed effects model is the logistic growth model, which was studied for example by [START_REF] Pinheiro | Mixed-Effects Models in S and S-PLUS[END_REF] in their well known example of orange trees growth. In this model, a logistic curve is used to model the growth of each individual in the population as a nonlinear function of three individual parameters. Denoting by y ij the variable measured for individual i at age x j (e.g. the trunk circumference in the orange trees example), for each individual i these three parameters are:

the asymptotic value of the trunk circumference ϕ i1 , the age at which the individual reaches half its asymptotic value, ϕ i2 , and the growth scale ϕ i3 .

More precisely, we have, for 1 ≤ i ≤ N, 1 ≤ j ≤ J:

y ij = ϕ i1 1 + exp - x j -ϕ i2 ϕ i3 + ε ij , with ϕ i = β + b i , where β = (β 1 , β 2 , β 3 ) ∈ R 3 , b i ∼ N (0, Γ) and ε ij ∼ N (0, σ 2 ).

Variance components testing

Description of the testing procedure

Let r ∈ {1, . . . , p}. We consider general test hypotheses of the following form, to test the nullity of r variances and of the corresponding covariances:

H 0 : θ ∈ Θ 0 against H 1 : θ ∈ Θ, (1)
where Θ 0 ⊂ Θ. Up to permutations of rows and columns of the covariance matrix Γ, we can assume that we are testing the nullity of the last r variances.

We write Γ in blocks as follows:

Γ =   Γ 1 Γ t 12 Γ 12 Γ 2   , with Γ 1 a (p -r) × (p -r) matrix, Γ 2 a r × r matrix, Γ 12 a r × (p -r) matrix
and where A t denotes the transposition of matrix A, for any matrix A.

The spaces associated to the null and alternative hypotheses are then:

Θ 0 = {θ ∈ R q | β ∈ R p , Γ 1 ∈ S p-r + , Γ 12 = 0, Γ 2 = 0, Σ ∈ S J + } Θ = {θ ∈ R q | β ∈ R p , Γ ∈ S p + , Σ ∈ S J + }.
We emphasize that the parameter space Θ is not open, and that the tested parameter values are on the boundary of Θ.

We recall below the likelihood ratio test procedure. Let us denote by y N 1 the joint vector of a N -sample (y 1 , . . . , y N ), and by L(y N 1 ; θ) the joint likelihood. We then define the likelihood ratio test statistics by:

LRT N = -2 log sup θ∈Θ 0 L(y N 1 ; θ) sup θ∈Θ L(y N 1 ; θ) . (2)
For a nominal level 0 < α < 1, the rejection region R α is defined by

R α = {LRT N ≥ q α }
, where q α is the (1 -α) quantile of the distribution of LRT N under the null hypothesis. However in practice the finite sample distribution of LRT N is generally intractable in the case of nonlinear mixed effects models. Therefore we focus on its asymptotic distribution.

Asymptotic property of the likelihood ratio test

Let us denote by θ * the true value of the parameters. We assume that the following condition (C1) is satisfied:

(i) The value θ * is in Θ 0 , i.e. θ * is of the form θ * = (β * , Γ * , Σ * ), with Γ * =   Γ * 1 0 (p-r)×r 0 r×(p-r) 0 r×(r)   where Γ * 1 is a (p -r) × (p -r) matrix.
(ii) The matrices Γ * 1 and Σ * are positive definite. In particular, we assume that the variances that are not being tested are strictly positive.

To derive the asymptotic distribution of the likelihood ratio test statistics under the null hypothesis, we need to ensure the consistency of the maximum likelihood estimate (MLE). Therefore, we assume that the following general condition (C2) is fulfilled [START_REF] Silvapulle | Constrained statistical inference: Order, inequality, and shape constraints[END_REF]: Remark 1. Note that the consistency and asymptotic normality of the MLE models in the specific context of nonlinear mixed effects have been studied in [START_REF] Nie | Strong consistency of the maximum likelihood estimator in generalized linear and nonlinear mixed-effects models[END_REF][START_REF] Nie | Convergence rate of MLE in generalized linear and nonlinear mixed-effects models: theory and applications[END_REF]. He exhibited specific assumptions that ensure these theoretical results. However, these assumptions might be difficult to verify in practice.

(i) The function L is injective in θ (to
Before stating the expression of the asymptotic distribution of the likelihood ratio test statistics, we recall the definition of the chi-bar-square distribution (for more details, see [START_REF] Shapiro | Asymptotic distribution of test statistics in the analysis of moment structures under inequality constraints[END_REF]; [START_REF] Silvapulle | Constrained statistical inference: Order, inequality, and shape constraints[END_REF]).

Definition 1. Let C be a closed convex cone of R q , V a positive definite matrix of size q × q and Z ∼ N (0, V ). The distribution of the random variable defined by

χ2 (V, C) = Z t V -1 Z -min θ∈C (Z -θ) t V -1 (Z -θ)
is called a chi-bar-square distribution. It is equal to a mixture of chi-square distributions with different degrees of freedom as follows:

∀t ≥ 0 P ( χ2 (V, C) ≤ t) = q i=0 w i (q, V, C)P (χ 2 i ≤ t),
where the weights (w i (q, V, C)) 0≤i≤q are non-negative numbers summing up to one, and where χ 2 i is a random variable following the chi-square distribution with i degrees of freedom, with the convention that χ 2 0 ≡ 0.

We can now establish the asymptotic property of the likelihood ratio test statistics.

Theorem 1. Assume that conditions (C1) and (C2) are fulfilled. Consider the test defined in (1). We denote by I * the Fisher information matrix evaluated at the true value θ * ∈ Θ 0 . Then:

LRT N ---→ N →∞ χ2 (I -1 * , T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ ), (3)
where T (Θ, θ) is the tangent cone to Θ at θ, and S ⊥ is the orthogonal complement of S, for any subset S of R q .

Proof. First we apply Theorem 3 of [START_REF] Self | Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions[END_REF] to derive a general expression of the asymptotic distribution of the likelihood ratio test statistics.

Indeed under conditions (C1) and (C2), the conditions stated in Section 1 of [START_REF] Self | Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions[END_REF] are fulfilled. In particular, condition (C2) ensures the consistency of the maximum likelihood estimates on Θ 0 and Θ. Moreover we prove that both Θ and Θ 0 can be approximated by cones, on which the consistency of the MLE also holds. We recall that the sets Θ 0 and Θ are defined by:

Θ 0 = {θ ∈ R q | β ∈ R p , Γ 1 ∈ S p-r + , Γ 12 = 0, Γ 2 = 0, Σ ∈ S J + }, Θ = {θ ∈ R q | β ∈ R p , Γ ∈ S p + , Σ ∈ S J + },
where the constraints are that Γ 1 , of size (p -r) × (p -r) and Γ, of size p × p, are positive semi-definite. Following Sylvester's criterion (see [START_REF] Gilbert | Positive definite matrices and Sylvester's criterion[END_REF]), a matrix is positive semi-definite if and only if all its leading principal minors are positive. Therefore both Θ 0 and Θ can be written as a set of polynomial equalities and inequalities corresponding to the different determinants. In particular, it means that both Θ 0 and Θ are semi-algebraic sets. They are therefore Chernoff-regular, i.e. they admit approximating cones at every point, and in particular at every point θ * ∈ Θ 0 [START_REF] Drton | Likelihood ratio tests and singularities[END_REF].

Therefore by applying Theorem 3 of [START_REF] Self | Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions[END_REF], we get that the asymptotic distribution of the likelihood ratio test is the same as the one of the likelihood ratio test from the null hypothesis "θ ∈ A(Θ 0 , θ * )" versus the alternative one "θ ∈ A(Θ, θ * )", where A(Θ 0 , θ * ) and A(Θ, θ * ) are respectively the approximation cones of Θ 0 and Θ at θ * :

LRT N → inf θ∈A(Θ 0 ,θ * )\θ * Z -θ 2 I * - inf θ∈A(Θ,θ * )\θ * Z -θ 2 I * , ( 4
)
where Z is a random variable with multivariate Gaussian distribution with mean θ and covariance matrix I -1 * . We use now the link between approximating and tangent cones, due to [START_REF] Geyer | On the asymptotics of constrained m-estimation[END_REF]. More precisely, if we denote by T (Θ, θ * ) the tangent cone of Θ at θ * and by A(Θ, θ * ) the approximating cone of Θ at θ * , then A(Θ, θ * ) = θ * + T (Θ, θ * ). Thus, we obtain:

LRT N → inf θ∈T (Θ 0 ,θ * ) Z -θ 2 I * -inf θ∈T (Θ,θ * ) Z -θ 2 I * . (5)
Finally, applying Theorem 3.7.1, page 84 of [START_REF] Silvapulle | Constrained statistical inference: Order, inequality, and shape constraints[END_REF], we derive that:

LRT N ---→ N →∞ χ2 (I -1 * , T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ ), (6)
which concludes the proof.

Analytical expression of the asymptotic distribution

We emphasize that the cone T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ always admits an analytical expression in our context. We consider in the following three common cases for the correlation structure of the p-dimensional random effects, namely independent, block-correlated and full correlated, corresponding respectively to a covariance matrix Γ diagonal, block-diagonal and full. The following proposition details the expressions of the cone T (Θ, θ * )∩T (Θ 0 , θ * ) ⊥ in these three cases.

Proposition 1.

(i) Assume that Θ = {θ ∈ R q | β ∈ R p , Γ ∈ S p + , Γ diagonal, Σ ∈ S J + }. Then T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ = {0} p × {0} p-r × R r + × {0} J(J+1)/2 . (ii) Assume that Θ = {θ ∈ R q | β ∈ R p , Γ ∈ S p + , Γ block-diagonal, Σ ∈ S J + } with Γ = diag(Γ 1 , . . . , Γ K )
, where, for k = 1, . . . , K, Γ k is a full covariance matrix of size r k ×r k , associated with the k-th sub-group of correlated random effects. Assume here that we want to test that the K-th block Γ K of variances is null and consider the corresponding subset Θ 0 . Then

T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ = {0} p × K-1 k=1 {0} r k (r k +1)/2 × S r K + × {0} J(J+1)/2 . (iii) Assume that Θ = {θ ∈ R q | β ∈ R p , Γ ∈ S p + , Γ full, Σ ∈ S J + }. Then T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ = {0} p × {0} (p-r)(p-r+1)/2 × R r(p-r) × S r + × {0} J(J+1)/2 .
The proof of Proposition 1 relies on technical elements from convex analysis (see Hiriart-Urruty and Lemarechal (1996); Hiriart-Urruty and Malick ( 2012)) and can be easily adapted to general covariance matrix structure using similar tools. It is postponed to the Appendix.

Moreover, thanks to the expressions of the cone T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ established in Proposition 1, we can deduce that several weights involved in the chi-bar-square distribution defined in (3) are equal to 0. The following corollary details this result for the two common cases of a diagonal or full covariance matrix Γ, described in Proposition 1.

Corollary 1.

(i) Assume that Θ = {θ ∈ R q | β ∈ R p , Γ ∈ S p + , Γ diagonal, Σ ∈ S J + }. Then the distribution of the random variable χ2 (I -1 * , T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥
) is a mixture of (r + 1) chi-square distributions with degrees of freedom between 0 and r.

(ii) Assume that Θ = {θ ∈ R q | β ∈ R p , Γ ∈ S p + , Γ full, Σ ∈ S J + }. Then the distribution of the random variable χ2 (I -1 * , T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥
) is a mixture of (r(r + 1)/2 + 1) chi-square distributions with degrees of freedom between r(p -r) and (r(p -r) + r(r + 1)/2).

The proof is postponed to the Appendix.

Theoretical result for non identically distributed variables

The theoretical results above can be extended in a natural way to non identically distributed observations, in particular to models involving covariates as defined in Section 2, or to models involving a specific number of observations J i per individual possibly different for each i.

However Theorem 3 of [START_REF] Self | Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions[END_REF], which leads to the asymptotic representation of the likelihood ratio test, can not be applied for non identically distributed random variables. In particular, the Fisher information matrix of N independent non identically distributed observations is no more defined as N times the Fisher information matrix of one observation as for a N -sample. Nevertheless, the asymptotic distribution of the LRT can be identified provided that suitable assumptions are fulfilled. We consider here the same conditions as those proposed in [START_REF] Silvapulle | Constrained statistical inference: Order, inequality, and shape constraints[END_REF] p156. Therefore we define the condition (C2)(iv) for all θ, there exists some positive definite matrix ν(θ) such that

N -1/2 ∇ θ log L(y N 1 ; θ) → N (0, ν(θ)) and N -1 ∇ 2 θ log L(y N 1 ; θ) → -ν(θ) a.s.,
and the function θ → ν(θ) is continuous. Moreover we define the condition (C3) the maximum likelihood estimates on Θ 0 and Θ are consistent. General assumptions ensuring this condition in mixed effects models are provided in [START_REF] Nie | Strong consistency of the maximum likelihood estimator in generalized linear and nonlinear mixed-effects models[END_REF]. We can then state the following result:

Theorem 2. Assume that conditions (C1), (C2)(i)(ii)(iv) and (C3) are fulfilled. Consider the test defined in (1). Then:

LRT N ---→ N →∞ χ2 (ν(θ * ) -1 , T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ ), (7)
where T (Θ, θ) is the tangent cone to Θ at θ, and S ⊥ is the orthogonal complement of S, for any subset S of R q .

Proof. We adapt the proof of our theoretical result stated in Subsection 3.2 to the case of non identically distributed variables using suited tools presented in [START_REF] Silvapulle | Constrained statistical inference: Order, inequality, and shape constraints[END_REF]. We detail all the steps below. Assuming (C2)(i)(ii)(iv), we establish the following fundamental quadratic approximation using standard tools of asymptotic statistics:

log L(y N 1 ; θ) = log L(y N 1 ; θ * ) + 1 2N ∇ θ log L(y N 1 ; θ * ) t ν(θ * ) -1 ∇ θ log L(y N 1 ; θ * ) - 1 2 (Z N -u) t ν(θ * )(Z N -u) + δ N (u)
where

Z N = N -1/2 ν(θ * ) -1 ∇ θ log L(y N 1 ; θ) and sup u ≤K |δ N (u)
| converges toward zero in probability for K > 0 given. Then, applying Lemma 4.2.3 in [START_REF] Silvapulle | Constrained statistical inference: Order, inequality, and shape constraints[END_REF], we obtain the √ N -consistency of the maximum likelihood estimates on Θ 0 and on Θ (see [START_REF] Andrews | Estimation when a parameter is on a boundary[END_REF] for the proof of this lemma in the non identically distributed setting). Following the lines of the proof of Theorem 1, we prove that both sets Θ and Θ 0 are Chernoff-regular and can be approximated by cones. By Corollary 4.7.5 in [START_REF] Silvapulle | Constrained statistical inference: Order, inequality, and shape constraints[END_REF], we get that the consistency of the maximum likelihood estimates also holds on the corresponding approximation cones. Then applying Proposition 4.8.2 in [START_REF] Silvapulle | Constrained statistical inference: Order, inequality, and shape constraints[END_REF], we get

LRT N → inf θ∈A(Θ 0 ,θ * )\θ * Z -θ 2 ν(θ * ) - inf θ∈A(Θ,θ * )\θ * Z -θ 2 ν(θ * ) ,
where Z is a random variable with multivariate Gaussian distribution with mean θ and covariance matrix ν(θ * ) -1 . Finally, we conclude the proof following the same lines as in the proof of Theorem 1.

Extension of the test and discussion

Testing simultaneously fixed effects and variance components.

Θ = R b × R p + × S J , Θ 0 = s 1 {β * j } × R b-s × {0} r × (R * + ) p-r × S J , and
T (Θ, θ * ) = R b × R p + × S J , T (Θ 0 , θ * ) = s 1 {0} × R b-s × {0} r × R p-r × S J ,
leading finally to

T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ = R s × {0} b-s × R r + × {0} p-r × {0} J(J+1)/2 .
Choosing s = r = 1, we recover in particular the asymptotic distribution 0.5χ 2 1 + 0.5χ 2 2 given in the case 6 of [START_REF] Self | Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions[END_REF].

Effects of nuisance parameters. Assumption (C1) ii) ensures that the nuisance variance parameters are not on the boundary of the parameter space.

Indeed those nuisance parameters have no action on the asymptotic distribution whatever their number and their type (mean or variance components).

Their corresponding components in the subsets Θ and Θ 0 are the same, and consequently also in the tangent cones T (Θ, θ * ) and T (Θ 0 , θ * ). This leads to null components in the cone T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ . This appears also in the asymptotic distributions obtained in cases 5 and 6 of [START_REF] Self | Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions[END_REF], which do not depend on the number of nuisance parameters. However as illustrated by the case 8 in [START_REF] Self | Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions[END_REF], the case of a test involving nuisance parameters on the boundary of the parameter space is more intricate and requires further studies.

Practical implementation

Computation of the likelihood ratio test statistic

The computation of the likelihood ratio test requires the computation of the maximum likelihood values under the null and alternative hypotheses, denoted respectively by θ0 and θ1 , as well as the values of the likelihood at these two points, L(y N 1 ; θ0 ) and L(y N 1 ; θ1 ).

However, in the context of nonlinear mixed effects models the likelihood is not available in a closed form, and we need to resort to stochastic variants of the Expectation-Maximization (EM) algorithm, such as the Stochastic Approximation EM algorithm for example [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF], to compute θ0 and θ1 . For the same reason, L(y N 1 ; θ0 ) and L(y N 1 ; θ1 ) cannot be computed explicitly, and should be approximated using appropriate methods such as numerical or stochastic integration.

Since the decision to reject the null hypothesis relies on the value of the test statistics, the approximation of L(y N 1 ; θ) must be computed precisely. Let us denote by L(y i ; θ) the marginal likelihood of the i-th individual, and by (y N 1 ; θ) the joint log-likelihood. Then:

(y N 1 ; θ) = log N i=1 L(y i ; θ) = N i=1 log R p f (y i | ϕ i ; θ)p(ϕ i ; θ)dϕ i ,
where f (• | ϕ i ; θ) is the conditional probability density function of y i given the random effect ϕ i , and p(•; θ) is the probability density function of the random effect ϕ i . This quantity can be approximated using classical methods for integral approximations. However in the case of high dimensional random effects, stochastic integration is preferred over numerical approximations, allowing a better approximation for comparable computation times.

In practice, each L(y i ; θ) can be approximated independently from the others using Monte Carlo methods, and this can be done in parallel to optimize the execution time. Computing the joint log-likelihood using the sum of the marginal log-likelihoods, instead of taking the logarithm of their product can also result in less numerical issues. To further reduce the variability of the LRT statistics estimate, we can compute directly the Monte Carlo esti-mate of LRT N based on the same sample of size M rather than using two estimates of (y N 1 ; θ0 ) and (y N 1 ; θ1 ) based on two different samples. Thus, let us consider:

LRT N,M = -2 N i=1 log M m=1 f (y i | φm i,0 ; θ0 ) M m=1 f (y i | φm i,1 ; θ1 )
,

where φm i,0 = β0 + Γ1/2 0 Z m i , φm i,1 = β1 + Γ1/2 1 Z m i and Z m i ∼ N (0, I p ).

Computation of chi-bar-square weights when Γ is diagonal

In general, the weights involved in the definition of the chi-bar-square distribution are not available in a tractable form. However, in the special case of a diagonal covariance matrix, the cone C involved in the chi-barsquare distribution is polyhedral of dimension r (see Proposition 1). Indeed, in this case the cone can be written as

C = {θ ∈ R q | Rθ ≥ 0}, with R = 0 r×(p+p-r) | I r | 0 r× J(J+1)

2

, a full-rank matrix of dimension r × q. For polyhedral cones of this type, [START_REF] Shapiro | Asymptotic distribution of test statistics in the analysis of moment structures under inequality constraints[END_REF] provided the exact weights expressions for 1 ≤ r ≤ 3. The case r = 2 is also treated by [START_REF] Self | Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions[END_REF]. Following the notation of [START_REF] Shapiro | Asymptotic distribution of test statistics in the analysis of moment structures under inequality constraints[END_REF], we denote by

ρ ij = v ij /(v ii v jj ) 1/2 and ρ ij.k = (ρ ij -ρ ik ρ jk )/((1 -ρ 2 ik )(1 -ρ 2 jk )) 1/2
, respectively the correlation coefficient, and the partial correlation coefficient associated with the covariance matrix RI -1 * R t , where v ij stands for the element in row i and column j of matrix RI -1 * R t . Using Proposition 3.6.1 of Silvapulle and Sen (2011), we have w i (q, I -1 * , C) = w i (r, RI -1 * R t , R r + ) and denoting by w i,r = w i (r, RI -1 * R t , R r + ), we have the following expressions: (i) for r = 1, we have w 0,1 = w 1,1 = 1/2, (ii) for r = 2, we have w 0,2 = 1/2 π -1 cos -1 (ρ 12 ), w 1,2 = 1/2, and w 2,2 = 1/2 -1/2 π -1 cos -1 (ρ 12 ), (iii) for r = 3, we have: w 1,3 = 1/2 -w 3,3 , andw 0,3 = 1/2w 2,3 .

w 3,3 = (4π) -1 (2π -cos -1 (ρ 12 ) -cos -1 (ρ 13 ) -cos -1 (ρ 23 )), w 2,3 = (4π) -1 (3π - cos -1 (ρ 12.3 ) -cos -1 (ρ 13.2 ) -cos -1 (ρ 23.1 )),
For r > 3, and in more general settings, e.g. when C is not a polyhedral cone, one has to either approximate the weights through numerical integration or Monte Carlo simulations, or to directly compute the tail probability of the chi-bar-square distribution (see [START_REF] Silvapulle | Constrained statistical inference: Order, inequality, and shape constraints[END_REF], page 78).

Experiments

Simulation study

Simulation settings. We consider a linear mixed effects model and the logistic mixed effects model described in Section 2.2 with two and three random effects to explore different settings. We denote by θ * the true parameter value used to generate the data under H 0 .

Let us denote by M 1 the linear model with three individual parameters where we set g(ϕ i , x j ) = ϕ i1 + ϕ i2 x j + ϕ i3 x 2 j and ϕ i = β * + b i . We choose β * = (0, 7, 2) t , γ * 1 = 1.3, γ * 2 = 1 and γ * 12 = 1.04, corresponding to a correlation coefficient of 0.8 between b i1 and b i2 . We consider the null hypothesis H 0 defined by γ * 3 = γ * 13 = γ * 23 = 0. In the sub-model with two individual parameters, we set g(ϕ i , x j ) = ϕ i1 + ϕ i2 x j , where ϕ i = β * + b i , β * = (0, 7) t and γ * 1 = 1.3. In this case we consider H 0 defined by γ * 2 = γ * 12 = 0. In each simulation settings, x j = j, J = 20 and σ * = 1.5.

Let us denote by M 2 the logistic model with three random effects where we set g(ϕ i , x j ) = ϕ i1 /(1 + exp (-(x j -ϕ i2 )/ϕ i3 )), where ϕ i = β * + b i . We set β * = (200, 500, 150) t , γ * 2 = 50, γ * 3 = 15 and γ * 23 = 375, corresponding to a correlation coefficient of 0.5 between ϕ i2 and ϕ i3 . We consider here the null using the SAEM algorithm implemented in the R package saemix [START_REF] Comets | SAEMIX, an R version of the SAEM algorithm[END_REF], for the nonlinear model. Others parts of the codes are also developed in R and are available on request.

Empirical level. For each test case, to evaluate the level of the test we generate K datasets D 0,1 , . . . , D 0,K under the null hypothesis, and we denote by θ0,k (resp. θ1,k ) the maximum likelihood estimates of θ * using dataset D 0,k under H 0 (resp. H 1 ). The likelihood ratio test statistics estimate is denoted by LRT k . Then, the empirical level of the test for a sample size K is equal

to αK = 1 K K k=1 1 LRT k >cα
, where c α is the (1 -α) quantile of the limiting distribution of the LRT statistics. In practice, c α is not always available in a closed form and may be estimated as mentioned in Section 4.

Note that in the linear mixed model case, the Fisher information matrix I * is known and is given by:

(I * ) i,j = ∂Xβ ∂θ i t θ=θ * (Ω * ) -1 ∂Xβ ∂θ j θ=θ * + 1 2 tr (Ω * ) -1 ∂Ω ∂θ i θ=θ * (Ω * ) -1 ∂Ω ∂θ j θ=θ * ,
where Ω * = ZΓ * Z t + (σ * ) 2 I J , θ i is the i-th element of vector θ, where tr(A) denotes the trace of A for any matrix A and where for a matrix A of size m × n, ∂A/∂x is the matrix of size m × n whose element (i, j) is given by

(∂A/∂x) i,j = ∂A ij /∂x.
We first analyze the finite sample size properties of the LRT statistics when performing the test in the linear model M 1 involving two or three random effects, with and without correlations between the random effects. We start by testing that the variance of one random effect is zero. We compute the empirical level as detailed above for nominal level α in {0.01, 0.05, 0.10} and for a sample size N varying in {100, 500, 800}. Results are presented q q q q q q q q q 0.025 Full covariance structure q q q q q q q q q 0.025 in Figures 1 and2. We observe that the empirical levels are closer to the nominal ones when N grows, for random effects involving two and three components. We also observe that the empirical levels are lower than the nominal ones, particularly when N is smaller. This is a known property of chi-bar-square distributions (see for example [START_REF] Fitzmaurice | A note on permutation tests for variance components in multilevel generalized linear mixed models[END_REF] and [START_REF] Drikvandi | Testing multiple variance components in linear mixed-effects models[END_REF]). However, the theoretical results developed in this paper are asymptotic ones. For small values of N , non-asymptotic test procedures such as permutation tests are more appropriate (see [START_REF] Drikvandi | Testing multiple variance components in linear mixed-effects models[END_REF]).

Let us now highlight that one can be led to false conclusions when performing the LRT in a model without taking into account the presence of q q q q q q q q q 0.025 Full covariance structure q q q q q q q q q 0.025 correlations between random effects. For example, let us consider Case 2, corresponding to a limiting distribution 0.5χ 2 1 + 0.5χ 2 2 . If we do not assume a correlation structure between the two random effects although there exists one, then the corresponding limiting distribution is 0.5χ 2 0 + 0.5χ 2 1 . We compute the corresponding empirical quantiles. Results are displayed in Table 1. We observe that the empirical levels in column 3 are too large, leading to possibly wrong conclusions. This emphasizes that the presence of correlations between the random effects in the model plays a crucial role when performing a test on variance components in mixed effects models.

We then evaluate the finite sample size properties of the LRT statistics when performing the test of two variances equal to zero in the linear model M 1 involving three independent random effects. Results are detailed in Figure 2. As previously, we observe that the empirical levels converge to the nominal ones when N grows. However the asymptotic effect seems to occur slower when testing that two variances are equal to zero. This may be explained by the fact that more parameters have to be estimated than in the case where one variance is tested equal to zero.

We also study the properties of the LRT statistics when testing that the two variances are equal to zero in the linear model M 1 involving two independent random effects, which amounts to testing the absence of random effects in the model. Results are presented in Table 2. We observe the convergence of the empirical levels towards the nominal ones as N increases.

Let us now focus on the nonlinear model M 2 . We analyze the finite sample size properties of the LRT statistics when performing the test in the 3. We observe that in the case of this nonlinear mixed effects model, the empirical levels are not as close to the nominal ones as in the linear case and are always lower than the nominal levels in our example. We study the asymptotic behaviour of the empirical levels and notice that they are not improved when N is increased to 1000 (results not presented). These numerical results might be explained by the numerical integrations which have to be performed to compute the LRT statistics in nonlinear mixed effects models, which is not the case in the linear setting. To illustrate this, we use the linear model M 1 where the likelihood can be computed exactly. We compare the exact value (given by the lmer function in R) with the approximation obtained by an importance sampling scheme (given by the saemix function in R), for K = 1000 datasets. It appears on this specific example that the approximation procedure of saemix has a tendency to under-estimate the likelihood value, especially under the alternative hypothesis. Consequently, the LRT statistics is in general lower than the exact value when computed by saemix, which implies lower empirical levels (see Table 4).

When we consider correlations between the random effects in model M 2 , the empirical levels are globally slightly lower than in the models without correlations. This may be explained again by the fact that the models with correlations involved more parameters to estimate. We also observe that the empirical levels in model M 2 with 3 random effects are lower than those in model M 2 involving only 2 random effects. The same argument can be retained in this case. q q q q q q q q q 0.025 Empirical power. We assess the empirical power of the procedure in both models M 1 and M 2 , when performing the test of one variance equal to zero in a model with two random effects. Results are given in Tables 5 and6 for the linear model, and in Table 7 for the nonlinear model. We observe that the empirical power rises rapidly as the true value of the variance being tested γ * 2 increases. As β * 2 = 7, it corresponds to a relative standard deviation of approximately 1.4%. In the model with correlated random effects, the empirical powers increase with the correlation coefficient between the random effects.

Misspecification of the random effects distribution.. We assess the effect of a misspecification of the random effects distribution on the empirical level.

More precisely, we generate K = 10000 datasets in the linear model M 1 with two independent random effects, following non-Gaussian distributions.

Then, we compute the maximum likelihood estimators under H 0 and H 1 and q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 2 the corresponding test statistics for each dataset, assuming that the random effects are Gaussian. We consider three types of distribution: a Student, a log-normal and a mixture of two Gaussian distributions. Hence, we cover three different departures from the Gaussian distribution: higher tails with the Student distribution, asymmetry with the log-normal distribution, and bimodality with the mixture of Gaussian distributions. Data are generated under the null hypothesis, with γ * 2 = 0, and ϕ i1 generated according to one of the three types of distribution mentioned above. For the Student distribution, we choose a degree of freedom of 4.89, leading to a standard deviation of 1.3. For the log-normal distribution, we choose the parameters such that the variable log ϕ i1 is centered and its standard deviation takes two different values, namely 0.25 and 1.5, corresponding respectively to a distribution which is slightly asymmetric with a mode around 0.95, and to a distribution which is highly asymmetric, with a mode around 0.10. For the mixture of Gaussian distributions, we consider two cases: (i) a symmetric one where the expectations of the Gaussian distributions are set respectively to -1 and 1, both standard deviations to 1, and both weights to 0.5, and (ii) an asymmetric one where the expectations are set respectively to -3 and 1, both standard deviations to 1, and the weights to 0.25 and 0.75. Results for random effects following the Student distribution t 4.89 , the log-normal distribution such that log(ϕ i1 ) ∼ N (0, 1.5 2 ) and the Gaussian mixture 0.5N (-1, 1)+0.5N (1, 1) are presented in Figure 5. Results for the two other distributions are very similar and therefore not presented here. We observe on these specific examples that the empirical levels are globally closed to the nominal ones. Moreover, except in the case of the highly asymmetric log-normal distribution, the empirical levels are closer to the nominal ones than the empirical levels obtained in the case where the random effects are Gaussian. This numerical study gives some insight on the effects of misspecification of the random effects distribution on the proposed test procedure. Yet, it can not lead to a general conclusion about the robustness of the proposed test procedure to misspecification of the random effects distribution. However let us mention that there exist diagnostic tools for the normality assumption of the random effects distribution in mixed effects models. Let us quote for example those proposed recently by [START_REF] Drikvandi | Diagnosing misspecification of the random-effects distribution in mixed models[END_REF] and Drikvandi (2017).

q q q q q q q q q 0.025 0.050 0.075 0.100 200 400 600 800

Sample size N Empirical level α q q q 0.01 0.05 0.1 (a) Student q q q q q q q q q 0.025 0.050 0.075 0.100 200 400 600 800

Sample size N Empirical level α q q q 0.01 0.05 0.1 (b) log-normal q q q q q q q q q 0.025 The model can be written in the following form :

y i = X i β + Z i b i + ε i X i =         1 x i t 1 x i t 1 1 x i t 2 x i t 2 . . . . . . . . . . . . 1 x i t J x i t J         , β =         α 0 α 1 β 0 β 1         , Z i =         1 t 1 1 t 2 . . . . . . 1 t J        
.

First, we assess if the random effects are correlated or not using a LRT test, and we obtain a test statistics of 0.3027. In this case, the limiting distribution is a chi-square with one degree of freedom. Thus, at the level of 5%, the null hypothesis is not rejected. Therefore in the sequel we consider a diagonal covariance matrix Γ = diag(γ 2 1 , γ 2 2 ). We first test H 0 : {γ 2 1 = 0, γ 2 2 = 0} against H 1 : {γ 2 1 ≥ 0, γ 2 2 ≥ 0}, i.e. we test for the absence of randomness in the model. We compute the likelihood ratio test statistics and compare it to the rejection threshold q d α associated with the limiting distribution d and corresponding to a level α. The test statistics is equal to 50.13, the limiting distribution is the mixture w 0 χ 2 0 + 0.5χ 2 1 + (1 -w 0 )χ 2 2 , with w 0 = 0.37, and the rejection threshold is q 0.37χ 2 0 +0.5χ 2 1 +0.63χ 2 2 0.05 = 3.61. Therefore, at the level of 5% we reject the null hypothesis that there is no random effects in the model.

Next, we want to test H 0 : {γ 2 1 ≥ 0, γ 2 2 = 0} against H 1 : {γ 2 1 ≥ 0, γ 2 2 ≥ 0}, i.e. we want to test if there is a random slope. In this case, the test statistics is equal to 0.5304, the rejection threshold is equal to q 0.5χ 2 0 +0.5χ 2 1 0.05 = 2.706 and the p-value is equal to 0.233. Therefore the null hypothesis is not rejected, which suggests that there is no randomness in the slope. These results are consistent with previous results obtained on this dataset, in particular in [START_REF] Drikvandi | Testing variance components in balanced linear growth curve models[END_REF].

0.139χ 2 0 + 0.5χ 2 1 + 0.361χ 2 2 . The likelihood ratio test statistics is equal to LRT = 3.119, and the rejection threshold is equal to q 0.139χ 2 0 +0.5χ 2 1 +0.361χ 2 2 0.05 = 4.682. The corresponding p-value is evaluated at 0.114 on this dataset. In other words, assuming a diagonal covariance matrix Γ, we do not reject the null hypothesis that both the inflection point and the growth rates are fixed effects and do not vary among individuals, at the asymptotic level of 5%.

Discussion

Several perspectives of this work are of great interest, both from a theoretical and a practical point of view. For a wide-range application of these results, adapted tools to compute precisely the test statistics in nonlinear mixed effects models have to be developed. Indeed the reliability of the test is linked to the precise evaluation of the test statistics. Moreover, the procedure opens very promising perspectives in applications with high dimensional random effects. For example, models of plant growth have raised expectations to help improving the understanding of gene by environment interactions by developing a predictive capacity that scales from genotype to phenotype [START_REF] Letort | Quantitative genetics and functional-structural plant growth models: sim-ulation of quantitative trait loci detection for model parameters and application to potential yield optimization[END_REF]. In plant ecophysiological models, one genotype should be represented by one unique set of parameters, and reversely, two different genotypes should potentially be characterized by two different sets of parameters [START_REF] Tardieu | Virtual plants: modelling as a tool for the genomics of tolerance to water deficit[END_REF]. Such models are often descriptive ones, involving mechanistic parameters. Therefore, considering these parameters as random effects is relevant in order to understand how they vary within a given population [START_REF] Baey | A nonlinear mixed effects model of plant growth and estimation via stochastic variants of the EM algorithm[END_REF]. However this leads to many computational issues, since the number of parameters is high. Therefore, the issue of reducing the model size in a mixed effects model has gained an increasing interest in plant growth modelling. Also more advanced works on the computational methods are still necessary for the computation of the Fisher information matrix. Besides, from a theoretical point of view, it would be interesting to establish results when considering nuisance parameters lying on the boundary of the parameter space, as well as results for the test power.

Finally, since the asymptotic regime is not always reached in practice, it would be of a great interest to develop a finite sample-size procedure using for example bootstrap methods or permutation tests in the spirit of the ones developed in the context of linear mixed effects models [START_REF] Fitzmaurice | A note on permutation tests for variance components in multilevel generalized linear mixed models[END_REF][START_REF] Drikvandi | Testing multiple variance components in linear mixed-effects models[END_REF].

Appendix

Proof of Proposition 1

The proof relies on technical elements from convex analysis and can be easily adapted to general covariance matrix structure using similar tools.

To calculate the tangent cones to Θ 0 and Θ at θ * , we can use general results from Hiriart-Urruty and [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms I: Fundamentals[END_REF] on the definition of tangent cones, and more recent results by Hiriart-Urruty and Malick (2012) on the tangent cone of the set of symmetric positive semi-definite matrices.

The proof is carried out in the cases where the covariance matrix Γ is diagonal and full. Similar tools can be used in other cases where a more sparse structure is assumed for Γ.

Similarly, we have:

T (S J + , Σ * ) = M ∈ S J | ∀u ∈ ker Σ * , M u, u ≥ 0 = M ∈ S J | ∀u ∈ {0}, M u, u ≥ 0 = S J .
In the end, we get:

T (Θ, θ * ) = R p × R (p-r)(p-r+1)/2 × R r(p-r) × S r + × S J .

In the case where Γ is diagonal, i.e. when the effects are assumed to be independent, we have:

T (Θ, θ * ) = R p × R p-r × R r + × S J .
Tangent cone of Θ 0 . We recall that Θ 0 is defined as: Θ 0 ={θ ∈ R q | β ∈ R p , Γ 1 ∈ S p-r + , Γ 12 = 0, Γ 2 = 0, Σ ∈ S J + } =R p × S p-r + × {0} r(p-r) × {0} r(r+1)/2 × S J + .

Using similar tools than for Θ, we get:

T (Θ 0 , θ * ) = R p × S p-r × {0} r(p-r) × {0} r(r+1)/2 × S J = R p × R (p-r)(p-r+1)/2 × {0} r(p-r)+r(r+1)/2 × R J(J+1)/2 .

Note that in the particular case where Γ is diagonal, i.e. when the effects are supposed to be independent, the parameter space Θ 0 and hence its tangent cone can be simplified, and we have:

T (Θ 0 , θ * ) = R p × R p-r × {0} r × R J(J+1)/2 .
Finally, we derive the expressions stated in (i) and (iii).

  ensure the identifiability of the model), (ii) The first three derivatives of the log-likelihood w.r.t. θ exist and are bounded by a function whose expectation exists, (iii) The Fisher information matrix is finite and positive definite.

  The same kind of theoretical results can be established when testing simultaneously fixed effects values and variance components. Indeed, the components of the null hypothesis subset Θ 0 corresponding to the fixed effects values tested are modified as follow: they are equal to the corresponding value of β * in place of R. Based on this new expression, we calculate the corresponding tangent cone to be able to derive finally the expression of the asymptotic distribution of the LRT. Let us detail the case of a mixed model with p independent random effects where we are interested in testing simultaneously that 0 < s ≤ b fixed effects are equal to the corresponding values of β * and that 0 < r ≤ p variances are equal to zero. We obtain:

Figure 1 :

 1 Figure 1: Empirical level and Monte-Carlo error in the linear model for a theoretical level α ∈ {0.01, 0.05, 0.10} when testing that one variance is null in a model with 2 random effects (Cases 1-2), evaluated on K = 10000 datasets of size N ∈ {100, 500, 800}

Figure 2 :

 2 Figure 2: Empirical level and Monte-Carlo error in the linear model for a theoretical level α ∈ {0.01, 0.05, 0.10} when testing that one variance is null in a model with 3 random effects (Cases 3-4), evaluated on K = 10000 datasets of size N ∈ {100, 500, 800}

Figure 3 :

 3 Figure 3: Empirical level and Monte-Carlo error in the linear model for a theoretical level α ∈ {0.01, 0.05, 0.10}, when testing that the two variances γ 2 1 and γ 2 2 are null (Case 5), evaluated on K = 10000 datasets of size N ∈ {100, 500, 800}

Table 4 :

 4 Comparison of the empirical level and Monte Carlo standard error in parenthesis computed using the exact value of the LRT statistics and an approximation of the LRT statistics using an importance sampling (IS) scheme, in the linear model, when testing that one variance is null, and evaluated on K = 1000 datasets of size N = 500 α = 0.01 α = 0.05 α = 0.10 Exact value of the LRT 0Empirical power and Monte Carlo standard error in parenthesis when testing that one variance is null in the linear model with two independent random effects, evaluated on K = 10000 datasets of size N . α = 0.01 α = 0.05 α = 0.10 γ 2 N = 100 N = 500 N = 100 N = 500 N = 100 N = 500 0.01 0

Figure 4 :

 4 Figure4: Comparison of the results based on the exact evaluation of the likelihood (comby the lmer function) and on its approximation using importance sampling (computed by the saemix function), evaluated on 1000 datasets. Figure4a: difference between the exact value of the likelihood and its approximation, under each hypotheses. Figure4b: LRT statistics based on the approximated values of the likelihood, as a function of its exact value.

Figure 5 :

 5 Figure 5: Empirical level and Monte-Carlo error in the linear model for a theoretical level α ∈ {0.01, 0.05, 0.10}, evaluated on K = 10000 datasets of size N ∈ {100, 500, 800}, when testing that γ 2 2 is null, under different misspecification of the random effects distribution: (a) ϕ i1 ∼ t 4.89 , (b) log(ϕ i1 ) ∼ N (0, 1.5 2 ) , (c) ϕ i1 ∼ 0.5N (-1, 1) + 0.5N (1, 1)

Table 1 :

 1 Empirical level and Monte Carlo standard errors in parenthesis in the linear model for a theoretical level α ∈ {0.01, 0.05, 0.10}, when testing that one variance is null in a model with two correlated random effects, using the theoretical quantiles of the

	limiting distribution 0.5χ 2 1 + 0.5χ 2 2 (column 2), and using the quantiles of the limiting
	distribution 0.5χ 2 0 + 0.5χ 2 1 obtained when considering a model with uncorrelated random
	effects (column 3), for N = 500
	α	α0.5χ 2 1 +0.5χ 2 2	α0.5χ 2 0 +0.5χ 2 1
	0.01	0.009	0.050
		(0.0009)	(0.0019)
	0.05	0.045	0.174
		(0.0021)	(0.0037)
	0.10	0.092	0.311
		(0.0029)	(0.0046)

Table 2 :

 2 Empirical level and Monte-Carlo standard errors in parenthesis in the linear

	model for a theoretical level α ∈ {0.01, 0.05, 0.10}, when testing that two variances are
	null in a model with two independent random effects (Case 6), evaluated on K = 10000
	datasets of size N ∈ {100, 200, 500}
	α	N = 100 N = 500
	0.01	0.008	0.009
		(0.0009) (0.0009)
	0.05	0.041	0.046
		(0.0020) (0.0021)
	0.10	0.082	0.088
		(0.0027) (0.0028)
	model involving two and three random effects, with and without correlations
	between random effects. Results are detailed in Table

Table 3 :

 3 Empirical level in the nonlinear model and Monte Carlo standard errors in

	parenthesis for a theoretical level α ∈ {0.01, 0.05, 0.10} when testing that one variance is
	null (Cases 1-4), evaluated on K = 1000 datasets of size N = 500
		2 random effects	3 random effects
	α	Γ diagonal (Case 1) Γ full (Case 2) Γ diagonal (Case 3) Γ full (Case 4)
	0.01	0.003	0.007	0	0.003
		(0.0017)	(0.0020)	(0)	(0.0033)
	0.05	0.038	0.033	0.040	0.033
		(0.0060)	(0.0048)	(0.0113)	(0.0104)
	0.10	0.082	0.073	0.077	0.073
		(0.0087)	(0.0068)	(0.0154)	(0.0151)

Table 6 :

 6 Empirical power and Monte-Carlo standard error in parenthesis when testing that one variance is null in the linear model with two correlated random effects, evaluated on K = 10000 datasets of size N .

			α = 0.01	α = 0.05	α = 0.10
	γ 2	ρ 12 N = 100 N = 500 N = 100 N = 500 N = 100 N = 500
	0.01 0.1	0.0135	0.0371	0.0624	0.1250	0.1193	0.2051
			(0.0005) (0.0019) (0.0021) (0.0033) (0.0035) (0.0040)
		0.5	0.0385	0.2674	0.1280	0.4957	0.2152	0.6236
			(0.0019) (0.0044) (0.0033) (0.0050) (0.0041) (0.0048)
		0.75 0.0737	0.5617	0.2059	0.7827	0.3103	0.8640
			(0.0026) (0.0049) (0.0040) (0.0041) (0.0046) (0.0034)
		0.9	0.1010	0.7265	0.2547	0.8857	0.3705	0.9365
			(0.0030) (0.0044) (0.0043) (0.0032) (0.0048) (0.0024)
	0.02 0.1	0.0694	0.4182	0.1869	0.6585	0.2825	0.7699
			(0.0025) (0.0049) (0.0039) (0.0047) (0.0045) (0.0042)
		0.5	0.2556	0.9652	0.4728	0.9912	0.5993	0.9952
			(0.0044) (0.0018) (0.0050) (0.0009) (0.0049) (0.0007)
		0.75 0.4825	0.9991	0.7008	0.9999	0.8022	1
			(0.0050) (0.0003) (0.0046) (0.0001) (0.0040)	(0)
		0.9	0.6025	1	0.8058	1	0.8756	1
			(0.0049)	(0)	(0.0039)	(0)	(0.0033)	(0)
	0.05 0.1	0.9644	1	0.9925	1	0.9964	1
			(0.0018)	(0)	(0.0009)	(0)	(0.0006)	(0)

Table 7 :

 7 Empirical power Monte-Carlo standard error in parenthesis when testing that one variance is null in the nonlinear model with two independent random effects, evaluated on K = 1000 datasets of size N .

			α = 0.01	α = 0.05	α = 0.10
	γ 2 ρ N = 100 N = 500 N = 100 N = 500 N = 100 N = 500
	2	0	0.30	0.98	0.46	0.995	0.55	0.999
			(0.0145) (0.0044) (0.0158) (0.0022) (0.0157) (0.001)
	2 0.2	0.182	0.945	0.310	0.985	0.400	0.991
			(0.0122) (0.0072) (0.0146) (0.0038) (0.0155) (0.0030)
	2 0.5	0.448	1	0.639	1	0.717	1
			(0.0157)	(0)	(0.0152)	(0)	(0.0142)	(0)
	2 0.8	0.880	1	0.953	1	0.970	1
			(0.0103)	(0)	(0.0067)	(0)	(0.0054)	(0)

hypothesis H 0 defined by γ * 1 = γ * 12 = γ * 13 = 0. In the sub-model with two random effects, we set β * = (200, 500) t and γ * 2 = 50. In this case we consider H 0 defined by γ * 1 = γ * 12 = 0. The vector of observation times (x j ) is the same for all the individuals, and is defined as a vector of 20 equally spaced values between 50 and 1000, plus 5 equally spaced values between 1100 and 1500. In each simulation settings, σ * = 10. When only two random effects are considered in the model, β 3 is fixed to 150 and not estimated by the algorithm.

For the two models M 1 and M 2 , we consider several test cases, according to the number of variances being tested, to the total number of random effects in the model, and to the presence of correlations between the random effects. More precisely, we consider the following cases: (i) testing that one variance is null in a model with two independent (Case 1) or non independent (Case 2) random effects, and in a model with three independent (Case 3) or non independent (Case 4) random effects, (ii) testing that two variances are null in a model with three independent (Case 5) random effects, and (iii) testing that two variances are null in a model with two independent random effects (Case 6), i.e. testing the absence of random effects. The limiting distribution of the test statistics in each of these cases is respectively 0.5χ 2 0 + 0.5χ 2 1 (Case 1), 0.5χ 2 1 + 0.5χ 2 2 (Case 2), 0.5χ 2 0 + 0.5χ 2 1 (Case 3), 0.5χ 2 2 + 0.5χ 2 3 (Case 4), and w 0,2 χ 2 0 + 0.5χ 2 1 + (0.5 -w 0,2 )χ 2 2 (Case 5 and Case 6) (see Section 4.2). We emphasize here that the limiting distribution strongly depends on the correlation structure of the random effects.

Parameter estimation is performed either using the lmer function implemented in the R package lme4 [START_REF] Bates | Fitting linear mixedeffects models using lme4[END_REF], for the linear model, or

Real data analysis

The method is illustrated on two sets of real data. The first one is the famous dental growth dataset from [START_REF] Potthoff | A generalized multivariate analysis of variance model useful especially for growth curve problems[END_REF], in which the distance from the center of the pituitary gland to the pteryomaxillary fissure was measured at 4 different ages for 27 children (16 boys and 11 girls). This dataset is available in the R package mice. The second dataset comes from a study of coucal growth rates, available as a Dryad package [START_REF] Goymann | Data from: Sex roles, parental care and offspring growth in two contrasting coucal species[END_REF]. Body masses of 678 nestlings from two species (white-browed coucals and black coucals) were recorded every two days from their hatching date until they left the nests. In this paper, we only consider data from the white-browed coucals species, corresponding to the highest sample size (385 individuals).

A linear model is fitted to the dental growth data using the lme4 package, with two random effects as described in Section 5.1 (model M 1 ). More precisely, if we denote by y ij , 1 ≤ i ≤ 27, 1 ≤ j ≤ 4, the dental measurement of child i of sex x i , at age t j , the following model is considered, with a random slope and a random intercept:

A nonlinear model is fitted to the white-browed coucals data using the saemix package, with three random effects as described in Section 5.1 (model M 2 ). More precisely, if we denote by y ij , 1 ≤ i ≤ 385, 1 ≤ j ≤ n i the body mass of nestling i at age x j , we consider the following model:

where ϕ i1 is the asymptotic body mass of individual i, ϕ i2 the age in days at which individual i reaches half its asymptotic body mass and ϕ i3 the growth rate of individual i.

We test whether the variances of the inflexion point and the growth rate are equal to 0, and thus we consider the Case 4 described previously. We consider a diagonal covariance matrix Γ = diag(γ 2 1 , γ 2 2 , γ 2 3 ) and we test H 0 :

The limiting distribution of the LRT statistics is the mixture w 0,2 χ 2 0 + 0.5χ 2 1 + (0.5w 0,2 )χ 2 2 , where w 0,2 is defined in Section 4 and can be computed from the correlation coefficient between parameters γ 2 2 and γ 2 3 obtained from the Fisher information matrix.

The estimated Fisher information matrix Î is obtained as an output of the saemix package when the full model corresponding to H 1 is fitted to the data. From this matrix we can easily compute the covariance matrix

), and the corresponding correlation matrix Ĉ = diag( V ) -1/2 V diag( V ) -1/2 . Then, the correlation coefficient ρ12 needed to compute w 0,2 , is the element (1, 2) of Ĉ.

In our case, ρ12 = -0.644, leading to the three following weights: w 0,2 = 0.139, w 1,2 = 0.5 and w 2,2 = 0.361, and thus to the limiting distribution Tangent cone of Θ. We recall that Θ is defined as:

Now, since each term in the above product is a convex cone, we have that (Hiriart-Urruty and Lemarechal, 1996, Proposition 5.3.1.)). Therefore, the tangent cone of Θ at θ * is given by:

where T (S p + , Γ * ) is the tangent cone of the set of symmetric positive semidefinite matrices of size p × p at Γ * , and T (S J + , Σ * ) the tangent cone of the set of symmetric positive semi-definite matrices of size J × J at Σ * .

To identify T (S p + , Γ * ) and T (S J + , Σ * ), we can use the result established by Hiriart- [START_REF] Hiriart-Urruty | A Fresh Variational-Analysis Look at the Positive Semidefinite Matrices World[END_REF]. According to the authors, the tangent cone of S p + at A ∈ S p + is given by T A = {M ∈ S p | M u, u ≥ 0 for all u ∈ ker A}, where S p is the set of symmetric matrices of size p × p.

In our case, since Γ * = Γ * 1 0 0 0 , we have:

Proof of Corollary 1

Let V be a positive-definite matrix and C a closed convex cone of R q . We denote by C o = {x ∈ R q | x t y ≤ 0, ∀y ∈ C} its polar cone. We recall the following properties for the weights of the chi-bar-square distribution χ2 (V, C) [START_REF] Shapiro | Asymptotic distribution of test statistics in the analysis of moment structures under inequality constraints[END_REF][START_REF] Shapiro | Towards a Unified Theory of Inequality Constrained Testing in Multivariate Analysis[END_REF]:

then the first k weights {w i (q, V, C o ), i = 0, . . . , k -1} are zero, 3. if C contains a linear space of dimension l, for 1 ≤ l ≤ q, then the last l weights {w i (q, V, C o ), i = q -l + 1, . . . , q} are zero.

In our case, C = T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ and V = I -1 * , and we have for both cases mentioned in the corollary: (i) C = {0} p × {0} (p-r)(p-r+1)/2 × R r(p-r) × S r + × {0} J(J+1)/2 which is included in R r(p-r)+r(r+1)/2 , i.e. in a linear space of dimension q -(p + (p -r)(p -r + 1)/2 + J(J + 1)/2). Therefore, using properties (i) and (ii) above, the weights w i (q, I -1 * , C), for i = r(pr) + r(r + 1)/2 + 1, . . . , q are zero. Moreover, C contains R r(p-r) , i.e. a linear space of dimension r(p -r), which means using properties (i) and (iii) above, that the weights w i (q, I -1 * , C), for i = 0, . . . , r(p -r) -1 are zero. (ii) C = {0} p × {0} p-r × R r + × {0} J(J+1)/2 which is included in R r , a linear space of dimension q -p -(p -r) -J(J + 1)/2. It follows that the weights w i (q, I -1 * , C), for i = r + 1, . . . , q are zero. Then, since C does not contain any linear space of dimension k > 0, all the other weigths are non-zero.