# Estimation and variable selection in high dimension in nonlinear mixed-effects models.

Antoine Caillebotte, Estelle Kuhn, Sarah Lemler

<figure><img src="image_1.png" /></figure>

## To cite this version:

Antoine Caillebotte, Estelle Kuhn, Sarah Lemler. Estimation and variable selection in high dimension in nonlinear mixed-effects models.. 2025. hal-05005304v3

HAL Id: hal-05005304

https://hal.inrae.fr/hal-05005304v3

Preprint submitted on 1 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ESTIMATION AND VARIABLE SELECTION IN HIGH DIMENSION IN
NONLINEAR MIXED-EFFECTS MODELS.

Antoine Caillebotte ¹,² & Estelle Kuhn ² & Sarah Lemler ³

¹ Université Paris-Saclay, INRAE, UMR GQE-Moulon, France, caillebotte.antoine@inrae.fr,

² Université Paris-Saclay, INRAE, UR MAIAGE, France, estelle.kuhn@inrae.fr,

³ Université Paris-Saclay, CentraleSupélec, Laboratoire MICS, France,
sarah.lemler@centralesupelec.fr

**Abstract.** We consider nonlinear mixed effects models including high-dimensional covariates to model individual parameters variability. The objective is to identify relevant covariates among a large set under sparsity assumption and to estimate model parameters. To face the high dimensional setting we consider a regularized estimator namely the maximum likelihood estimator penalized with the $l_1$-penalty. We rely on the use of the eBIC model choice criteria to select an optimal reduced model. Then we estimate the parameters by maximizing the likelihood of the reduced model. We calculate in practice the maximum likelihood estimator penalized with the $l_1$-penalty though a weighted proximal stochastic gradient descent algorithm with an adaptive learning rate. This choice allows us to consider very general models, in particular models that do not belong to the curved exponential family. We demonstrate first in a simple linear toy model through a simulation study the good convergence properties of this optimization algorithm. We compare then the performance of the proposed methodology with those of the GLMM-LASSO procedure in a linear mixed effects model in a simulation study. We illustrate also its performance in a nonlinear mixed-effects logistic growth model through simulation. We finally highlight the benefit of the proposed procedure relying on an integrated single step approach regarding two others two steps approaches for variable selection objective.

**Keywords.** nonlinear mixed effects model, high dimension, variable selection, LASSO penalty, stochastic gradient descent, weighted proximal.
# 1 Introduction

Mixed effects models are very fine and very useful statistical modeling tools for analyzing data with hierarchical structures and repeated measurements (see Pinheiro and Bates [2006]). In particular it is possible to account for several levels of variability of a phenomenon observed within a population of individuals, each individual being observed several times in different conditions. They are widely used in many applied fields such as agronomy, pharmacology, and even economics. Mixed effects models are composed of two nested levels of modeling: on the one hand a common structural modeling of the phenomenon of interest parameterized for each individual in the population by specific individual parameters, on the other hand a modeling of these individual parameters as random variables which accounts for their variability within the population. The first modeled the intra individual variability among the repeated measurements of each individual, the latter stands for the inter individual variability between individuals in the population.

The structural modeling of the phenomenon of interest can be linear or non-linear in the individual parameters. Non-linear type modeling can in particular account for complex phenomena modeled mechanistically by models which integrate, for example, physical or biological knowledge. The parameters of such models are often of strong practical interest because they are interpretable from an applied point of view. In such setting a central objective is to characterize and explain the variation of these individual parameters within the population. Therefore the modeling of individual parameters is done through random effects at the individual level and can also integrate descriptive covariates of the individuals.

Depending on the context, these descriptive covariates can be of high dimension and the objective is then to identify those which are the most relevant to explain the variabilities observed within the population by estimating the associated vector of regression parameters. Let us consider the field of plant ecophysiology. Mechanistic models have been proposed to describe plant development processes. These models integrate descriptive variables of the environment as covariates acting directly on the plant development process. The parameters of these models are often physical quantities such as leaf appearance speeds or light interception capacities. These parameters may vary when considering a population of plants from different varieties each characterized by its genotype (see Baey et al. [2018]). Their variations can be modeled into a mixed effects model and also integrate large genetic markers characterizing the genotypic variability within the population. In this context, identifying the relevant covariates among a set of high-dimensional covariates amounts to identifying the genetic markers that influence the phenomenon of interest (see Bhatnagar et al. [2020]).

From the point of view of inference in mixed effects models including high-dimensional covariates, the objective is on the one hand to select the relevant covariates from a set of high-dimensional covariates in order to identify a parsimonious model with a reduced number of covariates and on the other hand to estimate the parameters in the reduced model. There are two main difficulties for the inference task. The first one is the high dimension of the covariates. Usually the selection of relevant covariates can be done via a regularization approach. Furthermore, in the context of mixed effects models, an additional difficulty appears due to the presence of random effects which are not observed. This is a classic context of latent variable models. Inference is complex to carry out in this framework due to the latent structure and involves the use of efficient numerical methods. Parameter inference can be done for example by maximum likelihood via Expectation Maximization (EM) or stochastic gradient type algorithms. In the context of exponential family models, EM type algorithms are easy to implement and have good theoretical properties (see Dempster et al. [1977]; Delyon et al. [1999]). On the other hand, to our knowledge, there are no theoretical results of convergence outside this framework. In addition, the implementation of these algorithms outside the exponential family is more complex. To get around this limitation, a trick called exponentialization trick is sometimes used in practice. However, its limits have been highlighted in Debavelaere and Allassonnière [2021].
In particular, this procedure can generate significant estimation biases due to the fact that the inference of the parameters is carried out in an extended model different from the initial one. Stochastic gradient methods can be applied in more generic models, in particular outside the exponential family. Theoretical guarantees of convergence towards an extremum of the target function have been established under general assumptions. However, these methods are quite sensible in practice to the tuning of the learning rate, in particular in high-dimensional parame- ter spaces due to the heterogeneity of the different components of the gradient. Adaptive choices of the learning rate and procedures based on gradient preconditioning have been proposed in generic contexts to overcome this computational difficulty [Duchi et al., 2011; Kingma and Ba, 2017]. More recently in the context of maximum likelihood estimation in general latent variable models, a stochastic gradient algorithm integrating a gradient preconditioning step based on an estimator of the Fisher information matrix obtained as a by-product of the algorithm has been proposed, opening new possibilities for maximum likelihood inference (Baey et al. [2023]).

In order to achieve the objective of selecting relevant variables from a set of high-dimensional variables in mixed effects models, several regularization approaches have been developed. Jürg Schell-dorfer and Bühlmann [2014] proposed a maximum likelihood estimator with a LASSO penalty (Tibshirani [1996]) in the context of linear mixed effects models with a Gaussian error term and developed a R package (see Schelldorfer et al. [2014]). Fort et al. [2019] and Ollier [2022] proposed estimators with more general penalties in the case of nonlinear mixed effects models belonging to the curved exponential family. Bertrand and Balding [2013] have also proposed an stochastic penalized versions of the EM (Delyon et al. [1999]) algorithm to take into account multiple parameters in pharmacokinetic models. Bayesian approaches based on “spike and slab” distributions have also been developed by Heuclin et al. [2020] in the case of linear mixed models and by Naveau et al. [2024] in the non linear case. On the other hand, to our knowledge, there are no high-dimensional variable selection methods for general nonlinear mixed-effects models, in particular outside the exponential family.

In this contribution, we propose a variable selection and re-estimation procedure in a gen-
eral mixed effects model with high dimensional covariates explaining individual parameters vari-
ability. We deal with the variable selection in high dimension through a maximum likelihood
estimate regularized via a $\ell_1$-penalty term. To calculate this penalized estimator, we use an
adaptive weighted proximal stochastic gradient algorithm to simultaneously handle the presence
of latent variables and the non differentiability of the penalty term. The latter allows to handle
very general models, in particular our procedure does not require that the models belong to the
curved exponential family. This paper is organised as follows. The second section introduces the
mixed-effects models with high-dimensional covariates and presents two classical examples. The
third section is devoted to the description of the proposed variable selection procedure in the
high-dimensional setting and the parameter’s re-estimation in the reduced model. The fourth
section gathers practical details of the numerical methodology. Finally, we present a simulation
study and discuss the potential of the proposed method with respect of existing ones. The paper
ends with some conclusion and perspectives.

2 Nonlinear mixed effects model integrating high dimensional covariates to model individual variability

2.1 Model description

Let $N$ be a positive integer. We consider $J$ repeated measurements for each individual $i$.
Therefore we have $J$ observations per individual. Let us denoted by $Y_{i,j}$ the $j$-th observation
of the $i$-th individual for $1 \le i \le N$ and $1 \le j \le J$. We assume that $Y_{i,j}$ takes value in $\mathbb{R}^d$. We
model this observation with a nonlinear mixed effects model (see Pinheiro and Bates [2006] and
Davidian [2017]):

$$ 
\left\{ 
\begin{array}{l}
Y_{i,j} = m(\alpha, t_{i,j}, Z_i) + \varepsilon_{i,j}, \\
Z_i \stackrel{i.}{\sim} \mathcal{N}(X_i\beta, \Omega), \\
\varepsilon_{i,j} \stackrel{i.i.d}{\sim} \mathcal{N}(0, \Sigma).
\end{array}
\right. 
$$

where $m$ is a function depending on population parameter $\alpha$ taking values in $\mathbb{R}^a$, individual covariates $(t_{i,j})$ and individual parameters modeled by the latent variable $Z_i$. We consider functions that can be either linear or nonlinear with respect to the latent variable. Note that in the case of longitudinal data, the covariates $t_{i,j}$ stands usually for the $j$-th observation time of individual $i$. Individual parameters $(Z_i)$ describes the inter-individual variability of the population. For the $i$-th subject, the quantity $Z_i$ is a $q$-dimensional random vector, independent of $(\varepsilon_{i,j})_j$ and assumed to be distributed as a Gaussian distribution with expectation $X_i\beta$ and covariance matrix $\Omega$, where $X_i$ is a matrix of size $q \times p$ and $\beta$ is a vector of size $p$. The noise term $\varepsilon_{i,j}$ is usually assumed centered Gaussian with unknown covariance matrix $\Sigma$. The unknown parameters of the nonlinear mixed-effects model are therefore $\theta = (\alpha, \beta, \Omega, \Sigma) \in \mathbb{R}^a \times \mathbb{R}^p \times \mathbb{S}_q^{++} \times \mathbb{S}_d^{++}$ where $\mathbb{S}_d^{++}$ stands for the set of symmetric positive definite matrices of size $d \times d$. We emphasize that our model formulation is very general and allows us to consider model function $m$ depending on both population parameters and individual parameters. We consider the setting of high-dimensional covariates where $p$ can be much larger than $N$ and assume that only a small subset of covariates are relevant to be included in the model.

**Remark 2.1.** More general distributions than the Gaussian one can be chosen for the noise term as well as more general settings with various numbers of observations per individual can be considered.

**Remark 2.2.** It may be useful in practice to reparametrize the latent variables $(Z_i)$ and/or the parameters. For example, consider as individual parameter $\log(Z_i)$ rather than $Z_i$ if individual parameters are assumed to be positive. We refer for more details to the reparametrization cookbook (Leger [2023]).

## 2.2 Practical examples

We present in this section two classical nonlinear mixed effects models, one modeling general growing process and one from the pharmacological field.

### 2.2.1 Logistic growth curves

We consider as first example the specific case of the logistic curve model used to model growth dynamic, which is commonly used in the nonlinear mixed-effect models' community and presented in Pinheiro and Bates [2006]. The observation $Y_{i,j}$ is the circumference of the tree $i$ at time $t_{ij}$ and is modeled through the logistic model with the function $m$ given by:

$$ m(\alpha, t_{ij}, Z_i) = \frac{Z_{i1}}{1 + \exp\left(-\frac{t_{ij} - Z_{i2}}{\alpha}\right)}, \quad (1) $$

where the individual parameters $Z_{i1}$ represents the asymptotical maximum value of the circumference, $Z_{i2}$ represents the value of the sigmoid's midpoint, and $\alpha$ represents the logistic growth rate.

### 2.2.2 Pharmacodynamic model

We consider as a second example the two compartments pharmacodynamic model used by Pinheiro and Bates [2006]. The observation $Y_{i,j}$ is the serum concentration measured at time
$t_{i,j}$. The model is given by:

$$m(t_{ij}, Z_i) = \frac{D_i k_{ai}}{V_i(k_{ai} - C l_i / V_i)} \left( \exp(-k_{ai} t_{ij}) - \exp\left(-\frac{C l_i}{V_i} t_{ij}\right) \right), \quad (2)$$

where the individuals parameters $V_i$ represents the distribution volume, $k_{ai}$ the absorption rate, $C l_i$ the clearance and $C l_i / V_i$ the elimination rate; $D_i$ is the known dose of the drug receive by individual $i$.

# 3 Variable selection and estimation procedure

In this section, we present the proposed method called NLMEM-LASSO for both selecting the relevant covariates and estimating the model parameters.

## 3.1 Estimation in latent variable model

We consider the maximum likelihood estimator to infer the non-linear mixed-effects model's parameters. In the context of latent variable models, the marginal likelihood, denoted by $g$, is obtained by integrating the complete likelihood over the latent variables, which are not observed.

$$\begin{aligned} g(\theta; Y) &= \prod_{i=1}^{N} \int f(\theta; Y_i, Z_i) dZ_i = \prod_{i=1}^{N} \int p_{\theta}(Y_i|Z_i) p_{\theta}(Z_i) dZ_i \\ &= \prod_{i=1}^{N} \int \left\{ \prod_{j=1}^{J} p_{\theta}(Y_{i,j}|Z_i) \right\} p_{\theta}(Z_i) dZ_i \end{aligned} \quad (3)$$

where $f(\theta; Y, Z)$, $p_{\theta}(Y|Z)$, $p_{\theta}(Z)$ are respectively the density of the pair $(Y, Z)$, the density of $Y$ conditionally to $Z$, and the density of $Z$.

One can usually estimate the model parameters by maximizing the logarithm of the marginal likelihood using the maximum likelihood estimator written as follows:

$$\hat{\theta}^{\text{MLE}} = \arg \max_{\theta \in \Theta} \log g(\theta; Y), \quad (4)$$

where $\Theta$ denotes the parameter space. From a practical point of view, this estimate can often not be calculated explicitly. To deal with latent variables, classical methods used to infer the unknown parameters are Expectation Maximization like algorithms (see Ng et al. [2012]). However two major limits of these procedures is that i) theoretical convergence properties are established only for models belonging to the curved exponential family; ii) their implementation is easy for models of the curved exponential family but much tricky outside this setting. However one can use the exponentiation trick by defining certain parameters as Gaussian variables and study an augmented model that belongs to the exponential family. Nevertheless, its limits have been highlighted in Debavelaere and Allassonnière [2021]. Another type of numerical approaches to calculate the MLE are stochastic gradient descent algorithms (see Cappé et al. [2005]). Recently Baey et al. [2023] have presented a preconditioned stochastic gradient descent for estimation in a latent variable model adapted to general latent variables models, in particular theoretical guarantees are obtained for general models, without assuming that the models belong to the curved exponential family.

## 3.2 Penalized likelihood in high-dimensional setting

In our context, we must deal with the high dimension of the covariates, therefore we introduce a regularization term and consider a penalized maximum likelihood estimator. We aim to select
relevant variables among the covariates and follow the idea of the LASSO (Least Absolute Shrinkage and Selection Operator) procedure, which was initially developed for linear regression models in Tibshirani [1996]. This method enables us to handle high-dimensional covariates and select a subset of explanatory covariates from a large collection. Therefore we introduce a $\ell_1$ penalty term, which only depends on the parameter $\beta$:

$$ \text{pen}_{\lambda}(\theta) = \lambda \|\beta\|_1 = \lambda \sum_{k=1}^{p} |\beta_k|, \qquad (5) $$

where $\lambda$ is a positive real called regularization parameter. Our goal is then to maximize the penalized criterion defined as the difference of the logarithm of the marginal likelihood and of the penalty term. Let us define the penalized maximum likelihood estimator by:

$$ \hat{\theta}_{\lambda}^{\text{PEN}} = \arg\max_{\theta \in \Theta} \{\log g(\theta; Y) - \text{pen}_{\lambda}(\theta)\}, \qquad (6) $$

where $\Theta$ denotes the parameter space and where $\lambda$ is a positive parameter. The larger the value of $\lambda$, the more $\beta$ will be constrained to have zero components. Conversely, the smaller the value of $\lambda$, the freer the components of $\beta$ will be. It is customary to determine the value of $\lambda$ using cross-validation or using model criterion (see Tibshirani [1996]). We will consider the eBIC model choice criterion well adapted to the high dimensional setting to find on optimal regularization value (see Chen and Chen [2008]).

### 3.3 Regularization path procedure

To select only the most explanatory covariates, it's important to choose a well-balanced value for the regularization parameter $\lambda$. We choose an optimal parameter by minimizing the extended Bayesian Information Criterion (eBIC) (see Chen and Chen [2008]). Guided by the intuition given by Delattre et al. [2014], we penalize the number of degrees of freedom by the log of the total number of observations $NJ$. We consider a grid $\Lambda$ of values for the regularization parameter $\lambda$. We conduct then the following inference methodology :

i) For all $\lambda \in \Lambda$ repeat the following steps:

* calculate $\hat{\theta}_{\lambda}^{\text{PEN}} = \arg\max_{\theta \in \Theta} \{\log g(\theta; Y) - \text{pen}_{\lambda}(\theta)\}$.

* deduce the associated support $\hat{S}_{\lambda} = \{j \in \{1, \dots, p\}, \hat{\beta}_{\lambda,j}^{\text{PEN}} \neq 0\}$

* calculate $\hat{\theta}_{\lambda}^{\text{MLE}} = \arg\max_{\theta \in \Theta} \{\log g_{\lambda}(\theta; Y)\}$, where $g_{\lambda}$ is the marginal likelihood in the model restricted to the support $\hat{S}_{\lambda}$.

* calculate the eBIC criterion:

$$ \text{eBIC}(\lambda) = -2 \log g_{\lambda}(\hat{\theta}_{\lambda}^{\text{MLE}}; Y) + |\hat{S}_{\lambda}| \log(NJ) + 2 \log \left( \binom{p}{|\hat{S}_{\lambda}|} \right) \qquad (7) $$

ii) Select the parameter $\lambda$ that minimizes the eBIC:

$$ \hat{\lambda} = \arg\min_{\lambda \in \Lambda} \text{eBIC}(\lambda), $$

and consider the final estimator defined by $\hat{\theta}_{\lambda}^{\text{MLE}}$.
# 4 Numerical methodology in practice

In order to carry out in practice our variable selection and re-estimation procedure, we need to be able to efficiently and accurately calculate both the maximum likelihood estimate in a reduced model and the penalized estimate for numerous different values of the regularization parameter. For both objectives, we have to deal with the integral on the latent variables, which means that there are no explicit expression for the criteria to be optimised in most cases. Therefore, we will consider numerical methods to solve these maximization problems. We propose to use an adaptive stochastic gradient algorithm for the first objective and a proximal weighted extended version for the second objective.

## 4.1 Adaptive Stochastic Gradient Descent Algorithm

Adaptive algorithms such as AdaGrad [Duchi et al., 2011], RMSPprop [Tieleman, 2012] and Adam [Kingma and Ba, 2017] have proved their worth in many settings. The objective of these strategies is to automatically rescale through the adaptive learning rate each gradient descent direction, in order to scale the different components of the gradient, homogenizing the evolution of the algorithm. We choose to benefit from the advantages offered by the adaptive algorithm AdaGrad. The main advantage is the adaptive property of the proposed learning rate, since the convergence behavior of stochastic gradient algorithms depends strongly of the choice of the learning rate. The second advantage relies on the component wise independence property of the Adagrad learning rate which will be of great importance when dealing with the non-differentiable penalty term of the criterion in the next section. The algorithm is divided into two steps. The first step consists in the simulation of a realization of the latent variables, through a direct sampling from the posterior distribution using the current value of the parameter or a step of a Metropolis-Hastings sampler (see Marin and Robert [2007]). The second step is the update of the parameter estimate using the adaptive preconditioned gradient descent on the approximate complete likelihood, More precisely, given the realization of the latent variable, we evaluate the gradient of the complete log-likelihood and calculate for each component of the parameter the cumulative sums of the square of these gradients, and define the preconditioning diagonal matrix proportional to the inverse of the square root of these values.

## 4.2 Adaptive Weighted Proximal Stochastic Gradient Descent Algorithm

When evaluating the penalized estimate in practice, we have to face to the non-differentiability of the $l_1$-penalty term and consider therefore a proximal algorithm as presented for example by Fort et al. [2019]. Note that in our context of high-dimensional setting, it is difficult to manage a full matrix as preconditioner from a computational point of view. Therefore the adaptive vectorial learning rate detailed in the previous section is particularly well-adapted since the preconditioning matrix is diagonal and thus allows to evaluate the proximal step component by component. The main idea is to combine the stochastic gradient descent algorithm with a Proximal Forward-Backward algorithm (Chen and Rockafellar [1997]; Tseng [2000]). The algorithm is divided into three steps; a realization of the latent variables is sampled with a first step called *Simulation*, as in the previous case. The second step is the adaptive preconditioned gradient descent on the approximate complete likelihood, called the *Forward* step. The last step, called *Backward*, deals with the penalty term. We apply the proximal operator (Moreau [1962]; Rockafellar [1976]; Chouzenoux et al. [2014]) detailed in the following paragraf.

We first describe the weighted proximal operator. First of all, for any symmetric positive-definite matrix $A \in S_d^{++}$, let's denote by $\|\cdot\|_A$ the A-weighted norm defined with the A-weighted scalar product by $\forall (x, y) \in \mathbb{R}^d, \langle x, y \rangle_A = \langle x, Ax \rangle$ and $\|x\|_A^2 = \langle x, x \rangle_A$.

**Definition 4.1.** Assume that $A \in S_d^{++}$. Let $\psi: \Theta \to \mathbb{R}$ be a proper, lower semicontinuous, convex function and $\theta \in \Theta$. The weighted proximal operator of $\psi$ in $\theta$ relative to the metric
induced by A is the unique minimizer of $\psi(\theta') + \frac{1}{2} \|\cdot - \theta\|_A^2$:

$$ \operatorname{Prox}_{A,\psi}(\theta) = \arg \min_{\theta' \in \mathbb{R}^d} \left( \psi(\theta') + \frac{1}{2} \|\theta' - \theta\|_A^2 \right) \quad (8) $$

The general resulting update step of the proximal gradient descent algorithm targeting the minimum of a given function f can then be written as:

$$ \theta_{k+1} = \operatorname{Prox}_{A_k, \operatorname{pen}_{\lambda}} (\theta_k - A_k^{-1} \nabla_{\theta} f(\theta_k)) $$

where $(A_k)_{k \ge 1}$ is a sequence of symmetric positive-definite preconditioning matrices.

**Remark 4.2.** We emphasize that when $A_k = \gamma_k^{-1} \mathrm{Id}$, the operator $\operatorname{Prox}_{\gamma^{-1} \mathrm{Id}, \psi}$ is the classical isotropic proximal operator defined by Moreau [1962] and we obtain the standard proximal stochastic gradient descent where $(\gamma_k)_{k \ge 1}$ is the common step size.

In general, $\operatorname{Prox}_{A,\mathrm{pen}}$ is well defined but may not be quickly computable. Becker and Fadili state several assumptions on A and the function pen that allow the proximal operator to be computable as for example, rank restriction and pen is separable, i.e. $\mathrm{pen}(x) = \sum_{l=1}^p \mathrm{pen}_l(x_l)$. If no assumptions are made, the proximal operator is often not explicit and will require additional steps. Since we choose the AdaGrad type learning rate and the $\ell_1$-penalty, the calculation of the proximal operator is separable and can be perfomed component per component independently and has an explicit form: for a given $\lambda > 0$, we get for all $s > 0$ and $x \in \mathbb{R}$

$$ \operatorname{Prox}_{s, \lambda|\cdot|}(x) = \begin{cases} 0 & \text{if } |x| < \lambda/s \\ x - \lambda/s & \text{if } x \ge \lambda/s \\ x + \lambda/s & \text{if } x \le -\lambda/s \end{cases} \quad (9) $$

We provide here in detail the steps of the adaptive weighted proximal stochastic gradient descent algorithm.

<table><caption>**Algorithm 1:** Adaptive Weighted Proximal Stochastic Gradient (AWPSG)</caption><tr><td colspan="2">Require: $\gamma_0 > 0$, $\epsilon > 0$, $K_{max} \in \mathbb{N}$</td></tr><tr><td>1</td><td>Initialize Initialize starting point $\theta_0 \in \mathbb{R}^d$; adaptive stepsize vector $s_0 = 0$; iteration $k = 0$.</td></tr><tr><td>2</td><td>while ($\theta_k$) not converged or $k \le K_{max}$ do</td></tr><tr><td>3</td><td>&nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow k + 1$</td></tr><tr><td>4</td><td>&nbsp;&nbsp;&nbsp;&nbsp;• Simulation step :</td></tr><tr><td>5</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Draw $Z^{(k)}$ from the posterior distribution $p(\cdot|\theta_k)$ or using a single step of a Hastings Metropolis procedure having as stationary distribution the posterior.</td></tr><tr><td>6</td><td>&nbsp;&nbsp;&nbsp;&nbsp;• Gradient computation : $v_k \leftarrow \frac{1}{N} \sum_{i=1}^N \nabla \log f_{\theta_k}(Y_i, Z_i^{(k)})$</td></tr><tr><td>7</td><td>&nbsp;&nbsp;&nbsp;&nbsp;• Adaptive stepsize : $s_{k,l}^2 \leftarrow s_{k-1,l}^2 + v_{k,l}^2$ for $l \in \{1, \dots, p\}$</td></tr><tr><td>8</td><td>&nbsp;&nbsp;&nbsp;&nbsp;• Update parameters :</td></tr><tr><td>9</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• Forward step : $\omega_{k,l} \leftarrow \theta_{k-1,l} + \gamma_0 v_{k,l} / (s_{k,l} + \epsilon)$ for $l \in \{1, \dots, p\}$</td></tr><tr><td>10</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;• Backward step : $\theta_{k,l} \leftarrow \operatorname{Prox}_{(s_{k,l}+\epsilon)/\gamma_0, \lambda|\cdot|(\omega_{k,l})}$ for $l \in \{1, \dots, p\}$</td></tr><tr><td>11</td><td>end</td></tr><tr><td>12</td><td>Return $\hat{\theta} = \theta_k$</td></tr></table>

**Remark 4.3.** Regarding the simulation of the realization of the latent variable $Z^{(k)}$ at iteration $k$, it is often the case that direct sampling from the posterior distribution is not possible and one must require to Monte Carlo Markov Chain procedure such as the Hastings Metropolis algorithm or the Gibbs algorithm (see Marin and Robert [2007]).
## 4.3 Details on the regularization path

In this section, we give more detail about the procedure used to choose a regularization parameter $\lambda$ (see section 3.3). The minimization of the eBIC criterion (7) calculated on a proposal grid $\Lambda$ can be represented by the regularization path which can be plotted to visualize the evolution of the support of $\beta$ with respect to the regularization parameter $\lambda$ (e.g. figure 7). This type of graph shows that the larger the regularization parameter, the more the component of the vector $\beta$ is restricted to zero. The smaller $\lambda$ is, the more $\beta$ components are free. In practice, the proposal grid $\Lambda$ is built on a logarithmic scale. In fact, it's preferable to explore the smaller values a little more, in order to capture rapid changes in the support of $\beta$. This way, we don't spend too much time on ranges of regularization values that would have returned the same support.

Note also that to evaluate the eBIC criterion defined in (7) one must evaluate the quantity $g_{\lambda}(\hat{\theta}_{\lambda}^{MLE}; Y)$ which does not have an explicit form due to the presence of the latent variables $Z$ and the non-linearity of the model. Therefore we use a Monte Carlo procedure to calculate an approximation of $g_{\lambda}(\hat{\theta}_{\lambda}^{MLE}; Y)$.

# 5 Numerical experiments

In this section, we study the performance of the proposed procedure through simulations. The numerical study is divided into three parts. First, we demonstrate the good behavior of the optimization algorithm presented in section 4.2. Then we compare the performance of the variable selection procedure with re-estimation in a linear mixed-effects model to those of the existing method GLMMCLASSO. We assess the robustness of our method when dealing with a nonlinear mixed-effect model in different configurations. We emphasize that to the best of our knowledge, there exists actually no other methodology to deal with this issue in general nonlinear mixed-effect model. Finally, we highlight the benefits of the proposed integrated procedure, which combines the mixed population model and the variable selection procedure in a single step, compared with two other two-step approaches.

## 5.1 Convergence analysis of the adaptive proximal weighted gradient stochastic descent algorithm

In this section, we demonstrate through simulation the convergence properties of the proposed adaptive proximal weighted gradient stochastic descent algorithm to optimize for each value of the regularization parameter $\lambda$ in a grid $\Lambda$ the penalized log-likelihood criterion defined in Equation (6). To that purpose we consider a simple toy model where the value of the parameter $\theta$ that maximizes the targeted criterion can be calculated exactly though an explicit expression in order to highlight the accuracy of the proposed algorithm. Let us therefore consider for $N \in \mathbb{N}^*$ and $J \in \mathbb{N}$, $J \ge 2$ the following linear mixed effects model where for all $i \in \{1, ..., N\}$

$$ \left\{ \begin{array}{l} Y_i = X_i\beta + WZ_i + \varepsilon_i \\ Z_i \stackrel{i.i.d}{\sim} \mathcal{N}(0, \gamma^2) \\ \varepsilon_i \stackrel{i.i.d}{\sim} \mathcal{N}(0, \sigma^2 I_J) \end{array} \right. \qquad (10) $$

where $Y_i = (Y_{i,1}, ..., Y_{i,J}) \in \mathbb{R}^J$, $X_i = (X_{i,1}, ..., X_{i,J}) \in \mathcal{M}_{J,p}(\mathbb{R})$, $\varepsilon_i = (\varepsilon_{i,1}, ..., \varepsilon_{i,J}) \in \mathbb{R}^J$ and $W$ is the vector of $\mathbb{R}^J$ with all components equal to 1. We assume that for each $i$ the vectors $Z_i$ and $\varepsilon_i$ are independent. Thus the marginal distribution of $Y_i$ is Gaussian with expectation $X_i\beta$ and covariance matrix $\Gamma = \sigma^2 I_J + \gamma^2 WW^t$. We assume also that the random vectors $(Z_i)$ are independent., as well as the random vectors $(\varepsilon_i)$.
Note that thanks to the Woodbury formulae, we get an explicit expression for the inverse of $\Gamma$, given by:

$$ \Gamma^{-1} = \frac{1}{\sigma^2} I_J - \frac{\gamma^2}{\sigma^2(J\gamma^2 + \sigma^2)} WW^t \quad (11) $$

We assume that the variance parameters $\gamma^2$ and $\sigma^2$ are known and that parameters $\beta \in \mathbb{R}^p$ are estimated. Thus the penalized criterion writes, up to an additive constant:

$$ C_{\lambda}(\beta) = \log g(\beta; Y) - \lambda \sum_{k=1}^{p} |\beta_k| \quad (12) $$

$$ = -\frac{1}{2} \sum_{i=1}^{N} (Y_i - X_i \beta)^t \Gamma^{-1} (Y_i - X_i \beta) - \lambda \sum_{k=1}^{p} |\beta_k| \quad (13) $$

Using the expression given in (11), we get:

$$ C_{\lambda}(\beta) = -\frac{1}{2\sigma^2} \sum_{i=1}^{N} \|Y_i - X_i\beta\|^2 - \frac{\gamma^2}{2\sigma^2(J\gamma^2 + \sigma^2)} \sum_{i=1}^{N} \left( \sum_{j=1}^{J} Y_{ij} - \left(\sum_{j=1}^{J} X_{ij}\right)\beta \right)^2 - \lambda \sum_{k=1}^{p} |\beta_k| $$

where $X_{ij}$ is the j-th line of matrix $X_i$. Let us define the assumption (H1): $\forall 1 \le i \le N$ and $\forall 1 \le l \le p \sum_{j=1}^{J} X_{ijl} = 0$. If the design matrix $X$ satisfies (H1), then the associated penalized criterion writes, up to an additive constant:

$$ C_{\lambda}(\beta) = -\frac{1}{2\sigma^2} \sum_{i=1}^{N} \|Y_i - X_i\beta\|^2 - \lambda \sum_{k=1}^{p} |\beta_k| $$

Note that this expression of the criterion corresponds exactly to the LASSO criterion expression in the classical homoscedastic linear regression. If we assume additionally that the design $X$ is orthogonal, meaning that $X^t X = I_p$, we get an explicit expression of the value $\hat{\beta}_{\lambda}^{\text{PEN}}$ of $\beta$ that maximizes this criterion, depending on the expression of the ordinary least square estimate denoted by $\hat{\beta}^{\text{OLS}}$ (see Tibshirani [1996]):

$$ \hat{\beta}_{\lambda}^{\text{PEN}} = \operatorname{sgn}(\hat{\beta}^{\text{OLS}})(|\hat{\beta}^{\text{OLS}}| - \lambda)|_{+}, \quad (14) $$

with $\operatorname{sgn}(x)$ equals 1 if $x > 0$, -1 if $x < 0$ and 0 if $x = 0$ and where the penalized estimate is defined as

$$ \hat{\beta}_{\lambda}^{\text{PEN}} = \operatorname*{arg\,max}_{\beta} \{\log g(\beta; Y) - \operatorname{pen}_{\lambda}(\beta)\} \quad (15) $$

and where

$$ \hat{\beta}^{\text{OLS}} = (X^t X)^{-1} X^t Y = X^t Y. \quad (16) $$

We detail here how to build in practice an orthogonal design matrix $X$ satisfying (H1). First we sample for $1 \le i \le N$ a matrix $U_i$ of size $J \times p$ with independent standard Gaussian components. We then centered these matrices by column and define the matrices $(V_i)$ of size $J \times p$ by $V_i = U_i - Wm_i$ where $m_i$ is the vector of size $p$ such that $m_{il} = \sum_{j=1}^{J} U_{ijl}/J$ for each $1 \le l \le p$. We denote by $T$ the matrix of size $p \times p$ defined by $T = V^t V$ and consider its Cholesky decomposition $T = LL^t$ where $L$ is a lower triangular matrix. Finally we define the design $X$ by $X = V(L^t)^{-1}$. It is easy to check that this design is orthogonal by construction and satisfies (H1) by construction.

Given this design $X$, we generate one dataset according to the model presented in equation (10) with $N = 100$, $J = 5$, $\gamma^2 = 4^2$, $\sigma^2 = 2^2$, $p = 200$. We choose for the components of $\beta$ the values 4 and -3 for the first and second components respectively, whereas the others are
fixed equal to zero. We run the proposed algorithm for three different values of the regularization parameters $\lambda$, namely one leading to select more variables than in the true model, one leading to select the same variables as in the true model, and one leading to select less variables than in the true model. For each configuration, we run the algorithm several times with different initializations chosen randomly. We present five trajectories of the estimates obtained with the algorithm for the fourth first components of $\beta$ to highlight the convergence toward the target in different settings of regularization: settings with more selected variables than those in the true support corresponding to over selection case on Figure 1, settings with the same selected variables than the variables in the true support corresponding to exact selection case on Figure 2, settings with less selected variables than the variables in the true support corresponding to under selection case on Figure 3. We emphasize again that in this example the true values of the penalized estimates for $\beta$ are calculated explicitely though Equation (14). These exact values are represented in dotted horizontal lines on the graphs. We observe in all settings that the algorithm's trajectories converge with high accuracy toward the targeted values, which differ depending on the regularization parameter value.

<figure><img src="image_3.png" /><figcaption>Figure 1: Trajectories of 5 runs of the algorithm performed on the same dataset for the four first components of $\beta$ with randomly chosen initializations and regularization parameter $\lambda$ equal to 0.6 leading to over selection of the support.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Trajectories of 5 runs of the algorithm performed on the same dataset for the four first components of $\beta$ with randomly chosen initializations and regularization parameter $\lambda$ equal to 1 leading to exact selection of the support.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: Trajectories of 5 runs of the algorithm performed on the same dataset for the four first components of $\beta$ with randomly chosen initializations and regularization parameter $\lambda$ equal to 2 leading to under selection of the support.</figcaption></figure>

## 5.2 Variable selection in high dimension and re-estimation analysis in a linear mixed effects model

For this section and the next one, we generate $n_{run}$ independent data sets, and we fit the corresponding model using the routine described in section 3.3 and get an estimate $\hat{\theta}_{i}^{MLE}$ for each run $i \in \{1, ..., n_{run}\}$. To compare the results across different scenarios, we calculate several metrics. We evaluate the sensitivity, specificity, and accuracy to study the selection capacity of the method. Sensitivity (Se) measures the proportion of true positives (TP) correctly identified, while Specificity (Sp) quantifies the proportion of true negatives (TN) correctly identified. Accuracy (Ac) represents the overall proportion of correctly classified instances, including both TP and TN. We abbreviated false negatives and false positives respectively by FN and FP.

$$ \text{Se} = \frac{TP}{TP + FN} \quad \text{Sp} = \frac{TN}{TN + FP} \quad \text{Ac} = \frac{TP + TN}{P + N} $$

We also evaluate the Relative Root Mean Square Errors (RRMSE) to measure the estimation quality of the method defined as:

$$ \text{RRMSE}(\hat{\theta}_{k}^{MLE}) = \sqrt{\frac{1}{n_{\text{run}}} \sum_{i=1}^{n_{\text{run}}} \frac{(\hat{\theta}_{i,k}^{\text{MLE}} - \theta_{k})^{2}}{\theta_{k}^{2}}} \quad \text{mse}(\hat{\theta}_{k}^{\text{MLE}}) = \frac{1}{n_{\text{run}}} \sum_{i=1}^{n_{\text{run}}} (\hat{\theta}_{i,k}^{\text{MLE}} - \theta_{k})^{2} $$

where $\hat{\theta}_{i,k}^{MLE}$ is the k-th component of $\hat{\theta}_{i}^{MLE}$.

In this section, we compare our method NLMEM-LASSO to the R package GLMMCLASSO (Groll and Tutz [2014]). We use the same procedure to select the regularization parameter. We consider the following linear model :

$$ \left\{ \begin{array}{l} Y_{i,j} = Z_{i,j1} + Z_{i2}t_j + \varepsilon_{i,j} \\ Z_{i,j,1} \stackrel{i.}{\sim} \mathcal{N}(\mu_1 + X_{i,j}\beta, \gamma_1^2) \\ Z_{i,2} \stackrel{i.i.d}{\sim} \mathcal{N}(\mu_2, \gamma_2^2) \\ \varepsilon_{i,j} \stackrel{i.i.d}{\sim} \mathcal{N}(0, \sigma^2) \end{array} \right. \qquad (17) $$

The model parameters are $\beta \in \mathbb{R}^p$, $\mu = (\mu_1, \mu_2) \in \mathbb{R}^2$, and $\gamma_1^2, \gamma_2^2, \sigma^2 \in (\mathbb{R}^*_{+})^3$. In this model, high-dimensional covariates model variability between individuals. We generated 100 data sets independently according to equation (17) for several scenarios. It is assumed that each individual is observed $J = 10$ times, at the same instants spread uniformly over a range between 0 and
1. We use the following values for the parameters: $\mu = (2, 5)$, $\gamma_1 = 1$, $\gamma_2 = 2$ and $\sigma = 1$. For each different value of $p$, we choose the vector $\beta$ such that the first three components are equal to $(\beta_1, \beta_2, \beta_3) = (8, -10, 20)$ and the rest are equal to zero. Additionally, we generate the matrix of covariates $X$ following a uniform distribution $X_{ijk} \sim U([-1, 1])$; $\forall i \in \{1, ..., N\}, \forall j \in \{1, ..., J\}, \forall k \in \{1, ..., p\}$.

By testing scenarios with increasing numbers of individuals $N \in \{100, 200\}$, we want to assess how both methods estimate more or less accurately when the sample size increases. In parallel we show different scenarios where the number of covariates increases $p \in \{200, 1000\}$. We assess how both methods, despite the high-dimensional context, manage to select the most explanatory variables. To evaluate the selection variable capacity of both methods, we present selection scores in Table 1 and proportions of well-selected variables regarding the support on Figure 4. Whatever the configuration, the proposed procedure NLMEM-LASSO and the procedure GLMMLASSO perform very similar regarding the considered indicators of average sensitivity, specificity and accuracy. Let's take a closer look at the performance of the two methods in terms of support selection: both generally find the right support. However, the procedure NLMEM-LASSO seems to over-select support less often than the procedure GLMMLASSO. Tables 2 and 6 (in Appendix) show the mean value over the 100 data sets of the Relative Root Mean Square Errors (RRMSE) for each model parameter, along with their estimates obtained after variable selection in the reduced model. Regarding the parameter's re-estimation in the reduced model, we emphasize that both methods give similar good results regarding estimation of regression parameters and of the residual variance; nevertheless the procedure NLMEM-LASSO also adequately estimates the variances of the random effects, unlike the procedure GLMMLASSO. This might be explained by the tendency of the GLMMLASSO method to over-select the support.

<table><caption>Table 1: Average of sensitivity (Se), specificity (Sp), accuracy (Ac), mean square error (mse) over 100 repetitions within the LMEM using the NLMEM-LASSO and GLMMLASSO methods.</caption><thead><tr><th rowspan="3"></th><th colspan="4">N = 100</th><th colspan="4">N = 200</th></tr><tr><th colspan="2">P = 200</th><th colspan="2">P = 1000</th><th colspan="2">P = 200</th><th colspan="2">P = 1000</th></tr><tr><th>NLMEM-LASSO</th><th>GLMMLASSO</th><th>NLMEM-LASSO</th><th>GLMMLASSO</th><th>NLMEM-LASSO</th><th>GLMMLASSO</th><th>NLMEM-LASSO</th><th>GLMMLASSO</th></tr></thead><tbody><tr><th>Ac</th><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><th>Se</th><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><th>Sp</th><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><th>mse</th><td>0.005</td><td>0.005</td><td>0.001</td><td>0.001</td><td>0.003</td><td>0.002</td><td>0.001</td><td>0.001</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 4: Proportion of correct support selection (solid bars) and support overselection (striped bars) using both procedures NLMEM-LASSO and GLMMLASSO for $N \in \{100, 200\}$ and for $P \in \{200, 1000\}$.</figcaption></figure>
<table><caption>Table 2: Average Relative Root Mean Square Errors (RRMSE) and parameter estimates over 100 repetitions for the Linear Mixed Effects Model (LMEM) where $\hat{\theta}_{\text{NLMEM-LASSO}}$ and $\hat{\theta}_{\text{GLMMLASSO}}$ are the estimates in the reduced model after variable selection using the NLMEM-LASSO and GLMMLASSO methods, respectively.</caption><thead><tr><th rowspan="3"></th><th rowspan="3">θ*</th><th colspan="4">P = 200</th><th colspan="4">P = 1000</th></tr><tr><th colspan="2">NLMEM-LASSO</th><th colspan="2">glmmLasso</th><th colspan="2">NLMEM-LASSO</th><th colspan="2">glmmLasso</th></tr><tr><th>θ̂</th><th>RRMSE</th><th>θ̂</th><th>RRMSE</th><th>θ̂</th><th>RRMSE</th><th>θ̂</th><th>RRMSE</th></tr></thead><tbody><tr><td rowspan="8">N = 100</td><td>μ<sub>1</sub></td><td>2.00</td><td>1.99</td><td>3.15</td><td>1.98</td><td>6.13</td><td>1.99</td><td>3.07</td><td>1.98</td><td>6.14</td></tr><tr><td>μ<sub>2</sub></td><td>5.00</td><td>5.02</td><td>1.92</td><td>5.03</td><td>3.70</td><td>5.02</td><td>1.91</td><td>5.03</td><td>3.68</td></tr><tr><td>γ<sub>1</sub><sup>2</sup></td><td>1.00</td><td>0.99</td><td>19.96</td><td>0.61</td><td>39.74</td><td>0.98</td><td>19.28</td><td>0.61</td><td>39.40</td></tr><tr><td>γ<sub>2</sub><sup>2</sup></td><td>4.00</td><td>3.93</td><td>17.10</td><td>0.88</td><td>78.10</td><td>3.90</td><td>17.01</td><td>0.88</td><td>78.12</td></tr><tr><td>σ<sup>2</sup></td><td>1.00</td><td>1.00</td><td>5.28</td><td>0.98</td><td>5.07</td><td>1.00</td><td>5.29</td><td>0.98</td><td>5.15</td></tr><tr><td>β<sub>1</sub></td><td>8.00</td><td>8.03</td><td>7.20</td><td>8.05</td><td>7.20</td><td>8.05</td><td>7.23</td><td>8.06</td><td>7.74</td></tr><tr><td>β<sub>2</sub></td><td>-10.00</td><td>-9.95</td><td>5.49</td><td>-9.97</td><td>5.81</td><td>-9.87</td><td>6.04</td><td>-9.92</td><td>6.30</td></tr><tr><td>β<sub>3</sub></td><td>20.00</td><td>19.99</td><td>2.74</td><td>19.98</td><td>2.89</td><td>19.91</td><td>2.65</td><td>19.86</td><td>2.80</td></tr><tr><td rowspan="8">N = 200</td><td>μ<sub>1</sub></td><td>2.00</td><td>2.00</td><td>1.95</td><td>1.99</td><td>3.65</td><td>2.00</td><td>1.83</td><td>2.00</td><td>3.65</td></tr><tr><td>μ<sub>2</sub></td><td>5.00</td><td>5.00</td><td>1.85</td><td>5.01</td><td>3.49</td><td>5.00</td><td>1.84</td><td>5.01</td><td>3.49</td></tr><tr><td>γ<sub>1</sub><sup>2</sup></td><td>1.00</td><td>0.99</td><td>14.63</td><td>0.61</td><td>39.12</td><td>0.98</td><td>14.09</td><td>0.61</td><td>38.97</td></tr><tr><td>γ<sub>2</sub><sup>2</sup></td><td>4.00</td><td>4.06</td><td>12.77</td><td>0.89</td><td>77.83</td><td>4.04</td><td>12.41</td><td>0.89</td><td>77.80</td></tr><tr><td>σ<sup>2</sup></td><td>1.00</td><td>1.00</td><td>4.00</td><td>0.98</td><td>4.07</td><td>1.00</td><td>4.09</td><td>0.98</td><td>3.94</td></tr><tr><td>β<sub>1</sub></td><td>8.00</td><td>8.03</td><td>4.52</td><td>8.02</td><td>4.59</td><td>7.91</td><td>5.01</td><td>7.92</td><td>5.21</td></tr><tr><td>β<sub>2</sub></td><td>-10.00</td><td>-9.94</td><td>3.70</td><td>-9.96</td><td>3.85</td><td>-9.88</td><td>4.12</td><td>-9.90</td><td>4.11</td></tr><tr><td>β<sub>3</sub></td><td>20.00</td><td>20.01</td><td>1.93</td><td>20.01</td><td>1.92</td><td>20.02</td><td>2.13</td><td>20.02</td><td>2.13</td></tr></tbody></table>

## 5.3 Variable selection in high-dimension and re-estimation analysis in a logistic non-linear mixed-effects model

We study in this section the model presented in the first example in Equation (1), where we explain a part of the variability of the logistic’s midpoint by the high-dimensional covariates. The model can be written as follows:

$$ 
\left\{ 
\begin{aligned}
Y_{i,j} &= \frac{Z_{i1}}{1 + \exp\left(-\frac{t_{ij} - Z_{i2}}{\alpha}\right)} + \varepsilon_{i,j} \\
Z_{i,1} &\stackrel{i.i.d}{\sim} \mathcal{N}(\mu_1, \gamma_1^2) \\
Z_{i,2} &\stackrel{i.}{\sim} \mathcal{N}(\mu_2 + X_i\beta, \gamma_2^2) \\
\varepsilon_{i,j} &\stackrel{i.i.d}{\sim} \mathcal{N}(0, \sigma^2)
\end{aligned}
\right. 
\qquad (18) 
$$

The model parameters are $\alpha \in \mathbb{R}$, $\beta \in \mathbb{R}^p$, $\mu = (\mu_1, \mu_2) \in \mathbb{R}^2$, and $\gamma_1^2, \gamma_2^2, \sigma^2 \in (\mathbb{R}_+^*)^3$. The quantity $Z_{i1}$ represents the asymptotic maximum value of the curve, $Z_{i2}$ represents the value of the logistic’s midpoint, and $\alpha$ represents the logistic growth rate. The two first parameters are modeled at the individual level of the mixed model, whereas the third one is modeled at the population level. Note that due to the presence of the fixed effect $\alpha$, the joint density of $(Y_i, Z_i)$ does not belong to the curved exponential family as defined in Delyon et al. [1999]. We generate 100 data sets independently according to equation (18) for several scenarios. It is assumed that each individual was observed $J = 15$ times, at the same instants equally spaced over a range between 150 and 3000. We use the following values for the parameters: $\mu = (200, 1200)$, $\gamma_1^2 = 7^2$, $\gamma_2^2 = 30^2$, $\alpha = 300$ and $\sigma^2 = 30$. For each different value of $p$, we choose the vector $\beta$ such that the first three components are equal to $(100, 200, -300)$ and the rest equal to zero. Furthermore, we generate the matrix of covariates $X$ with $N$ rows and $p$ columns, following a uniform distribution $X_{i,k} \sim U([-1, 1])$; $\forall i \in \{1, ..., N\}, \forall k \in \{1, ..., p\}$.

As in the previous analysis, we present scenarios with an increasing number of individuals $N \in \{100, 200\}$ and different covariates size $p \in \{200, 500, 1000\}$. We first assess the selection variable capacity of our method, by presenting selection scores in Table 3 and proportions of well-selected variables regarding the support on Figure 5. We also indicate over-selection cases.
The later shows that as the number of variables *p* increases, support selection performance deteriorates; in particular, a higher proportion of support over-selection is observed. However, for a fixed number of variables *p*, the support selection performance improves as the number of individuals *N* increases.

Table 4 summarizes the results of the estimation of the model parameters in the reduced model after variable selection step, for *N* = 100 and *N* = 200. We observe that the estimates of the parameters are accurate in all configurations. As expected, the relative root mean square errors decrease with *N* for a given covariate size *p*. For *N* = 100, we observe that the relative root mean square errors generally increase with *p*, which can be explained by the fact that increasing the size of the covariates increases the difficulty of the selection procedure. The phenomenon is less pronounced as the sample size increases, reflecting the asymptotically expected behavior.

<table><caption>Table 3: Average sensitivity (Se), specificity (Sp), accuracy (Ac), and mean square error (mse) using the procedure NLMEM-LASSO over 100 repetitions for the logistic non-linear mixed effects model</caption><thead><tr><th rowspan="2"></th><th colspan="3">N = 100</th><th colspan="3">N = 200</th></tr><tr><th>P = 200</th><th>P = 500</th><th>P = 1000</th><th>P = 200</th><th>P = 500</th><th>P = 1000</th></tr></thead><tbody><tr><td>Ac</td><td>0.997</td><td>0.999</td><td>0.999</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><td>Se</td><td>0.980</td><td>0.927</td><td>0.933</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><td>Sp</td><td>0.997</td><td>0.999</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><td>mse</td><td>0.819</td><td>0.851</td><td>0.387</td><td>0.127</td><td>0.057</td><td>0.026</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 5: Proportion of correct support selection (solid bars) and overselection (striped bars) in the logistic non-linear mixed effects model using the procedure NLMEM-LASSO for $N \in \{100, 200\}$ and for $P \in \{200, 500, 1000\}$.</figcaption></figure>

As an illustration of the behavior of the optimization algorithm in the reduced model, Figure 6 displays the estimated parameter as a function of iterations during the execution of the adaptive stochastic gradient algorithm. We can observe the good convergence of the algorithm toward a limit which is very closed to the true parameter value, as expected for the maximum likelihood estimate.
<table><caption>Table 4: Average Relative Root Mean Square Errors (RRMSE) and parameter estimates over 100 repetitions for the logistic non-linear mixed effects model where $\hat{\theta}$ is the estimate in the reduced model after variable selection using the NLMEM-LASSO.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">θ*</th><th colspan="2">P = 200</th><th colspan="2">P = 500</th><th colspan="2">P = 1000</th></tr><tr><th>θ̂</th><th>RRMSE</th><th>θ̂</th><th>RRMSE</th><th>θ̂</th><th>RRMSE</th></tr></thead><tbody><tr><td rowspan="7">N = 100</td><td>μ<sub>1</sub></td><td>200.00</td><td>199.90</td><td>0.38</td><td>199.90</td><td>0.40</td><td>199.93</td><td>0.39</td></tr><tr><td>μ<sub>2</sub></td><td>1200.00</td><td>1200.02</td><td>0.35</td><td>1199.87</td><td>0.38</td><td>1200.20</td><td>0.37</td></tr><tr><td>γ<sub>1</sub><sup>2</sup></td><td>49.00</td><td>48.46</td><td>17.97</td><td>48.20</td><td>17.68</td><td>48.31</td><td>17.82</td></tr><tr><td>γ<sub>2</sub><sup>2</sup></td><td>900.00</td><td>861.29</td><td>38.85</td><td>1062.15</td><td>63.75</td><td>1057.36</td><td>67.47</td></tr><tr><td>α</td><td>300.00</td><td>299.86</td><td>0.77</td><td>299.99</td><td>0.71</td><td>300.16</td><td>0.76</td></tr><tr><td>σ<sup>2</sup></td><td>30.00</td><td>30.06</td><td>5.74</td><td>30.09</td><td>5.82</td><td>30.09</td><td>5.77</td></tr><tr><td>β<sub>0</sub></td><td>120.00</td><td>119.90</td><td>3.47</td><td>119.38</td><td>3.85</td><td>118.44</td><td>3.92</td></tr><tr><td rowspan="5">N = 200</td><td>β<sub>1</sub></td><td>70.00</td><td>69.74</td><td>5.86</td><td>69.27</td><td>6.80</td><td>69.68</td><td>6.75</td></tr><tr><td>β<sub>2</sub></td><td>40.00</td><td>37.81</td><td>26.31</td><td>31.32</td><td>47.85</td><td>31.56</td><td>45.43</td></tr><tr><td>μ<sub>1</sub></td><td>200.00</td><td>199.97</td><td>0.26</td><td>199.99</td><td>0.25</td><td>199.99</td><td>0.26</td></tr><tr><td>μ<sub>2</sub></td><td>1200.00</td><td>1199.75</td><td>0.23</td><td>1199.94</td><td>0.22</td><td>1199.96</td><td>0.23</td></tr><tr><td>γ<sub>1</sub><sup>2</sup></td><td>49.00</td><td>48.76</td><td>11.83</td><td>48.74</td><td>12.06</td><td>48.71</td><td>11.95</td></tr><tr><td></td><td>γ<sub>2</sub><sup>2</sup></td><td>900.00</td><td>894.39</td><td>14.65</td><td>891.41</td><td>14.39</td><td>890.88</td><td>14.16</td></tr><tr><td></td><td>α</td><td>300.00</td><td>299.82</td><td>0.50</td><td>299.91</td><td>0.50</td><td>299.91</td><td>0.51</td></tr><tr><td></td><td>σ<sup>2</sup></td><td>30.00</td><td>30.10</td><td>3.76</td><td>30.10</td><td>3.97</td><td>30.12</td><td>3.86</td></tr><tr><td></td><td>β<sub>0</sub></td><td>120.00</td><td>119.87</td><td>2.35</td><td>119.88</td><td>2.53</td><td>119.35</td><td>2.42</td></tr><tr><td></td><td>β<sub>1</sub></td><td>70.00</td><td>69.97</td><td>4.50</td><td>70.37</td><td>4.50</td><td>69.13</td><td>4.36</td></tr><tr><td></td><td>β<sub>2</sub></td><td>40.00</td><td>40.53</td><td>6.89</td><td>39.39</td><td>7.11</td><td>40.29</td><td>6.59</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 6: Parameter estimates in the reduced model across the iterations of the adaptive stochastic gradient algorithm within the logistic mixed model. Dotted lines: true values</figcaption></figure>

## 5.4 Comparison of integrated and two steps approaches for variable selection in a non-linear mixed effects model

In this section, following the study performed in [Naveau et al., 2024], we compare the performance for variable selection of three procedures. The first procedure consists in two steps, namely first in treating the observations of each individual separately by adjusting one non-linear
<figure><img src="image_9.png" /><figcaption>Figure 7: Left: Regularization path, i.e. values of estimates of $\beta$, in solid line; Right: values of the eBIC criterion in dotted line; the dotted vertical line represents the chosen regularization value.</figcaption></figure>

regression by individual to get estimates of its parameters and second performing a high dimensional linear regression between these estimates and the covariates to select the more relevant. The second procedure consists also in two steps, namely first adjusting a non linear mixed model to estimate the individual parameters of all individuals simultaneously, and second performing a high dimensional linear regression between these estimates and the covariates to select the more relevant. The third procedure is the one proposed in section 3.3 and integrates both steps in one single step.

We carry out the first procedure called NLM-GLMNET using the method of least squares implemented in the NLM R function [Dennis and Schnabel, 1996] to estimate individual parameters and then perform variable selection with the LASSO procedure using the GLMNET package [Friedman et al., 2010]. This last step is performed using cross-validation in order to select an optimal regularization parameter chosen such that the mean squared error is within one standard deviation of the minimum mse as recommended by Hastie et al. [2009]. The second procedure called SAEMIX-GLMNET uses the Stochastic Approximation Expectation Maximization (SAEM) algorithm [Kuhn and Lavielle, 2004] implemented in the R package SAEMIX [Comets et al., 2017] to adjust the population mixed model and obtain the individual parameters estimates. The variable selection is then performed as in the first procedure. The third procedure called NLMEM-LASSO is performed as detailed in section 3.3.

For this comparison, we use the pharmacodynamic non linear mixed effect model, presented in the second example in Equation (2). We generate 100 independent datasets according to the following model with $N = 200$ individuals and $p = 500$ covariates :

$$ 
\left\{ 
\begin{aligned}
Y_{ij} &= \frac{DZ_{i1}}{V(Z_{i1} - Z_{i2})} \left( \exp\left(-\frac{Z_{i2}}{V}t_{ij}\right) - \exp(-Z_{i1}t_{ij}) \right) + \varepsilon_{ij}, \\
Z_i &= \mu + X_i\beta + \xi_i, \quad \xi_i \stackrel{\text{i.i.d}}{\sim} N_2\left(0, \begin{pmatrix} \gamma_1^2 & 0 \\ 0 & \gamma_2^2 \end{pmatrix}\right), \quad \varepsilon_{ij} \stackrel{\text{i.i.d}}{\sim} N(0, \sigma^2),
\end{aligned}
\right. 
\quad (19) $$

The model parameters are $\mu = (\mu_1, \mu_2) \in \mathbb{R}^2$, $\beta \in \mathcal{M}_{p,2}$ and $(\gamma_1^2, \gamma_2^2, \sigma^2) \in (\mathbb{R}_+^*)^3$. The quantity $V$ is the distribution volume, whereas $D$ is the administered dose. These two parameters are supposed to be known. The quantity $Z_{i1}$ represents the absorption rate, $Z_{i2}$ represents the elimination rate. We consider the following observations times [0.05, 0.15, 0.25, 0.4, 0.5, 0.8, 1, 2, 7, 12, 24, 40]. These times are chosen to cover a wide range of the pharmacokinetic profile.

We generate datasets with a censored rate $\rho_{censored}$ varying in $\{0\%, 10\%, 20\%, 30\%, 40\%\}$. More precisely a proportion $\rho_{censored}$ of the individuals are observed at only the first three times
and supposed to leave the study after that. The other individuals are observed at all times. We use the following values for the parameters: $\mu = (6, 8)$, $\gamma_1 = 0.2$, $\gamma_2 = 0.1$ and $\sigma = 10^{-3}$. The high-dimensional parameter $\beta$ is chosen such that:

$$ \beta = \begin{pmatrix} 3 & 2 & 1 & 0 & 0 & 0 & \dots & 0 \\ 0 & 0 & 3 & 2 & 1 & 0 & \dots & 0 \end{pmatrix}^T $$

The individual covariates $X_i$ are generated from a binomial distribution with a probability of success of 0.2 and are standardized. We emphasize here that the second procedure SAEMIX-GLMNET is carried out on the following mixed model without high-dimensional covariates presented in Equation (20) to estimate the individual parameter $Z_i$ with the package SAEMIX in the first step.

$$ \left\{ \begin{aligned} Y_{ij} &= \frac{DZ_{i1}}{V(Z_{i1} - Z_{i2})} \left( \exp\left(-\frac{Z_{i2}}{V}t_{ij}\right) - \exp(-Z_{i1}t_{ij}) \right) + \varepsilon_{ij}, \\ Z_i &= \mu + \xi_i, \quad \xi_i \stackrel{i.i.d}{\sim} N_2 \left(0, \begin{pmatrix} \gamma_1^2 & 0 \\ 0 & \gamma_2^2 \end{pmatrix}\right), \quad \varepsilon_{ij} \stackrel{i.i.d}{\sim} N(0, \sigma^2), \end{aligned} \right. \quad (20) $$

Table 5 shows the mean estimations errors (mee) in the first step of NLM-GLMNET and SAEMIX-GLMNET methods over the 100 datasets. As expected, both methods estimate the first parameter $Z_1$ with a smaller error than the second parameter $Z_2$. Moreover, both struggle to estimate the second parameter $Z_2$ when the censoring rate increases. However the mixed-effect approach handle the missing data better than the individual by individual approach as the error of the second parameter $Z_2$ is smaller for SAEMIX-GLMNET.

<table><caption>Table 5: Mean estimation errors (mee) of the parameters $Z_1$ and $Z_2$ for the two-step approaches NLM-GLMNET and SAEMIX-GLMNET.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Method</th><th colspan="5">Censored data rate</th></tr><tr><th>0%</th><th>10%</th><th>20%</th><th>30%</th><th>40%</th></tr></thead><tbody><tr><td rowspan="2">mee(Z<sub>1</sub>)</td><td>NLM-GLMNET</td><td>0.088</td><td>0.096</td><td>0.104</td><td>0.112</td><td>0.120</td></tr><tr><td>SAEMIX-GLMNET</td><td>0.102</td><td>0.107</td><td>0.111</td><td>0.116</td><td>0.120</td></tr><tr><td rowspan="2">mee(Z<sub>2</sub>)</td><td>NLM-GLMNET</td><td>0.143</td><td>0.657</td><td>1.188</td><td>1.701</td><td>2.222</td></tr><tr><td>SAEMIX-GLMNET</td><td>0.150</td><td>0.407</td><td>0.663</td><td>0.905</td><td>1.151</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>Figure 8: Proportion of correct support selection (solid bars) and overselection (striped bars) for the NLMEM-LASSO, SAEMIX-GLMNET and NLM-GLMNET methods under increasing censoring rates.</figcaption></figure>
Figure 8 shows the comparison of variable selection performance between the proposed NLMEM-LASSO method and the two-step approaches NLM-GLMNET and SAEMIX-GLMNET. On one hand variable selection performs well on the first parameter $Z_1$ and their associated covariates, and is not impacted by censoring. On the other hand, variable selection on the second parameter is deteriorated by censorship in all procedures. Figure 8 highlights the interests of using procedures based on mixed effects models and integrated approach. The two-step approach SAEMIX-GLMNET is more successful than the two-step approach NLM-GLMNET in finding exact support for different censoring levels. This can be explained by the lower estimation errors obtained by SAEMIX-GLMNET observed in Table 5. However, unlike NLMEM-LASSO, it overselects explanatory variables. In all cases, the NLMEM-LASSO method outperforms the two-step approaches.

# 6 Conclusion and perspectives

In this work, we jointly address variable selection in high-dimension and parameter estimation in general mixed-effect model. We emphasize that the proposed procedure called NLMEM-LASSO can handle well linear and non-linear mixed effects models, without requiring them to belong to the curved exponential family. We perform an adaptive weighted proximal stochastic gradient algorithm to deal simultaneously with the latent variables and the LASSO penalty when maximizing the penalized likelihood. We highlight the well behavior of this optimization algorithm in various settings through simulation. We use the eBIC criterion to perform the choice of the regularization parameter leading to the selection of the variables corresponding to the reduced model. Our integrated methodology for variable selection has been thoroughly evaluated in simulation study to demonstrate its performance in linear and non-linear mixed models. Its performance has also been compared to the ones of the procedure GLMMLASSO in the case of linear mixed models. Moreover, we emphasize that the proposed procedure can handle a wide range of modeling choices for the random effects and the error term, including in particular correlation structure for the random effects.

There are several interesting perspectives to this work. The first one consists in studying the effect of the presence of correlation between the variables involved in the selection which is known to deteriorate the accuracy of the selection procedure. The second one consists in studying theoretical properties of the proposed procedure, in particular convergence of the optimization algorithm combining weighted proximal operator and adaptive preconditioned gradient descent and support consistency guarantees. Finally, it would be of great interest to assort the re-estimation results in the reduced model with confidence region, by studying theoretical properties of the estimator and prediction after the variable selection step as done in post model selection inference.

# Funding and Acknowledgements

This work was funded by the https://stat4plant.mathnum.inrae.fr/(Stat4Plant) project ANR-20-CE45-0012.

# References

* Charlotte Baey, Amélie Mathieu, Alexandra Jullien, Samis Trevezas, and Paul-Henry Cournède. Mixed-effects estimation in dynamic models of plant growth for the assessment of inter-individual variability. *Journal of agricultural, biological and environmental statistics*, 23: 208–232, 2018.
Charlotte Baey, Maud Delattre, Estelle Kuhn, Jean-Benoist Leger, and Sarah Lemler. Efficient preconditioned stochastic gradient descent for estimation in latent variable models. In *Proceedings of the 40th International Conference on Machine Learning*, ICML’23. JMLR.org, 2023.

Stephen Becker and Jalal Fadili. A quasi-newton proximal splitting method. In F. Pereira, C.J. Burges, L. Bottou, and K.Q. Weinberger, editors, *Advances in Neural Information Processing Systems*, volume 25. Curran Associates, Inc., 2012. URL https://proceedings.neurips.cc/paper_files/paper/2012/file/e034fb6b66aacc1d48f445ddfb08da98-Paper.pdf.

Julie Bertrand and David J Balding. Multiple single nucleotide polymorphism analysis using penalized regression in nonlinear mixed-effect pharmacokinetic models. *Pharmacogenetics and genomics*, 23(3):167–174, 2013.

Sahir R Bhatnagar, Yi Yang, Tianyuan Lu, Erwin Schurr, JC Loredo-Osti, Marie Forest, Karim Oualkacha, and Celia MT Greenwood. Simultaneous snp selection and adjustment for population structure in high dimensional prediction models. *PLoS genetics*, 16(5):e1008766, 2020.

Olivier Cappé, Eric Moulines, and Tobias Rydén. Springer series in statistics, 2005.

George H.-G. Chen and R. Tyrrell Rockafellar. Convergence rates in forward-backward splitting. *SIAM J. Optim.*, 7:421–444, 1997. URL https://api.semanticscholar.org/CorpusID:7104716.

Jiahua Chen and Zehua Chen. Extended bayesian information criteria for model selection with large model spaces. *Biometrika*, 95(3):759–771, 2008. ISSN 00063444, 14643510. URL http://www.jstor.org/stable/20441500.

Emilie Chouzenoux, Jean-Christophe Pesquet, and Audrey Repetti. Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function. *J. Optim. Theory Appl.*, 162(1):107–132, July 2014. ISSN 0022-3239. doi: 10.1007/s10957-013-0465-7. URL https://doi.org/10.1007/s10957-013-0465-7.

Emmanuelle Comets, Audrey Lavenu, and Marc Lavielle. Parameter estimation in nonlinear mixed effect models using saemix, an r implementation of the saem algorithm. *Journal of Statistical Software*, 80(3):1–41, 2017. doi: 10.18637/jss.v080.i03. URL https://www.jstatsoft.org/index.php/jss/article/view/v080i03.

Marie Davidian. *Nonlinear models for repeated measurement data*. Routledge, 2017.

Vianney Debavelaere and Stéphanie Allassonnière. On the curved exponential family in the stochastic approximation expectation maximization algorithm. *ESAIM: Probability & Statistics*, 25, 2021.

Maud Delattre, Marc Lavielle, and Marie-Anne Poursat. A note on BIC in mixed-effects models. *Electronic Journal of Statistics*, 8(1):456–475, 2014. doi: 10.1214/14-EJS890. URL https://doi.org/10.1214/14-EJS890.

Bernard Delyon, Marc Lavielle, and Eric Moulines. Convergence of a stochastic approximation version of the EM algorithm. *The Annals of Statistics*, 27(1):94–128, 1999. ISSN 00905364. URL http://www.jstor.org/stable/120120. Publisher: Institute of Mathematical Statistics.

Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the em algorithm. *Journal of the royal statistical society: series B (methodological)*, 39(1):1–22, 1977.
J. E. Dennis and Robert B. Schnabel. *Numerical Methods for Unconstrained Optimization and Nonlinear Equations*. Society for Industrial and Applied Mathematics, 1996. doi: 10.1137/1.9781611971200. URL https://epubs.siam.org/doi/abs/10.1137/1.9781611971200.

John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. *Journal of Machine Learning Research*, 12(61):2121–2159, 2011. URL http://jmlr.org/papers/v12/duchi11a.html.

Gersende Fort, Edouard Ollier, and Adeline Samson. Stochastic proximal-gradient algorithms for penalized mixed models. *Statistics and Computing*, 29:231–253, 2019.

Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Regularization paths for generalized linear models via coordinate descent. *Journal of Statistical Software*, 33(1):1–22, 2010. doi: 10.18637/jss.v033.i01.

Andreas Groll and Gerhard Tutz. Variable selection for generalized linear mixed models by l1-penalized estimation. *Statistics and Computing*, 24(2):137–154, 2014. ISSN 1573-1375. doi: 10.1007/s11222-012-9359-z. URL https://doi.org/10.1007/s11222-012-9359-z.

Trevor Hastie, Robert Tibshirani, and Jerome Friedman. *The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition (Springer Series in Statistics)*. 02 2009. ISBN 0387848576.

Benjamin Heuclin, Frédéric Mortier, Catherine Trottier, and Marie Denis. Bayesian varying coefficient model with selection: An application to functional mapping. *Journal of the Royal Statistical Society: Series C (Applied Statistics)*, 70, 11 2020. doi: 10.1111/rssc.12447.

Lukas Meier Jürg Schelldorfer and Peter Bühlmann. Glmmlasso: An algorithm for high-dimensional generalized linear mixed models using l1-penalization. *Journal of Computational and Graphical Statistics*, 23(2):460–477, 2014. doi: 10.1080/10618600.2013.773239. URL https://doi.org/10.1080/10618600.2013.773239.

Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization, 2017. URL https://arxiv.org/abs/1412.6980.

Estelle Kuhn and Marc Lavielle. Coupling a stochastic approximation version of em with an mcmc procedure. *ESAIM: PS*, 8:115–131, 2004. doi: 10.1051/ps:2004007. URL https://doi.org/10.1051/ps:2004007.

Jean-Benoist Leger. Parametrization cookbook: A set of bijective parametrizations for using machine learning methods in statistical inference. *arXiv preprint arXiv:2301.08297*, 2023.

Jean-Michel Marin and Christian P Robert. *Bayesian core: a practical approach to computational Bayesian statistics*. Springer, 2007.

Jean Jacques Moreau. Fonctions convexes duales et points proximaux dans un espace hilbertien. *Comptes rendus hebdomadaires des séances de l'Académie des sciences*, 255:2897–2899, 1962.

Marion Naveau, Guillaume Kon Kam King, Renaud Rincent, Laure Sansonnet, and Maud Delattre. Bayesian high-dimensional covariate selection in non-linear mixed-effects models using the saem algorithm. *Statistics and Computing*, 34(1):53, 2024.

Shu Kay Ng, Thriyambakam Krishnan, and Geoffrey J McLachlan. The em algorithm. *Handbook of computational statistics: concepts and methods*, pages 139–172, 2012.

Edouard Ollier. Fast selection of nonlinear mixed effect models using penalized likelihood. *Computational Statistics & Data Analysis*, 167:107373, 2022.
José Pinheiro and Douglas Bates. *Mixed-effects models in S and S-PLUS*. Springer science & business media, 2006.

R. Tyrrell Rockafellar. Monotone operators and the proximal point algorithm. *SIAM Journal on Control and Optimization*, 14(5):877-898, 1976. ISSN 0363-0129, 1095-7138. doi: 10.1137/0314056. URL http://epubs.siam.org/doi/10.1137/0314056.

Jürg Schelldorfer, Lukas Meier, and Peter Bühlmann. Glmmlasso: an algorithm for high-dimensional generalized linear mixed models using $l_1$-penalization. *Journal of Computational and Graphical Statistics*, 23(2):460-477, 2014.

Robert Tibshirani. Regression shrinkage and selection via the lasso. *Journal of the Royal Statistical Society: Series B (Methodological)*, 58(1):267-288, 1996. ISSN 00359246. doi: 10.1111/j.2517-6161.1996.tb02080.x. URL https://onlinelibrary.wiley.com/doi/10.1111/j.2517-6161.1996.tb02080.x.

T. Tieleman. Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude, 2012. URL https://cir.nii.ac.jp/crid/1370017282431050757.

P. Tseng. A modified forward-backward splitting method for maximal monotone mapping. *Siam Journal on Control and Optimization - SIAM*, 38, 01 2000.
# Appendix

## Appendix A. Simulation study results for LMEM

<table><caption>Table 6: Average Relative Root Mean Square Errors (RRMSE) and parameter estimates over 100 repetitions for the Linear Mixed Effects Model (LMEM) with *p* = 500 where $\hat{\theta}_{\text{NMEM-LASSO}}$ and $\hat{\theta}_{\text{GLMMLASSO}}$ are the estimates in the reduced model after variable selection using the NMEM-LASSO and GLMMLASSO methods, respectively.</caption><thead><tr><th rowspan="3"></th><th rowspan="3">θ*</th><th colspan="8">P = 500</th></tr><tr><th colspan="4">N = 100</th><th colspan="4">N = 200</th></tr><tr><th colspan="2">NMEM-LASSO</th><th colspan="2">GLMMLASSO</th><th colspan="2">NMEM-LASSO</th><th colspan="2">GLMMLASSO</th></tr><tr><th></th><th></th><th>θ̂<sub>NMEM-LASSO</sub></th><th>RRMSE</th><th>θ̂<sub>GLMMLASSO</sub></th><th>RRMSE</th><th>θ̂<sub>NMEM-LASSO</sub></th><th>RRMSE</th><th>θ̂<sub>GLMMLASSO</sub></th><th>RRMSE</th></tr></thead><tbody><tr><td>μ<sub>1</sub></td><td>2.00</td><td>1.98</td><td>6.36</td><td>1.98</td><td>6.15</td><td>1.98</td><td>3.97</td><td>2.00</td><td>3.64</td></tr><tr><td>μ<sub>2</sub></td><td>5.00</td><td>5.03</td><td>4.03</td><td>5.03</td><td>3.69</td><td>5.01</td><td>2.70</td><td>5.01</td><td>3.49</td></tr><tr><td>γ<sub>1</sub><sup>2</sup></td><td>1.00</td><td>1.00</td><td>20.55</td><td>0.61</td><td>39.38</td><td>1.01</td><td>13.92</td><td>0.61</td><td>39.53</td></tr><tr><td>γ<sub>2</sub><sup>2</sup></td><td>4.00</td><td>3.93</td><td>17.32</td><td>0.88</td><td>78.11</td><td>4.03</td><td>12.72</td><td>0.88</td><td>77.92</td></tr><tr><td>σ<sup>2</sup></td><td>1.00</td><td>0.99</td><td>4.99</td><td>0.98</td><td>5.00</td><td>1.01</td><td>3.91</td><td>0.98</td><td>3.97</td></tr><tr><td>β<sub>1</sub></td><td>8.00</td><td>8.14</td><td>7.33</td><td>8.13</td><td>7.45</td><td>8.08</td><td>5.35</td><td>7.98</td><td>5.45</td></tr><tr><td>β<sub>2</sub></td><td>-10.00</td><td>-9.98</td><td>6.05</td><td>-10.00</td><td>6.63</td><td>-9.99</td><td>3.91</td><td>-9.96</td><td>4.22</td></tr><tr><td>β<sub>3</sub></td><td>20.00</td><td>19.98</td><td>3.30</td><td>19.97</td><td>3.42</td><td>19.98</td><td>1.91</td><td>19.97</td><td>2.14</td></tr></tbody></table>