nlmixr2 — un modèle, plusieurs estimateurs

Les méthodes d'estimation de population présentées ici traitent la même intégrale de différentes façons. Pour chaque sujet, la vraisemblance des observations exige de moyenner sur tous les $\eta$ compatibles avec la population, et cette intégrale n'a pas de forme close dès que le modèle est non linéaire en $\eta$ — situation fréquente en PK. Le modèle ne dit rien de cette intégrale : il fournit seulement de quoi l'écrire. Ce que est = choisit, c'est la manière de l'approximer ou de l'échantillonner ; posthoc est le cas distinct où les paramètres de population restent fixés.

FOCEI utilise une approximation locale. Il remplace le modèle par sa tangente au mode individuel $\hat{\eta}_i$, ce qui rend l'intégrande gaussien et l'intégrale analytique. Le prix : il faut retrouver $\hat{\eta}_i$ pour chaque sujet à chaque itération de population, et il faut les dérivées $\partial f / \partial \eta$. Une boucle d'optimisation dans une boucle d'optimisation.

SAEM ne repose pas sur la même approximation locale déterministe. Il traite les paramètres individuels comme des données manquantes et les simule par MCMC au lieu de rechercher leur mode conditionnel à chaque itération. Sa boucle interne propose un $\eta$, évalue le modèle, puis accepte ou rejette la proposition. Il n'utilise donc ni les mêmes dérivées ni la même optimisation interne que FOCEI.

nlme utilise une autre approximation locale, par l'algorithme alterné de Lindstrom et Bates, qui enchaîne un pas de moindres carrés non linéaires pénalisés et un pas de modèle linéaire mixte. Il aboutit en pratique à une approximation très proche de FOCE sans interaction.

posthoc ne l'aborde pas du tout. Les paramètres de population sont fixés ; il ne reste qu'à trouver les $\hat{\eta}_i$. Il n'y a plus de vraisemblance de population à maximiser.

Le modèle, écrit une fois, ne définit que deux densités : celle des observations sachant les paramètres individuels, $p(y_i \mid \psi_i)$, et celle des paramètres individuels dans la population, $p(\psi_i \mid \theta)$. La vraisemblance de population en découle mécaniquement, et elle est la même pour tout le monde :

Relisez maintenant la structure d'un modèle nlmixr2 avec cette formule en tête. Le bloc ini({}) fixe le point de départ de $\theta$ et la structure de $p(\psi_i \mid \theta)$ ; le bloc model({}) définit $p(y_i \mid \psi_i)$. Ni l'un ni l'autre ne dit un mot de l'intégrale. Le signe $\int$ n'appartient à aucun des deux blocs — il appartient au moteur. C'est toute la justification de l'argument est =.

est =	Traitement de l'intégrale	Point d'appui	Besoin de $\partial f / \partial \eta$
"focei"	Taylor d'ordre 1 + interaction $\eta$–$\varepsilon$	$\hat{\eta}_i$	oui
"saem"	échantillonnage MCMC, aucune linéarisation	—	non
"nlme"	Lindstrom-Bates (PNLS/LME alternés), $\approx$ FOCE sans interaction	$\hat{\eta}_i$	oui
"posthoc"	pas d'intégrale : $\theta$ et $\Omega$ fixés, on ne cherche que les $\hat{\eta}_i$	—	oui

Le modèle, puis les quatre appels :

mod <- function() {
  ini({
    tka <- log(1.1);  tcl <- log(0.135);  tv <- log(7.8)
    eta.cl ~ 0.09
    eta.v  ~ 0.04
    prop.err <- 0.12
  })
  model({
    ka <- exp(tka)
    cl <- exp(tcl + eta.cl)
    v  <- exp(tv  + eta.v)
    d/dt(depot) = -ka * depot
    d/dt(centr) =  ka * depot - (cl / v) * centr
    cp = centr / v
    cp ~ prop(prop.err)
  })
}

f1 <- nlmixr2(mod, dat, est = "focei",
              control = foceiControl(maxOuterIterations = 5000, covMethod = "r,s"))
f2 <- nlmixr2(mod, dat, est = "saem",
              control = saemControl(nBurn = 500, nEm = 300, nmc = 5, seed = 20260716))
f3 <- nlmixr2(mod, dat, est = "nlme")
f4 <- nlmixr2(mod, dat, est = "posthoc")

mod n'a pas bougé. Chaque moteur a son objet de contrôle — foceiControl(), saemControl(), nlmeControl() — et n'accepte que le sien : les réglages ne se transposent pas, parce que les algorithmes n'ont pas les mêmes leviers.

Ce que rxode2 compile. Entre f1 et f2, ce qui change n'est pas le modèle : c'est son consommateur. rxode2 traduit le bloc model({}) en C, le compile en bibliothèque partagée (.dll sous Windows, .so ailleurs) et la charge dans la session R. Le coût est payé une fois ; les dizaines de milliers de résolutions qui suivent sont des appels à du code machine, pas à de l'interprété. C'est aussi pourquoi un compilateur est requis — Rtools sous Windows — et pourquoi le premier ajustement d'un modèle marque un temps d'arrêt qui n'a rien à voir avec l'estimation.

Le SAEM utilise des résolutions en avant. FOCEI réclame en plus $\partial f / \partial \eta$ ; rxode2 dérive alors symboliquement le système et compile les équations de sensibilité en avant. Le fichier source est identique, mais l'objet compilé peut différer. Un est = "focei" qui suit un est = "saem" peut donc nécessiter une compilation adaptée aux sensibilités requises.

Ce que deviennent les métriques. SAEM estime les paramètres sans fournir, à lui seul, une valeur unique directement comparable à toutes les sorties d'autres moteurs. Selon les options et la version, nlmixr2 peut ajouter plusieurs évaluations après l'estimation. Avec logLik = FALSE, une OBJF peut notamment être évaluée par l'approximation FOCEi aux estimations SAEM ; avec logLik = TRUE, une log-vraisemblance peut être évaluée par quadrature de Gauss-Hermite.

# l'OFV rendu apres un run SAEM : par defaut, l'objectif FOCEi evalue aux estimations SAEM
f2 <- nlmixr2(mod, dat, est = "saem", control = saemControl(logLik = FALSE))

# ou une vraie quadrature de Gauss-Hermite : nnodesGq = 1 donne Laplace
f2q <- nlmixr2(mod, dat, est = "saem",
               control = saemControl(logLik = TRUE, nnodesGq = 3, nsdGq = 1.6))

Les noms d'options sont dépendants de la version : les exemples utilisent la convention actuelle nnodesGq / nsdGq. Vérifiez la documentation de votre installation avant de recycler un ancien script.

OBJF et logLik ne sont donc pas des synonymes. AIC et BIC dérivent d'une log-vraisemblance donnée, avec une convention de comptage des paramètres et des observations ; ils ne décrivent pas automatiquement la même évaluation que l'OBJF affichée. La mention « SAEM OBJF by FOCEi approximation » décrit uniquement l'évaluateur de cette OBJF aux estimations SAEM : elle ne transforme pas toutes les vraisemblances présentes dans l'objet en évaluations FOCEi. Avant toute comparaison, consignez estimateur, évaluateur, métrique, convention de signe et de constante, options et version.

Sortie après est = "saem"	Ce qu'elle représente	Interprétation correcte
paramètres SAEM	estimations de $\theta$, $\Omega$ et de l'erreur obtenues par approximation stochastique EM	résultat de l'estimation SAEM
OBJF avec logLik = FALSE	approximation FOCEi évaluée aux paramètres SAEM dans la configuration décrite ici	diagnostic post-estimation, pas « vraisemblance SAEM » ni preuve de comparabilité LRT
logLik, AIC et BIC	quantités liées à l'évaluation de vraisemblance retenue par l'objet et aux conventions de calcul	identifier explicitement l'évaluateur ; ne pas les déduire du seul libellé OBJF
saemControl(logLik = TRUE)	évaluation numérique postérieure à l'estimation par quadrature de Gauss-Hermite (nnodesGq, nsdGq)	approximation distincte, dont le coût croît rapidement avec la dimension
est = "imp", "impmap" ou "qrpem"	estimateurs par importance sampling disponibles dans la série 7 et encore signalés comme méthodes de recherche	ne pas les confondre avec une simple OFV calculée après SAEM ; documenter méthode, version et réglages

Les exemples ont été revérifiés avec la documentation nlmixr2 7.0.1 / nlmixr2est 7.1. Les valeurs par défaut FOCEI et SAEM ont changé dans la série 7 ; reproduire un ancien run exige donc de consigner la version et les objets de contrôle, pas seulement la graine.

Scénario pédagogique, pas benchmark reproduit. Imaginez le même modèle oral à 1 compartiment ajusté successivement par FOCEI, SAEM, nlme et posthoc. Les estimations peuvent être proches sur des données informatives, mais les temps d'exécution dépendent du matériel et des réglages, et les nombres appelés OFV ne deviennent pas comparables parce qu'ils figurent dans le même tableau.

• FOCEI optimise une approximation conditionnelle de la vraisemblance.
• SAEM estime les paramètres par approximation stochastique ; une OBJF FOCEi éventuellement évaluée après coup n'est pas l'objectif que SAEM a optimisé.
• nlme utilise l'approximation de Lindstrom-Bates, sans interaction.
• posthoc fixe les paramètres de population et n'est donc pas un nouvel ajustement de population.

Une différence entre ces sorties peut refléter un changement d'approximation ou de convention, pas une amélioration du modèle. Seules des exécutions documentées, avec jeu de données, graine, versions et contrôles, autoriseraient à présenter des résultats numériques comme un benchmark.

Deuxième temps : là où les moteurs peuvent cesser d'être interchangeables. Imaginez maintenant un protocole épars et une élimination de Michaelis-Menten.

  model({
    ka <- exp(tka)
    vm <- exp(tvm + eta.vm)
    km <- exp(tkm)
    v  <- exp(tv  + eta.v)
    d/dt(depot) = -ka * depot
    d/dt(centr) =  ka * depot - (vm * (centr / v)) / (km + (centr / v))
    cp = centr / v
    cp ~ prop(prop.err)
  })

Dans ce type de situation, FOCEI peut atteindre maxInnerIterations, produire une covariance instable ou devenir sensible aux valeurs initiales, car il recherche un mode conditionnel pour chaque sujet. SAEM peut être plus robuste sur cette surface grâce à son exploration stochastique, sans que cela garantisse une convergence correcte ni une information suffisante dans les données.

Le piège de nlmixr2 n'est pas qu'un moteur soit mauvais. C'est que changer de moteur coûte un mot, et qu'aucune trace ne subsiste dans le nombre qui en sort.

Le scénario est banal, et c'est ce qui le rend redoutable. Un modèle de base ne démarre pas en FOCEI ; on l'ajuste avec SAEM, puis on utilise ses paramètres comme valeurs initiales d'un modèle avec covariable ajusté par FOCEI. Soustraire l'OBJF associée au premier ajustement de celle du second mélange alors deux chaînes d'estimation.

Il faut réajuster le modèle de base et le modèle avec covariable avec le même FOCEI, les mêmes conventions et des contrôles documentés. Seule cette différence constitue une candidate au LRT, sous ses autres hypothèses. L'écart qui disparaît après ce réalignement ne venait pas de la covariable, mais du changement d'estimateur ou d'évaluation.

Sous NONMEM, changer d'estimateur oblige à éditer $ESTIMATION : un acte délibéré, sur une ligne visible, archivée avec le run. Sous nlmixr2, c'est un argument d'appel, et l'objet mod — celui que vous relirez dans six mois — n'en garde aucune trace. Six mois plus tard, fit1 et fit2 sont deux objets R avec deux OFV, et votre tableau de ΔOFV ne se souvient plus que l'un était SAEM.

La discipline est simple, et elle est entièrement à votre charge : un seul moteur par chaîne de comparaison, décidé avant de commencer, et le est = écrit à côté de chaque OFV dans le journal de run — au même titre que le numéro du modèle. Le nombre, seul, ne le porte pas.

• Le trait de conception de nlmixr2 : le modèle est un objet R, l'estimateur est un argument. ini({}) et model({}) décrivent les densités ; le signe $\int$ appartient au moteur, d'où est =.
• Quatre workflows représentatifs : "focei" approche l'intégrale par tangente en $\hat{\eta}_i$ avec interaction, "saem" l'échantillonne, "nlme" l'approche autrement (Lindstrom-Bates, $\approx$ FOCE sans interaction), et "posthoc" fixe $\theta$ pour ne chercher que les $\hat{\eta}_i$. D'autres estimateurs existent dans nlmixr2.
• rxode2 compile le modèle en C, une fois, et toutes les résolutions suivantes sont du code machine. FOCEI a besoin des sensibilités $\partial f / \partial \eta$, que rxode2 dérive et compile en plus : même source, objet compilé différent — d'où la recompilation en passant de SAEM à FOCEI.
• SAEM peut être plus robuste sur certaines surfaces conditionnelles difficiles, avec des données éparses ou des valeurs initiales sensibles ; cela ne garantit ni convergence ni identifiabilité. La raideur des EDO relève séparément du choix et du réglage du solveur rxode2, et ne justifie pas à elle seule SAEM. FOCEI est un choix cohérent lorsqu'une chaîne de ΔOFV fondée sur FOCEI est planifiée.
• OBJF, logLik, AIC et BIC ne sont pas interchangeables. Après SAEM, nlmixr2 peut ajouter plusieurs évaluations postérieures ; rapportez toujours l'évaluateur, la métrique, la convention, les options et la version.
• Même évaluateur ne veut pas dire LRT valide. Une OBJF FOCEi calculée aux estimations SAEM et une OBJF obtenue à l'optimum FOCEI correspondent à des paramètres issus de deux procédures différentes. Pour comparer des modèles emboîtés, réajustez-les avec la même chaîne d'estimation et d'évaluation, puis vérifiez aussi les hypothèses du LRT.
• Un seul moteur par chaîne de comparaison, choisi d'avance, et le est = consigné à côté de chaque OFV. La facilité du changement est exactement ce qui rend l'erreur facile.