dataesr/scientific-paragraphs-categorization huggingface dataset:
60 GHz processing speed and 16 GB memory. In order to implement the different algorithms described in this work, we used the following MATLAB functions: expm: it allows to calculate the exponential of a square matrix. This function is based on a scaling and squaring algorithm with a Padé approximation  [24]. lyap: it allows to solve Sylvester or Lyapunov matrix equations. For our purposes, the instruction lyap(A,B,-C) delivers the matrix X solution of the algebraic Sylvester equation A X + X B = C. integral: it allows to calculate numerically an integral, using the arguments "ArrayValued" and "true". Furthermore, we precise that when the constant solution or integral formula methods are combined with the block Arnoldi process to obtain an approximate solution to the differential equation, the iterations are stopped as soon as the dimension of the Krylov subspace generated by the block Arnoldi process reaches a maximum value m = M max = 110 or as soon as the maximal norm r max computed by the algorithm is lower than 10 -10 µ where µ = ∥A∥ + ∥B∥ + ∥C∥ in the full rank case and µ = ∥A∥ + ∥B∥ + ∥E∥ ∥F∥ in the low rank case.
Experiment 1 In this first example, the numerical tests are done with moderate size matrices A and B. We compare the solution provided by our proposed constant solution method implemented via Algorithm 1 with the one obtained using the integral formula (5) as well as with the solutions given by some classical ODE's solvers from Matlab. The solvers ode15s, ode23s, ode23t and ode23tb are usually used for stiff ODE's, while the other solvers ode45, ode23 and ode113 are used for non stiff ODE's. Note that since some ODE solvers behave similarly and in order not to overload the plots, we only give the results obtained with the four methods ode15s, ode23s, ode23tb and ode45. In the following two experiments, we consider the time intervals [0, T ] with T is either T = 1 with the number of nodes is N = 10 or T = 10 with the number of nodes is N = 50 which means that the step time is δ T = 0.
1 when T = 1 while δ T = 0. 2 when T = 10. Here, we consider the matrices A 0 =gallery('leslie',n 0 ), B 0 = gallery('minij',s 0 ) with n 0 = 50 and s 0 = 10 and the coefficient matrices A, B of the differential Sylvester equations are generated by  (36) , as explained in the benchmark example. The parameters α, β are equal to -2 and -1 respectively. As the matrices K, R are those given at the beginning of section 4, the size of the matrices A, B are now n = 150 and s = 30 respectively. Here, we point out that the solution computed by Algorithm 1 and those computed by the Algorithm based on integral formula or issued by the Matlab ODE solvers are compared to the exact one given in  (37)  which is considered as the reference solution X ref. Thus, in the plots, we represent the behaviour of the norm of the relative error t k → ∥X k -X ref k ∥ ∥X ref k ∥ as a function of t k where t k = k δ T.
This gives A ∈ R n×n , B ∈ R s×s with n = n 2 0 and s = s 2 0. We examine the performances of CS-BA and IF-BA for four choices of n 0 and s 0 which are (n 0 , s 0 ) = (30, 3), (n 0 , s 0 ) = (50, 3), (n 0 , s 0 ) = (30, 5) and (n 0 , s 0 ) = (50, 5). The considered time intervals are [0, T ] where T = 1 and N = 10 or T = 2 and N = 20. This means that the step time is always δ T = 0. 1. In Table  2 , we reported the time ratio (TR) and the relative difference norm (RDN) between the CPU-time of CS-BA and IF-BA. Experiment 2. 2. In this test, we took A 0 =gallery('hanowa',1500,-5) and B 0 = gallery('leslie',6) from the Matlab gallery and transform them into A et B of sizes n = 4500 and s = 18 respectively by using  (36)  in which we took α = -7 and β = -5. The obtained results for different time intervals which are summarized in Table  3  include the time ratio TR and the relative error norms REN CS-BA , REN IF-BA between the approximate solutions X CS-BA , X IF-BA given by CS-BA and IF-BA respectively and the X Exact the exact solution computed by  (37).
Declarations • Ethical Approval. "Not Applicable". • Availability of supporting data. "No data sets were generated or analyzed during the current work" • Competing interests. "The authors declare no competing interests. " • Funding "No funding was provided to support this work • Authors' contributions. "The three authors contributed equally to this work" • Acknowledgments. We thank the anonymous reviewers and the editor for any criticism, comments, suggestions, or remarks that will help us improve this manuscript.
Acknowledgements Part of this paper was elaborated during a research period of LB at the  Université Savoie Mont Blanc. The  Université Savoie Mont Blanc  and its facilities are kindly acknowledged. LB is member of the  Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA)  of the  Istituto Nazionale di Alta Matematica (INdAM). Data sharing not applicable to this article as no datasets were generated or analysed during the current study.
The issue of variance components testing arises naturally when building mixed-effects models, to decide which effects should be modeled as fixed or random. 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 any subset of the variances and covariances corresponding to any 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. Theoretical properties of the used likelihood ratio test are recalled, numerical methods used to implement the test procedure are detailed and examples based on different real datasets using different mixed models are provided.
Introduction Mixed-effects models are widely used in many fields of applications, ranging from agronomy, social science, medical science, or 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 softwares are available to fit these types of model, for example the SAS procedures MIXED and NLMIXED for linear, generalized linear and nonlinear mixed models, or the dedicated softwares NONMEN and Monolix for nonlinear mixed models, with an emphasis on pharmacodynamic and pharmacokinetic models for the former. In R, three main packages are available to deal with these models: lme4  (Bates, Mächler, Bolker, and Walker 2015) , nlme  (Pinheiro, Bates, DebRoy, Sarkar, and R Core Team 2019)  and saemix  (Comets, Lavenu, and Lavielle 2017).
g. , in R. Among others, we can cite for example the glmmLasso  (Groll 2017 ) package, computing an 1penalized maximum likelihood estimator in the context of generalized linear mixed models, or the cAIC4 (Saefken and Ruegamer 2018) 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  (Kuznetsova, Brockhoff, and Christensen 2017)  provides enhanced versions of functions anova() and summary() of the lme4 package by computing pvalues for the corresponding 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.
It also provides a step() function, which simplify the random effects structure of the model 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 chi-square one, 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  (Halekoh and Højsgaard 2014)  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 implements the tests proposed by  Crainiceanu and Ruppert (2004) , 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 and from function gamm of package mgcv. 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 all types of mixed-effects models, 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 dedicated to test that the variances of any subset of the random effects are equal to zero in LMMs, GLMMs and NLMMs, 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 returns the p-value based on the asymptotic distribution of the LRT statistic associated with the two nested models corresponding to the null and alternative hypotheses. The theory behind this test has been established in  Baey et al. (2019). The purpose of the current paper is to provide hints for practical implementation and a complete description of the different numerical methods used for the computation of each quantity involved.
This approach is more time consuming but more accurate, in particular for small samples. 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 (ϕ b, * i , y b, * i ), 1 ≤ i ≤ n, 1 ≤ b ≤ 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 θ * b. The bootstrap estimate of the FIM is then given by the empirical covariance matrix between the bootstrap estimates: Î * = 1 B B b=1 θ * b ( θ * b ) -θ * ( θ * ) , ( 5 ) where θ * = 1 B B b=1 θ * b. The second option consists in using the FIM estimate as it is computed by the package that was used to fit the models, when possible.
Finally, the third option is to provide a user-defined positive semidefinite matrix which should be an estimate of the FIM. It can be for example computed by using the R package PFIM 4. 0  (Dumont, Lestini, Le Nagard, Mentré, Comets, Nguyen et al. (2018) ) based on Monte Carlo tools or by implementing the method presented in  (Delattre and Kuhn (2019) ) based on stochastic approximation tools.
Computation of the tangent cone In  Baey et al. (2019) , it was highlighted that the limiting distribution of the likelihood ratio test statistic depends on the structure of Γ and its expression is exhibited in two common cases, when Γ is 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       , ( 6 ) 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 = p. 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). It the appendix of this book, it is mentioned that more general structures for the covariance matrix Γ could be implemented using modularized functions of the lme4 package.
Nevertheless these functionalities are not treated in the varTestnlme package at the moment. 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 block-diagonal case with full blocks. This restriction allows for a simple block-diagonal structure of the covariance matrix and eases the computation. In this context of a block-diagonal covariance matrix Γ, three cases can arise when testing that a subset of the 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 into which covariances are tested equal to 0 without testing that the corresponding variances are equal to 0 (only testing non-diagonal 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.
Orthodontal data The first dataset comes from a study on dental growth  (Potthoff and Roy 1964) , where the distance between the pituitary and the 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 , ( 18 ) ε ij ∼ N (0, σ 2 ), (b i1 , b i2 ) ∼ N (0, Γ). ( 19 ) 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)  (Lesnoff, Laval, Bonnet, Abdicho, Workalemahu, Kifle, Peyraud, Lancelot, and Thiaucourt 2004) , where the number of new serological cases occuring during a given time period was recorded at 4 occasions, on 15 herds (see Figure  2 ). 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  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 R> data("cbpp", package = "nlme") This data can be fitted using a generalized linear mixed-effects model and a logistic regression.
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 the i-th tree at age x j , by Asym i the asymptote for the i-th tree, lrc i the logarithm of the growth rate for the i-th tree and R 0,i the height of tree i at age 0. We consider the following model: y ij = Asym i + (R 0,i -Asym i ) exp(-e lrc i x j ) + ε ij , ε ij ∼ N (0, σ 2 ) (22) (Asym i , R 0,i , lrc i ) = β + b i , b i ∼ N 3 (0, Γ) (23).
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 one of the following three packages: nlme, lme4 and 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 (18). 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: 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. lme4 <-lmer(distance ~1 + Sex + age + age * Sex + + (1 + age || Subject), data = Orthodont, REML = FALSE) R> lm2.
Nonlinear model Let us consider the nonlinear model defined in (  22 ). 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    , σ > 0} and Θ = {β ∈ R 3 , Γ =    γ 2 1 0 0 0 γ 2 2 0 0 0 γ 2 3    , σ > 0}, i. e. , that only the asymptote is random. The corresponding code using nlme and lme4 is: Using saemix, the syntax is given by: R> modelSSasymp <-function(psi, id, xidep){ Asym <-psi[id,1] R0 <-psi[id,2] lrc <-psi[id,3] age <-xidep ypred <-Asym + (R0 -Asym) * exp( -exp(lrc) * age) return(ypred)} R> nlm1. saemix <-saemixModel(model = modelSSasymp, description = "Asymptotic regression", psi0 = matrix(start, ncol = 3, byrow = TRUE, dimnames = list(NULL, c("Asym", "R0", "lrc"))), transform. par = c(0,0,0), fixed.
Variance component testing Once the models have been fitted using one of the three packages lme4, nlme or saemix, the test of the variance components can be run using the varCompTest() function of the varTestnlme package. varCompTest(m1, m0, control = list(M=5000, parallel=F, nbcores=1, B=1000), fim = "extract", pval. comp = "bounds") 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 size of the Monte Carlo sample for the computation of the weights of the chi-bar-square distribution (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 Fisher Information Matrix should be extracted from the fitted model m1; "compute" if it should be computed using parametric bootstrap ; or FIM a user-defined matrix to be used as the Fisher Information Matrix, • pval.
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 any subset of random effects is null, taking into account correlation between them. 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. Additional tools are provided to approximate the Fisher Information Matrix when it is needed to compute the limiting distribution of the likelihood ratio test statistic.
Acknowledgments The results of this work and its companion [1] were first announced in 2016. We warmly thank the institutions which have welcomed us during these years: the institutions of  Rio de Janeiro PUC ,  IMPA ,  UFF ,  UFRJ ,  UERJ , the  Universidad de la República of Montevideo , and the  Institut de Mathématiques de Bourgogne of Dijon. P. G. B. was supported by  Ministerio de Educación, Cultura y Deporte  (Grant No. MTM2017-87697-P ) and  Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). D. F. was partially supported by the  RFBR  project  20-01-00420. V. K. was partially supported by the  Réseau France-Brésil en Mathématiques , by the project  ANR Gromeov  ( ANR-19-CE40-0007 ), and by  Centre Henri Lebesgue  ( ANR-11-LABX-0020-01 ). C. M. was partially supported by  CAPES postdoc program ,  Brazil  ( INCTMat  &  PNPD 2015-2016 ), the  Ministerio de Economía y Competitividad (Spain) , Grant  MTM2014-56950-P , Spain ( 2015-2017 ), the  Programa Cientista do Estado do Rio de Janeiro, FAPERJ, Brazil  ( 2015-2018 ),  CNPq  research grant  310915/2019-8  (Brazil) and  MathAmSud:  " Rigidity and Geometric Structures in Dynamics "  2019-2020  ( CAPES-Brazil ).
He also received support from the  CSIC , via the  Programa de Movilidad e Intercambios Académicos. M. T. was partially supported by  PEPS -Jeunes Chercheur -e-s -2017  and  2019  ( CNRS ), the  Réseau France-Brésil en Mathématiques ,  MathAmSud RGSD (Rigidity and Geometric Structures in Dynamics)   19-MATH-04 , the project  ANER Agroupes  ( AAP 2019 Région Bourgogne-Franche-Comté ), and by the project  ANR Gromeov  ( ANR-19-CE40-0007 ). His host department IMB receives support from the  EIPHI Graduate School  ( ANR-17-EURE-0002 ).
Numerical results Before presenting numerical results, let us comment briefly on some practical details concerning the effective implementation. Our code in based on Matlab. The resolution of the nonlinear system (  18 )-(  19 ), in its compact form (39) is achieved thanks to Newton's method: (69) J(ρ n, )δρ n, = -H n (ρ n, ), ρ n, +1 = ρ n, + δρ n, , with J standing for the Jacobian matrix of H n. Note that ρ n σ , σ ∈ E ext is not considered as an unknown and is deduced from the cell values thanks to  (20). We initialize (69) by setting ρ n,0 = ρ n-1 and than iterate until δρ n, ∞ / ρ n, +1 ∞ ≤ 10 -12. Then we set ρ n = ρ n, +1. 5. 1. Numerical evidence of the convergence. The first numerical test we propose aims at confirming our intuition concerning the second order accuracy in space of the scheme sketched in Section 2.
Acknowledgements. This project has received funding from the  European Union 's  Horizon 2020 research and innovation programme under grant agreement No 847593 (EURAD program ,  WP DONUT) , and was further supported by  Labex CEMPI  ( ANR-11-LABX-0007-01 ) and the  Fédération de Recherche Mathématique des Hauts-de-France  ( FR CNRS 2037 ). C. Cancès also acknowledges support from the  COMODO  ( ANR-19-CE46-0002 ) and  MICMOV  ( ANR-19-CE40-0012 ) projects. J. Venel  warmly thanks the  Inria research center of the University of Lille  for its hospitality and support, and the authors further thank  Claire Chainais  for stimulating discussions.
Study design There are three major steps involved in this study (Figure  1 ). First, we collected data and methods for calculating reproduction numbers. Then, we described the epidemiological characteristics of the outbreaks in four regions, H, P, X, and Z. Finally, we obtained reproduction numbers of COVID-19 in those regions through three different methods, which are the definition method, the method based on the next-generation matrix, and method based on SI and epidemic curve, to evaluate transmissibility of the diseases under various circumstances.