Monolix — le modèle d'erreur résiduelle et la censure

Le modèle d'erreur répond à une seule question : à quelle distance de sa prédiction une observation devient-elle surprenante ?

Un dosage réel a deux régimes de bruit, et ils n'ont rien à voir l'un avec l'autre :

• un plancher, en mg/L, indifférent à la concentration — bruit de fond, ligne de base, tout ce qui subsiste quand il n'y a presque plus rien à mesurer. C'est constant(a) ;
• un pourcentage, qui grandit avec la concentration — dilutions, pipetage, calibration. C'est proportional(b).

Un jeu de données de PK couvre couramment deux ou trois ordres de grandeur, du pic au dernier creux. Il traverse donc les deux régimes. Aucun modèle à un seul terme n'est juste sur toute la gamme : c'est très exactement pourquoi les modèles combinés existent, et pourquoi ils sont un candidat fréquent.

Monolix écrit toutes ses observations sous une seule forme :

Une ligne de DEFINITION: ne fait que remplir les deux trous de cette formule :

• distribution= choisit $h$, la transformation de l'observation ;
• errorModel= choisit $g$, la fonction d'écart-type.

Ce sont deux réglages indépendants, et la moitié des malentendus vient de les confondre.

Les quatre errorModel.

errorModel	$g(f)$	ce que ça dit
constant(a)	$a$	plancher seul, en mg/L
proportional(b)	$b \cdot f$	pourcentage seul
combined1(a, b)	$a + b \cdot f$	somme des écarts-types
combined2(a, b)	$\sqrt{a^2 + (b \cdot f)^2}$
c34,79.3,68.167,158.7,102.5,238c34.3,79.3,51.8,119.3,52.5,120
c340,-704.7,510.7,-1060.3,512,-1067
l0 -0
c4.7,-7.3,11,-11,19,-11
H40000v40H1012.3
s-271.3,567,-271.3,567c-38.7,80.7,-84,175,-136,283c-52,108,-89.167,185.3,-111.5,232
c-22.3,46.7,-33.8,70.3,-34.5,71c-4.7,4.7,-12.3,7,-23,7s-12,-1,-12,-1
s-109,-253,-109,-253c-72.7,-168,-109.3,-252,-110,-252c-10.7,8,-22,16.7,-34,26
c-22,17.3,-33.3,26,-34,26s-26,-26,-26,-26s76,-59,76,-59s76,-60,76,-60z
M1001 80h400000v40h-400000z"/>​	somme des variances

DEFINITION:
y1 = {distribution=normal, prediction=Cc, errorModel=constant(a)}
y1 = {distribution=normal, prediction=Cc, errorModel=proportional(b)}
y1 = {distribution=normal, prediction=Cc, errorModel=combined1(a, b)}
y1 = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}

combined1 contre combined2. Ce n'est pas une question de goût. combined2 est le modèle que l'on obtient en supposant deux sources de bruit indépendantes — un plancher d'écart-type $a$ et un terme proportionnel d'écart-type $b f$. Indépendance, donc les variances s'additionnent : $g^2 = a^2 + (bf)^2$. C'est la forme qui a une justification statistique. combined1 additionne les écarts-types : aucune paire de sources indépendantes ne produit cela, c'est une paramétrisation linéaire commode de l'écart-type, rien de plus.

Leur écart se calcule exactement. En posant $r = bf/a$ le rapport des deux régimes :

Ce rapport vaut 1 aux deux extrémités ($r \to 0$ : les deux valent $a$ ; $r \to \infty$ : les deux valent $bf$) et il est maximal en $r = 1$, c'est-à-dire exactement à la concentration $f = a/b$ où les deux régimes s'équilibrent. Il y vaut $\sqrt{2}$. Donc : combined1 n'est jamais plus étroit que combined2, et il est au plus 41 % plus large, dans une bande étroite de concentrations au milieu de la gamme.

Deux conséquences pratiques, opposées :

1. Les deux modèles ajustent presque pareil. Leur différence est confinée au milieu de la gamme. Les départager sur un $\Delta$OFV de 2 points, c'est départager du bruit.
2. Leurs paramètres ne sont pas interchangeables. Un $a$ et un $b$ estimés sous combined1 ne décrivent pas la même erreur une fois recopiés dans combined2.

distribution : normal contre logNormal.

• normal : $h$ = identité, donc $y = f + g\varepsilon$. Le support est $\mathbb{R}$ tout entier : le modèle accorde une probabilité non nulle à des concentrations négatives. Anodin tant que $f \gg g$, absurde dès que $f$ approche la LOQ — et cela se voit, la bande basse de la VPC passe sous zéro.
• logNormal : $h = \log$, donc $\log(y) = \log(f) + g\varepsilon$, soit $y = f \cdot e^{g\varepsilon}$. Support strictement positif, distribution asymétrique à droite.

Avec logNormal et constant(a), on retrouve exactement le modèle d'erreur exponentiel : $y = f \cdot e^{a\varepsilon}$. Comme $e^{a\varepsilon} \approx 1 + a\varepsilon$ pour $a$ petit, l'écart-type sur l'échelle naturelle vaut environ $a \cdot f$ : le modèle se comporte comme un proportionnel, mais sans jamais rendre de valeur négative.

Corollaire : sous logNormal, constant(a) produit déjà une erreur proportionnelle à la prédiction. Empiler un proportional(b) par-dessus compte deux fois le même effet. La combinaison est légale syntaxiquement ; elle est rarement ce que vous vouliez dire.

La censure (BLQ). Une observation « < LOQ » n'est ni une valeur manquante, ni un nombre : c'est une inégalité. Elle affirme que la vraie concentration est quelque part dans $(0, LOQ)$ — et c'est une information réelle, souvent la seule dont vous disposiez sur la phase terminale.

Monolix la prend au mot, et le fait dans le jeu de données, pas dans le modèle. Vous déclarez une colonne de type censored (usuellement cens) : cens = 1 marque un point censuré à gauche, et la colonne d'observation porte alors la LOQ elle-même, pas une valeur imputée. Une colonne limit facultative fournit l'autre borne de l'intervalle.

Ce qui change dans la vraisemblance tient en un mot. Une observation ordinaire contribue par une densité :

une observation censurée contribue par une probabilité :

C'est la méthode M3 de Beal, et c'est tout : vous cochez un type de colonne, Monolix écrit la vraisemblance. Là où NONMEM demande un F_FLAG, un PHI à la main et un LAPLACIAN, la censure est ici un attribut de la donnée. Ce n'est pas un détail de confort : le coût d'entrée de M3 étant nul, il n'y a plus aucune excuse pour un LOQ/2.

Notez que $g$ apparaît aussi dans la contribution censurée. Le modèle d'erreur et la censure ne sont pas deux choix indépendants : c'est le terme additif $a$ qui gouverne la vitesse à laquelle $\Phi((LOQ-f)/g)$ sature quand $f$ passe sous la LOQ.

Un médicament oral, profils riches sur 24 h, concentrations de 4 mg/L au pic jusqu'à la LOQ à 0,05 mg/L — près de deux ordres de grandeur. Le SAEM tourne en combined2 et rend $a = 0{,}02$ mg/L et $b = 0{,}12$.

Les deux régimes, en chiffres. Le basculement se situe en $f = a/b = 0{,}167$ mg/L :

prédiction $f$	$g$ sous combined2	$g$ sous combined1	écart
0,05 (la LOQ)	0,021	0,026	+24 %
0,167 ($=a/b$)	0,028	0,040	+41 %
1,0	0,122	0,140	+15 %
4,0 (le pic)	0,480	0,500	+4 %

On lit directement le résultat de la section précédente : l'écart culmine à $\sqrt{2}$ au point de bascule et s'efface aux deux bouts.

Pourquoi le plancher existe. Prenons un creux tardif : observation $y = 0{,}09$ mg/L, prédiction individuelle $f = 0{,}06$ mg/L. L'écart est de 0,03 mg/L — 30 ng/mL, soit à peu près le bruit du dosage à ce niveau. Un point sans histoire.

• Sous combined2(0,02 ; 0,12) : $g = \sqrt{0{,}02^2 + 0{,}0072^2} = 0{,}0213$, donc $\text{IWRES} = 0{,}03/0{,}0213 = 1{,}41$. Rien à signaler, et c'est correct.
• Sous proportional(0,12) seul : $g = 0{,}12 \times 0{,}06 = 0{,}0072$ mg/L, soit 7,2 ng/mL. Le modèle vient d'affirmer qu'il résout cette concentration à 7 ng/mL près — sept fois mieux que la LOQ du dosage. $\text{IWRES} = 0{,}03/0{,}0072 = 4{,}17$. Un aberrant à 4 sigma, entièrement fabriqué par le modèle d'erreur.

Ce n'est pas cosmétique. À elle seule, cette observation pèse environ 13 unités de $-2LL$ de plus sous proportional que sous combined2. Le SAEM fait ce qu'on lui demande : il déformera $Cl$ et $V$ sur toute la population pour aller chercher ce creux. Le terme additif n'est pas un paramètre de nuisance — c'est l'affirmation que le dosage a un plancher de bruit, et c'est lui qui empêche la vraisemblance de prendre au sérieux des prédictions proches de zéro.

Comment on le voit. L'IWRES est le résidu divisé par l'écart-type que le modèle d'erreur revendique :

Cette division est tout l'intérêt du diagnostic : si $g$ est juste, les IWRES forment un nuage normal centré sur 0, d'écart-type 1, et de même largeur partout. D'où le graphique qui juge le modèle d'erreur — les IWRES contre les prédictions, pas contre le temps — et la règle de lecture : on lit la largeur, pas le centre.

• Nuage qui s'ouvre vers la droite (étroit en bas, large en haut) : $g$ est trop petit aux fortes concentrations. Vous êtes en constant, il manque le pourcentage.
• Nuage qui se resserre vers la droite (large en bas, étroit en haut) : $g$ est trop petit aux faibles concentrations. Vous êtes en proportional, il manque le plancher — c'est le creux à 4 sigma ci-dessus, répété sur tous les sujets.
• Nuage courbé, ou décentré de 0 : cela ne parle plus du modèle d'erreur. Un U contre le temps accuse le modèle structural. Aucun modèle d'erreur ne répare une faute de forme : il ne fait qu'élargir la bande jusqu'à ce que le défaut cesse d'être signalé.

L'erreur évidente — laisser constant sur un jeu qui couvre deux ordres de grandeur — est bruyante : le nuage des IWRES s'ouvre en trompette, et vous le voyez en trois secondes. Bruyante, donc bon marché.

Deux réflexes en découlent.

Lisez $\mathrm{SD}(\text{IWRES})$ avant la forme du nuage. C'est une ligne sur les résidus exportés par Monolix. Loin de 1, le graphique n'a plus de puissance de détection, et la platitude du nuage ne prouve rien sur $g$. La séquence est : d'abord le chiffre, ensuite seulement la forme.

Puis changez de diagnostic. L'$\varepsilon$-shrinkage frappe tout ce qui conditionne aux paramètres individuels. Passez aux diagnostics par simulation — VPC, NPDE : ils sont construits en simulant depuis le modèle de population, donc ils voient le modèle d'erreur que vous avez réellement déclaré, et non celui qui survit une fois chaque sujet ajusté à ses propres points.

• Le modèle d'erreur ne décrit pas du bruit : $g(f)$ est le poids de chaque point dans la vraisemblance. Il arbitre quelle partie du profil le SAEM a le droit de rater.
• Une ligne de DEFINITION: remplit deux trous indépendants de $h(y) = h(f) + g\varepsilon$ : distribution= choisit $h$, errorModel= choisit $g$.
• Les quatre $g$ : constant(a) $= a$ ; proportional(b) $= bf$ ; combined1(a,b) $= a+bf$ ; combined2(a,b) $= \sqrt{a^2+(bf)^2}$.
• combined2 additionne les variances (deux sources indépendantes) et c'est l'équivalent du Y = F + F*EPS(1) + EPS(2) de NONMEM ; combined1 additionne les écarts-types. combined1 n'est jamais plus étroit, au plus $\sqrt{2}$ fois plus large, exactement en $f = a/b$. Les $(a, b)$ ne se transposent pas de l'un à l'autre.
• logNormal + constant(a) = erreur exponentielle $y = f e^{a\varepsilon}$ : support positif, et a sur l'échelle log — un CV, pas des mg/L.
• BLQ : colonne de type censored, la LOQ dans la colonne d'observation. Le point contribue par $\Phi((LOQ-f)/g)$ — le M3 de Beal, natif, sans code, et donc plus aucune raison d'imputer LOQ/2.
• Jugez $g$ sur les IWRES contre prédictions, en lisant la largeur : nuage qui s'ouvre = pourcentage manquant ; qui se resserre = plancher manquant ; courbé = structural, pas résiduel.
• $\mathrm{SD}(\text{IWRES})$ avant la forme. $\varepsilon\text{-shrinkage} = 1 - \mathrm{SD}(\text{IWRES})$ : quand il est élevé, un beau nuage d'IWRES ne prouve rien — passez à la VPC et aux NPDE.