Ana Arribas-Gil

Karine Bertin
email: karine.bertin@uv.cl

Cristian Meza
email: cristian.meza@uv.cl

Vincent Rivoirard
email: vincent.rivoirard@dauphine.fr

LASSO-type estimators for semiparametric nonlinear mixed-effects models estimation

Keywords: LASSO, Nonlinear mixed-effects model, On-line auction, SAEM algorithm, Semiparametric estimation

published or not. The documents may come

Introduction

We consider the semiparametric nonlinear mixed effects model (SNMM) as defined by [START_REF] Ke | Semiparametric nonlinear mixedeffects models and their applications (with discussion)[END_REF] in which we have n individuals and we observe:

y i j = g(x i j , φ i , f ) + ε i j , ε i j ∼ N (0, σ 2 ) i.i.d., (1)
i = 1, . . . , N, j = 1, . . . , n i , where y i j ∈ R is the jth observation in the ith individual, x i j ∈ R d is a known regression variable, g is a common known function governing within-individual behaviour and f is an unknown nonparametric function to be estimated. The random effects φ i ∈ R p satisfy

φ i = A i β + η i , η i ∼ N (0,Γ ) i.i.d.
where A i ∈ M p,q are known design matrices, β ∈ R q is the unknown vector of fixed effects and we suppose that ε i j and η i are mutually independent. The parameter of the model is (θ , f ), where θ = (β ,Γ , σ 2 ) belongs to a finite dimensional space whereas f belongs to an infinite dimensional space of functions denoted H . [START_REF] Ke | Semiparametric nonlinear mixedeffects models and their applications (with discussion)[END_REF] consider the most common type of SNMM in practice, in which g is linear in f conditionally on φ i , g(x i j , φ i , f ) = a(φ i ; x i j ) + b(φ i ; x i j ) f (c(φ i ; x i j )),

(2)

where a, b and c are known functions which may depend on i.

Different formulations of SNMM's have been recently used to model HIV dynamics [START_REF] Wu | The study of longterm HIV dynamics using semi-parametric non-linear mixed-effects models[END_REF]Liu andWu, 2007, 2008), time course microarray gene expression data [START_REF] Luan | Model-based methods for identifying periodically expressed genes based on time course microarray gene expression data[END_REF], circadian rhythms [START_REF] Wang | A flexible model for human circadian rhythms[END_REF][START_REF] Wang | Shape-invariant modeling of circadian rhythms with random effects and smoothing spline anova decompositions[END_REF], as in the following example, or to fit pharmacokinetic and pharmacodynamic models [START_REF] Wang | Semiparametric mixed-effects analysis of PKPD models using differential equations[END_REF], among many other applications.

Example 1 The following model was proposed by [START_REF] Wang | A flexible model for human circadian rhythms[END_REF] to fit human circadian rhythms:

y i j = µ + η 1i + exp(η 2i ) f x i j - exp(η 3i ) 1 + exp(η 3i ) + ε i j , ε i j ∼ N (0, σ 2 ) i.i.d.
η i ∼ N (0,Γ ) i.i.d.

for i = 1, . . . , N, j = 1, . . . , n i , where y i j is the physiological response of individual i at the jth time point x i j . This model can be written in the general form (1) as: y i j = g(x i j , φ i , f ) + ε i j , ε i j ∼ N (0, σ 2 ) i.i.d.,

g(x i j , φ i , f ) = φ 1i + exp(φ 2i ) f x i j -exp(φ 3i ) 1 + exp(φ 3i ) φ i = (1, 0, 0) µ + η i , η i ∼ N (0,Γ ) i.i.d.

where φ i = (φ 1i , φ 2i , φ 3i ) and η i = (η 1i , η 2i , η 3i ) . In this example f represents the common shape of the observed curves, and φ 1i , exp(φ 2i ), and exp(φ 3i )/(1 + exp(φ 3i )) stand for the individual vertical shift, individual amplitude and individual horizontal shift respectively. Here d = 1, p = 3, q = 1 and the parameter of the model is (µ,Γ , σ 2 , f ). This model was also used by [START_REF] Ke | Semiparametric nonlinear mixedeffects models and their applications (with discussion)[END_REF] for modeling Canadian temperatures at different weather stations.

Let us introduce the following vectorial notation: y i = (y i1 , . . . , y in i ) , y = (y 1 , . . . , y N ) , φ = (φ 1 , . . . , φ N ) , η = (η 1 , . . . , η N ) , g i (φ i , f ) = (g(x i1 , φ i , f ), . . . , g(x in i , φ i , f )) , g(φ , f ) = (g 1 (φ 1 , f ) , . . . , g N (φ n , f ) ) , A = (A 1 , . . . , A N ) , Γ = diag(Γ , . . . ,Γ ) and n = ∑ N i=1 n i . Then, model (1) can be written as:

y|φ ∼ N (g(φ , f ), σ 2 I n ) (3) φ ∼ N (Aβ , Γ )
where I n is the identity matrix of dimension n, and the likelihood of observations y is:

p(y; (θ , f )) = p(y|φ ; (θ , f ))p(φ ; (θ , f ))dφ = 1 (2πσ 2 ) n 2 exp -1 2σ 2 y -g(φ , f ) 2 × 1 (2π) N p 2 |Γ | N 2 exp -1 2 Γ -1/2 (φ -Aβ ) 2 dφ = 1 (2π) n+N p 2 (σ 2 ) n 2 |Γ | N 2 × exp -1 2 1 σ 2 y -g(φ , f ) 2 + Γ -1/2 (φ -Aβ ) 2 dφ , ( 4
)
where • is the L 2 norm. In their seminal paper, Ke and Wang consider a penalized maximum likelihood approach for the estimation of (θ , f ). That is, they propose to solve max

θ , f { (y; (θ , f )) -nλ J( f )} (5)
where (y; (θ , f )) is the marginal log-likelihood, J( f ) is some roughness penalty and λ is a smoothing parameter. Moreover, they assume that f belongs to some reproducing kernel Hilbert space (RKHS) H = H 1 ⊕ H 2 , where H 1 is a finite dimensional space of functions, H 1 = span{ψ 1 , . . . , ψ M }, and H 2 is a RKHS itself. Since the nonlinear function f interacts in a complicated way with the random effects and the integral in (4) is intractable, they replace (y; (θ , f )) by a first-order linearization of the likelihood with respect to the random effects. Then, they propose to estimate (θ , f ) by iterating the following two steps: i) given an estimate of f , get estimates of θ and φ by fitting the resultant nonlinear mixed model by linearizing the log-likelihood (replacing by ˜ ). In practice they use the R-function nlme [START_REF] Pinheiro | Mixed-Effects Models in S and S-PLUS[END_REF] to solve this step. ii) given an estimate of θ , θ , estimate f as the solution to max

f ∈H ˜ (y; ( θ , f , φ )) -nλ J( f ) .
Since in ii) the approximated log-likelihood involves a bounded linear functional, the maximizer in H of ˜ (y; ( θ , f , φ ))-Nλ J( f ) given θ and φ belongs to a finite dimensional space and it is estimated as a linear combination of functions from H 1 and H 2 . Conceptually, the whole approach is equivalent to solving (5) not on H but on a finite-dimensional approximation space of H at each iteration. As it is discussed in that article, despite of the lack of an exact solution, the spline smoothing method provides good results and its use in this framework is largely justified. However, the method relies on prior knowledge of the nonlinear function f and provides better results when this kind of information is available.

In practice, the Ke and Wang's method is implemented in the R package assist [START_REF] Wang | Assist: A suite of s functions implementing spline smoothing techniques[END_REF]) and in particular in the snm function which is directly related with the nlme function. As for the parametric estimation, it is important to point out some drawbacks of the approximated methods based on linearization of the log-likelihood, such as the first-order linearization conditional estimates (FOCE) algorithm used in the snm funtion [START_REF] Wang | Assist: A suite of s functions implementing spline smoothing techniques[END_REF]). It has been shown that they can produce inconsistent estimates of the fixed effects, in particular when the number of measurements per subject is not large enough [START_REF] Ramos | Estimation of nonlinear random coefficient models[END_REF]; [START_REF] Vonesh | A note on the use of Laplace's approximation for nonlinear mixed-effects models[END_REF]; [START_REF] Ge | An approximate likelihood approach to nonlinear mixed effects models via spline approximation[END_REF]). Furthermore, simulation studies have shown unexpected increases in the type I error of the likelihood ratio and Wald tests based on these linearization methods [START_REF] Ding | Assessing antiviral potency of anti-HIV therapies in vivo by comparing viral decay rates in viral dynamic models[END_REF]). In addition, from of statistical point of view, the theoretical basis of this linearizationbased method is weak.

Since estimation in SNMMs is an important problem and a difficult task from which many challenging aspects arise, in this paper we propose an alternative estimation procedure to tackle some of these points. On the one hand, for the parametric step we will focus on the maximization of the exact likelihood. We propose to use a stochastic version of the EM algorithm, the so-called SAEM algorithm introduced by [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF] and extended by [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] for nonlinear mixed models, to estimate θ without any approximation or linearization. This stochastic EM algorithm replaces the usual E step of EM algorithm [START_REF] Dempster | Maximumlikelihood from incomplete data via the EM algorithm[END_REF] by a simulation step and a stochastic procedure, and converges to a local maximum of the likelihood. The SAEM has been proved to be computationally much more efficient than other stochastic algorithms as for example the classical Monte Carlo EM (MCEM) algorithm [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithm[END_REF] thanks to a recycling of the simulated variables from one iteration to the next (see [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF]). Indeed, previous attempts to perform exact ML estimation in SNMMs have been discarded because of the computational problems related to the use of an MCEM algorithm (see [START_REF] Liu | Simultaneous inference for semiparametric nonlinear mixed-effects models with covariate measurement errors and missing responses[END_REF], 2008, 2009)). Moreover we use a Restricted Maximum Likelihood (REML) version of the SAEM algorithm to correct bias estimation problems of the variance parameters following the same strategy as [START_REF] Meza | Estimation in the probit normal model for binary outcomes using the saem algorithm[END_REF]. On the other hand, for the nonparametric step we will propose a LASSO-type method for the estimation of f . The popular LASSO estimator (least absolute shrinkage and selection operator, [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF]) based on 1 penalized least squares, has been extended in the last years to nonparametric regression (see for instance Bickel et al (2009)).

It has been also used by [START_REF] Schelldorfer | Estimation for high-dimensional linear mixed-effects models using l1-penalization[END_REF] in highdimensional linear mixed-effects models. In the nonparametric context, the idea is to reconstruct a sparse approximation of f with linear combinations of elements of a given set of functions { f 1 , . . . , f M }, called dictionary. That is, we are implicitly assuming that f can be well approximated with a small number of those functions. In practice, for the nonparametric regression problem, the dictionary can be a collection of basis functions from different bases (splines with fixed knots, wavelets, Fourier, etc.). The difference between this approach and the smoothing splines, is that the selection of the approximation function space is done automatically and based on data among a large collection of possible spaces spanned by very different functions. This is particularly important in situations in which little knowledge about f is available. This approach allows us to construct a good approximation of the nonparametric function which is sparse thanks to the large dictionary. The sparsity of the approximation gives a model more interpretable and since few coefficients have to be estimated, this minimizes the estimation error. The LASSO algorithm allows to use the dictionary approach to select a sparse approximation, unlike to wavelet thresholding or 0 -penalization. Moreover the LASSO algorithm has a low computational cost since it is based on a convex penalty.

We can summarize our iterative estimation procedure as: i) given f , an estimate of f , get estimates of θ and φ by fitting the resulting NLME with the SAEM algorithm (using either ML or REML). ii) given estimates of θ and φ , solve the resulting nonparametric regression problem using a LASSO-type method.

The rest of the article is organized as follows. In Section 2.1 we describe the SAEM algorithm and its REML version in the framework of SNMMs. In Section 3 we propose a LASSO-type method for the estimation of f in the resulting nonparametric regression problem after estimation of θ and φ . Oracle inequalities and subset selection properties for the proposed estimator are provided in the Supplementary Material. In Section 4, we describe the algorithm that combines both procedures to perform joint estimation of (θ , f ) in the SNMM. Finally, in Section 5, we illustrate our method through a simulation study and the analysis of price dynamics in on-line auction data. We conclude the article in Section 6.

2 Estimation of the finite-dimensional parameters

SAEM estimation of θ and φ

In this subsection we consider that we have an estimate of f , f , obtained in the previous estimation step that does not change during the estimation of θ . Thus, we can proceed as if f was a known nonlinear function and we fall into the SAEM estimation of nonlinear mixed-effects model framework (see [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF]). In this setting, convergence of the algorithm to a local maximum of the likelihood is guaranteed. In fact, note that since the estimation of f is performed by solving a nonparametric regression problem with regression variables c( φ i ; x i j ), i = 1, . . . , N, j = 1, . . . , n i (see Section 3), it will depend on the estimated value of φ at the precedent iteration. Then, we will note fthe current estimated function.

The complete likelihood for model (1) is:

p(y, φ ; θ ) = p(y|φ ; θ )p(φ ; θ ) = 1 (2π) n+N p 2 (σ 2 ) n 2 |Γ | N 2 exp -1 2 1 σ 2 y -g(φ , f-) 2 + Γ -1/2 (φ -Aβ ) 2
where n = ∑ N i=1 n i . Then, the complete log-likelihood is:

log p(y, φ ; θ ) = - 1 2 C + n log σ 2 + N log |Γ | + 1 σ 2 y -g(φ , f-) 2 + N ∑ i=1 (φ i -A i β ) Γ -1 (φ i -A i β ) (6)
where C is a constant that does not depend on θ .

The distribution of the complete-data model belongs to the exponential family, that is log p(y, φ ; θ ) = -Ψ (θ ) + S(y, φ ), Φ(θ ) , where •, • stands for the scalar product and S(y, φ ) is the sufficient statistics. The EM algorithm in this framework would involve the computation of E[S(y, φ )|y; θ (k) ] in the E step, which in our case is intractable. The SAEM algorithm replaces, at each iteration, the step E by a simulation step (S) of the missing data (φ ) and an approximation step (A). Then, iteration k of the SAEM algorithm writes:

-S step: simulate m values of the random effects, φ (k+1,1) , . . . , φ (k+1,m) , from the conditional law p(•|y; θ (k) ). -A step: update s k+1 according to:

s k+1 = s k + χ k 1 m m ∑ l=1 S(y, φ (k+1,l) ) -s k .
-M step: update the value of θ :

θ (k+1) = arg max θ {-Ψ (θ ) + s k+1 , Φ(θ ) } (7)
where (s k ) k is initialized at s 0 and (χ k ) k is a decreasing sequence of positive numbers which accelerates the convergence [START_REF] Kuhn | Coupling a stochastic approximation version of EM with an MCMC procedure[END_REF]. The role of the sequence (χ k ) k is crucial in the SAEM algorithm since it performs a smoothing of the calculated likelihood values from one iteration to another. In practice, this smoothing parameter is defined as follows. During the first L iterations, χ k = 1, and from iteration (L + 1) the smoothing parameter starts to decrease in order to stabilize the estimates and provide a faster convergence towards the true ML estimates. For example, [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] recommend to take χ k = (k -L) -1 for k ≥ (L + 1). The choices of the total number of iterations, K, and of L are then crucial. In order to define these constants, following [START_REF] Jank | Implementing and diagnosing the stochastic approximation em algorithm[END_REF] and [START_REF] Meza | Reml estimation of variance parameters in nonlinear mixed effects models using the SAEM algorithm[END_REF], we may use a graphical approach based on the likelihood difference from one iteration to the next one and monitor SAEM by estimating its progress towards θ ML by using the property of increasing likelihood of the EM algorithm (see for more details [START_REF] Meza | Reml estimation of variance parameters in nonlinear mixed effects models using the SAEM algorithm[END_REF]). Then, the total number of iterations can be fixed and the smoothing step can be defined. However, it is important to note that this procedure implies to run the SAEM algorithm twice. Furthermore, as all EM-type algorithms, SAEM is sensitive to the choice of the initial values. From (6), the sufficient statistics for the complete model are given by k+1) is obtained in the maximization step as follows:

s 1,i,k+1 = s 1,i,k + χ k 1 m m ∑ l=1 φ (k+1,l) i -s 1,i,k , i = 1, . . . , N s 2,k+1 = s 2,k + χ k 1 m m ∑ l=1 N ∑ i=1 φ (k+1,l) i φ (k+1,l) i -s 2,k s 3,k+1 = s 3,k + χ k 1 m m ∑ l=1 y -g(φ (k+1,l) , f-) 2 -s 3,k . Now, θ (
β (k+1) = N ∑ i=1 A i Γ (k) -1 A i -1 N ∑ i=1 A i Γ (k) -1 s 1,i,k+1 Γ (k+1) = 1 N s 2,k+1 - N ∑ i=1 A i β (k+1) s 1,i,k+1 - N ∑ i=1 s 1,i,k+1 A i β (k+1) + N ∑ i=1 A i β (k+1) A i β (k+1) σ 2 (k+1) = s 3,k+1 n .
When the simulation step cannot be directly performed, [START_REF] Kuhn | Coupling a stochastic approximation version of EM with an MCMC procedure[END_REF] propose to combine this algorithm with a Markov Chain Monte Carlo (MCMC) procedure. Then, the simulation step becomes:

-S step: using φ (k,l) , draw φ (k+1,l) with transition probability Π θ (k) (•|φ (k,l) ), l = 1, . . . , m, that is, (φ (k+1,1) ), . . . , (φ (k+1,m) ) are m Markov chains with transition kernels Π θ (k) . In practice, these Markov chains are generated using a Hastings-Metropolis algorithm (see [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] for details).

With respect to the number of chains, the convergence of the whole algorithm to a local maximum of the likelihood is granted even for m = 1. Greater values of m can accelerate the convergence, but in practice m is always lower than 10. This is the main difference with the MCEM algorithm, in which very large samples of the random effects have to be generated to obtain convergence of the algorithm.

REML estimation of variance components

It is well known that the maximum likelihood estimator of variance components in mixed effects models can be biased downwards because it does not adjust for the loss of degrees of freedom caused by the estimation of the fixed effects. This is also true in the context of SNMMs as [START_REF] Ke | Semiparametric nonlinear mixedeffects models and their applications (with discussion)[END_REF] point out in their paper.

To overcome this problem we consider restricted maximum likelihood (REML) estimation. REML, as originally formulated by [START_REF] Patterson | Recovery of inter-block information when block sizes are unequal[END_REF] in the context of linear models, is a method that corrects this problem by maximizing the likelihood of a set of linear functions of the observed data that contain none of the fixed effects of the model. But this formulation does not directly extend beyond linear models, where in general it is not possible to construct linear functions of the observed data that do not contain any of the fixed effects. However, in the case of nonlinear models, other alternative formulations of REML have been proposed. Here, we will consider the approach of [START_REF] Harville | Bayesian inference for variance components using only error contrasts[END_REF], that consists in the maximization of the likelihood after integrating out the fixed effects. To perform this integration we follow [START_REF] Foulley | Heterogeneous variances in gaussian linear mixed models[END_REF] and consider the fixed effects as random with a flat prior. The combination of this REML approach with the SAEM algorithm in the context of nonlinear mixed effects models has been studied recently by [START_REF] Meza | Estimation in the probit normal model for binary outcomes using the saem algorithm[END_REF]. The authors showed the efficiency of the method against purely ML estimation performed by SAEM and against REML estimation based on likelihood approximation methods. Following these ideas we note z = (φ , β ) the random effects and θ = (Γ , σ 2 ) the new parameter of the model. As in the general case, the simulation step is performed through an MCMC procedure. Here, since we have to draw values from the joint distribution of (φ , β )|y; θ (k) , we use a Gibbs scheme, i.e., we iteratively draw values from the conditional distributions of φ |y, β (k) ; θ (k) and β |y, φ (k) ; θ (k) . Then, we use again a Hastings-Metropolis algorithm to obtain approximations of these conditional distributions. Finally, iteration k of the SAEM-REML algorithm for model (3) writes:

-S step: using z (k,l) = (φ (k,l) , β (k,l) ), simulate z (k+1,l) = (φ (k+1,l) , β (k+1,l) ), l = 1, . . . , m with a Metropolis-within-Gibbs scheme. -A step: update sk+1 by sk+1 = sk +

χ k 1 m m ∑ l=1
S(y, z (k+1, j) ) -sk , namely:

s1,k+1 = s1,k + χ k 1 m m ∑ l=1 N ∑ i=1 η (k+1,l) i η (k+1,l) i -s1,k s2,k+1 = s2,k + χ k 1 m m ∑ l=1 y -g(z (k+1,l) , f-) 2 -s2,k
where η ,l) . -M step: update θ by θ (k+1) = arg max θ {-Ψ ( θ ) + sk+1 , Φ( θ ) }, namely:

(k+1,l) i = φ (k+1,l) i -A i β (k+1
Γ (k+1) = s1,k+1 N and σ 2 (k+1) = s2,k+1 n . (8)
In many situations, it is important to obtain inference on the fixed effects in the context of REML estimation of variance components. Following [START_REF] Meza | Estimation in the probit normal model for binary outcomes using the saem algorithm[END_REF], estimation of fixed effects can be directly obtained as a by-product of the SAEM-REML algorithm via the expectation of the conditional distribution of the fixed effects given the observed data, the estimate, f , of the unknown function f and the REML estimates of the variance-covariance components. This estimator makes sense in an Empirical Bayes framework.

Estimation of the function f using a LASSO-type method

In this part, our objective is to estimate f in the model (1) using the observations y i, j and assuming that for i = 1, . . . , N we have φ i = φ i and σ 2 = σ 2 where the estimates φ i and σ 2 have been obtained in the precedent SAEM step. Since g satisfies (2), model (1) can be rewritten as

ỹi j = b(φ i ; x i j ) f (x i j ) + ε i j , i = 1 . . . , N, j = 1, . . . , n i with ỹi j = y i j -a(φ i ; x i j ) and xi j = c(φ i ; x i j ).
Of course, since the φ i 's and σ 2 depend on the observations, the distribution of σ -1 ỹi j is no longer Gaussian and the ε i j 's are not i.i.d. but dependent. But in the sequel, to be able to derive theoretical results, we still assume that

ε i j iid ∼ N (0, σ 2 ), (9)
where the value of σ 2 is given by σ 2 . Simulation studies of Section 5 show that this assumption is reasonable. However, note that ( 9) is true at the price of splitting the data set into two parts: the first part for estimating θ and φ , the second part for estimating f . Now, reordering the observations, it is equivalent to observing (y 1 , . . . , y n ) with n = ∑ N i=1 n i , such that

y i = b i f (x i ) + ε i , ε i ∼ N (0, σ 2 ) i.i.d. ( 10
)
where the b i 's and the design (x i ) i=1,...,n are known and depend on the estimators of the precedent SAEM step and the ε i 's are random variables with variance σ 2 estimated by σ 2 . Note that the notation y i , i = 1, . . . , n, does not correspond to the original observations in the SNMM or to any of the values introduced in the previous sections, and it is used in this section for the sake of simplicity. Without loss of generality, we suppose that b i = 0 for all i = 1, . . . , n.

In the sequel, our objective is then to estimate f nonparametrically in model ( 10). A classical method would consist in decomposing f on an orthonormal basis (Fourier basis, wavelets,...) and then to use a standard nonparametric procedure to estimate the coefficients of f associated with this basis ( 0 -penalization, wavelet thresholding,...). In the same spirit as Bertin et al (2011) who investigated the problem of density estimation, we wish to combine a more general dictionary approach with an estimation procedure leading to fast algorithms. The dictionary approach consists in proposing estimates that are linear combinations of various types of functions. Typically, the dictionary is built by gathering together atoms of various classical orthonormal bases. This approach offers two advantages. First, with a more wealthy dictionary than a classical orthonormal basis, we aim at obtaining sparse estimates leading to few estimation errors of the coefficients. Secondly, if the estimator is sparse enough, interesting interpretations of the results are possible by using the set of the non-zero coefficients, which corresponds to the set of functions of the dictionary "selected" by the procedure. For instance, we can point out the frequency of periodic components of the signal if trigonometric functions are selected or local peaks if some wavelets are chosen by the algorithm. Both aspects are illustrated in the next sections. 0 -penalization or thresholding cannot be combined with a dictionary approach if we wish to obtain fast and good algorithms. But LASSO-type estimators based on 1penalization, leading to minimization of convex criteria, constitute a natural tool for the dictionary approach. Furthermore, unlike ridge penalization or more generally p -penalization with p > 1, 1 -penalization leads to sparse solutions for the minimization problem, in the sense that if the tuning parameter is large enough some coefficients are exactly equal to 0 (see [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF]).

There is now huge literature on LASSO-type procedures. From the theoretical point of view and in the specific context of the regression model close to (10), we mention that LASSO procedures have already been studied by Bunea et al (2006), Bunea et al (2007a), Bunea et al (2007b), Bunea (2008), Bickel et al (2009), van de Geer (2010), and Bühlmann and van de Geer (2011) among others.

In our setting, the proposed procedure is the following. For M ∈ N * , we consider a set of functions {ϕ 1 , . . . , ϕ M }, called the dictionary. We denote for λ ∈ R M ,

f λ = M ∑ j=1 λ j ϕ j .
Our objective is to find good candidates for estimating f which are linear combinations of functions of the dictionary, i.e. of the form f λ . We consider, for

λ ∈ R M crit(λ ) = 1 n n ∑ i=1 (y i -b i f λ (x i )) 2 + 2 M ∑ j=1 r n, j |λ j |,
where r n, j = σ ϕ j n

τ log M n
with τ > 0 and for a function h

h 2 n = 1 n n ∑ i=1 b 2 i h 2 (x i ).
We call the LASSO estimator λ the minimizer of λ -→ crit(λ ) for λ ∈ R M and we denote f = f λ .

The function λ -→ crit(λ ) is the sum of two terms: the first one is a goodness-of-fit criterion based on the 2 -loss and the second one is a penalty term that can be viewed as the weighted 1 -norm of λ .

Before going further, let us discuss the important issue of tuning. In our context, the tuning parameter is the constant τ. From a theoretical point of view (see Theorem 1 in the supplementary material), the benchmark value for τ is 2. In the sequel, τ will be chosen satisfying two criteria: to be as close as possible to this benchmark value and allowing the stability of the SAEM algorithm. In Section 5, we will see that sometimes we choose values of τ smaller than 2 but relatively close of it, in particular to obtain the convergence of the variance components estimates, which is always challenging in NLME models.

Once we have chosen a value for τ satisfying these two criteria, the numerical scheme of the nonparametric step is the following:

-Using the estimates of the φ i 's and of σ 2 obtained in the previous iteration of SAEM, compute for i = 1, . . . , n, the observations y i , the constants b i and the design x i . -Evaluate the dictionary {ϕ 1 , . . . , ϕ M } at the design and calculate r n, j . -Obtain the LASSO estimates λ and f λ .

In practice, there exist many efficient algorithms to tackle this third point, namely, the minimization on λ of crit(λ ). For the implementation of our estimation procedure we have considered the approach used by Bertin et al (2011) which consists in using the LARS algorithm.

Numerical results of our procedure are presented in next sections but we also validate our approach from a theoretical point of view. Theoretical results are presented in the supplementary material. We prove oracle inequalities and properties of support for sparse functions under the mild assumption log(M) = o(n). Oracle inequalities ensure that the LASSO estimator of f behaves as well as the best linear combination of functions of the dictionnary. Moreover, we obtain that if the function f is a sparse linear combination of functions from the dictionnary, then the support of the LASSO estimator (functions of the dictionary selected in the LASSO estimator) is included in the support of the function f . These results are generalizations of the results of Bunea et al (2006), Bunea et al (2007a), Bunea et al (2007b), van de Geer (2010) and Bunea (2008) and they are obtained under more general assumptions on the dictionnary. In particular, in our results, the functions of the dictionary do not need to be bounded independently of n and M, which allow us to take wavelet functions.

Estimation algorithm and inferences

We propose the following estimation procedure for semiparametric estimation of (θ , f ) in model ( 3), combining the algorithms described in sections 2.1 and 3:

Estimation Algorithm -ML version: at iteration k, -Given the current estimate of θ , θ (k) = (β (k) ,Γ (k) , σ 2(k)
), and m sampled values of the random effects φ (k,l) , l = 1, . . . , m, update the estimates of f , f (k,l) , l = 1, . . . , m, with the algorithm described in Section 3. -Given the current estimates of f , f (k,l) , l = 1, . . . , m, sample m values of the random effects φ (k,l) , l = 1, . . . , m, and update the value of θ , θ (k+1) = (β (k+1) ,Γ (k+1) , σ 2(k+1) ) with algorithm (7). ( 11)

Estimation Algorithm -REML version: at iteration k, -Given the current estimate of θ , θ (k) = (Γ (k) , σ 2(k)
), and m sampled values of the missing data z (k,l) = (φ (k,l) , β (k,l) ), l = 1, . . . , m, update the estimates of f , f (k,l) , l = 1, . . . , m, with the algorithm described in Section 3.

-Given the current estimates of f , f (k,l) , l = 1, . . . , m, sample m values of the missing data z (k+1,l) = (φ (k+1,l) , β (k+1,l) ), l = 1, . . . , m, and update the value of θ , θ (k+1) = (Γ (k+1) , σ 2(k+1) ) with algorithm (8). ( 12)

As it is explained in Section 2.1, for parametric estimation (SAEM or SAEM-REML algorithms alone) the number of chains, m, can be set to 1, which still guarantees the convergence towards a local maximum of the log-likelihood.

Higher values of m, may accelerate the convergence of the algorithms (but in practice, m is always lower than 10).

For the global semiparametric estimation procedure, we extend this idea of "parallel chains" of values to the estimation of f . Indeed, at iteration k, the estimation of f depends on the value of the missing data, and thus, from m sampled values z (k,1) , . . . , z (k,m) we obtain m estimates of f , f (k,1) , . . . , f (k,m) (see Section 3). Then, in the second step, we use each one of these different estimates of f in parallel to perform parametric estimation (using f (k,l) to sample z (k+1,l) and replacing fby f (k,l) in ( 8) for the estimation of θ ). This is in the case of the REML version of the algorithm, but the same idea underlies the ML version.

Inferences on model and individual parameters, β ,Γ , σ 2 and φ , are performed as in NLMEs (see [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF] and [START_REF] Meza | Estimation in the probit normal model for binary outcomes using the saem algorithm[END_REF]). For inferences on the nonlinear function f , we propose an empirical approach based on the fact that our algorithm automatically provides large samples of estimates of f . Indeed, at each iteration of algorithms ( 11) and ( 12) we obtain m estimates of f . The last iterations of the algorithms typically correspond to small values of χ k in algorithms ( 7) and ( 8), see Section 5 for the details. This can be seen as a phase in which the estimates of parameters are stabilized since we assume that convergence has been reached. Let us note by K and L < K the total number of iterations and the number of iterations in the "stabilization phase" of the algorithm. Then, by considering the last L 0 < L iterations of the algorithm, we get a large sample of estimates of f : f (k,l) , l = 1, . . . , m, k = K -L 0 + 1, . . . , K. These m × L 0 estimates of f are obtained conditionally on values of θ which are supposed to be close to the corresponding ML or REML estimates. Then, we obtain a point estimate for f as:

f = 1 m × L 0 K ∑ k=K-L 0 +1 m ∑ l=1 f (k,l) . ( 13
)
We think that it will be interesting to study how to exploit the estimates f (k,l) to obtain pointwise confidence intervals for f (x). An intuitive empirical pointwise (1-α)100% confidence interval for f (x) could be dfined as follows:

  f (x) -z α 2 S 2 f (x) m × L 0 , f (x) + z α 2 S 2 f (x) m × L 0   . ( 14
)
where

S 2 f (x) = 1 m × L 0 -1 ∑ K k=K-L 0 +1 ∑ m l=1 ( f (k,l) (x)-f (x)) 2 and z α 2 is the 1 -α
2 percentile of a standard normal distribution. This interval is of course not a true (1 -α)100% confidence interval for f (x) but constitutes an approximation of it. It provides a starting point for further research on how function samples generated by semiparametric stochastic approximation algorithms, such us saem-lasso, can be used for inference.

Application to synthetic and real data

Since our procedure consists in the combination of a parametric and a nonparametric estimation algorithm, one may be interested in evaluating the performance of both components separately. In Section 5.1 we provide a simulation study to compare only the parametric versions of our method and Ke and Wang's procedure. In Section 5.2 we compare both procedures in the whole semiparametric setting.

Simulation study: parametric estimation

As a first step, we want to validate through simulation our parametric estimation strategy alone, based on the SAEM algorithm, and to compare it, in the framework of SNMMs, to the FOCE method implemented in [START_REF] Ke | Semiparametric nonlinear mixedeffects models and their applications (with discussion)[END_REF] via the nlme function. In order to be able to assess only the differences induced by the use of different parametric estimation algorithms, we will use the same nonparametric estimation algorithm for the estimation of f , namely the procedure proposed by [START_REF] Ke | Semiparametric nonlinear mixedeffects models and their applications (with discussion)[END_REF]. In Section 5.2, we compare the whole versions, including nonparametric estimation, of both approaches. To this end, we performed the following simulation study based in [START_REF] Ke | Semiparametric nonlinear mixedeffects models and their applications (with discussion)[END_REF] where data were generated from the model:

y i j = φ 1i + exp(φ 2i )2 f j N - exp(φ 3i ) 1 + exp(φ 3i ) + ε i j , i = 1 . . . , N, j = 1, . . . , J,
where ε i j ∼ N (0, σ 2 ) and φ i = (φ 1i , φ 2i , φ 3i ) ∼ N (µ,Γ ) with µ = (µ 1 , µ 2 , µ 3 ) . The nonlinear function was set to f (t) = sin(2πt). As in the original setting, we choose a complex scenario with small sizes of individuals and observations and with high variance values: N = J = 10, µ = (1, 0, 0) , σ 2 = 1 and Γ is diagonal with diag(Γ ) = (1, 0.25, 0.16).

These data were analyzed using two semiparametric procedures: our SAEM based method combined with the nonparametric algorithm of Ke and Wang's (called semi-SAEM) and Ke and Wang's procedure for semiparametric models (called snm). For the SAEM algorithm, we used 80 iterations and the following sequence (χ k ): χ k = 1 for 1 ≤ k ≤ 50 and χ k = 1/(k -50) for 51 ≤ k ≤ 80. We also considered m = 5 chains in each iteration. For the nonparametric estimation algorithm common to both procedures, following [START_REF] Ke | Semiparametric nonlinear mixedeffects models and their applications (with discussion)[END_REF] we considered that f is periodic with period equal to 1 and 1 0 f = 0, i.e. f ∈ W 0 2 (per) = W 2 (per) span{1} where W 2 (per) is the periodic Sobolev space of order 2 in L 2 and span{1} represents the set of constant functions. The same initial values were used for both methods:

µ 0 = (1, 0, 0), σ 2 0 = 2 and diag(Γ 0 ) = (γ 0 1 , γ 0 2 , γ 0 3 ) = (1, 0.3, 0.1).
Tables 1 and2 summarize the performance of both methods over 100 simulated data sets. For each parameter we show the sample mean, the mean squared error (MSE( θ ) = 1 100 ∑ 100 i=1 (θ -θi ) 2 ), and a 95% confidence interval computed over the total number of simulations.

We also compared the REML estimates obtained with our method and with snm (using the REML version of nlme) for the same simulated data sets. The results are summarized in Tables 3 and4. It can be seen that the mean values for the REML estimates obtained with both procedures were closer to the simulated values, especially for the parameter γ 1 . Moreover, the individual confidence intervals of REML estimates of this parameter, at a 95% level, include the true value for these parameters on the contrary to the ML estimates, showing that REML versions of the algorithms were able to correct the bias observed with ML. If we compare our method and snm, for both procedures ML and REML, we obtained results that are similar but it seems that our REML estimates are closer to the simulated values than those obtained with Ke and Wang's method. Furthermore, we can observe that our REML version, in comparison with our ML method, allows to reduce the bias of estimation of variance components in a better way. For instance, in Tables 2 and4, we see that, for γ 1 , we reduce the bias in almost 93% with our REML method whereas with Ke and Wang's REML method this reduction is only of 27%. Finally, let us point out that fixed effects estimates are more accurate with our REML method than with Ke and Wang's one. Let us remind that for SAEM-REML these estimates are the expectation of the conditional distribution of fixed effects given the observed data and the REML estimates of the variancecovariance parameters.

An important issue to discuss is the convergence of estimates with this kind of iterative maximization algorithms. It is well known that approximate methods for maximum likelihood estimation often present numerical problems and even fail to converge in the framework of NLME estimation (see [START_REF] Hartford | Consequences of misspecifying assumptions in nonlinear mixed effects models[END_REF] for instance). An advantage of the exact likelihood method is exactly to avoid those convergence problems as it was established by [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF]. In this simulation study, we have to say that both semi-SAEM and snm achieved convergence for all the data sets. However, we also tried to fit a nonlinear mixed effects model to the simulated data, that is, assuming that f was known and estimating only the fixed and random effects with SAEM and nlme, and in that case the second algorithm failed to converge for several data sets. It seems that in this case the combination of nlme with a nonparametric algorithm to perform semiparametric estimation solves the numerical problems encountered by nlme on its own. However, this is not true in general as we will see in the next simulation study.

Simulation study: semiparametric estimation

In order to test our LASSO-based estimator we consider the same general model of the previous section

y i j = φ 1i + exp(φ 2i )2 f j N - exp(φ 3i ) 1 + exp(φ 3i ) + ε i j , i = 1 . . . , N, j = 1, . . . , J,
where ε i j ∼ N (0, σ 2 ) and φ i = (φ 1i , φ 2i , φ 3i ) ∼ N (µ,Γ ) with µ = (µ 1 , µ 2 , µ 3 ) . Now, f (•) is supposed to be unknown and must be estimated. It is generated as a mixture of one trigonometric function and two Laplace densities (see Figure 1).

f (t) = 0.6 × sin(2πt) +0.2× e -40|t-0.75| 2 × 1 0 e -40|t-0.75| +0.2× e -40|t-0.8| 2 × 1 0 e -40|t-0.80|
. Data were simulated using the following parameters: N = 10, J = 20, µ = (1, 0, 0) , σ 2 = 0.4 and Γ is diagonal with diag(Γ ) = (0.25, 0.16, 0.04).

The chosen function exhibits two sharp peaks that can not be clearly distinguished by only looking at the resulting data (Figure 2). We propose this setting in order to compare the performance of our method and snm in a situation in which the underlying function is not smooth. Indeed, the definition of Ke and Wang's method guarantees that it will achieve very good results if the function to be estimated is well approximated by combinations of spline functions. However, there might be practical situations in which assessing the smoothness of the underlying function might not be easy. It is then interesting to investigate the performance of both methods in such cases. Data were analyzed using the two following semiparametric procedures: our SAEM and LASSO based method (called LASSO-SAEM) and Ke and Wang's procedure for semiparametric models, still denoted snm. For both methods we obtained the REML estimates of parameters.

It is necessary to specify several values in order to run our algorithm, such as the choice of the LASSO's tuning parameter τ and the inputs of the SAEM algorithm (initial values, step sizes χ k , number of chains in the MCMC step, number of burn-in iterations, and total number of iterations). For the latter, we used again 80 iterations with χ k = 1 for 1 ≤ k ≤ 50 and χ k = 1/(k -50) for 51 ≤ k ≤ 80, and we considered m = 5 chains in each iteration. The initial values, which were also used with snm, were: µ 0 = (1, 0, 0), σ 2 0 = 2 and diag(Γ 0 ) = (γ 0 1 , γ 0 2 , γ 0 3 ) = (1, 0.3, 0.1). The nonparametric LASSO step has been performed with τ = 1/3. For some datasets, larger values of τ did not lead to the stabilization of the convergence of some parameters, in particular the variance γ 2 , and smaller values of τ provided similar results to the one presented here. The dictionary chosen combined very different orthonormal families, namely Fourier functions with Haar wavelets, which ensured a sufficiently incoherent design in the spirit of Section 3. More precisely, our dictionary was composed by the following Fourier functions {t →1,t →cos(πt),t →sin(πt),t →cos(2π jt), t → sin(2π jt), j = 1, • • • , 5} and by the Haar wavelet basis with resolution between 2 4 and 2 7 , with a total size of 245 functions. Note that the data xi j = c(φ i ; x i j ) belongs approximately to [-0.4, 1.6]. For snm, we took f ∈ W 0 2 (per). Of course, the true function does not belong to that space and a partial spline model with possible change points would be more appropriate for modeling it. However, we want to reflect the fact that in a real situation the only information available is the one provided by the observed data set. In this case the simulated data exhibit a clear periodic structure which we try to capture with a function in W 0 2 (per). In Figures 1 and2, we can see the estimates of f compared with the true function and the fitted data with the two meth-ods for a specific simulated data set. Results for REML estimates obtained with LASSO-SAEM and snm for 100 simulated data sets are summarized in Tables 5 and6. We can see that the means of the estimates obtained with our method are close to their real values except for the variance of the error, σ 2 , since our method tends to overestimate that parameter. However, we get overall better results than using the snm methodology (except for γ 1 ).

An important issue for this kind of problem is the estimation of the nonlinear function f . Ke and Wang's method based on splines works very well for regular functions. So, it is interesting to study its performance on less smooth functions, which is typically the case with the function f considered here. Then, to evaluate the accuracy of the estimation, we calculated the Integrated Square Error (ISE) of f for each simulated data set. Figure 3 provides a summary of estimates of f using LASSO-SAEM and snm. We computed the ISE for each estimate of f and plotted the estimates corresponding to (a) the minimum, (b) 1/4 quantile, (c) median, (d) 3/4 quantile and (e) maximum ISEs. We can see that our method outperforms snm in the estimation of f , in the sense that our estimates are able to detect the presence of the peaks in the original function.

As for the functions of the dictionary selected with our LASSO method, it is interesting to note that the 100 linear combinations of functions of the dictionary obtained for each one of the 100 data sets have a length which varies between 10 and 32 functions, with an average length equal to 20. Furthermore, in 98% of the cases, the method selects the function sin(2πt) with the highest coefficient. For the remaining two data sets, the functions sin(6πt) and sin(10πt) are selected. For all the replicates, in addition to these sine functions, the rest of the selected functions are related to the Haar wavelets with smaller coefficients. So, our method is quite robust.

It is important to point out that the results obtained with snm are based only on 51 data sets since the function did not reach convergence in 46 data sets and in other 3 data sets we obtained incoherent estimation of the nonlinear function, when using the default setup of the snm algorithm (REML estimation and Generalized Cross Validation for the choice of the penalized parameter). By contrast, our method achieved convergence for all simulated data sets with the specific setup used here (choice of τ, initial values, number of chains, step sizes χ k , number of iterations, etc . . .).

To assess the robustness of the LASSO procedure, we have also performed an analysis of these data sets with a dictionary that is composed by the union of the dictionary defined above (the 245 functions) and the dictionary used in Section 5.3 (the 64 functions). The results obtained are very similar to those presented in Tables 5 and6, so we have not included them here. Moreover, the estimates of f are also very similar. In particular, for 50% of the data sets, the estimates of f select only components in the old dictionary (with Fourier and wavelet functions) and for all the datasets, only 7% of the selected functions belongs to the dictionary defined in Section 5.3. Additionally, the function sin(2πt) is selected with the highest coefficient in 90% of the cases.

Finally, we compute the confidence intervals defined in ( 14) with L 0 = 20. We obtained very thin confidence intervals and a poor coverage (less to 40%) with these datasets. We think that it is a challenging issue to study if an appropriate choice of L 0 in ( 14) may yield to more robust intervals.

Application to on-line auction data

Modeling of price paths in on-line auction data has received a lot of attention in the last years [START_REF] Shmueli | Visualizing online auctions[END_REF][START_REF] Jank | Functional data analysis in electronic commerce research[END_REF][START_REF] Shmueli | The BARISTA: a model for bid arrivals in online auctions[END_REF][START_REF] Liu | Functional data analysis for sparse auction data[END_REF]. One of the reasons is the availability of huge amounts of data made public by the on-line auction and shopping website eBay.com, which has become a global market place in which millions of people worldwide buy and sell products. The price evolution during an auction can be thought as a continuous process which is observed discretely and sparsely only at the instants in which bids are placed. In fact, bids tend to concentrate at the beginning and at the end of the auction, responding to two typically observed phenomena, "early bidding" and "bid sniping" (a situation in which "snipers" place their bids at the very last moment). To our knowledge, [START_REF] Reithinger | Modelling price paths in on-line auctions: smoothing sparse and unevenly sampled curves by using semiparametric mixed models[END_REF] provide the first attempt to model price paths taking into account the dependence among different auctions. This is an important consideration, since in practice bidders can participate in different auctions that take place simultaneously. They propose a semiparametric additive mixed model with a boosting estimation approach. In the same line, but considering a more complex interaction of the random effects and the unknown nonlinear function, we propose the following shapeinvariant model for the price paths:

y i j = φ 1i + exp(φ 2i ) f (t i j -φ 3i ) + ε i j , i = 1, • • • , N, j = 1, • • • , n i , where ε i j ∼ N (0, σ 2 ) and φ i = (φ 1i , φ 2i , φ 3i ) ∼ N (µ,Γ ) with µ = (µ 1 , µ 2 , µ 3 ) .
We introduce an individual random horizontal shift, φ 3i , to model the possible delay of the price dynamics in some auctions with respect to the rest. We analyzed a set of 183 eBay auctions for Palm M515 Personal Digital Assistants (PDA), of a fixed duration of seven days, that took place between March and May, 2003. This is the data set used in [START_REF] Reithinger | Modelling price paths in on-line auctions: smoothing sparse and unevenly sampled curves by using semiparametric mixed models[END_REF] and it is publicly available at http://www.rhsmith.umd.edu/digits/ statistics/data.aspx. We were interested in modeling the live bids, that is, the actual prices that are shown by eBay during the live auction. Note that these are different from the bids placed by bidders during the auction, which are the prices recorded in the bid history published by eBay after the auction closes. Then, a transformation on the bid records is required to recover the live bids (see [START_REF] Shmueli | Visualizing online auctions[END_REF] for details). The live bids range from $0.01 to $300 and form a sequence of non decreasing prices for each auction. We typically observe between 10 and 30 bids per auction, although there are auctions with only two bids. We have a total of 3280 bids for the 183 auctions. Following [START_REF] Reithinger | Modelling price paths in on-line auctions: smoothing sparse and unevenly sampled curves by using semiparametric mixed models[END_REF], we considered the square root of live bids to reduce the price variability. We run the REML version of our LASSO-SAEM algorithm, of which we performed 100 iterations with the following sequence of decreasing steps (χ k ) k : χ k = 1 for 1 ≤ k ≤ 60 and χ k = 1/(k -60) for 61 ≤ k ≤ 100. We also considered m = 3 chains in each iteration. The dictionary for nonparametric estimation was composed by a combination of B-splines of degrees three and four, with 17 knots unequally spaced so that most of the knots were in those places with more data observed (at the beginning, at the end and at the middle of the interval), 10 power functions, 10 exponential functions and 5 logit functions, with a total size of 64 functions. The estimate of f is monotone, as expected by the nature of the data, and presents two steepest parts at the beginning and at the end of the interval. At each iteration of the algorithm the estimated function at the nonparametric step is a sparse combination of the functions of the dictionary. In fact, the set of functions selected by the LASSO method at the last iterations of the algorithm is almost constant, containing mainly two functions, ϕ(x) = x 0.35 and ϕ(x) = exp(0.9 x), and in some iterations a small component of a cubic B-spline around the middle of the interval. In Figure 4 we present the last 24 estimates f (k,l) from which we have obtained f as in ( 13), and f , together with a 95% pointwise confidence band. These results have been obtained with τ = 2 as the value for the tuning parameter in the LASSO estimation step. The estimates for µ and Γ are presented in Table 7.

To assess the robustness of the LASSO procedure, we have also performed an analysis of this data set with a dictionary that is composed by the union of the dictionary defined above and the dictionary used in Section 5.2 to analyze the simulated data. That is, we have added Fourier and Haar Fig. 5 Left: Estimates of f obtained with snm (dashed line) and saem-lasso with the large dictionary (solid lines) and 95% confidence bands (gray shadows). Right: Last 24 LASSO estimates in saem-lasso.

Time

Square root of live bids Auction 177

Fig. 6 Observed live bids (circles) and fitted price curves for a subset of 18 auctions obtained with snm (dahsed lines) and saem-lasso (solid lines) with the large dictionary. wavelets bases to the dictionary initially chosen. The results are very similar to those obtained with the original dictionary. They are shown in Figure 5. In particular, the estimates of f are almost identical. Among the last 24 estimates of f , f (k,l) , obtained with this new dictionary, only two estimates contain a significant component of functions not included in the original dictionary. To compare our method to Ke and Wang's, in Figures 5 and6 we also present the results of the analysis of this data set with snm. We have to mention that we have performed this analysis with five different function models for f and two different criteria for the estimation of the smoothing parameter, namely, general cross validation (GCV) and generalized maximum likelihood (GML). So, we ran snm with ten different specifications, among which we got convergence for only six specifications. None of the six estimates of f is strictly monotone and five of them are extremely rough.

In Figure 5 we present the smoothest snm f -estimate, which is obtained by modeling f with cubic splines and by using the GLM criterion, together with the saem-LASSO estimate obtained with the largest dictionary. In Figure 6 we present the observed live bids and the model fits for 18 chosen auctions with different price profiles. We can appreciate how the fitted models provide in general an accurate fit of the final price, even in the cases when bid sniping is present.

There are some differences between the two fits, mostly at the beginning of each auction, although the fitted curves are in general similar with the two methods. For the rest of the combinations of a function model and a smoothing estimation criterion used with snm, the fits of the data are suboptimal. Indeed, the fitted price curves produce almost perfect interpolation of the data.

As for the computation time, saem-lasso took 300 seconds to run on these data on a 2.5 GHz Mac OS X whereas the average time for snm over the six runs was about six hours on the same computer.

Conclusions and discussion

Semiparametric nonlinear mixed effects models cover a wide range of situations and generalize a large class of models, such as nonlinear mixed effects models or self-modelling nonlinear regression models among others. We have proposed a new method for estimation in SNMMs combining an exact likelihood estimation algorithm with a LASSOtype procedure. Our strategy relies on an iterative procedure to estimate θ conditioned on f and vice versa, which allow us to tackle the parametric and the nonparametric problem independently. This makes possible the use of fast algorithms providing an accurate and computationally efficient estimation method. Concerning parametric estimation, our simulation results illustrate our method and point out some important advantages of using an exact likelihood estimation algorithm instead of likelihood approximation methods, such as convergence of the estimates. The REML version of our algorithm, corrects the estimation of variance components accounting for the loss of degrees of freedom from estimating the fixed effects and provide satisfactory results. However, as it was already pointed out in the comments to [START_REF] Ke | Semiparametric nonlinear mixedeffects models and their applications (with discussion)[END_REF], it will be important to define a REML estimator that can also take into account the loss of degrees of freedom from estimating the nonlinear function. As for computational aspects, we have to mention that the SAEM algorithm avoids the convergence problems encountered by nlme based routines.

For nonparametric aspects, the dictionary approach based on LASSO algorithms shows, in some situations, some improvements when compared with Ke and Wangs' methodology. This is the case for instance for spiky or non-continuous functions to be estimated. Our dictionary method can adapt to different features of signals for wealthy enough dictionaries. Furthermore, our methodology allows us to obtain interesting interpretation with respect to the functions of the dictionary selected by the procedure. For instance, we can detect trends, frequencies of sinusoids or location and heights of peaks of the common shape represented by the estimated function f . We have observed that our LASSO estimate achieves good theoretical and numerical results if the dictionary is wealthy and incoherent enough. From the theoretical point of view, incoherence is expressed, in this paper, by Assumption A1(s) or by the quantity ρ(S * ) defined in the Supplementary Material. These incoherence assumptions are hard to check in practice and we do not know if they can be relaxed in our setting.

We mention that our method can be non robust if the dictionary is not wealthy enough. That is, if the function to be estimated cannot be well approximated by linear combinations of the functions of the dictionary, the functions that are selected can vary from one simulation to another, which may lead to different estimates. However, it the main features of a signal (periodicity, smoothness, peaks,...) are included in the dictionary, our method is very robust to the enlarging of the dictionary with additional functions, as seen in Sections 5.2 and 5.3.

In Section 3, the particular structure of the observations (where we have n i observations for each individual i) is not used for applying the standard LASSO-procedure. But a natural and possible extension of this work would be to take into account this structure and then to apply a more sophisticated LASSO-type procedure inspired, for instance, by the group-LASSO proposed by [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] to achieve better results. This is a challenging research axis we wish to investigate from a theoretical and practical point of view. The LASSO is a very popular algorithm, but Hybrid Adaptive Spline, MARS or BSML (see [START_REF] Sklar | Non-parametric regression with basis selection from multiple libraries[END_REF]) could also be combined with the dictionary approach proposed in this paper. Since results of our paper show that the dictionary approach seems promising, results of our paper could be extended by using algorithms mentioned previously from both theoretical and practical points of view.

Among other possible extensions of this work, a very promising one would be the use of the nonparametric techniques herein described for density estimation (in the spirit of (Bertin et al, 2011)) of the random errors, assuming that they do not need to be normal. Indeed, the recent work of [START_REF] Comte | Nonparametric estimation of random effects densities in linear mixed-effects model[END_REF] deals with this problem in the case of a linear mixed effects model. Its generalization to NLMEs or even SNMMs is a real challenge.

Supplementary Material for "LASSO-type estimators for Semiparametric Nonlinear Mixed-Effects Models Estimation" Ana Arribas-Gil • Karine Bertin • Cristian Meza • Vincent Rivoirard the date of receipt and acceptance should be inserted later 1 Theoretical results for the LASSO-type estimator

Assumptions

As usual, assumptions on the dictionary are necessary to obtain oracle results for LASSOtype procedures. We refer the reader to van de Geer and [START_REF] Van De Geer | On the conditions used to prove oracle results for the Lasso[END_REF] for a good review of different assumptions considered in the literature for LASSO-type estimators and connections between them. The dictionary approach aims at extending results for orthonormal bases. Actually, our assumptions express the relaxation of the orthonormality property. To describe them, we introduce the following notation. For l ∈ N, we denote

ν min (l) = min |J|≤l min λ ∈R M λ J ̸ =0 || f λ J || 2 n ||λ J || 2 ℓ 2 and ν max (l) = max |J|≤l max λ ∈R M λ J ̸ =0 || f λ J || 2 n ||λ J || 2 ℓ 2
, where || • || ℓ 2 is the l 2 norm in R M . The notation λ J means that for any k ∈ {1, . . . , M}, (λ J ) k = λ k if k ∈ J and (λ J ) k = 0 otherwise. Previous quantities correspond to the "restricted" eigenvalues of the Gram matrix G = (G j, j ′ ) with coefficients

G j, j ′ = 1 n n ∑ i=1 b 2 i φ j (x i )φ j ′ (x i ).
Assuming that ν min (l) and ν max (l) are close to 1 means that every set of columns of G with cardinality less than l behaves like an orthonormal system. We also consider the restricted correlations

δ l,l ′ = max |J|≤l |J ′ |≤l ′ J∩J ′ = / 0 max λ ,λ ′ ∈R M λ J ̸ =0,λ ′ J ′ ̸ =0 ⟨ f λ J , f λ ′ J ′ ⟩ ||λ J || ℓ 2 ||λ ′ J ′ || ℓ 2 , where ⟨ f , g⟩ = 1 n ∑ n i=1 b 2 i f (x i )g(x i
). Small values of δ l,l ′ means that two disjoint sets of columns of G with cardinality less than l and l ′ span nearly orthogonal spaces. We will use the following assumption considered in Bickel et al (2009).

Assumption 1 For some integer 1 ≤ s ≤ M/2, we have ν min (2s) > δ s,2s .

(A1(s))

Oracle inequalities of the Dantzig selector were established under this assumption in the parametric linear model by [START_REF] Candès | The Dantzig selector: statistical estimation when p is much larger than n[END_REF] and for density estimation by Bertin et al (2011). It was also considered by Bickel et al (2009) for nonparametric regression and for the LASSO estimate.

Let us denote 2s) .

κ s = √ ν min (2s) ( 1 - δ s,2s ν min (2s) ) > 0, µ s = δ s,2s √ ν min (
We will say that λ ∈ R M satisfies the Dantzig constraints if for all j = 1, . . . , M

(Gλ ) j -β j ≤ r n, j , (1)
where

β j = 1 n n ∑ i=1 b i φ j (x i )Y i .
We denote D the set of λ that satisfies (1). The classical use of Karush-Kuhn-Tucker conditions shows that the LASSO estimator λ ∈ D, so it satisfies the Dantzig constraint. Finally, we assume in the sequel M ≤ exp(n δ ), for δ < 1. Therefore, if ∥φ j ∥ n is bounded by a constant independent of n and M, then r n, j = o(1) and oracle inequalities established below are meaningful.

Oracle inequalities

We obtain the following oracle inequalities.

Theorem 1 Let τ > 2. With probability at least 1 -M 1-τ/2 , for any integer s < n/2 such that (A1(s)) holds, we have for any α > 0,

|| f -f || 2 n ≤ inf λ ∈R M inf J 0 ⊂{1,...,M} |J 0 |=s { || f λ -f || 2 n + α ( 1 + 2µ s κ s ) 2 Λ (λ , J c 0 ) 2 s + 16s ( 1 α + 1 κ 2 s ) r 2 n } (2) where r n = sup j=1,...,M r n, j , Λ (λ , J c 0 ) = ||λ J C 0 || ℓ 1 + ( || λ || ℓ 1 -||λ || ℓ 1 ) + 2 , for any x ∈ R x + := max(x, 0) and || • || ℓ 1 is the l 1 norm in R M .
Theorem 2 Let τ > 2. With probability at least 1 -M 1-τ/2 , for any integer s < n/2 such that (A1(s)) holds, we have for any α > 0,

|| f -f || 2 n ≤ inf λ ∈D inf J 0 ⊂{1,...,M} |J 0 |=s    || f λ -f || 2 n + α ( 1 + 2µ s κ s ) 2 ||λ J C 0 || ℓ 1 + || λJ C 0 || ℓ 1 s + 32s ( 1 α + 1 κ 2 s ) r 2 n    . (3)
Similar oracle inequalities were established by Bunea et al (2006), Bunea et al (2007a), Bunea et al (2007b), or van de Geer (2010). But in these works, the functions of the dictionary are assumed to be bounded by a constant independent of M and n. Let us comment the right-hand side of inequalities ( 2) and (3) of Theorems 1 and 2. The first term is an approximation term which measures the closeness between f and f λ and that can vanish if f is a linear combination of the functions of the dictionary. The second term can be considered as a bias term. In both theorems, the term ||λ J C 0 || ℓ 1 corresponds to the cost of having λ with a support different of J 0 . For a given λ , this term can be minimized by choosing J 0 as the set of largest coordinates of λ . Note that if the function f has a sparse expansion on the dictionary, that is f = f λ where λ is a vector with s non-zero coordinates, then by choosing J 0 as the set of the s non-zero coordinates, the approximation term and the term ||λ J C 0 || ℓ 1 vanish.

In Theorem 1, the term

( || λ || ℓ 1 -||λ || ℓ 1 ) +
will be smaller as the ℓ 1 -norm of the LASSO estimator is small and this term is equal to 0 if || λ || ℓ 1 ≤ ||λ || ℓ 1 , which is frequently the case.

In Theorem 2, given a vector λ such that f λ approximates well f , the term || λJ C 0 || ℓ 1 will be small if the LASSO estimator selects the largest coordinates of λ . The last term can be viewed as a variance term corresponding to the estimation of f as linear combination of s functions of the dictionary (see Bertin et al (2011) for more details). Finally, the parameter α calibrates the weights given for the bias and variance terms.

The following section deals with estimation of sparse functions.

The support property of the LASSO estimate

Let τ > 2. In this section, we apply the LASSO procedure with rn, j instead of r n, j , with

rn, j = σ ∥φ j ∥ n √ τ log M n , τ > τ.
We assume that the regression function f can be decomposed on the dictionary: there exists

λ * ∈ R M such that f = M ∑ j=1 λ * j φ j .
We denote S * the support of λ * :

S * = { j ∈ {1, . . . , M} : λ * j ̸ = 0 } ,
and by s * the cardinal of S * . We still consider the LASSO estimate λ and, similarly, we denote Ŝ the support of λ :

Ŝ = { j ∈ {1, . . . , M} : λ j ̸ = 0 } .
One goal of this section is to show that with high probability, we have: Ŝ ⊂ S * .

We have the following result.

Theorem 3 We define

ρ(S * ) = max k∈S * max j̸ =k | < φ j , φ k > | ∥φ j ∥ n ∥φ k
∥ n and we assume that there exists c ∈ (0, 1/3) such that s * ρ(S * ) ≤ c.

If we have

√ τ + √ τ √ τ - √ τ ≤ 1 -c 2c , then P { Ŝ ⊂ S * } ≥ 1 -2M 1-τ/2
. A similar result was established by Bunea (2008) in a slightly less general model. However, her result is based on strong assumptions on the dictionary, namely each function is bounded by a constant L (see Assumption (A2)(a) in Bunea (2008)). This assumption is mild when considering dictionaries only based on Fourier bases. It is no longer the case when wavelets are considered and Bunea's assumption is satisfied only in the case where L depends on M and n on the one hand and is very large on the other hand. Since L plays a main role in the definition of the tuning parameters of the method, with too rough values for L, the procedure cannot achieve satisfying numerical results for moderate values of n even if asymptotic theoretical results of the procedure are good. In the setting of this paper, where we aim at providing calibrated statistical procedures, we avoid such assumptions.

Finally, we have the following corollary.

Corollary 1 We suppose that A1(s * ) is satisfied and that there exists c ∈ (0, 1/3) such that s * ρ(S * ) ≤ c.

If we have

√ τ + √ τ √ τ - √ τ ≤ 1 -c 2c , then, with probability at least 1 -4M 1-τ/2 , || f -f || 2 n ≤ 32s * r2 n κ s * , where rn = sup j=1,...,M rn, j .
This corollary is a simple consequence of Theorem 2 with λ = λ * and J 0 = S * . Taking λ = λ * implies that the approximation term vanishes. Taking J 0 = S * implies that the bias term vanishes since the support of the LASSO estimator is included in the the support of λ * . In this case, assuming that sup j ∥φ j ∥ n < ∞, the rate of convergence is the classical rate

s * log M n .
2 The proofs

Preliminary lemma

Lemma 1 For 1 ≤ j ≤ M, we consider the event

A j = { |V j | < r n, j } where V j = 1 n ∑ n i=1 b i φ j (x i )ε i . Then, P (A j ) ≥ 1 -M -τ/2 .
Proof of Lemma 1: We have

P ( A c j ) ≤ P (√ n|V j |/(σ ∥φ j ∥ n ) ≥ √ nr n, j /(σ ∥φ j ∥ n ) ) ≤ P ( |Z| ≥ √ τ log M ) ≤ M -τ/2
where Z is a standard normal variable.

Proof of Theorem 1

Let λ ∈ R M and J 0 such that |J 0 | = s. We have

∥ f λ -f ∥ 2 n = ∥ f -f ∥ 2 n + ∥ f λ -f ∥ 2 n + 2 n n ∑ i=1 b 2 i ( f (x i ) -f (x i ) ) ( f λ (x i ) -f (x i ) ) . We have ∥ f λ -f ∥ 2 n = ∥ f ∆ ∥ 2 n where ∆ = λ -λ . Moreover A = 2 n n ∑ i=1 b 2 i ( f (x i ) -f (x i ) ) ( f λ (x i ) -f (x i ) ) = 2 M ∑ j=1 (λ j -λ j ) [ (G λ ) j -β j ] ,
where

β j = 1 n n ∑ i=1 b 2 i φ j (x i ) f (x i ).
Since λ satisfies the Dantzig constraint, we have with probability at least 1 -M 1-τ/2 , for any j ∈ {1, . . . , M},

|(G λ ) j -β j | ≤ |(G λ ) j -β j | + | β j -β j | ≤ 2r n, j
and |A| ≤ 4r n ∥∆ ∥ 1 . This implies that

∥ f -f ∥ 2 n ≤ ∥ f λ -f ∥ 2 n + 4r n ∥∆ ∥ 1 -∥ f ∆ ∥ 2 n .
Moreover using Lemma 1 and Proposition 1 of Bertin et al (2011) (where the norm ∥ • ∥ 2 is replaced by ∥ • ∥ n ), we obtain that

( ||∆ J C 0 || ℓ 1 -||∆ J 0 || ℓ 1 ) + ≤ 2||λ J C 0 || ℓ 1 + ( || λ || ℓ 1 -||λ || ℓ 1 ) + (4)
and

|| f ∆ || n ≥ κ s ||∆ J 0 || ℓ 2 - µ s √ |J 0 | ( ||∆ J C 0 || ℓ 1 -||∆ J 0 || ℓ 1 ) + ≥ κ s ||∆ J 0 || ℓ 2 -2 µ s √ |J 0 | Λ (λ , J c 0 ).
Note that Proposition 1 of Bertin et al (2011) is obtained using Lemma 2 and Lemma 3 of Bertin et al (2011). In our context, Lemma 2 and Lemma 3 can be proved in the same way by replacing the norm ∥ • ∥ 2 by ∥ • ∥ n and by considering P J 01 as the projector on the linear space spanned by (φ j (x 1 ), . . . , φ j (x n )) j∈J 01 . Now following the same lines as Theorem 2 of Bertin et al (2011), replacing κ J 0 by κ s and µ J 0 by µ s , we obtain the result of the theorem.

Proof of Theorem 2

We consider λ D defined by λ D = argmin λ ∈R M ||λ || ℓ 1 such that λ satisfies the Dantzig constraint (1).

Denote by f D the estimator f λ D . Following the same lines as in the proof of Theorem 1, it can be obtained that, with probability at least 1 -M 1-τ/2 , for any integer s < n/2 such that (A1(s)) holds, we have for any α > 0,

|| f D -f || 2 n ≤ inf λ ∈R M inf J 0 ⊂{1,...,M} |J 0 |=s { || f λ -f || 2 n + α ( 1 + 2µ s κ s ) 2 Λ (λ , J c 0 ) 2 s + 16s ( 1 α + 1 κ 2 s ) r 2 n } ,
where here

Λ (λ , J c 0 ) = ||λ J C 0 || ℓ 1 + ( || λ D || ℓ 1 -||λ || ℓ 1 ) + 2 .
If the infimum is only taken over the vectors λ that satisfy the Dantzig constraint, then, with the same probability we have

|| f D -f || 2 n ≤ inf λ ∈D inf J 0 ⊂{1,...,M} |J 0 |=s { || f λ -f || 2 n + α ( 1 + 2µ s κ s ) 2 ∥λ J C 0 ∥ 2 l 1 s + 16s ( 1 α + 1 κ 2 s ) r 2 n } .
(5) Following the same lines as the proof of Theorem 1, replacing λ by λ D , we obtain, with probability at least 1 -M 1-τ/2 ,

∥ f -f ∥ 2 n ≤ ∥ f D -f ∥ 2 n + 4r n ∥∆ ∥ 1 -∥ f ∆ ∥ 2 n ,
with ∆ = λλ D . Applying (4) where λ plays the role of λ and λ D the role of λ , the vector ∆ satisfies

( ||∆ J C 0 || ℓ 1 -||∆ J 0 || ℓ 1 ) + ≤ 2|| λJ C 0 || ℓ 1 .
Following the same lines as in the proof of Theorem 1, we obtain that for each J 0 ⊂ {1, . . . , M} such that

|J 0 | = s || f -f || 2 n ≤    || f D -f || 2 n + α ( 1 + 2µ s κ s ) 2 ∥ λJ C 0 ∥ 2 l 1 s + 16s ( 1 α + 1 κ 2 s ) r 2 n    . (6)
Finally, ( 5) and ( 6) imply the theorem. Using notation of Lemma 3, we have:

∥ f λ * -f μ ∥ 2 n = ∥ ∑ k∈S * (λ * k -μk )φ k ∥ 2 n = ∑ k∈S * (λ * k -μk ) 2 ∥φ k ∥ 2 n + ∑ k∈S * ∑ j∈S * , j̸ =k (λ * k -μk )(λ * j -μ j ) < φ j , φ k >,
and

∑ k∈S * (λ * k -μk ) 2 ∥φ k ∥ 2 n ≤ ∥ f λ * -f μ ∥ 2 n + ρ(S * ) ∑ k∈S * ∑ j∈S * , j̸ =k |λ * k -μk |∥φ k ∥ n × |λ * j -μ j |∥φ j ∥ n ≤ ∥ f λ * -f μ ∥ 2 n + ρ(S * ) ( ∑ k∈S * |λ * k -μk |∥φ k ∥ n ) 2
.

Finally,

( ∑ k∈S * |λ * k -μk |∥φ k ∥ n ) 2 ≤ s * ∑ k∈S * (λ * k -μk ) 2 ∥φ k ∥ 2 n ≤ s *   ∥ f λ * -f μ ∥ 2 n + ρ(S * ) ( ∑ k∈S * |λ * k -μk |∥φ k ∥ n ) 2   ,
which shows that

( ∑ k∈S * |λ * k -μk |∥φ k ∥ n ) 2 ≤ s * 1 -ρ(S * )s * ∥ f λ * -f μ ∥ 2 n . Now, 1 n n ∑ i=1 ( y i -b i ∑ j∈S * μ j φ j (x i ) ) 2 + 2 ∑ j∈S * rn, j | μ j | ≤ 1 n n ∑ i=1 ( y i -b i ∑ j∈S * λ * j φ j (x i ) ) 2 + 2 ∑ j∈S * rn, j |λ * j |. So, ∥ ∑ j∈S * μ j φ j ∥ 2 n - 2 n n ∑ i=1 b i y i ∑ j∈S * μ j φ j (x i ) + 2 ∑ j∈S * rn, j | μ j | ≤ ∥ ∑ j∈S * λ * j φ j ∥ 2 n - 2 n n ∑ i=1 b i y i ∑ j∈S * λ * j φ j (x i ) + 2 ∑ j∈S *
rn, j |λ * j |, and using previous notation,

∥ f μ ∥ 2 n - 2 n n ∑ i=1 b i y i ∑ j∈S * μ j φ j (x i ) + 2 ∑ j∈S * rn, j | μ j | ≤ ∥ f λ * ∥ 2 n - 2 n n ∑ i=1 b i y i ∑ j∈S * λ * j φ j (x i ) + 2 ∑ j∈S * rn, j |λ * j |.
Therefore,

∥ f λ * -f μ ∥ 2 n = ∥ f μ ∥ 2 n + ∥ f λ * ∥ 2 n -2 < f μ , f λ * > ≤ 2∥ f λ * ∥ 2 n -2 < f μ , f λ * > + 2 n n ∑ i=1 b i y i ∑ j∈S * ( μ j -λ * j )φ j (x i ) + 2 ∑ j∈S * rn, j (|λ * j | -| μ j |) = 2 n n ∑ i=1 b i y i ( f μ (x i ) -f λ * (x i )) - 2 n n ∑ i=1 b 2 i f λ * (x i )( f μ (x i ) -f λ * (x i )) +2 ∑ j∈S * rn, j (|λ * j | -| μ j |) = 2 n n ∑ i=1 b i (y i -E(y i ))( f μ (x i ) -f λ * (x i )) + 2 ∑ j∈S * rn, j (|λ * j | -| μ j |) = 2 n n ∑ i=1 b i ε i ( f μ (x i ) -f λ * (x i )) + 2 ∑ j∈S * rn, j (|λ * j | -| μ j |) = 2 M ∑ j=1 V j ( μ j -λ * j ) + 2 ∑ j∈S * rn, j (|λ * j | -| μ j |).
Now let us assume that for any j ∈ S * , V j < r n, j . Then,

∥ f λ * -f μ ∥ 2 n < 2 ∑ j∈S * (r n, j + rn, j )| μ j -λ * j | < 2σ √ log M n ( √ τ + √ τ) ∑ j∈S * ∥φ j ∥ n | μ j -λ * j |.

Fig. 1

 1 Fig.1True function f (solid line) and its estimates obtained with LASSO-SAEM (dashed line) and snm (dash-dotted line) for a particular data set in the semiparametric simulation study.

Fig. 2

 2 Fig.2Simulated data and fitted curves obtained with LASSO-SAEM (solid line) and snm (dashed line) for a particular data set in the semiparametric simulation study.

Fig. 3

 3 Fig. 3 Estimated functions corresponding to the five quantiles of ISE ((a) minimum, (b) 1/4 quantile, (c) median, (d) 3/4 quantile and (e) maximum) obtained with LASSO-SAEM (dashed line) and snm (dash-dotted line) compared to the true function f (solid line) for the total of the 100 simulated data sets in the semiparametric simulation study.

Fig. 4

 4 Fig. 4 Left: Estimated nonlinear function f (solid line) and 95% confidence band (gray shadow) in the on-line auction data set. Right: Last 24 LASSO estimates whose empirical mean provides f .

<|

  i b i φ j (x i )) -∑ k∈S * μk < φ j , φ k > ≥ rn, j -r n, j φ j , f λ * > -∑ k∈S * μk < φ j , φ k > ≥ rn, j -r n, j k -μk ) < φ j , φ k > ≥ rn, j -r n, j * )∥φ j ∥ n ∑ k∈S * |λ * k -μk | ∥φ k ∥ n ≥ rn, j -r n, j < φ j , φ k > | ∥φ j ∥ n ∥φ k ∥ n .

  (S * )s * and for any j / ∈ S * ,ρ(S * )∥φ j ∥ n ∑ k∈S * |λ * k -μk |∥φ k ∥ n < 2σ √ log M n ∥φ j ∥ n ( √ τ + √ τ) ρ(S * )s * 1 -ρ(S * )s * ≥ r n, j } and using Lemma 1, since P {S c } ≤ A + B, P {S } ≥ 1 -2M 1-τ/2 .

Table 1

 1 ML procedure: Mean, MSE and 95% confidence interval of mean components.

	True Value	Method	µ 1 1	µ 2 0	µ 3 0
	Mean MSE 95 % C.I.	semi-SAEM snm semi-SAEM snm semi-SAEM [0.99;1.12] [0.27;0.36] 1.06 0.31 1.05 0.26 0.12 0.16 0.12 0.11 snm [0.99;1.12] [0.22;0.30] [-0.02;0.01] 0.27 -0.01 0.10 0.01 [0.23;0.30]

Table 2

 2 ML procedure: Mean, MSE and 95% confidence interval of variance components obtained with semi-SAEM and snm.

	True Value	Method	γ 1 1	γ 2 0.25	γ 3 0.16	σ 2 1
	Mean MSE 95 % C.I.	semi-SAEM snm semi-SAEM snm semi-SAEM [0.77;0.95] [0.21;0.27] [0.14;0.17] [0.92;0.98] 0.86 0.24 0.16 0.95 0.89 0.19 0.14 0.99 0.22 0.02 0.01 0.03 0.22 0.02 0.01 0.03 snm [0.80;0.98] [0.17;0.21] [0.13;0.16] [0.96;1.02]

Table 3

 3 REML procedure: Mean, MSE and 95% confidence interval of mean components.

	True Value	Method	µ 1 1	µ 2 0	µ 3 0
	Mean MSE 95 % C.I.	semi-SAEM snm semi-SAEM snm semi-SAEM [1.01;1.07] [-0.03;0.02] [-0.02;0.01] 1.04 -0.01 -0.01 1.05 0.26 -0.01 0.03 0.02 0.01 0.12 0.11 0.01 snm [0.99;1.12] [0.22;0.30] [-0.02;0.01]

Table 4

 4 REML procedure: Mean, MSE and 95% confidence interval of variance components obtained with semi-SAEM and snm.

	True Value	Method	γ 1 1	γ 2 0.25	γ 3 0.16	σ 2 1
	Mean MSE 95 % C.I.	semi-SAEM snm semi-SAEM snm semi-SAEM [0.89;1.08] [0.22;0.28] [0.14;0.18] [0.92;0.98] 0.99 0.25 0.16 0.95 0.92 0.19 0.15 1.02 0.21 0.03 0.01 0.03 0.23 0.02 0.01 0.03 snm [0.83;1.02] [0.17;0.22] [0.13;0.17] [0.98;1.05]

Table 5

 5 REML procedure: Mean, MSE and 95% confidence interval of means components obtained with LASSO-SAEM and snm.

	True Value	Method	µ 1 1	µ 2 0	µ 3 0
	Mean MSE 95 % C.I.	LASSO-SAEM snm LASSO-SAEM snm LASSO-SAEM snm	0.97 1.09 0.009 0.019 [0.949;0.984] [1.057;1.119 ] [1.293;1.482] [-0.025;0.015] 0.02 0.01 1.39 -0.01 0.009 0.003 2.035 0.005 [0.005;0.041] [-0.006;0.014]

Table 6

 6 REML procedure: Mean, MSE and 95% confidence interval of variance components obtained with LASSO-SAEM and snm.

	True Value	Method	γ 1 0.25	γ 2 0.16	γ 3 0.04	σ 2 0.4
	Mean MSE 95 % C.I.	LASSO-SAEM snm LASSO-SAEM snm LASSO-SAEM [0.16;0.20] [0.12;0.15] [0.030;0.037] [0.66;0.73] 0.18 0.14 0.03 0.69 0.21 0.11 0.03 0.90 0.01 0.01 4.0e-4 0.12 0.02 0.01 5.9e-4 0.27 snm [0.18;0.25] [0.09;0.14] [0.028;0.042] [0.86;0.94]

Table 7

 7 Estimated mean vector and covariance matrix of the random effects and estimated error variance in the on-line auction data set.

	φ 1 1.04 Correlation 1 (7.68) Mean Matrix -0.02 (variances) 0.41 σ 2	φ 2 0.18 -0.02 1 (0.19) 0.37 1.93	φ 3 -0.06 0.41 0.37 1 (0.23) φ 3 φ 1 φ 2

Acknowledgements

The authors would like to thank the anonymous Associate Editor and two referees for valuable comments and suggestions.

Suppementary Material

In the supplement available on-line we provide theoretical results and proofs for the LASSO-type estimator of Section 3.

is supported by projects FONDECYT 1090285 and ECOS/CONICYT C10E03 2010, Chile. Cristian Meza is supported by project FONDECYT 11090024, Chile.

Proof of Theorem 3

We first state the following lemma.

Lemma 2 We have for any u ∈ R M ,

Proof of Lemma 2: Since for any λ ,

Since λ minimizes λ -→ crit(λ ), we have for any k,

where |s( λk )| ≤ 1 and s( λk

, which proves the result. Now, still with s * = card(S * ), we consider for

and μ = arg min µ∈R s * critS * (µ). Then we set

and we state the following lemma.

Lemma 3 On the set S , the non-zero coordinates of λ are included into S * .

Proof of Lemma 3: Recall that λ is a minimizer of λ -→ crit(λ ). Using standard convex analysis arguments, this is equivalent to say that for any 1

Similarly, on S , we have

So, on S , the vector μ such μ j = μ j if j ∈ S * and μ j = 0 if j / ∈ S * is also a minimizer of λ -→ crit(λ ). Using Lemma 2, we have for any 1 ≤ i ≤ n:

Therefore, on S , the non-zero coordinates of λ are included into S * .

Lemma 3 shows that we just need to prove that

(see Lemma 1) and

Proof of Corollary 1

First note that λ * satisfies the Dantzig constraint (1) where r n, j is replaced by rn, j with probability larger than 1 -M 1-τ/2 . On the event Ŝ ⊂ S * , we have λ * (S * ) C = λ(S * ) C = 0, then applying Theorem 2, we obtain that for any α > 0

which implies the result of the theorem.