Charlotte Baey

Estelle Kuhn

An R Package for Variance Components Testing in Linear and Nonlinear Mixed-Effects Models

Keywords: generalized mixed-effects models, nonlinear mixed-effects models, variance components, likelihood ratio test, random effects, R

The issue of variance components testing arises naturally when building mixed-effects models, to decide which effects should be modeled as fixed or random or to build parsimonious models. While tests for fixed effects are available in R for models fitted with lme4, tools are missing when it comes to random effects. The varTestnlme package for R aims at filling this gap. It allows to test whether a subset of the variances and covariances corresponding to a subset of the random effects, are equal to zero using asymptotic property of the likelihood ratio test statistic. It also offers the possibility to test simultaneously for fixed effects and variance components. It can be used for linear, generalized linear or nonlinear mixed-effects models fitted via lme4, nlme or saemix. Numerical methods used to implement the test procedure are detailed and examples based on different real datasets using different mixed models are provided. Theoretical properties of the used likelihood ratio test are recalled.

Introduction

Mixed-effects models are widely used in many fields of applications, ranging from agronomy, social science, medical science, to biology. One can distinguish three main types of mixed models: linear mixed models (LMMs), generalized linear mixed models (GLMMs) and nonlinear mixed models (NLMMs). Several packages and software are available to fit these types of model, for example the SAS (SAS Institute Inc. 2020) procedures MIXED and NLMIXED for linear, generalized linear and nonlinear mixed models, or the dedicated software NONMEM (ICON plc 2022) and Monolix (Lixoft 2021) for nonlinear mixed models, with an emphasis on pharmacodynamic and pharmacokinetic models for the former. In R (R Core Team 2023), three main packages are available to deal with these models: lme4 (Bates, Mächler, Bolker, and Walker 2015, for LMMs, GLMMs and NLMMs), nlme (Pinheiro, Bates, DebRoy, Sarkar, and R Core Team 2022, for LMMs and NLMMs) and saemix (Comets, Lavenu, and Lavielle 2017, for LMMS and NLMMs). The way they handle nonlinear mixed models differs: both lme4 and nlme rely on an approximation of the likelihood function based on a linearization of the model, while saemix uses an EM-type algorithm to derive the maximum likelihood estimate, without resorting to any approximation or linearization of the model.

When building a mixed-effects model, one is facing the delicate choice of which effects should be modeled as random. The objective can be twofold: on the one hand, being able to compare a simpler model with less random effects to a larger one can help choosing a parsimonious model, on the other hand, it can help identifying the sources of variability. This is the case for example in the context of nonlinear mixed-effect models applied to mechanistic models, where the model parameters are considered as random variables. In plant growth modeling for example (see [START_REF] Baey | Mixed-Effects Estimation in Dynamic Models of Plant Growth for the Assessment of Inter-Individual Variability[END_REF], identifying parameters that vary across individuals or genotypes and those that are constant allows to identify key biological processes exhibiting genotypic variability, the other being considered as speciesspecific. Understanding this variability is crucial to enhance our knowledge of genotype by environment interactions and then, to better understand plant adaptation to environmental constraints.

This issue has been studied from both a theoretical and a practical point of view. On the theoretical side, two main approaches have been explored. The first one follows the idea of variable selection. [START_REF] Chen | Random Effects Selection in Linear Mixed Models[END_REF] focused on linear mixed models and proposed a Bayesian formulation using mixture priors with a point mass at zero for the random effects. More recently, [START_REF] Ibrahim | Fixed and Random Effects Selection in Mixed Effects Models[END_REF]; [START_REF] Groll | Variable Selection for Generalized Linear Mixed Models by ℓ 1 -Penalized Estimation[END_REF] used a penalized likelihood approach to simultaneously select fixed and random effects. Selection criteria adapted to the context of mixed-effects models were also suggested by [START_REF] Vaida | Conditional Akaike Information for Mixed-Effects Models[END_REF]; [START_REF] Gurka | Selecting the Best Linear Mixed Model under REML[END_REF] and [START_REF] Delattre | A Note on BIC in Mixed-Effects Models[END_REF]. The second approach involves hypothesis testing based on variance components, with the seminal work of Stram andLee (1994, 1995) in linear mixed models (see also the review by [START_REF] Molenberghs | Likelihood Ratio, Score, and Wald Tests in a Constrained Parameter Space[END_REF]. Recently, [START_REF] Baey | Asymptotic Distribution of Likelihood Ratio Test Statistics for Variance Components in Nonlinear Mixed Effects Models[END_REF] exhibited the limiting distribution of the likelihood ratio test (LRT) statistic, for testing that the variances of any subset of the random effects are equal to zero.

On the practical side, several of the aforementioned methods have been implemented, e.g., in R. Among others, we can cite for example the glmmLasso [START_REF] Groll | glmmLasso: Variable Selection for Generalized Linear Mixed Models by ℓ 1 -Penalized Estimation[END_REF]) package, computing an ℓ 1 -penalized maximum likelihood estimator in the context of generalized linear mixed models, or the cAIC4 (Säfken, Rügamer, Kneib, and Greven 2021) package, implementing the conditional AIC criterion. As far as hypothesis testing is concerned, most of the available tools are designed for linear or generalized linear mixed-effects models, and are mostly adapted to outputs from lme4 package. The lmerTest package [START_REF] Kuznetsova | lmerTest Package: Tests in Linear Mixed Effects Models[END_REF] provides enhanced versions of functions anova() and summary() from lme4 package by computing p values associated to the tests for fixed effects. It also implements both the Satterthwaite's and the Kenward-Roger's approximations to correct p values in unbalanced or small sample size situations. Finally, it provides a step() function to simplify the random effects structure using a step-down approach. Each step is based on the computation of p values from the LRT for testing that the variance of each random effect is equal to zero. However, the asymptotic distribution used to compute these p values is a χ 2 , with the degree of freedom computed as the difference between the number of parameters under the alternative hypothesis and the number of parameters under the null hypothesis. In the context of mean components testing in linear and generalized linear models, the pbkrtest package [START_REF] Halekoh | A Kenward-Roger Approximation and Parametric Bootstrap Methods for Tests in Linear Mixed Models -The R Package pbkrtest[END_REF] provides two alternatives to the LRT in small sample size situations, based either on the Kenward-Roger approximation or on parametric bootstrap. The two aforementioned packages can only be used with models fitted with the lme4 package. The RLRSim package [START_REF] Scheipl | Size and Power of Tests for a Zero Random Effect Variance or Polynomial Regression in Additive and Linear Mixed Models[END_REF] implements the tests proposed by [START_REF] Crainiceanu | Likelihood Ratio Tests in Linear Mixed Models with One Variance Component[END_REF], based on the exact finite sample distribution of the LRT statistic. It can be used with a larger variety of packages and can handle not only models fitted from packages lme4 and nlme, but also from package gamm4 (Wood and Scheipl 2020) and from function gamm of package mgcv [START_REF] Wood | Generalized Additive Models: An Introduction with R. 2nd edition[END_REF]. However, only linear mixed models with one single random effect are currently covered by this package. In SAS, the option COVTEST can be used for linear mixed models in PROC MIXED and provides Wald tests for the variance components. More recently, a SAS macro %COVTEST performing a bootstrap parametric test has been proposed for linear mixed models.

To the best of our knowledge, there is no software or R package implementing the LRT for variance components in LMMS, GLMMs and NLMMs, and for an arbitrary subset of variance components corresponding to any subset of the random effects.

In this paper, we present the varTestnlme R package [START_REF] Baey | varTestnlme: Variance Components Testing for Linear and Nonlinear Mixed Effects Models[END_REF] dedicated to test that the variances of a subset of the random effects are equal to zero in LMMs, GLMMs and NLMMs with one level of random effects, using models that were fitted using either lme4, nlme or saemix. Moreover it is possible to test simultaneously mean and variance components. The varTestnlme package takes as inputs the two competitor models, fitted with the same package, and computes the asymptotic distribution of the LRT statistic associated with the two nested models corresponding to the null and alternative hypotheses. It provides the test statistic value, the limiting distribution and the associated p value. The theory behind this test has been established in [START_REF] Baey | Asymptotic Distribution of Likelihood Ratio Test Statistics for Variance Components in Nonlinear Mixed Effects Models[END_REF]. The purpose of the current paper is to provide hints for a practical implementation and a complete description of the different numerical methods used for the computation of each quantity involved. In particular, we detailed the algorithms used to sample from the limiting χ2 distribution and to estimate the weights of each component of the χ2 distribution. We also extend results of [START_REF] Baey | Asymptotic Distribution of Likelihood Ratio Test Statistics for Variance Components in Nonlinear Mixed Effects Models[END_REF], where only the cases of a diagonal or full covariance matrix for the random effects were exemplified, to the general case of a block-diagonal covariance matrix, covering a large range of covariance structures. The package varTestnlme is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=varTestnlme.

Mixed-effects models and the LRT procedure are recalled in Section 2. Practical implementations of the package are presented through various examples detailed in Section 3. Technical tools for the identification of the limiting distribution of the LRT statistic are detailed in Section 4. This section can be skipped if one is interested in the practical implementation of the functions of the package. Finally, current limitations and possible extensions are discussed in Section 5.

Mixed-effects models and likelihood ratio test procedure

Description of the mixed-effects models

We consider the following nonlinear mixed-effects model (Davidian and Giltinan 1995, p. 98;Pinheiro and Bates 2000, p. 306;Lavielle 2014, p. 24):

y i = g(φ i , x i ) + ε i , ( 1
)
where y i denotes the vector of observations of size m for individual i, 1 ≤ i ≤ n, g a nonlinear function, φ i the vector of individual parameters, x i the 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, ( 2
)
where β is the vector of fixed effects taking values in R p , U i and V i are covariates matrices of appropriate dimensions and b i is the centered vector of random effects. We assumed that the random vectors (b i ) 1≤i≤n are independent and identically distributed centered Gaussian vectors with covariance matrix Γ of size q × q , and that (ε i ) 1≤i≤n are independent and identically distributed centered Gaussian vectors with covariance matrix σ 2 I m , where I m is the identity matrix of size m×m. Note that such models can handle non identically distributed observations (y i ), for example due to the presence of covariates as illustrated in Section 3. However, multilevel random effects are not yet covered by the current version of the package. Finally, we assume that the sequences (ε i ) and (b i ) are mutually independent. We denote by c the total number of parameters and by θ = (β, Γ, σ 2 ) the vector of parameters taking value in the set Θ = {θ ∈ R c | β ∈ R p , Γ ∈ S q + , σ 2 > 0} where S q + denotes the set of positive definite matrices of size q.

Testing variance components

We consider general test hypotheses of the following form, to test for the nullity of r variances among q:

H 0 : θ ∈ Θ 0 against H 1 : θ ∈ Θ, (3)
where Θ 0 ⊂ Θ. Up to permutations of rows and columns of Γ, we can assume that we are testing for the nullity of the last r variances. We then write Γ in blocks as follows:

Γ = Γ 1 Γ ⊤ 12 Γ 12 Γ 2 ,
with Γ 1 a (q -r) × (q -r) matrix, Γ 2 a r × r matrix, Γ 12 a r × (q -r) matrix and where A ⊤ denotes the transposition of matrix A, for any matrix A.

We thus have:

Θ 0 = {θ ∈ R c | β ∈ R p , Γ 1 ∈ S q-r + , Γ 12 = 0 r,q-r , Γ 2 = 0 r,r , σ 2 > 0} Θ = {θ ∈ R c | β ∈ R p , Γ ∈ S q + , σ 2 > 0}.
(4)

The likelihood ratio test (LRT) statistic is defined by:

LRT n = -2 log sup θ∈Θ 0 L(y 1 , . . . , y n ; θ) sup θ∈Θ L(y 1 , . . . , y n ; θ) ,

where L(y 1 , . . . , y n ; θ) is the likelihood function.

As shown in [START_REF] Baey | Asymptotic Distribution of Likelihood Ratio Test Statistics for Variance Components in Nonlinear Mixed Effects Models[END_REF], assuming conditions (C1) and (C2) required for Theorem 1 to hold are fulfilled, the limiting distribution of the LRT statistic associated with the two hypotheses H 0 and H 1 defined in Equation 3 is:

LRT n ---→ n→∞ χ2 (I -1 * , T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ ), (5)
where I * is the Fisher information matrix, T (A, θ) is the tangent cone to space A at point θ, and where χ2 (V, C) denotes the χ2 distribution parametrized by the positive definite matrix V and the cone C. The χ2 distribution is a mixture of χ 2 distributions, where the degrees of freedom and the weights involved in the mixture depend both on V and C. More details about the computation of the χ2 distribution parameters are given in Section 4.

Testing simultaneously fixed effects and variance components

Results of the previous section can be extended to the case where one is testing simultaneously that a subset of the mean parameters and a subset of the covariance parameters are null. In this case, if we denote by I f := {k 1 , . . . , k r f } the set of indices corresponding to the fixed effects which means are tested equal to 0, the parameter space Θ 0 is slightly modified as follows:

Θ 0 = {θ ∈ R c | β k = 0, k ∈ I f , β k ∈ R, k / ∈ I f , Γ 1 ∈ S q-r + , Γ 12 = 0 r,q-r , Γ 2 = 0 r,r , σ 2 > 0} (6)
The theoretical result also holds in this test setting, leading to a χ2 distribution.

Practical implementation through examples

The varTestnlme package provides a unified framework for likelihood ratio tests of fixed and random effects of a linear, generalized linear or nonlinear mixed-effects models fitted either with the nlme, lme4 or saemix package. The main function varCompTest takes, in its simplest form, two arguments: m1 the fitted model under the alternative hypothesis H 1 and m0 the fitted model under the null hypothesis H 0 . It gives as outputs the test statistic value, the limiting distribution (i.e. the weights and degrees of freedom associated to each χ 2 component of the χ2 distribution), as well as the p value of the test.

In this section, we give an overview of how the package works, through three examples involving linear, generalized linear and nonlinear mixed-effects models. These examples gather cases where all the quantities characterising the limiting distribution can be obtained explicitly, as well as cases where the Fisher information matrix (FIM), but also the weights of the χ2 distribution should be approximated. In the following, we describe the datasets used to illustrate the package (Section 3.1), the different test hypotheses and the implementation of the models under each hypothesis (Section 3.2), and finally the results obtained using the varCompTest() function (Section 3.3). In Section 3.4, we illustrate the added-value of the FIM estimation proposed by the package, on two examples based on a linear and a nonlinear mixed models.

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 q q q q q q q q q q q q q q q q q q q q

F06 F01 F05 F07 F02 F08 F03 F04 F11 M14 M09 M15 M06 M04 M01 M10 F10 F09 M16 M05 M02 M11 M07 M08 M03 M12 M13 8

Data

For illustrative purposes, three datasets will be used throughout the paper, each of them covering one of the three types of models that can be treated with the varTestnlme package. These datasets are available in packages nlme, lme4 or datasets.

Orthodontal data

The first dataset comes from a study on dental growth [START_REF] Potthoff | A Generalized Multivariate Analysis of Variance Model Useful Especially for Growth Curve Problems[END_REF], where the distance between the pituitary and pterygomaxillary fissure was recorded every two years from the age of 8 to the age of 14, on 27 children, 16 boys and 11 girls (see Figure 1).

R> data("Orthodont", package = "nlme")

These data can be fitted using a linear mixed-effects model, with a random slope and a random intercept. Let us denote by y ij , 1 ≤ i ≤ 27, 1 ≤ j ≤ 4 the dental measurement of child i of sex x i at age t j . Then the model can be written as:

y ij = (β 1 + β 2 x i + b i1 ) + (β 3 + β 4 x i + b i2 )t j + ε ij , ε ij ∼ N (0, σ 2 ), (b i1 , b i2 ) ⊤ ∼ N (0, Γ) . ( 7
)
We denote by β = (β 1 , β 2 , β 3 , β 4 ) ⊤ the vector of fixed effects.

Bovine pleuropneumonia

The second dataset comes from a study on contagious bovine pleuropneumonia (cbpp) [START_REF] Lesnoff | Within-Herd Spread of Contagious Bovine Pleuropneumonia in Ethiopian Highlands[END_REF], where the number of new serological cases occurring during a given time period was recorded at 4 occasions, on 15 herds (see Figure 2). R> data("cbpp", package = "lme4")

This data can be fitted using a generalized linear mixed-effects model and a logistic regression.

Let us denote by y ij the number of serological cases in herd i and at time period t j . If we denote by g the logit function, then the model is given by:

E(y ij | b i ) = g -1 (b i + βt j ), b i ∼ N (0, Γ) (8)

Loblolly pine trees

The third dataset comes from a study on Loblolly pine trees [START_REF] Kung | Fitting Logistic Growth Curve with Predetermined Carrying Capacity[END_REF], where the growth of 14 trees was recorded on 6 occasions between the ages of 3 and 25 years (see Figure 3).

R> data("Loblolly", package = "datasets")

This data can be fitted using a nonlinear mixed-effects model. Let us denote by y ij , 1 ≤ i ≤ 14, 1 ≤ j ≤ 6 the height of tree i at age x j . We consider the following model:

y ij = Asym i + (R 0i -Asym i ) exp(-e lrc i x j ) + ε ij , ε ij ∼ N (0, σ 2 ) (Asym i , R 0i , lrc i ) ⊤ = β + b i , b i ∼ N 3 (0, Γ), (9)
where Asym i , lrc i and R 0i are the individual parameters corresponding respectively to the asymptote, the logarithm of the growth rate and the height at age 0 of tree i. 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 323 309 303 305

Preliminary step: Fitting the models under H 0 and H 1

The first step of the analysis consists in specifying H 0 and H 1 , the two hypotheses defining the test. Then, one needs to fit the model under both hypotheses, using nlme, lme4 or saemix.

The varTestnlme package automatically detects the structure of the models under H 0 and H 1 .

Note that both models m1 and m0 should be fitted with the same package, except when there is no random effects in m0 (see details below).

Linear model

We will consider three hypothesis testing configurations for the linear model presented in Equation 7. We detail below the null and alternative hypotheses in each case, and the code to fit the associated models.

1. Case 1: Testing that there is a random slope, in a model where the slope and the intercept are correlated:

H 0 : θ ∈ Θ 0 against H 1 : θ ∈ Θ, with Θ 0 = {β ∈ R 4 , Γ = γ 2 1 0 0 0 , σ 2 > 0} and Θ = {β ∈ R 4 , Γ = γ 2 1 γ 12 γ 12 γ 2 2 , σ 2 > 0}.
The syntax using lme4 and nlme packages is: 2. Case 2: Testing that there is a random slope, in a model where the slope and the intercept are independent: ), data = Orthodont, + method = "ML") R> lm2.h0.nlme <-lm1.h0.nlme 3. Case 3: Testing for the presence of random effects, in a model where the slope and the intercept are independent:

R> lm1.h1.
H 0 : θ ∈ Θ 0 against H 1 : θ ∈ Θ, with Θ 0 = {β ∈ R 4 , Γ = γ 2 1 0 0 0 , σ 2 > 0} and Θ = {β ∈ R 4 , Γ = γ 2 1 0 0 γ 2 2 , σ 2 > 0} R> lm2.h1.
H 0 : θ ∈ Θ 0 against H 1 : θ ∈ Θ, with Θ 0 = {β ∈ R 4 , Γ = 0 0 0 0 , σ 2 > 0} and Θ = {β ∈ R 4 , Γ = γ 2 1 0 0 γ 2 2 , σ 2 > 0} R> lm3.h1.lme4 <-lm2.h1.lme4 R> lm3.h1.nlme <-lm2.h1.nlme R> lm3.h0 <-lm(distance ~1 + Sex + age + age * Sex, data = Orthodont)

Generalized linear model

In the model considered in Equation 8, there is only one random effect. We will then consider the following test:

H 0 : θ ∈ Θ 0 against H 1 : θ ∈ Θ, with Θ 0 = {β ∈ R, Γ = 0} and Θ = {β ∈ R, Γ ≥ 0}.
The corresponding code using lme4 is:

R> glm1 <-glmer(cbind(incidence, size -incidence) ~period + (1 | herd), + family = binomial, data = cbpp) R> glm0 <-glm(cbind(incidence, size -incidence) ~period, + family = binomial, data = cbpp)

Nonlinear model

Let us consider the nonlinear model defined in Equation 9. Here we will carry out the following test:

H 0 : θ ∈ Θ 0 against H 1 : θ ∈ Θ, with Θ 0 = {β ∈ R 3 , Γ =    γ 2 1 0 0 0 0 0 0 0 0    , σ 2 > 0} and Θ = {β ∈ R 3 , Γ =    γ 2 1 0 0 0 γ 2 2 0 0 0 γ 2 3    , σ 2 > 0},
i.e., that only the asymptote is random.

The corresponding code using nlme and lme4 is: Using saemix, the syntax is given by: name.group = "Seed", name.predictors = "age", name.response = "height") R> ctrl <-saemixControl(print = FALSE, save = FALSE, save.graphs = FALSE) R> nlm1.saemix <-saemix(saemix.modelH1, saemix.data, ctrl) R> nlm0.saemix <-saemix(saemix.modelH0, saemix.data, ctrl)

R> start <-c(Asym = 103, R0 = -8.5, lrc = -3.2) R> nlm1.
R> modelSSasymp <-function(psi, id, xidep) { + Asym <-psi[id, 1] + R0 <-psi[id, 2] + lrc <-psi[id, 3] + age <-xidep[, 1] + ypred <-Asym + (R0 -

Variance component testing

Once the models have been fitted using lme4, nlme or saemix, the test can be run using the varCompTest() function of the varTestnlme package. The main function is:

R> varCompTest(m1, m0, + control = list(M = 5000, parallel = FALSE, nbcores = 1, B = 1000), + fim = "extract", pval.comp = "bounds", output = TRUE)

where:

• m1 is the model fitted under H 1 , i.e., an object of class 'merMod', 'glmerMod', 'lme', 'nlme' or 'SaemixObject'.

• m0 is the model fitted under H 0 using the same package as for m1, except when no random effects are present under H 0 (when testing for the absence of random effects).

In this case, m0 should be fitted using lm() for linear, glm() for generalized linear, and nls() for nonlinear models (these functions are available in the stats package).

• control is a list with: M, the Monte Carlo sample size for the computation of the χ2 weights (5000 by default, see Algorithm 1); parallel a boolean to specify whether the Monte Carlo sampling should be done in parallel (FALSE by default); nbcores the number of cores to be used with parallel computing (1 by default); and B the size of the bootstrap sample used to estimate the Fisher information matrix (1000 by default).

• fim could be either "extract" (the default) if the FIM should be extracted from the fitted model m1; "compute" if it should be computed using parametric bootstrap implemented in the package (see Section 4.1 for more details on the method); or FIM a user-defined matrix to be used as the Fisher information matrix.

• pval.comp specifies the way the p value of the test is computed, and could be either "bounds" (the default), in which case only bounds on the true p value are computed (see Section 4.5 for more details), "approx", in which case a Monte Carlo estimation of the exact p value is provided (see Section 4.4 for more details), or "both" for a combination of both approaches. In the case where the weights are known explicitly, no approximation is made and the exact weights are return.

• output is a boolean indicating if any outputs should be printed out in the console. The default value is TRUE, in which case it displays information about the fitted models and the results of the test.

Note that the χ2 weights approximation can be time consuming, especially when the Monte Carlo sampling size is high, or when the number of components in the χ2 distribution increases.

It is recommended to first run the function using the default setting, i.e., with pval.comp = "bounds", and to possibly re-run it with pval.comp = "approx" if more precision is needed.

The function returns an object of classes 'htest' and 'vctest' and as such, can be printed using the provided print() function of varTestnlme or the print.htest() function of package EnvStats [START_REF] Millard | EnvStats: An R Package for Environmental Statistics[END_REF]. In particular, it contains the following slots:

• null.value: The value of the tested parameters under the null hypothesis.

• alternative: The value of the parameters under the alternative hypothesis.

• statistic: The value of the LRT statistic.

• method: The name of the statistical test.

• parameters: A list with the χ2 distribution parameters: df the degrees of freedom, weights the weights associated with each component of the limiting χ2 distribution, sdWeights the standard deviation associated with the estimation of each weights (equal 0 if exact weights are available) and FIM the estimate of the Fisher information matrix.

• pvalue: A named numeric vector with four elements: pvalue.weights, the p value obtained using Equation 17(equals NA if weights were not computed, e.g., if option pval.comp was set to "bounds"), pvalue.sample: the p value obtained using Equation 18 (equals NA if weights were not computed or if exact weights were available and no sampling was done), pvalue.lowerbound and pvalue.upperbound the bounds on the p value, obtained from Equation 19.

We describe in the following subsections the practical implementation of the tests defined in Section 3.2 and detail the outputs obtained in each case.

Linear model

Case 1: Testing that the variance of age is null, in a model with two correlated random effects. We first run the function with the default arguments.

R> vt <-varCompTest(lm1.h1.lme4, lm1.h0.lme4) R> print(vt)

Variance components testing in mixed effects models Testing that the variance of the random effect associated to scale(age) is equal to 0 Likelihood ratio test statistic: LRT = 0.8331072 p-value from exact weights: 0.5103454

Using the summary() function, we get the following output:

R> summary(vt)

Variance components testing in mixed effects models Testing that: variance of the random effect associated to scale(age) is equal to 0 against the alternative that: variance of the random effect associated to scale(age) > 0

Likelihood ratio test statistic: LRT = 0.8331072

Limiting distribution: mixture of 2 chi-bar-square distributions with degrees of freedom 1 2 associated (exact) weights: 0.5 0.5 p-value of the test: from exact weights: 0.5103454

The LRT statistic is computed, and the asymptotic distribution is identified as a mixture between two χ 2 distributions with degrees of freedom 1 and 2. In this case we can compute the exact weights of the χ2 distribution and hence the exact p value. Using the EnvStats, it is possible to print the results of the test using function print.htest().

R> library("EnvStats") R> print.htestEnvStats(vt) If we re-run the function with the option pval.comp = "both", we get the same results since the weights are explicit in this example. This time, we run the function with fits from nlme to show that results are identical.

Results of Hypothesis

R> vt <-varCompTest(lm1.h1.nlme, lm1.h0.nlme) R> summary(vt)

Variance components testing in mixed effects models Testing that: variance of the random effect associated to age is equal to 0 against the alternative that: variance of the random effect associated to age > 0 Likelihood ratio test statistic: LRT = 0.8331072

Limiting distribution: mixture of 2 chi-bar-square distributions with degrees of freedom 1 2 associated (exact) weights: 0.5 0.5 p-value of the test: from exact weights: 0.5103454

Case 2: Testing that the variance of age is null in a model with two independent random effects. The number of components in the mixture is the same as in Case 1, but the degrees of freedom have been shifted downward. Note that the weights also have simple analytical expressions in this case, as well as the p value. Results are presented using the lme4 package, but similar results are obtained with the nlme package.

R> vt <-varCompTest(lm2.h1.lme4, lm2.h0.lme4) R> summary(vt)

Variance components testing in mixed effects models Testing that: variance of the random effect associated to age is equal to 0 against the alternative that: variance of the random effect associated to age > 0 Likelihood ratio test statistic: LRT = 0.5304106

Limiting distribution: mixture of 2 chi-bar-square distributions with degrees of freedom 0 1 associated (exact) weights: 0.5 0.5 p-value of the test: from exact weights: 0.2332171

Limiting distribution: mixture of 3 chi-bar-square distributions with degrees of freedom 0 1 2 p-value of the test: bounds on p-value: lower 0.05851384 upper 0.2049047

Depending on the desired level of the test, more precision may be needed on the p value.

For models fitted with nlme, this can be obtained using the options pval.comp = "both" or pval.comp = "approx" and fim = "extract" to use the FIM computed by nlme, or fim = "compute" to approximate the FIM using parametric bootstrap, with a default bootstrap sample size of 1000.

R> vt <-varCompTest(nlm1.nlme, nlm0.nlme, pval.comp = "both", + fim = "compute") R> summary(vt)

Variance components testing in mixed effects models Testing that: covariance matrix of R0 and lrc is equal to 0 against the alternative that: covariance matrix of R0 and lrc > 0 Likelihood ratio test statistic: LRT = 2.519869

Limiting distribution: mixture of 3 chi-bar-square distributions with degrees of freedom 0 1 2 associated (exact) weights: 0.2450473 0.5000000 0.2549527 p-value of the test: from exact weights: 0.128533 bounds on p-value: lower 0.05620995 upper 0.1980462

Here again, the weights are exact once the FIM is known. For information purposes, the above code took approximatively 175 seconds on a laptop with a 8 cores Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz processor.

For nonlinear models fitted with lme4, to the best of our knowledge no method is available to extract the FIM (although merDeriv, Wang and Merkle 2018, provides the FIM for linear and generalized linear models fitted via lme4). Only the option fim = "compute" can then be used. We get the following results, which are very similar to the ones obtained with the nlme package. However the code was a bit longer to run, and took 430 seconds.

R> vt nlm0.lme4,pval.comp A linear mixed-effects model with two correlated random effects was fitted to 500 simulated datasets of size n = 100.

y ij = β 1 + β 2 t j + b i1 + b i2 t j + ε ij , i = 1, . . . , n, j = 1, . . . , m b i1 b i2 ∼ N 0 0 , γ 2 1 γ 12 γ 12 γ 2 2 , ε ij ∼ N (0, σ 2 )
We took β 1 = 5, β 2 = 7, γ 1 = 0.8, γ 2 = 1, γ 12 = 0.4 and σ 2 = 1.2 2 . The FIM was extracted and estimated for each of the simulated datasets, with a bootstrap sample size B ∈ {100, 300}. Results are given in Table 1. Note that small discrepancies can be observed for the nlme package when running the replication material, due to different versions of the linear algebra library and the way it handles matrix multiplication. Results show that good results are obtained for small values of B. Interestingly, the bootstrap estimate lead to better coverages than the extracted estimate obtained using lme4, for variance components. This is crucial since the FIM components involved in the χ2 weights computation are those corresponding to the variance components.

Tools for the identification of the limiting distribution

In this section, we detail the technical tools used to identify the key quantities characterizing the limiting χ2 distribution, which is a mixture of χ 2 distributions. These quantities are the degrees of freedom and the associated weights of the χ 2 distributions in the mixture and depend on the FIM and on the tangent cone. Note that an exact expression of the tangent cone can always be derived in the general case of a block-diagonal covariance matrix Γ. However the FIM is not always available in a closed form, especially for nonlinear mixed-effects models. In Section 4.1, we present the parametric bootstrap procedure used to estimate the FIM, in Section 4.2, we give details on how to identify the tangent cone T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ , and in Section 4.3 we describe the Monte Carlo algorithm used to compute the weights of each χ 2 distribution. We also explain how the p value of the test can be computed, and how bounds can be obtained.

Estimation of the Fisher information matrix

In linear mixed-effects models, the FIM can be computed exactly. However this is not the case for nonlinear mixed-effects models. The existing R packages provide in some cases an approximation of the FIM based for example on a linearisation of the model, or on a block-diagonal approximation of the FIM where a zero correlation is assumed between fixed effects and variance components. Nevertheless, in general the FIM is not block-diagonal. Moreover these estimates are in general only valid for large samples, leading to possible wrong conclusions when considering small samples. Specific tools have to be developed for such settings.

In the varTestnlme package, we provide three options for the FIM. The first option is specifically adapted to small samples, and is more accurate in this setting albeit more time consuming. It consists in building an estimate of the FIM via parametric bootstrap. Let us denote by θ1 the estimate of θ obtained for the alternative model under H 1 , with lme4, nlme or saemix. We generate B bootstrap samples (φ u, * i , y u, * i ), 1 ≤ i ≤ n, 1 ≤ u ≤ B from the alternative model, using the parameter value θ1 and Equations 1-2. Then, for each of these bootstrap samples, we compute the associated estimate θ * u . The bootstrap estimate of the FIM is then given by the empirical covariance matrix between the bootstrap estimates:

Î * = 1 B B u=1 θ * u ( θ * u ) ⊤ -θ * ( θ * ) ⊤ , (10) where θ
* = 1 B B u=1
θ * u . The second option consists in using the FIM estimate provided by the package that was used to fit the models, when this estimate is available. Finally, the third option relies on a user-defined positive semi-definite matrix which should be an estimate of the FIM. It can be for example computed using the R package PFIM 4.0 [START_REF] Dumont | PFIM 4.0, An Extended R Program for Design Evaluation and Optimization in Nonlinear Mixed-Effect Models[END_REF] based on Monte Carlo tools or by implementing the method presented in [START_REF] Delattre | Estimating Fisher Information Matrix in Latent Variable Models Based on the Score Function[END_REF] based on stochastic approximation tools.

Computation of the tangent cone

The identification of the tangent cone is crucial since it allows to identify the number of components in the χ2 distribution, and the degrees of freedom associated to each of these components. Indeed, results from [START_REF] Shapiro | Asymptotic Distribution of Test Statistics in the Analysis of Moment Structures under Inequality Constraints[END_REF] state that, provided that one can found the largest linear subspace contained in the tangent cone, and the smallest linear subspace containing the tangent cone, then the degrees of freedom of the χ 2 components take values between the dimensions of these two subspaces. In this section, we give the exact expression of the tangent cone, and we identify the dimensions of the two aforementioned linear subspaces in the general case of a block-diagonal covariance matrix Γ.

In [START_REF] Baey | Asymptotic Distribution of Likelihood Ratio Test Statistics for Variance Components in Nonlinear Mixed Effects Models[END_REF], it was highlighted that the limiting distribution of the LRT statistic depends on the structure of Γ and its expression was exhibited in two common cases, when Γ was either diagonal or full. In the general case, Γ can always be written in the following block-diagonal form:

Γ =       Γ 1 0 . . . 0 0 Γ 2 . . . 0 . . . . . . . . . . . . 0 . . . 0 Γ K       , ( 11
)
where K ∈ N * , and for all 1 ≤ k ≤ K, Γ k is a positive semi-definite matrix of size r k × r k with K 1 r k = q. In the sequel, we assume that the blocks Γ 1 , . . . , Γ K are full i.e., all the covariances inside these blocks are non-null. This is how the covariance matrix of random effects is implemented in package lme4 (Bates et al. 2015, p. 7, Section 2.2). In the appendix of their paper, the authors mentioned that more general structures for the covariance matrix Γ could be implemented using modularized functions of the lme4 package. Nevertheless these functionalities are not currently treated in the varTestnlme package. In nlme, more sophisticated covariance structures can be used, and in saemix, no specific restriction is imposed on the covariance matrix structure. However in varTestnlme we restrict ourselves to the blockdiagonal case with full blocks. This restriction allows for a simple block-diagonal structure of the covariance matrix and eases the computation.

When the covariance matrix Γ is block-diagonal, three cases can arise when testing that a subset of variances are null: i) we are testing that R blocks among the K blocks Γ 1 , . . . , Γ K are null, ii) we are testing that R sub-blocks of the K blocks Γ 1 , . . . , Γ K are null, or iii) a mixture of i) and ii). Let us denote by ℓ 0 the number of blocks where covariances are tested equal to 0 without testing that the corresponding variances are equal to 0 (only testing nondiagonal elements), by ℓ 1 the number of blocks which are tested entirely equal to 0, and by ℓ 2 the number of blocks where sub-blocks are tested equal to 0. By sub-block, we mean that we are testing a sub-matrix which is strictly smaller than the block matrix from which it was extracted. We assume that 0

≤ ℓ 0 ≤ K, 0 ≤ ℓ 1 ≤ K, 0 ≤ ℓ 2 ≤ K and 1 ≤ ℓ 0 + ℓ 1 + ℓ 2 ≤ K.
Without loss of generality, and up to a permutation of rows and columns of Γ, we can assume that the blocks are in the following order: first, the blocks which are not tested, then the blocks where only covariances are tested, next the blocks which are tested entirely equal to zero, and finally the blocks where diagonal sub-blocks are tested equal to zero. Then, the null and alternative hypotheses in Equation 3 can be written as:

Θ 0 = {θ ∈ R c | β k = 0, k ∈ I f , β k ∈ R, k / ∈ I f , ∀k = 1, . . . , K -(ℓ 0 + ℓ 1 + ℓ 2 ), Γ k ∈ S r k + , ∀k = K -(ℓ 0 + ℓ 1 + ℓ 2 ) + 1, . . . , K -(ℓ 1 + ℓ 2 ), Γ k ∈ S r k + , with t k covariances null, ∀k = K -(ℓ 1 + ℓ 2 ) + 1, . . . , K -ℓ 2 , Γ k = 0, ∀k = K -ℓ 2 + 1, . . . , K, Γ k = Γk 0 0 0 , with Γk ∈ S r k -s k + , σ 2 > 0} Θ = {θ ∈ R c | β ∈ R p , Γ k ∈ S r k + , σ 2 > 0}
From the general expressions of Θ 0 and Θ in Equation 4, we can derive the expressions of T (Θ, θ * ) and T (Θ 0 , θ * ) using the results of Hiriart-Urruty and Malick (2012), and consequently the expression of the closed convex cone C := T (Θ, θ * ) ∩ T (Θ 0 , θ * ) ⊥ , involved in the χ2 distribution in Equation 5.

We obtain:

C = {0} a × R r f tested fixed effects × K-(ℓ 1 +ℓ 2 ) k=K-(ℓ 0 +ℓ 1 +ℓ 2 )+1 R t k covariances tested without testing the corresponding variances × K-ℓ 2 k=K-(ℓ 1 +ℓ 2 )+1 S r k + fully tested blocks × K k=K-ℓ 2 +1 {R s k (r k -s k ) × S s k + } sub-blocks tested , ( 12
)
where

a = p -r f non tested fixed effects + K-(ℓ 0 +ℓ 1 +ℓ 2 ) k=1 r k (r k -1) 2 non tested blocks + K k=K-ℓ 2 +1 (r k -s k )(r k -s k + 1) 2
untested sub-blocks in blocks where a sub-block is tested

+ 1 residual variance . ( 13
)
To identify the components of the χ2 mixture, i.e., the degrees of freedom and the weights of all the χ 2 distributions involved in the mixture, we can use the properties enounced by [START_REF] Shapiro | Asymptotic Distribution of Test Statistics in the Analysis of Moment Structures under Inequality Constraints[END_REF], stating that if C contains a linear space of dimension d 1 , the first d 1 weights of the mixture are null, and if C is included in a linear space of dimension d 2 , the last (c-d 2 ) weights of the mixture are null. The χ2 distribution χ2 (I -1 * , C) is then a mixture of (d 2 -d 1 + 1) χ 2 distributions with degrees of freedom varying between d 1 and d 2 .

According to the general formulation of the cone C given in Equation 12, we have:

d 1 = r f + K-(ℓ 1 +ℓ 2 ) k=K-(ℓ 0 +ℓ 1 +ℓ 2 )+1 t k + K k=K-ℓ 2 +1 s k (r k -s k ), d 2 = c -a. (14)
In particular, if only blocks of variances are tested (i.e., if r f = 0, ℓ 0 = 0 and ℓ 2 = 0), then d 1 = 0 and there is a Dirac mass at 0 in the mixture.

From Equations 13 and 14, we can see that the number of elements in the χ 2 mixture only depends on the number of variances being tested, and on the structure of the covariance matrix. When r f fixed effects are tested simultaneously with a set of variance components, the number of elements in the mixture is the same as in the case where only the set of variance components is tested, but the degrees of freedom of each element in the χ 2 mixture is shifted upward by r f .

As an example, let us consider a model with 3 random effects, with

φ i = β + b i , b i ∼ N (0, Γ), ε i ∼ N (0, σ 2 I m )
with Γ a positive definite covariance matrice and σ 2 > 0. Let us consider the following hypotheses:

H 0 : θ ∈ Θ 0 against H 1 : θ ∈ Θ, (15)
with

Θ 0 = {β ∈ R 3 , Γ =    Γ 2 1 0 0 0 0 0 0 0 0    , σ 2 > 0} Θ = {β ∈ R 3 , Γ =    Γ 2 1 Γ 12 Γ 13 Γ 12 Γ 2 2 Γ 23 Γ 13 Γ 23 Γ 2 3    , σ 2 > 0}.
Since no fixed effect is tested, we have r f = 0. We can consider Γ as a block-diagonal matrix with only K = 1 block. In this block of size 3 × 3, we are testing that a sub-block of size 2 × 2 is equal to 0. Using the above notations, we have ℓ 0 = 0, ℓ 1 = 0, ℓ 2 = 1, r 1 = 3, s 1 = 2 and c = 10. We thus have: a = 3 -0 + 0 + (3-2)(3-2+1) 2 + 1 = 5, d 1 = 0 + 1(3 -1) = 2 and d 2 = 10 -5 = 5. The limiting distribution of the LRT associated to the two above hypotheses is then a mixture of d 2 -d 1 + 1 = 4 χ 2 distributions with degrees of freedom 2, 3, 4 and 5.

Computation of the χ2 weights

In some simple cases, the χ2 weights have an analytical expression, which may depends on the FIM. In general however, the weights are non explicit and have to be estimated via Monte Carlo techniques as described in Algorithm 1.

In varTestnlme, when the weights of the χ2 distribution are not available in a closed form, they are estimated according to the procedure proposed by [START_REF] Silvapulle | Constrained Statistical Inference: Order, Inequality, and Shape Constraints[END_REF], using an estimate of the FIM. The general idea relies on the simulation of M i.i.d. realizations X 1 , . . . , X M from the limiting chi-bar-square distribution using its definition as the norm of the projection of a multivariate Gaussian random variable on a closed convex cone. More precisely, let C be a closed convex cone of R c , V a positive definite matrix of size c × c and Z ∼ N (0, V ). Then, the random variable X defined below follows a χ2 distribution with parameters V and C as detailed in [START_REF] Silvapulle | Constrained Statistical Inference: Order, Inequality, and Shape Constraints[END_REF]:

X = Z ⊤ V -1 Z -min θ∈C (Z -θ) ⊤ V -1 (Z -θ). ( 16
)
The algorithm first simulate M realizations of a multivariate Gaussian distributions, Z 1 , . . . , Z M , and then computes each term of Equation 16 using quadratic programming when necessary to get the projection on C. This gives M realizations X 1 , . . . , X M following the χ2 distribution χ2 (V, C). Since the number of components in the χ2 mixture and the associated degrees of freedom are known, the weights of each χ 2 component can be retrieved by defining a system of linear equations in w, the vector of weights. More details are given in Algorithm 1.

Computation of the p value

Let F d be the cumulative distribution function (CDF) of the χ 2 distribution with d degrees of freedom. The p value of the test can then be estimated in two different ways, both of them Algorithm 1 Estimation of the weights of the χ2 distribution using Monte Carlo simulations. Note that due to the (pseudo-)randomness of the sequence (Z i ), results will differ from one run of the algorithm to the other. For reproducible results, one might want to fix the seed beforehand.

1: Define Î-1 an estimate of the inverse of the FIM under H 1 , and let c d 1 +2 < • • • < c d 2 be a sequence of non-negative increasing numbers (see Figure 4).

2: for i = 1, . . . , M do 1. simulate Z i ∼ N (0, Î-1 ) 2. compute X i := Z ⊤ i Î-1 Z i -min θ∈C (Z i -θ) ⊤ Î-1 (Z i -θ),
using quadratic programming when C is the non-negative orthant R c ′ , for some c ′ ≤ c, and using general nonlinear optimization tools otherwise. In varTestnlme, we use respectively the quadprog (Turlach, Weingessel, and Moler 2019) and the alabama (Varadhan 2022) packages.

3: end for

4: Compute matrix A of size (d 2 -d 1 + 1) × (d 2 -d 1 + 1): A =         1 1 1 . . . 1 1 1 0 1 . . . 1 0 F d 1 (c d 1 +2 ) F d 1 +1 (c d 1 +2 ) F d 1 +2 (c d 1 +2 ) . . . F d 2 -1 (c d 1 +2 ) F d 2 (c d 1 +2 ) . . . . . . . . . • • • . . . . . . F d 1 (c d 2 ) F d 1 +1 (c d 2 ) F d 1 +2 (c d 2 ) . . . F d 2 -1 (c d 2 ) F d 2 (c d 2 )        
, where F i is the cumulative distribution function of a χ 2 distribution with i degrees of freedom.

N.B. The expression of A given above corresponds to the case where (d 2 -d 1 + 1) is even.

In the case where (d 2 -d 1 + 1) is odd, the last column has a '1' in the second row, instead of a '0'. 5: Compute vector κ of size (d 2 -d 1 + 1) :

κ =         1 1/2 1 M M i=1 1 X i ≤c d 1 +2
. . . being computed in the varTestnlme:

p1 = d 2 j=d 1 ŵj (1 -F j (LRT n )), (17) p2
= 1 M M i=1 1 X i ≥LRT n , ( 18
)
where ŵj is the estimated weight associated to the χ 2 distribution with j degrees of freedom, and where X 1 , . . . , X M are simulated according to the limiting χ2 distribution (see Algorithm 1 for more details on the notations).

Note that for very small values of the real p value, a very large sample size M would be needed in order to get a non-zero estimate p2 .

Bounds on p values

Since the simulation of X 1 , . . . , X M can be time consuming, and since F i (c) ≥ F j (c) for i < j, it is possible to compute bounds on the p value of the test. Indeed, since the sum of all the weights is equal to 1, and the sum of even (resp. odd) weights is equal to 1/2, natural bounds are given for the p value by [START_REF] Silvapulle | Constrained Statistical Inference: Order, Inequality, and Shape Constraints[END_REF]:

1 - F d 1 (LRT n ) + F d 1 +1 (LRT n ) 2 ≤ p ≤ 1 - F d 2 -1 (LRT n ) + F d 2 (LRT n ) 2 . ( 19
)
Even if those bounds can be crude, they are available in closed form and might bring enough information in some cases to support or reject the null hypothesis. By default, the varTestnlme package returns these bounds, and if more precision is needed by the user, it is possible to re-run the analysis in order to estimate the weights of the limiting distribution.

Summary and discussion

In this paper, we present the R package varTestnlme, which provides a unified framework for variance components testing in linear, generalized linear and nonlinear mixed-effects models fitted with nlme, lme4 and saemix. The main function varCompTest() performs the comparison of two nested models m1 and m0. It is thus possible to test that the variances of a subset of random effects is null, taking into account correlation between them, in a model with one level of random effects. It is also possible to test that one or several correlations are null, and to test simultaneously if some fixed effects are null. The tests are performed using likelihood ratio test and rely on the asymptotic properties of the LRT statistic. Additional tools are provided to approximate the Fisher information matrix when it is needed to compute the limiting distribution of the LRT statistic.

Possible future developments of the package include in particular the consideration of more than one level of random effects, to account for multiple nested random effects. Another perspective would be to allow for more general structures for the covariance matrix Γ of the random effects. This would encompass for example parametric structures such as compound symmetry, auto-regressive or with spatial correlation, as proposed in the nlme package.

Another interesting development would be the generalisation to models with non-Gaussian random effects, for example using heavy-tailed distributions. Specific tools would be needed to extend the theoretical results to other families of distributions. Finally, it would be useful to derive non-asymptotic results allowing to apply the test to small sample sizes.

Figure 1 :

 1 Figure1: Data from the orthodontal study: distance between the pituitary and the pterygomaxillary fissure as a function of age, by subject.

Figure 2 :

 2 Figure 2: Data from the cbpp study: number of new cases per period, by herd.

Figure 3 :

 3 Figure 3: Data from the loblolly study: height as a function of age, by tree.

  nlme <-nlme(height ~SSasymp(age, Asym, R0, lrc), + fixed = Asym + R0 + lrc ~1, random = pdDiag(Asym + R0 + lrc ~1), + start = start, data = Loblolly) R> nlm0.nlme <-nlme(height ~SSasymp(age, Asym, R0, lrc), + fixed = Asym + R0 + lrc ~1, random = pdDiag(Asym ~1), + start = start, data = Loblolly) R> nlm1.lme4 <-nlmer(height ~SSasymp(age, Asym, R0, lrc) + ~(0 + Asym + R0 + lrc || Seed), start = start, data = Loblolly) R> nlm0.lme4 <-nlmer(height ~SSasymp(age, Asym, R0, lrc) + ~(0 + Asym | Seed), start = start, data = Loblolly)

  Asym) * exp( -exp(lrc) * age) + return(ypred) + } R> psi0 <-matrix(start, ncol = 3, + dimnames = list(NULL, c("Asym", "R0", "lrc"))) R> saemix.modelH1 <-saemixModel(model = modelSSasymp, + description = "Asymptotic regression", psi0 = psi0, + transform.par = c(0, 0, 0), covariance.model = diag(3), + omega.init = 10 * diag(3), error.model = "constant") R> saemix.modelH0 <-saemixModel(model = modelSSasymp, + description = "Asymptotic regression", psi0 = psi0, + transform.par = c(0, 0, 0), covariance.model = diag(c(1, 0, 0)), + omega.init = 10 * diag(3), error.model = "constant") R> saemix.data <-saemixData(name.data = Loblolly, +

  = 0.5103454 pvalue.upperbound = 0.5103454

Figure 4 :

 4 Figure 4: Illustration of the choice of the non-increasing sequence c d 1 +2 < • • • < c d 2 used in Algorithm 1 to approximate the χ2 weights.

Table 1 :

 1 Empirical coverage of the 95% asymptotic confidence intervals according to the type of FIM estimate, in a linear mixed-effect model with two correlated random effects, computed on 500 repetitions.

	= "both",

Case 3: Testing the presence of randomness in the model, i.e., testing that the variances of age and of the intercept are null, in a model with two independent random effects.

R> vt lm3.h0) R> summary (vt) Variance components testing in mixed effects models Testing that: covariance matrix of Intercept and age is equal to 0 against the alternative that: covariance matrix of Intercept and age > 0 Likelihood ratio test statistic: LRT = 50.13311

Limiting distribution: mixture of 3 chi-bar-square distributions with degrees of freedom 0 1 2 p-value of the test: bounds on p-value: lower 7.18311e-13 upper 7.215163e-12

Bounds based on the smallest and biggest degrees of freedom of the chi-barsquare distribution components. Re-run with option 'pval.comp="both" or pval.comp="comp" to approximate the weights of each chi-bar-square component and the p-value.

A message is printed in addition to the results, to give some explanation about the bounds on the p value. In this case, using the provided upper bound we see that the exact p value is smaller than 10 -11 , which is enough in practice to reject the null hypothesis without computing the χ2 weights and the associated p value. It is also possible to compute the weights using the option pval.comp = "both". In that case, the weights are explicit but depend on the FIM. Since we are dealing with a linear mixed-effects model, the FIM has an analytical expression and can be extracted from nlme or lme4 packages using the option fim = "extract". This is the default behaviour of varCompTest function. Note that results can differ between packages since the methods used to compute the likelihood and hence the FIM is not the same in nlme and lme4.

Using nlme package: Variance components testing in mixed effects models Testing that: covariance matrix of (Intercept) and age is equal to 0 against the alternative that: covariance matrix of (Intercept) and age > 0 Likelihood ratio test statistic: LRT = 50.13311

Limiting distribution: mixture of 3 chi-bar-square distributions with degrees of freedom 0 1 2 associated (exact) weights: 0.3609427 0.5000000 0.1390573 p-value of the test: from exact weights: 2.525181e-12 bounds on p-value: lower 7.18311e-13 upper 7.215163e-12

Generalized linear model

Results from the test of the null hypothesis that there is no random effect associated with the herd are given by: R> vt <-varCompTest(glm1, glm0) R> summary(vt)

Variance components testing in mixed effects models Testing that: variance of the random effect associated to (Intercept) is equal to 0 against the alternative that: variance of the random effect associated to (Intercept) > 0 Likelihood ratio test statistic: LRT = 14.00527

Limiting distribution: mixture of 2 chi-bar-square distributions with degrees of freedom 0 1 associated (exact) weights: 0.5 0.5 p-value of the test: from exact weights: 9.114967e-05

Nonlinear model

Results from the test of the null hypothesis that only the asymptote is random against the alternative hypothesis of a diagonal covariance matrix between the random effects, using nlme and lme4, are given by: R> vt <-varCompTest(nlm1.nlme, nlm0.nlme) R> summary(vt)

Variance components testing in mixed effects models Testing that: covariance matrix of R0 and lrc is equal to 0 against the alternative that: covariance matrix of R0 and lrc > 0 Using the saemix package and the option fim = "extract", we get the following results, which are consistent with those obtained with lme4 and nlme (note that the likelihood value is different since saemix uses an important sampling approach to estimate the likelihood when lme4 and nlme rely on a linearization of the likelihood):

R> vt <-varCompTest(nlm1.saemix, nlm0.saemix, pval.comp = "both", + fim = "extract") R> summary(vt)

Variance components testing in mixed effects models Testing that: covariance matrix of omega2.R0 and omega2.lrc is equal to 0 against the alternative that: covariance matrix of omega2.R0 and omega2.lrc > 0 Likelihood ratio test statistic: LRT = 1.983014

Limiting distribution: mixture of 3 chi-bar-square distributions with degrees of freedom 0 1 2 associated (exact) weights: 0.2805023 0.5000000 0.2194977 p-value of the test: from exact weights: 0.160974 bounds on p-value: lower 0.07953661 upper 0.2650452

Comparison using different FIM estimates

To evaluate the performance of the FIM estimation via parametric bootstrap implemented in varTestnlme, we compared the empirical coverages of the nominal 95% confidence intervals based on the asymptotic normality property of the MLE obtained using the FIM extracted from packages nlme, lme4 and saemix, and our estimate.


Ahmad Alawad

Laure Latapie

David Evrard

Pierre Gros

Georges Istamboulie

Thierry Noguer

Carole Calas-Blanchard

SECM for Studying the Immobilization and Repartition of a Redox Anti-tetracycline Aptamer on Screen-printed Carbon Electrodes

SECM for studying the immobilization and repartition of a redox anti-tetracycline aptamer on screen-printed carbon electrodes

Ahmad Alawad [a,b] , Laure Latapie [c] , David Evrard [c] , Pierre Gros [c] , Georges Istamboulie [a,b] , Thierry Noguer [a,b] and Carole Calas-Blanchard* [a,b] [a] Université de Perpignan Via Domitia, Biocapteurs-Analyse-Environnement, 66860, Perpignan, France [b] Laboratoire de Biodiversité et Biotechnologies Microbiennes, USR 3579 Sorbonne Universités (UPMC) Paris 6 et CNRS Observatoire Océanologique, 66650, Banyuls-sur-Mer, France [c] Laboratoire de Génie Chimique, Université de Toulouse, CNRS, INPT, UPS Toulouse, France Abstract: Scanning electrochemical microscopy (SECM) is discussed as a versatile tool to provide a unique approach of localized electrochemical information in the context of biosensing research. The step-by-step immobilization of DNA aptamer with intrinsic redox activity on screen-printed carbon electrode (SPCE) was successfully monitored using SECM imaging tool. Control experiments were performed with a non-electroactive aptamer having a similar sequence length with the same proportion of A, T, G, and C bases. After immobilization of these aptamers on SPCE, SECM images showed the repartition of the electroactive anti-tetracycline aptamer when comparing with images produced for control and for all modification steps of SPCE. The possibility of tetracycline detection was also effectively proved by causing a decrease in recorded currents with SECM.

Keywords: Scanning electrochemical microscopy, electrochemical aptasensor, electroactive aptamer, tetracycline Biosensors technology began when analytical scientists faced the challenge in developing novel methodologies for the detection and quantification of biomedical analytes and simple/small molecules [1,2]. Over the past decade, many important technological advances have provided tools and materials required for the construction of various types of biosensor devices [3]. However, the selectivity and the sensitivity of a biosensor are primarily conditioned by its sensing layer [4]. Many different types of biological recognition elements (enzymes, antibodies, DNA/RNA…) have been employed for highly selective responses to particular analytes [4,5]. Among affinity elements, aptamers are DNA or RNA oligonucleotides of particular interest since they fold into three-dimensional conformations capable of binding strongly and selectively to target molecules [6]. Biosensors in which aptamers substituted the common biological recognition elements are called aptasensors. Among various aptasensors reported in literature, electrochemical aptasensors have been demonstrated as rapid, sensitive, selective and lowcost devices [7,8].

We have recently described the development of a label-free electrochemical aptasensor for the reagentless determination of tetracycline (TET) antibiotic in water. The sensing layer consisted in a specific 76-mer oligonucleotide [9] that was immobilized by carbodiimide coupling on a screen-printed carbon electrode (SPCE). Cyclic voltammetry experiments were carried out to characterize the sensor and revealed a very unusual and unexpected quasi-reversible signal, with oxidation and reduction peaks at +0.19 and -0.25 V vs silver / silver chloride (Ag/AgCl) reference electrode, respectively. Further experiments showed a specific decrease of this signal in presence of TET, demonstrating its proper binding by the aptamer. This remarkable behavior was exploited to detect TET at concentrations ranging from 0.11 to 45 nM, with a calculated limit of detection of 0.078 nM [10].

In the present work, we used scanning electrochemical microscopy (SECM) as an imaging technique to confirm the proper immobilization of intrinsically active anti-TET aptamers and to evaluate their repartition on the surface of a SPCE-based sensor. SECM is based on the current recorded at an ultramicroelectrode (UME) which is positioned close to a substrate to obtain its topography or local reactivity [11,12]. SECM has already shown promising issues for the study of grafted organic layers [13], characterization of bio-receptors immobilization on modified electrode surfaces [14] and as useful technique for mapping of surface-confined DNA and visualization of DNA hybridization thanks to an electroactive hybridization indicator [15]. Several studies confirmed that SECM is a suitable technique for imaging DNA-biosensors surface thanks to its good resolution and fast scanning rate [16][17][18]. SECM was used in this work in substrate generation/tip collection mode (SG/TC), where an electroactive species generated at the substrate surface (SPCE) is detected at the UME tip. As shown in Figure 1-a, the UME tip was fixed at a distance of 10 µm from the SPCE surface. The UME was moved along the x axis line at 10 µm/s and the current was recorded point by point. The measurements were carried out in a buffer used as electrolyte without adding redox species in solution. SECM imaging was performed after each step required for aptamer immobilization in a scanning area of 2000 x 2000 μm 2 , with a step of 40 μm (Figure 1-b). After immobilization of anti-TET aptamer by carbodiimide coupling, a new SECM experiment was performed. The resulting map shows a strong increase of the current up to a calculated average value of 50 pA (Figure 2-c) with an electrochemical reactivity less uniform than previously, where green and red spots corresponding to currents ranging from 100 to 280 pA reveal the presence of aptamer molecules on the surface. These results demonstrate that the studied aptamer can be detected on SPCE surface by SECM in SG/TC mode. Proper immobilization of aptamer was confirmed, allowing the global repartition of the aptamer on the electrode surface to be evaluated. In order to confirm this result, a similar sensor was elaborated with a random aptamer having the same size and nucleotide ratio than anti-TET aptamer. Figure 3 makes comparisons between the signals obtained with the anti-TET aptamer (Figure 3-a) and with the random one (Figure 3-b). The average current obtained in the presence of the anti-TET aptamer was 50 pA, while the average signal for random aptamer was 1 pA. Such a signal was almost the same as that observed before immobilization (0 pA) (Figure 2-b). These results confirmed that green and red spots in figure 2-c and 3-a were due to the proper immobilization of the electroactive anti-TET aptamer on the SPCE surface.

A previous study using the TET aptasensor demonstrated that the aptamer electrochemical response decreased in presence of TET [10], probably due to a conformational change upon target binding. In order to confirm this observation, we evaluated the effect of TET binding using SECM by scanning the sensing area after one hour of incubation with 45 nM of TET. As shown in figure 2-d, the binding of TET induced a decrease in the average current from 50 to 2 pA, thus completely shutting off the aptamer electrochemical signal. This observation confirms the previous results obtained by cyclic voltammetry displaying that the sensor operates as a signal-off device [10]. As discussed herein, we used SECM to characterize the anti-TET aptamer immobilization on a SPCE by analyzing 3D images (x-y-current). SECM surface imagery technique allowed us to evaluate the distribution of this electroactive aptamer on a functionalized electrode by observing the apparition of intense spots on the scanned area. These spots correspond to an increase in the current recorded at UME tip, generated thanks to the aptamer oxidation in absence of additional redox probes. These results were validated by a similar immobilization of a non-electroactive aptamer, having the same nucleotides as the previous one, but distributed randomly [10]. Based on the results obtained after incubation with tetracycline, SECM imagery could be considered as a possible technique for detecting this antibiotic by the anti-TET aptamer immobilized on the electrode surface. These preliminary results pave the way for the use of SECM as a detection technique for electrochemical aptasensors.

Experimental

The following chemicals were purchased from Sigma-Aldrich (France): sodium phosphate dibasic (Na2HPO4), potassium phosphate monobasic (KH2PO4), sodium chloride (NaCl), magnesium chloride (MgCl2), potassium chloride (KCl), calcium chloride (CaCl2), sulfuric acid (98%), sodium nitrite (NaNO2), N-hydroxysuccinimide (NHS), N-(3-dimethyl amino propyl)-N ′ -ethylcarbodiimide hydrochloride (EDC), 2-(N-morpholino)ethanesulfonic acid (MES), 4-amino benzoic acid (4-ABA), ethanolamine, Trizma ® hydrochloride and Tween 20. All solutions were prepared in deionized Milli-Q water (Millipore, Bedford, MA, USA). Anti-TET aptamer was a 76-mer oligonucleotide with the following sequence:5'-NH2-CGT ACG GAA TTC GCT AGC CCC CCG GCA GGC CAC GGC TTG GGT TGG TCC CAC TGC GCG TGG ATC CGA GCT CCA CGT G-3', as described in [9]. A 76-mer aptamer (random aptamer) which has the same proportion of each base than the anti-TET aptamer was also tested to investigate the redox signal: 5'-NH2-AAC GTT TCA GCG CAG GCG CTC ACG GTT GCG GTT GGT CCG CGT GGA TCC GAG CTC CAC GTC CGC ACG GCC GCA CTG C-3'. The oligonucleotides were purchased from Microsynth (Switzerland). Aptamer solutions were prepared in a 10 mM phosphate buffer pH7.4 containing KCl (10 mM) and MgCl2 (5 SPCE electrodes were fabricated using a DEK248 screen-printing system as previously described [10]. Prior to modification, the SPCE were electrochemically pre-treated by applying five cyclic voltammetric scans from 1.0 to -1.5 V vs Ag/AgCl in a sulfuric acid solution (H2SO4 0.5 M) containing 0.1 M KCl. A diazonium salt was synthesized in situ by mixing 2 mM NaNO2 with 2 mM of 4-ABA in 1 mL of 0.5 M HCl [8,10].

The mixture was left to react for 5 min at room temperature, and electrochemical generation of surface 4-carboxyphenyl (4-CP) was performed by linear sweep voltammetry, with potential scanning from 0.6 to -0.8 V vs Ag/AgCl [8,10]. The carboxylic groups of 4-CP-modified SPCE were then activated by incubating the working electrode for 60 min in 5 mL of 100 mM MES buffer containing 100 mM of EDC and 25mM of NHS. After washing the electrode surface with deionized water, 20 μL of 1μM aptamer solution were drop-casted onto the activated surface and allowed to react for 4 h under humid atmosphere. The electrodes were then rinsed with distilled water to remove unbound aptamers. Modified SPCE were finally incubated with 50 μL of 1M ethanolamine solution for 3 hours in order to block unreacted succinimide groups. The buffer used as electrolyte in all measurements is a binding buffer for aptamer/TET recognition, based on Tris-HCl (20 mM, pH 7.6) containing NaCl (100 mM), MgCl2 (2 mM), KCl (5 mM), CaCl2(1 mM) and Tween 20 (0.02 %).

All SECM measurements were performed with a 25 μm diameter gold UME electrode and using a Uniscan M370 bipotentiostat. A three-electrode setup was employed consisting in the gold UME as working electrode, carbon and Ag/AgCl from (SPE) electrodes as counter and reference electrodes respectively. The UME was polarized at 0.5 V vs Ag/AgCl, corresponding to the oxidation potential of anti-TET aptamer, while the working SPCE was polarized at -0.3 V vs Ag/AgCl, to ensure that the immobilized aptamer was present in its reduced form on the substrate. Preliminary experiments based on approach curves in a feedback mode allowed the determination of an optimal UME-substrate distance of 10 µm. For surface electrochemical characterization, a Substrate Generation / Tip Collection mode (SG/TC) was used.

Fig. 1 :

 1 Fig.1: a) Scheme of aptasensor preparation and probing mechanism, b) Scanning grid realized on SPCE working electrode

Figure 2 -

 2 Figure2shows the different images resulting from SECM measurements over the aptamer immobilization steps. For each scan, the average current was calculated from the digital data excluding those of the first measurement line, strongly disturbed. Figure2-a corresponds to the image obtained with a bare SPCE electrode immersed in binding buffer showing a homogenous repartition of current over the whole scanned area, with an average value of 15 pA. Figure2-bshows the mapping of the electrode surface after its modification by 4-carboxyphenyl layer (4-CP) following the reduction of the diazonium salt. This activation step led to a decrease of the current down to values close to 0 pA over the scanning area, resulting from the passivation of the carbon surface as already observed by electrochemical impedance spectroscopy experiments[8,10].

Fig. 2 .Fig. 3 .

 23 Fig.2. SECM images realized in binding buffer of: a) bare SPCE, b) 4-CP-modified SPCE, c) aptasensor for tetracycline detection, d) aptasensor response after incubation with 45nM of tetracycline. The scanned surface was 2000 x 2000 μm 2 with step size x =y=40 μm/point, z = 10 μm; EUME = 0.5 V vs Ag/AgCl; ESPCE = -0.3 V vs Ag/AgCl

Acknowledgements

Ahmad Alawad sincerely thanks Erasmus Mundus Peace 2 Programme and Programme national d'aide à l'Accueil en Urgence des Scientifiques en Exil (PAUSE, Collège de France) for their kind help in financing his doctoral work.