Monolix — le moteur SAEM

La vraisemblance de population exige d'intégrer les paramètres individuels, qu'on n'observe pas. Cette intégrale n'a pas de forme close dès que le modèle est non linéaire en ces paramètres — situation fréquente en PK/PD. Deux écoles s'affrontent depuis quarante ans.

Approximer localement le modèle ou l'intégrale. FO et FOCE utilisent une approximation de Taylor du modèle autour de $\eta = 0$ ou autour du mode individuel ; Laplace approxime l'intégrale autour de ce mode. On maximise ensuite exactement une fonction approchée.

Ne pas calculer l'intégrale du tout. C'est l'EM, et sa version stochastique le SAEM : les paramètres individuels sont traités comme des données manquantes, qu'on simule au lieu de les intégrer. Le modèle n'est jamais déformé — il est seulement évalué, en avant, pour des valeurs de paramètres tirées au sort.

La différence se voit mécaniquement, dans la boucle interne. À chaque itération, FOCE doit résoudre pour chaque sujet un problème d'optimisation — trouver le mode $\hat{\eta}_i$ — et a besoin des dérivées du modèle par rapport à $\eta$. Deux choses peuvent casser : l'optimisation interne peut ne pas converger, et les dérivées peuvent ne rien vouloir dire (une différence finie à travers une EDO raide ou le flanc quasi vertical d'un Emax ne mesure que du bruit d'intégrateur). La boucle interne du SAEM, elle, est un pas de Metropolis-Hastings : on propose un $\eta$, on évalue le modèle une fois, on accepte ou on rejette sur un rapport de vraisemblance. Aucune dérivée, aucun optimum interne.

Le prix est double et il faut le payer les yeux ouverts. La réponse est stochastique : deux runs identiques ne rendent pas exactement le même chiffre. Et comme l'algorithme ne calcule jamais la vraisemblance, il ne vous en rend aucune à l'arrivée.

La vraisemblance à maximiser intègre les paramètres individuels $\psi_i$ :

L'EM contourne l'intégrale en changeant de cible. Si on connaissait les $\psi_i$, maximiser serait trivial : ce ne seraient que des moyennes et des variances. On alterne donc entre estimer ce qu'on ne voit pas et maximiser comme si on le voyait :

Sauf que cette espérance est, elle aussi, une intégrale sans forme close. Le SAEM la remplace par une approximation stochastique, en trois sous-étapes par itération $k$ :

S — Simulation. Tirer $\psi^{(k)}$ dans la distribution conditionnelle $p(\psi \mid y, \theta_k)$, par MCMC. C'est le seul endroit où le modèle structural est appelé.

A — Approximation. Mettre à jour une moyenne mobile des statistiques exhaustives $S$ du modèle à données complètes :

M — Maximisation. Ré-estimer $\theta_{k+1}$ à partir de $s_{k+1}$, en forme close pour la famille exponentielle. C'est ce qui rend le M-step gratuit : il n'y a rien à optimiser numériquement.

Tout se joue dans $\gamma_k$. Delyon, Lavielle et Moulines (1999) donnent les deux conditions de convergence presque sûre : $\sum_k \gamma_k = \infty$ — le pas ne doit pas s'éteindre trop vite, sinon la suite se fige n'importe où — et $\sum_k \gamma_k^2 < \infty$ — il doit s'éteindre quand même, sinon le bruit ne s'efface jamais. Kuhn et Lavielle (2004) ont porté ce résultat aux modèles non linéaires à effets mixtes, en remplaçant le tirage exact du S-step par un noyau MCMC.

Monolix choisit $\gamma_k$ en deux phases, et c'est exactement la lecture du graphe de convergence :

Phase 1 — exploration. $\gamma_k = 1$. La récurrence se réduit à $s_{k+1} = S(y, \psi^{(k)})$ : aucune mémoire, chaque itération oublie la précédente. La suite $\theta_k$ ne converge pas — elle erre, à grands pas, dans la région de forte vraisemblance. C'est voulu : c'est ce qui rend le SAEM peu sensible aux valeurs initiales.

Phase 2 — lissage. $\gamma_k$ décroît, typiquement en $1/k^a$ avec $a$ entre $0{,}5$ et $1$. La mémoire s'allonge, chaque nouvelle simulation ne corrige plus $s_k$ que d'une fraction, le bruit de Monte-Carlo se moyenne et $\theta_k$ converge.

Dans le fichier projet, tout cela se règle explicitement :

<MONOLIX>

[TASKS]
populationParameters()
fim(method = StochasticApproximation)
logLikelihood(method = ImportanceSampling)   ; tache SEPAREE, a demander

[SETTINGS]
POPULATION:
exploratoryiterations = 500     ; phase 1 : pas constant, gamma = 1
smoothingiterations   = 200     ; phase 2 : pas decroissant, gamma ~ 1/k
exploratoryautostop   = yes     ; coupe la phase 1 sur critere
smoothingautostop     = yes
nbchains              = 5       ; chaines MCMC : peu de sujets => en mettre plus
simulatedannealing    = yes     ; freine la decroissance des variances en phase 1

LL:
nbfixediterations     = 10000   ; taille Monte-Carlo de l'echantillonnage d'importance

Deux réglages méritent un mot. nbchains : quand les sujets sont peu nombreux, une seule chaîne par sujet produit un bruit de simulation trop gros devant l'information des données ; Monolix duplique alors les sujets en plusieurs chaînes pour moyenner ce bruit. simulatedannealing : pendant l'exploration, les variances ($\omega^2$, les paramètres du modèle d'erreur) n'ont pas le droit de décroître plus vite qu'un coefficient imposé. Sans ce frein, le SAEM referme les variances autour des valeurs initiales dès les premières itérations — et il n'explore plus rien.

La vraisemblance, calculée à part. Relisez le M-step : il ne touche qu'à $s_{k+1}$, des statistiques du modèle à données complètes. À aucun moment l'algorithme n'évalue $L(\theta)$. Le SAEM maximise la vraisemblance sans jamais la calculer. À la fin du run vous avez $\hat{\theta}$ et rien à mettre en face — d'où la tâche logLikelihood, séparée.

Deux méthodes s'offrent alors, et le choix n'est pas neutre. Linearization linéarise le modèle autour des modes individuels : rapide, mais elle réintroduit l'approximation locale que le SAEM avait évitée pendant l'estimation. ImportanceSampling est l'option la plus cohérente avec ce moteur : elle réécrit l'intégrale du sujet $i$ comme une espérance sous une loi de proposition $h$ qu'on sait simuler,

et l'estime par la moyenne de $M$ tirages :

La proposition $h$ est centrée sur la distribution conditionnelle du sujet — celle-là même que le SAEM vient d'échantillonner, on l'a donc gratuitement — mais avec des queues lourdes (une loi de Student plutôt qu'une gaussienne), pour qu'aucun tirage égaré dans une zone où $h$ est minuscule ne reçoive un poids explosif.

Un 1 compartiment oral, 60 sujets, 10 prélèvements chacun, soit $n_{obs} = 600$ observations. Le modèle est assez sage pour que FOCE-I y soit quasi exact : c'est le cas le plus favorable qui soit à une comparaison.

• Monolix, SAEM puis échantillonnage d'importance : $-2LL = 2149{,}2$, erreur standard de Monte-Carlo $0{,}31$.
• NONMEM, FOCE-I, mêmes données, même modèle : $OFV = 1045{,}9$.

Un écart de 1 103 points. Dans cet exemple construit, l'essentiel de l'écart vient d'une constante : NONMEM omet de son objectif le terme $n_{obs}\log(2\pi)$, que Monolix inclut dans sa log-vraisemblance :

Il reste $2149{,}2 - 1045{,}9 - 1102{,}7 = 0{,}6$ point, soit environ deux erreurs de Monte-Carlo plus le résidu de l'approximation FOCE. Ici, les chiffres deviennent proches après harmonisation. Cela ne constitue pas une règle générale : sur un modèle moins sage, une autre définition de vraisemblance, une autre approximation ou d'autres conventions peuvent laisser un écart non négligeable.

Second temps, plus utile au quotidien. On ajoute le poids sur $Cl$ et on relance :

• sans covariable : $-2LL = 2149{,}2 \pm 0{,}31$
• avec le poids : $-2LL = 2145{,}1 \pm 0{,}31$

$\Delta = 4{,}1$, contre un seuil de $3{,}84$ pour un $\chi^2$ à 1 degré de liberté à 5 %. On garde la covariable ? Prudence : les deux $-2LL$ sont bruités, et le bruit de leur différence vaut $\sqrt{0{,}31^2 + 0{,}31^2} \approx 0{,}44$. Le seuil se trouve à $(4{,}1 - 3{,}84)/0{,}44 \approx 0{,}6$ écart-type de la valeur estimée. Cette décision n'est pas portée par les données : elle est portée par la graine du générateur aléatoire.

Le SAEM échoue autrement. Il fait les itérations que vous lui avez demandées, il s'arrête, il affiche parfois des paramètres, mais cela ne garantit pas que la convergence soit satisfaisante : la documentation Monolix prévoit explicitement des critères de convergence non atteints et l'arrêt au nombre maximal d'itérations. Le juge pratique reste le graphe de convergence — et c'est précisément là que le piège se referme.

Le diagnostic ne se réduit pas à une phase. On examine les trajectoires et l'indicateur de convergence pendant l'exploration, puis leur stabilisation pendant le lissage, les statuts d'arrêt et la cohérence des estimations. Une trajectoire encore monotone à la bascule demande de la prudence ; un $\omega$ qui glisse vers zéro sans se stabiliser peut signaler un effet aléatoire peu soutenu. La robustesse se vérifie aussi en répétant l'estimation avec plusieurs graines et valeurs initiales.

L'auto-stop utilise la stabilité de l'indicateur de convergence, des omégas et des paramètres d'erreur. Le désactiver et prolonger les phases peut constituer une analyse de sensibilité sur un modèle difficile ; ce n'est pas une règle universelle ni une garantie de meilleur optimum.

Et le contrôle qui tranche : relancez avec une autre graine et d'autres valeurs initiales. Si $\hat{\theta}$ bouge au-delà du bruit annoncé, ce n'est pas le SAEM qui est en cause. C'est votre vraisemblance qui est plate ou multimodale — et aucun réglage d'algorithme ne réparera un modèle non identifiable.

• Monolix ne propose pas le SAEM parmi d'autres méthodes : il est construit autour. Ses trois bizarreries apparentes en découlent.
• FO/FOCE utilisent une approximation locale du modèle, Laplace une approximation locale de l'intégrale ; le SAEM évite cette linéarisation pendant l'estimation et échantillonne les effets individuels. Pas de dérivée, pas d'optimisation interne par sujet : d'où la robustesse aux modèles raides et aux données éparses.
• Deux phases : exploration à pas constant ($\gamma_k = 1$, aucune mémoire, on erre exprès), puis lissage à pas décroissant ($\gamma_k \approx 1/k$, le bruit se moyenne, on converge). Fondement : Delyon–Lavielle–Moulines (1999), porté au non linéaire par Kuhn–Lavielle (2004).
• Le SAEM maximise la vraisemblance sans jamais l'évaluer : le $-2LL$ est une tâche séparée, par échantillonnage d'importance ou par linéarisation.
• Le $-2LL$ par échantillonnage d'importance est bruité : lisez son erreur de Monte-Carlo, et augmentez nbfixediterations avant de trancher un LRT serré.
• L'évaluation de convergence combine trajectoires d'exploration, phase de lissage, indicateurs et statuts, puis stabilité entre graines et valeurs initiales ; aucune phase ne suffit seule.
• Un $-2LL$ de Monolix et une OFV de NONMEM ne se comparent qu'après vérification de la définition de l'objectif, des données, de l'approximation, des constantes et des conventions ; l'exemple montre comment une constante peut expliquer une grande partie d'un écart.