pith. machine review for the scientific record. sign in

arxiv: 2605.02534 · v1 · submitted 2026-05-04 · 📊 stat.ME · stat.OT

Recognition: 2 theorem links

Conditional bootstrap for non-linear mixed effects models

Sofia Kaisaridi , Moreno Ursino , Emmanuelle Comets
Authors on Pith no claims yet

Pith reviewed 2026-05-08 18:53 UTC · model grok-4.3

classification 📊 stat.ME stat.OT
keywords non-linear mixed effects modelsbootstrapSAEM algorithmuncertainty estimationcoverage probabilitynon-parametric resamplingdata structure preservation
0
0 comments X

The pith

A new conditional non-parametric bootstrap preserves data structure while estimating uncertainty in non-linear mixed effects models.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

Non-linear mixed effects models need reliable uncertainty estimates for both fixed effects and variance components, yet the asymptotic Fisher information matrix often under-covers the variances. The paper introduces cNP, a non-parametric bootstrap that draws individual random effects from the conditional distributions produced as a by-product of the SAEM algorithm and draws residuals from their fitted distribution. This resampling step keeps the original number of observations per subject and the original covariate layout intact, eliminating any need for stratification. In simulations with sigmoid Emax models across rich, sparse, and unbalanced designs, cNP and the simpler case bootstrap delivered coverage rates closer to the nominal 95 percent than either the asymptotic method or the classical non-parametric residual bootstrap, although performance degraded when residual error was large.

Core claim

The paper claims that resampling interindividual random effects from the conditional distribution of the individual parameters (a direct output of the SAEM algorithm) together with residuals from their estimated distribution produces a bootstrap distribution whose coverage matches or exceeds that of existing methods while automatically retaining the exact structure of the original dataset.

What carries the argument

The cNP procedure, which conditions the bootstrap draws of random effects on the SAEM-estimated conditional distributions of the individual parameters.

If this is right

  • The asymptotic method produces lower-than-nominal coverage for variance parameters across the tested designs.
  • cNP maintains coverage without requiring stratification when the number of observations per subject or the covariate repartition must stay fixed.
  • Both cNP and the case bootstrap outperform the classical non-parametric residual bootstrap in the reported simulations.
  • All bootstrap methods lose coverage once residual variability becomes large.
  • The procedure is directly available in the saemix R package for immediate use with existing SAEM fits.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • Analysts working with sparse pharmacokinetic datasets could adopt cNP to avoid distorting subject-level observation counts when constructing confidence intervals.
  • The same conditional-resampling idea might be tested in other hierarchical models whose estimation routines already produce subject-specific posterior distributions.
  • Performance under non-normal random effects or more intricate covariate structures remains an open question that could be checked with targeted simulations.
  • If the underlying SAEM conditional distributions themselves are misspecified, the bootstrap intervals may inherit that misspecification rather than correct it.

Load-bearing premise

The conditional distributions of the individual parameters obtained from SAEM accurately reflect the true interindividual variability without systematic bias.

What would settle it

A simulation in which the true interindividual variances are known and the empirical coverage of the cNP bootstrap intervals falls substantially below the nominal level for an unbalanced design.

Figures

Figures reproduced from arXiv: 2605.02534 by Emmanuelle Comets, Moreno Ursino, Sofia Kaisaridi.

Figure 1
Figure 1. Figure 1: : Measurement times in the four unbalanced designs. The first 3 designs involved 100 subjects divided in 5 view at source ↗
Figure 2
Figure 2. Figure 2: Coverage rates of the 90% CI along with the errorbars of their MC uncertainty for the original designs. view at source ↗
Figure 3
Figure 3. Figure 3: Coverage rates of the 90% CI and errorbars of their MC uncertainty for the unbalanced designs for the view at source ↗
Figure 4
Figure 4. Figure 4: Coverage rates of the 90% CI and errorbars of their MC uncertainty for the original design and the two view at source ↗
Figure 5
Figure 5. Figure 5: Coverage rates of 95% CI and errorbars of their MC uncertainty for the original designs. Dashed lines indicate view at source ↗
Figure 6
Figure 6. Figure 6: Relative bias on parameters (left) and SE (right) for the original designs. Dashed lines delineate absolute view at source ↗
Figure 7
Figure 7. Figure 7: Barplots of predicted RSE by PFIM 4.0, across the unbalanced designs compared to the original rich designs view at source ↗
Figure 8
Figure 8. Figure 8: Coverage rates of 90% CI and errorbars of their MC uncertainty for the unbalanced designs for the view at source ↗
Figure 9
Figure 9. Figure 9: Coverage rates of 95% CI and errorbars of their MC uncertainty for the unbalanced designs for the view at source ↗
Figure 10
Figure 10. Figure 10: Coverage rates of 95% CI and errorbars of their MC uncertainty for the unbalanced designs for the view at source ↗
Figure 11
Figure 11. Figure 11: Relative bias on parameters (left) and SE (right) for the original dosing scheme and the four unbalanced view at source ↗
Figure 12
Figure 12. Figure 12: Relative bias on parameters (left) and SE (right) for the original dosing scheme and the four unbalanced view at source ↗
Figure 13
Figure 13. Figure 13: Barplots of emprical SE for the scenarios with increased coefficient compared across the four designs tested view at source ↗
Figure 14
Figure 14. Figure 14: Coverage rates of 90% CI and errorbars of their MC uncertainty for the original design and the two designs view at source ↗
Figure 15
Figure 15. Figure 15: Coverage rates of 95% CI and errorbars of their MC uncertainty for the original design and the two designs view at source ↗
Figure 16
Figure 16. Figure 16: Coverage rates of 95% CI and errorbars of their MC uncertainty for the unbalanced designs. Dashed lines view at source ↗
Figure 17
Figure 17. Figure 17: Relative bias on parameters (up) and SE (down) for the original dosing scheme and the two scenarios with view at source ↗
Figure 18
Figure 18. Figure 18: Relative bias on parameters (up) and SE (down) for the original dosing scheme and the two scenarios with view at source ↗
read the original abstract

Background and Objective: Uncertainty in non-linear mixed effect models is often assessed using the Fisher information matrix to derive the standard errors of estimation. The bootstrap is an alternative to the asymptotic method, with different approaches to handle the different levels of individual and population variabilities. The simplest method is the Case bootstrap where the entire vector of individuals is resampled, but this approach does not take into account the hierarchical nature of non-linear mixed effect models (NLMEM). Methods: We propose here a non-parametric bootstrap, cNP, to preserve the structure of the original data. We resample interindividual random effects from the conditional distribution of the individual parameters, obtained as a by-product of the SAEM algorithm, and residuals from their distribution. cNP was implemented in the saemix package for R along with the case, parametric (Par), and non-parametric (NP) residual bootstraps. Coverage rates were compared in a simulation study using sigmoid Emax models, with rich, sparse and unbalanced designs, and 3 levels of residual variability. Results: The asymptotic method tended to produce lower than theoretical coverages for the variance terms. Bootstraps provided more adequate coverage, but none of the approaches maintained coverage when the residual error increased. Overall, the new cNP and the Case provided better coverage than the classical NP. Conclusion: The new conditional non-parametric bootstrap can be used when it is important to preserve the structure of the original dataset, such as the number of observations or the repartition of covariates as it does not require stratification.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The manuscript proposes a conditional non-parametric bootstrap (cNP) for non-linear mixed effects models that resamples individual random effects from the SAEM-derived conditional distributions p(η_i | y_i, θ̂) and residuals from their estimated distribution while holding the original design matrix fixed. Implemented in the saemix R package, cNP is compared to case, parametric, and classical non-parametric residual bootstraps. A simulation study using sigmoid Emax models under rich, sparse, and unbalanced designs at three residual variability levels shows that asymptotic methods undercover variance parameters, bootstraps improve coverage overall, and cNP together with case bootstrap outperform classical NP, although all methods lose coverage at high residual error. The authors conclude that cNP is suitable when preserving data structure (e.g., number of observations or covariate repartition) is required.

Significance. If the central claim holds, the work supplies a practical, structure-preserving alternative to existing bootstrap methods for NLMEM uncertainty quantification, directly addressing a common need in pharmacometrics and biostatistics. The simulation design—multiple designs, three residual levels, and explicit comparison to asymptotic and other bootstrap approaches—provides concrete, reproducible evidence of improved coverage for variance components relative to Fisher-information-based standard errors. Credit is due for the open implementation in saemix and for evaluating the method on both rich and sparse designs.

major comments (2)
  1. [Methods] The cNP procedure (Methods) resamples η_i from p(η_i | y_i, θ̂) evaluated at the fixed point estimate θ̂. Because the conditional distributions are not marginalized over the uncertainty in the population parameters, the resulting bootstrap distribution may be under-dispersed relative to the true sampling distribution of the NLMEM estimator. This approximation is load-bearing for the validity claim yet is not accompanied by a theoretical argument or diagnostic showing that the marginal variability is recovered; the observed coverage collapse at high residual variability could partly originate from this source.
  2. [Results] The simulation study (Results) reports only qualitative statements that cNP and case bootstrap achieve “better coverage” and that coverage “collapses” at high residual error. No numerical coverage rates, number of Monte Carlo replicates, exact model parameter values for the sigmoid Emax, or per-parameter coverage tables are supplied. Without these quantities it is impossible to judge the magnitude of improvement or to verify that the reported superiority is not driven by a small number of replicates or particular design realizations.
minor comments (2)
  1. [Abstract] The abstract states that coverage rates were compared but supplies neither the number of replicates nor any quantitative coverage figures; adding these would make the results summary self-contained.
  2. [Methods] Notation for the conditional distribution p(η_i | y_i, θ̂) is introduced without an explicit statement of how the SAEM by-product is extracted or whether any smoothing or kernel density estimation is applied before resampling.

Simulated Author's Rebuttal

2 responses · 1 unresolved

We thank the referee for their insightful comments and the opportunity to improve our manuscript on the conditional non-parametric bootstrap for non-linear mixed effects models. We address each of the major comments in detail below, indicating the revisions we plan to make.

read point-by-point responses

  1. Referee: [Methods] The cNP procedure (Methods) resamples η_i from p(η_i | y_i, θ̂) evaluated at the fixed point estimate θ̂. Because the conditional distributions are not marginalized over the uncertainty in the population parameters, the resulting bootstrap distribution may be under-dispersed relative to the true sampling distribution of the NLMEM estimator. This approximation is load-bearing for the validity claim yet is not accompanied by a theoretical argument or diagnostic showing that the marginal variability is recovered; the observed coverage collapse at high residual variability could partly originate from this source.

    Authors: We recognize the validity of this observation regarding the conditional nature of our cNP procedure. By resampling from the conditional distributions p(η_i | y_i, θ̂) at the fixed θ̂, the method approximates the bootstrap distribution without fully accounting for uncertainty in the population parameters. This choice prioritizes computational efficiency and structure preservation, which is a key advantage for applications where maintaining the original data design is important. Although we do not provide a new theoretical proof of marginal coverage recovery in this work, the empirical results from our simulation study support the practical utility of cNP, showing superior coverage for variance parameters compared to asymptotic and classical NP methods in most scenarios. We agree that the coverage degradation at high residual variability may relate to this or other factors and will revise the manuscript to include an expanded discussion on the limitations of the conditional approximation, potential reasons for coverage loss, and recommendations for when cNP is most appropriate. Additional simulation diagnostics will be considered if they can be implemented without excessive computation. revision: partial

  2. Referee: [Results] The simulation study (Results) reports only qualitative statements that cNP and case bootstrap achieve “better coverage” and that coverage “collapses” at high residual error. No numerical coverage rates, number of Monte Carlo replicates, exact model parameter values for the sigmoid Emax, or per-parameter coverage tables are supplied. Without these quantities it is impossible to judge the magnitude of improvement or to verify that the reported superiority is not driven by a small number of replicates or particular design realizations.

    Authors: We appreciate this feedback on the presentation of our simulation results. The current manuscript emphasizes qualitative comparisons to highlight the overall trends. In the revised version, we will add detailed numerical results, including tables with coverage percentages for each model parameter across all designs and residual error levels. We will also report the number of Monte Carlo replicates performed, the precise values used for the sigmoid Emax model parameters, and per-parameter coverage breakdowns to facilitate a quantitative assessment of the improvements and to confirm the robustness of our findings. revision: yes

standing simulated objections not resolved
  • A complete theoretical argument demonstrating that the conditional bootstrap fully recovers the marginal sampling distribution of the NLMEM estimator.

Circularity Check

0 steps flagged

No circularity: cNP bootstrap defined independently and evaluated externally via simulation

full rationale

The paper first specifies the cNP procedure as resampling individual random effects from the conditional distribution obtained as a by-product of SAEM, plus residual resampling, while preserving the original design matrix. This algorithmic definition stands alone. The simulation study on sigmoid Emax models (rich/sparse/unbalanced designs, three residual-variability levels) then measures coverage rates as an external performance check against known true values; coverage outcomes are not fed back to define, tune, or justify the resampling rule. No self-citation is invoked to establish uniqueness or to smuggle an ansatz; comparisons to case, parametric, and NP bootstraps rely on standard implementations. The derivation chain therefore contains no self-definitional, fitted-input, or load-bearing self-citation reductions.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

The proposal rests on standard NLMEM assumptions and the accuracy of SAEM conditional distributions; no new free parameters or invented entities are introduced beyond the bootstrap resampling rules.

axioms (2)
  • domain assumption The SAEM algorithm produces reliable estimates of the conditional distribution of individual parameters given the data.
    Invoked when stating that random effects are resampled from the conditional distribution obtained as a by-product of SAEM.
  • domain assumption Residuals can be treated as exchangeable draws from their marginal distribution for non-parametric resampling.
    Underlying the non-parametric residual bootstrap component.

pith-pipeline@v0.9.0 · 5581 in / 1606 out tokens · 34233 ms · 2026-05-08T18:53:30.471839+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

20 extracted references

  1. [1]

    Generate a bootstrap sample by resampling from the data and/or from the estimated model

  2. [2]

    Compute the estimates of the parameters of the model for the bootstrap sample

  3. [3]

    Case bootstrap (Case):This method consists of resampling with replacement the entire subjects ( ξi,yi where yi = (yi1, yi2, ...., yini)′) from the original data before modelling

    Repeat steps 1-2Btimes to obtain the bootstrap distribution of the parameter estimates. Case bootstrap (Case):This method consists of resampling with replacement the entire subjects ( ξi,yi where yi = (yi1, yi2, ...., yini)′) from the original data before modelling. It is also called thepaired bootstrap. It is the most straightforward way to do bootstrapp...

  4. [5]

    draw a sample{η ∗ i }i=1,...,N from a normal distribution with mean zero and covariance matrix ˆΩ

  5. [6]

    draw a sample {ϵ∗ ij}j=1,...ni from a normal distribution with mean zero and covariance matrix Ini where I denotes the identity matrix

  6. [7]

    The bootstrap sample is obtained as follows:

    generate the bootstrap responses y∗ ij =f(x ij,ˆµ, η∗ i ) +g(x ij,ˆµ, η∗ i , σ)ϵ ∗ ij (where we explicitly replace ψi by ˆµandη∗ i to highlight the dependency on the resampled random effects) Non-parametric residual bootstrap (NP):This method resamples with replacement from the residuals obtained after model fitting for all subjects. The bootstrap sample ...

  7. [8]

    fit the model to the data to obtain ˆθ

  8. [9]

    estimate the individual parametersˆηi as the EBE, then center and normalise them as detailed below 15 arXivpreprintCONDITIONAL BOOTSTRAP FORNLMEM

  9. [10]

    draw a sampleη ∗ i with replacement fori= 1, ...Nfrom the resulting set

  10. [11]

    calculate the residualsˆϵij =y ij −f(x ij,ˆµ,ˆηi), then center and normalise them

  11. [12]

    draw a sample{ϵ ∗ ij}j=1,...ni with replacement globally from the resulting set

  12. [13]

    To correct for this, Carpenter et al

    generate the bootstrap responsesy ∗ ij =f(x ij,ˆµ, η∗ i ) +g(x ij,ˆµ, η∗ i , σ)ϵ ∗ ij The normalisation in step (b) and (d) is necessary because the EBE generally suffer from regression to the mean and their variance may be considerably smaller than the true variability in the population (Karlsson and Savic (2007)). To correct for this, Carpenter et al. (...

    2007

  13. [14]

    center the raw estimated random effects:˜ηi = ˆηi −¯ηi

  14. [15]

    obtain the EVD of the estimated variance-covariance matrix: ˆΩ =V Ω DΩV T Ω where DΩ is the diagonal matrix containing the eigenvalues of ˆΩ

  15. [16]

    obtain the EVD of S, the variance-covariance matrix of the centered random effects: S=V S DSV T S where DS is the diagonal matrix containing the eigenvalues ofS

  16. [17]

    calculate the correction matrixA η using these two decompositions as Aη =V S D−1/2 S VΩ D−1/2 Ω

  17. [18]

    transform the centered random effects using the ratioA η:ˆη ′ i = ˜ηi ×A η Similarly, the residuals are transformed using the empirical variance:

  18. [19]

    center the raw estimated residuals:˜ϵij = ˆϵij −¯ϵij

  19. [20]

    calculate the correction factor Aσ = 1/σemp where σemp is the empirical standard deviation of the centered residuals

  20. [21]

    (2017)) and is avail- able on CRAN

    transform the centered residuals using the ratioA σ:ˆϵ ′ ij = ˜ϵij ×A σ Implementation and runtimes The bootstrap methods were implemented for the saemix package in version 3.2 (Comets et al. (2017)) and is avail- able on CRAN. Development versions are downloadable on github https://github.com/saemixdevelopment/ saemixextension. In terms of runtimes, the ...

    2017