Some contributions on variable selection in nonlinear mixed-effects
models.

Maud Delattre
Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France.
https://madelattre.github.io/1Introduction and model formulation

2

Covariate selection
On the use of BIC in NLMEM
High-dimensional variable selection in NLMEM

Maud Delattre

2 / 33

1

Introduction and model formulation

2

Covariate selection
On the use of BIC in NLMEM
High-dimensional variable selection in NLMEM

Maud Delattre

3 / 33

Motivating example : pharmacokinetic (PK) studies

Clinical issues : define the optimal modalities of drug administration
Statistical issues :describe the time course of drug absorption, distribution, metabolism, and excretionin the bodyestimate key PK parameters (such as absorption or elimination rates)
evaluate the influence of covariates (such as age or weight)

Maud Delattre

4 / 33Other example : growth of soybean plants (nlme)

Biological issues : optimize crops
Statistical issues :describe the growth curves and estimate meaningful quantities (growth rate, growthplateau, etc.)
characterize the part due to the genotype and the part due to the year of culture in
the observed differences
be able to make predictions

Maud Delattre

5 / 33

Common features that require special consideration in the statistical models
repeated measurements from several subjects
similarly-shaped individual profiles
peak/maximum value, rise, decay vary from one individual to another

Maud Delattre

6 / 33

Hierarchical formulation of the mixed-effects model(Pinheiro J. C. and Bates D. M., 2000), (Lavielle M., 2014)
1)

Description of the intra-individual variability
yij ∼ f (tij , ϕi ) + g(tij , ϕi )eij , eij ∼ N (0, 1) , j = 1 : ni , i = 1 : N
ind.

i.i.d.

where
yij : jth observation for individual if , g : (nonlinear, often mechanistic) functions governing within-individual
behavior/residual variance
ϕi ∈ Rd : parameters for individual i
2)

Description of the inter-individual variability
ϕi = µ + βCi + ηi , ηi ∼ N (0, Ω)
i.i.d.

Ci ∈ Rq : vector of covariates for individual iµ ∈ Rd : intercept vector
β : d × q matrix of fixed-effects
θ = (µ, β, Ω)Maud Delattre

7 / 33

PK example
1)(Data)

Observations :yij = fij + fij eij , fij = f (tij , Di , ϕi ) , eij ∼ N (0, 1)
i.i.d.

f (t, D, ka , V , Cl) =
2)


D ka
e −(Cl/V ) t − e −ka t
V ka − Cl

Individual parameters :Maud Delattre

log ka,i=µka ,log Vi=µV + βV ci + ηV ,i , ηV ,i ∼ N (0, ωV2 )log Cli=2µCl + ηCl,i , ηCl,i ∼ N (0, ωCl)

i.i.d.

i.i.d.

8 / 33

Typical questions
1) Estimation of the population parameter θComputation ?
Theoretical properties ?
2) Computation of the Fisher information matrix, computation of the likelihood, . . .
3) Model building, . . .

Difficult questions because of non linearity + latent variable structure of the model

Maud Delattre

9 / 33

1

Introduction and model formulation

2

Covariate selection
On the use of BIC in NLMEM
High-dimensional variable selection in NLMEM

1Joint work with Marie-Anne Poursat and Marc Lavielle published in
A note on BIC in mixed effects models.. (2014), Electronic Journal of Statistics 8(1) p.
456-475

2

Joint work with Marie-Anne Poursat published in
An iterative algorithm for joint covariate and random effect selection in mixed effects
models. (2020), The International Journal of Biostatistics 16(2), 20190082

Maud Delattre

10 / 33

Purpose : Characterizing inter-individual variability

yij
ϕi

∼
ind.

=

f (tij , ϕi ) + g(tij , ϕi )eij , eij ∼ N (0, 1)
i.i.d.

µ + βCi + ηi , ηi ∼ N (0, Ω)
i.i.d.

Questionswhich covariates for which components of ϕi ?
,→ support of β (i.e. non nul coefficients in β) ?
Sβ = {(k, k ′ ) ∈ {1, . . . , d} × {1, . . . q} : βk,k ′ ̸= 0}

on which parameters to include random effects, with which covariance structure (i.e.
non nul elements in Ω) ?
,→ Related work : tests on variance components in NLMEM, see (Baey et al., 2019)

Maud Delattre

11 / 33

State of the art
Many methods for linear mixed-effects models
LASSO (Bondell H. D. et al., 2010), (Schelldorfer J. et al., 2011)...
. . .but no contribution in nonlinear models
In practice : BIC (Bayesian information criterion) is widely used (Schwarz G.,
1978)−2LL + dim(θ) × log(sample size)→ Sample size definition is unclear : = N or ntot =

N
P

ni ?

i=1

no consensus, both definitions are implemented in softwares
Purpose : propose a good definition of the BIC by revisiting the theoreticalfoundations of the criterion

Maud Delattre

12 / 33

Theoretical fundations of the criterion (Schwarz G., 1978)
Model collection : M = {M1 , . . . , MK }Bayesian model choice :M̂=arg max p(Mk |Y)where p(M|Y)∝π(M)p(Y|M)

and p(Y|M)

=

M∈M

Z
p(Y|θ, M)p(θ|M)dθ.

Laplace approximation :log p(Y|M) ≈ log p(Y|θ̂, M) −∂ 2 log p(Y|θ, M)
1
log det (Hθ̂ ) , Hθ̂ = −
|θ=θ̂
2
∂θ∂θ′

In models based on n independent observations, Hθ̂ ≈ nIθ̂ when n → ∞ anddim(θ)
1
1
log det (Hθ̂ ) ≈log n + log det (Iθ̂ ) ,2
2
2Maud Delattre

13 / 33

Parameter partitioning :



ϕF ,i
ϕR,i




=

µF
µR




+

βF
βR




Ci +

0



ηR,i

, ηR,i ∼ N (0, ΩR )i.i.d.

Asymptotic framework : N → ∞, mini ni → ∞ [Nie, 2007]

√

ntot
0

√0
N



cF − βF⋆βc
βR − βR⋆



 
−→ N

0,

I1
0

0
I2



Hθ̂ ≈ block diagonallog det (Hθ̂ ) ≈ log det(ntot I1 ) + log det(N I2 )ResultSelect model (i.e. covariate structure) minimizing :BICh = −2LL + dim(βR ) log N + dim(βF ) log ntot
Now implemented in the saemix R package
Maud Delattre

14 / 33

BICh = −2LL + dim(βR ) log N + dim(βF ) log ntot
Remarks
Covariate selection depends on Ω.
Hybrid version of the two simplest (empirical) BIC definitions :If dim(βR ) = 0,
BICntot = −2LL + dim(β) log ntot .
If dim(βF ) = 0,
BICN = −2LL + dim(β) log N.
Numerical experiments
BICh performs better than BICN and BICntot
Results evoke theoretical consistency for BICh (theoretical study under progress)

Maud Delattre

15 / 33

Algorithm for joint selection
Simultaneous selection of covariates and random effects
ϕi = µ + βCi + ηi , ηi ∼ N (0, Ω)
i.i.d.

→ Exhaustive search impossible
Ex. q = 6 covariates, d = 3 individual parameters ⇒ ≈ 2 millions possible combinations→ Incorrect specification of β ⇔ Incorrect specification of Ω
Iterative algorithm, iteration t
i)Selection of the random effect structure (BICΩ , pen = · · · log N)ii)

Selection of the covariate structure (stepwise, BICβ ,pen = · · · log N + · · · log ntot )Numerical experiments :Performance equivalent to exhaustive search or competitive methods
Substantial reduction of the calculation time

Maud Delattre

16 / 33

Pharmacokinetics example
Biological context : PK of an antibiotic administered by infusion (Burdet C. et al.,2015)

Model :
dA1 (t, ϕ) Q
Q + CL=A2 (t, ϕ) −A1 (t, ϕ)dt
V2
V1
dose 1+10≤t−tD ≤Tinf
Tinf V1dA2 (ϕ) Q
Q=
A1 (t, ϕ) −A2 (t, ϕ)dt
V1
V2

Yij = A1 (tij , ϕi ) + (a + bA1 (tij , ϕi ))ϵij , ϵij ∼ N (0, 1)
i.i.d.

⊤

ϕi = (log CLi , log V1,i , log Qi , log V2,i ) = µ + βCi + ηi , ηi ∼ N (0, Ω)
i.i.d.

Observations (yij ) : blood concentration of the antibiotic
N = 53 patients, ntot = 247
q = 12 covariates
Maud Delattre

17 / 33

Pharmacokinetics example
Model for the individual parameters :ϕi = (log CLi , log V1,i , log Qi , log V2,i )⊤= µ + βCi + ηi , ηi ∼ N (0, Ω)
i.i.d.

Result :covariates}

random effects

ϕ

{

CL

ClCr

age

⋆•

⋆

sex

⋆•
⋆

Q
⋆

P/F

η

⋆•
⋆•

V1

V2

weight

⋆•

⋆
•

BICN
1492.624
1643.676•(⋆) our method⋆•(•) (Burdet C. et al., 2015)→ best compromise between covariates and random effects→ best choice of dose according to the sex and age of the patients ?

Maud Delattre

18 / 33

1

Introduction and model formulation

2

Covariate selection
On the use of BIC in NLMEM
High-dimensional variable selection in NLMEM

Maud Delattre

19 / 33

First part of M. Naveau’s PhD
Preprint available at https://hal.inrae.fr/hal-03685060
M. Naveau et al. (2022) Bayesian high-dimensional covariate selection in non-linear
mixed-effects models using the SAEM algorithm.

Maud Delattre

20 / 33

Model
1)

Description of the intra-individual variability
yij ∼ f (tij , ϕi ) + g(tij , ϕi )eij , eij ∼ N (0, 1) , j = 1 : ni , i = 1 : N
ind.

2)

i.i.d.

Description of the inter-individual variability
ϕi = µ + βCi + ηi , ηi ∼ N (0, Ω)
i.i.d.

Issue : covariate selection in large dimension framework, ie q >> N,→ BIC can not be used anymore

Consider d = dim(ϕi ) = 1 (but extension to d > 1 should be straightforward)

Maud Delattre

21 / 33

Frequentist framework :Linear mixed-effects models : theoretical consistency results on Lasso-type methods(Schelldorfer J. et al., 2011).
Nonlinear mixed-effects models : computational aspects only (Bertrand J. andBalding D. J., 2013), (Ollier, 2022).
Bayesian framework :Linear (fixed-effects) models : use of sparsity priors with both computational andtheoretical contributions (see Tadesse and Vannucci (2021)).
Nonlinear mixed-effects models : implementation of MCMC methods (Lee, 2022).
Our approach : Combining spike-and-slab priors for variable selection and the
MCMC-SAEM algorithm.

Maud Delattre

22 / 33

Spike-and-slab prior (George and McCulloch, 1997)
For each 1 ≤ ℓ ≤ q : inclusion variable


δℓ =

1
0

if covariate ℓ is to be included in the model,
otherwise.

and distributionπ(β|δ) = Np (0, diag((1 − δℓ )ν0 + δℓ ν1 )) , 0 ≤ ν0 < ν1 fixed.
i.e. the βℓ ’s are independent and :βℓ |(δℓ = 0) ∼ N (0, ν0 ) : "spike" distribution, ν0 smallβℓ |(δℓ = 1) ∼ N (0, ν1 ) : "slab" distribution, ν1 large

Figure – Prior spike-and-slab. Source : Deshpande et al. (2019)
Maud Delattre

23 / 33

Hierarchical bayesian model

Maud Delattre

24 / 33

Proposed method : explore different levels of sparsity for β by trying different values for ν01

For ν0 ∈ ∆,▶ Compute the maximum a posteriori estimator of Θ with an MCMC-SAEM
algorithm :b MAPΘ= argmax π(Θ|y )ν0Θ∈Λ▶ Select covariates (Ročková and George, 2014) :


bν0 =
S

2



ℓ ∈ {1, . . . , p}|(βbνMAP)ℓ | ≥ sβ (ν0 , ν1 , α
bνMAP
)
0
0

bν0 )ν0 ∈∆ with eBIC (Chen and Chen, 2008) :Select the "best" model among (S




bν0 )
ν̂0 = argmin eBIC (Sν0 ∈∆Maud Delattre

25 / 33

Proposed method : explore different levels of sparsity for β by trying different values for ν01

For ν0 ∈ ∆,▶ Compute the maximum a posteriori estimator of Θ with an MCMC-SAEM
algorithm :b MAPΘ= argmax π(Θ|y )ν0Θ∈Λ▶ Select covariates (Ročková and George, 2014) :ℓ ∈ {1, . . . , p}

bν0 =
S

|(βbνMAP)ℓ | ≥ sβ (ν0 , ν1 , α
bνMAP
)
0
0

Threshold ? Estimate the a posteriori inclusion probabilities of each covariate :MAPδ̂ = argmax P(δ|Θ̂MAPν0 ) s.t. δ̂ℓ = 1 ⇐⇒ P(δℓ = 1|Θ̂ν0 ) ≥ 0.5δ2bν0 )ν0 ∈∆ with eBIC (Chen and Chen, 2008) :Select the "best" model among (S




bν0 )
ν̂0 = argmin eBIC (Sν0 ∈∆Maud Delattre

25 / 33

Regularization plot

Maud Delattre

26 / 33

Numerical experiments
Model
iid

yij = f (φi , ψ, tij ) + εij , εij ∼ N (0, σ 2 ) wheref (φi , ψ, tij ) =ψ 1
1 + exp −tij − φiψ2

ψ = (ψ1 , ψ2 ) fixed effects.
iid

φi = µ + t βVi + ξi , ξi ∼ N (0, Γ2 )Results
N = 100

Maud Delattre

N = 200

27 / 33

Application : plant senescence genetic marker identification

Observations : proportion of senesced surface of the canopy (yij ) for N = 220 wheat
varieties at n = 18 time points

Goal : identify genetic markers associated to senescence
26 189 SNPs in total
chromosome by chromosome analysis ⇒ q ≡ 2000 >> N for each chromosome
Maud Delattre

28 / 33

Model100
 + εij , εij i.i.d.
yij =∼ N (0, σ 2 ),
tij − φi
1 + exp −ψi
φi = µ + λvi + βCi + ξi ,
ψi = η + ωi ,

i.i.d.

ξi ∼ N (0, Γ2 ),i.i.d.

ωi ∼ N (0, Ω2 ),where
Ci : q-vector of genetic markers for variety i, 1 ≤ i ≤ N
vi : sub-populations covariates, not subject to selection

Maud Delattre

29 / 33

Results
Regularization plot for Chromosome 6A

Maud Delattre

30 / 33

Results

Maud Delattre

31 / 33

References I

Baey, C., Cournède, P.-H., and Kuhn, E. (2019). Asymptotic distribution of likelihood ratio test statistics for variancecomponents in nonlinear mixed effects models. Computational Statistics and Data Analysis, 135 :107–122.
Bertrand J. and Balding D. J. (2013). Multiple single nucleotide polymorphism analysis using penalized regression in nonlinear
mixed-effect pharmacokinetic models. Pharmacogenetics and Genomics, 23 :167–174.
Bondell H. D., Krishna A., and Ghosh S. K. (2010). Joint Variable Selection for Fixed and Random Effects in Linear
Mixed-Effects Models. Biometrics, 66(4) :1069–1077.
Burdet C., Pajot O., Couffignal C., Armand-Lefevre L., Foucrier A., LaouÃ©nan C., Wolff M., Massias M., and MentrÃ© F.
(2015). Population pharmacokinetics of single-dose amikacin in critically ill patients with suspected ventilator-associatedpneumonia. European Journal of Clinical Pharmacology, 71 :75–83.
Chen, J. and Chen, Z. (2008). Extended bayesian information criteria for model selection with large model spaces. Biometrika,
95(3) :759–771.
Deshpande, S. K., Ročková, V., and George, E. I. (2019). Simultaneous variable and covariance selection with the multivariate
spike-and-slab lasso. Journal of Computational and Graphical Statistics, 28(4) :921–931.
George, E. I. and McCulloch, R. E. (1997). Approaches for bayesian variable selection. Statistica sinica, pages 339–373.
Lavielle M. (2014). Mixed Effects Models for the Population Approach. Models, Tasks, Methods and Tools. Chapman andHall/CRC Biostatistics Series. Chapman and Hall/CRC.
Lee, S. Y. (2022). Bayesian Nonlinear Models for Repeated Measurement Data : An Overview, Implementation, and Applications.
Publisher : Preprints.
Ollier, E. (2022). Fast selection of nonlinear mixed effect models using penalized likelihood. Computational Statistics & Data
Analysis, 167.
Pinheiro J. C. and Bates D. M. (2000). Mixed-Effects Models in S and S-PLUS. Statistics and Computing. Springer-Verlag New
York Inc., 1st ed. edition.

Maud Delattre

32 / 33

References II

Ročková, V. and George, E. I. (2014). Emvs : The em approach to bayesian variable selection. Journal of the American Statistical
Association, 109(506) :828–846.
Schelldorfer J., Buhlmann P., and Van De Geer S. (2011). Estimation for High-Dimensional Linear Mixed-Effects Models Using
L1-Penalization. Scandinavian Journal of Statistics, 38 :197–214.
Schwarz G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2) :461–464.
Tadesse, M. G. and Vannucci, M. (2021). Handbook of bayesian variable selection.

Maud Delattre

33 / 33