nlmixr2 — le modèle d'erreur résiduelle et le BLQ

Lisez cp ~ add(add.err) comme une phrase française : « cp est observée avec une erreur additive, dont l'écart-type s'appelle add.err ». Le tilde est exactement celui de lm(y ~ x) : à gauche ce qu'on observe, à droite le modèle de cette observation. nlmixr2 emprunte la grammaire de R plutôt que d'en inventer une.

Cette lecture a une conséquence pratique immédiate. cp = centr/v est une affectation : après cette ligne, cp vaut quelque chose. cp ~ add(add.err) ne change la valeur de rien du tout — elle branche cp sur la colonne d'observations et dit à l'estimateur quelle barre d'erreur mettre autour. Une fonction model() sans ligne à tilde ne définit aucune vraisemblance ; ce n'est plus un modèle à estimer, c'est un simulateur.

Reste à répondre à la seule question que pose le modèle d'erreur : quelle largeur a la barre, et comment varie-t-elle avec la concentration ? Un dosage réel a deux régimes de bruit, sans rapport 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 add() ;
• un pourcentage, qui grandit avec la concentration — dilutions, pipetage, calibration. C'est prop().

Un profil de PK couvre couramment deux ou trois ordres de grandeur, du pic au dernier creux. Il traverse donc les deux régimes, et aucune forme à un seul terme n'est juste sur toute la gamme. C'est très exactement pourquoi add() + prop() est un candidat fréquent, à confirmer sur les résidus.

Toutes les formes de nlmixr2 se rangent sous une écriture unique :

La ligne à tilde ne fait que choisir $g$, la fonction d'écart-type.

déclaration dans model()	$g(f)$	ce que ça dit
cp ~ add(a)	$a$	plancher seul, en mg/L
cp ~ prop(b)	$b \cdot f$	pourcentage seul
cp ~ pow(b, c)	$b \cdot f^{c}$	puissance, exposant $c$ estimé
cp ~ add(a) + prop(b)	$\sqrt{a^2 + (bf)^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"/>​ ou $a + bf$	combiné — voir plus bas
cp ~ add(a) + pow(b, c)	$\sqrt{a^2 + (bf^{c})^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"/>​ ou $a + bf^{c}$	combiné à exposant libre
cp ~ lnorm(s)	erreur exponentielle	$\log y = \log f + s,\varepsilon$

prop(b) n'est que le cas $c = 1$ de pow(b, c). Laisser l'exposant libre revient à demander aux données où se situe le régime réel du dosage entre le plancher pur ($c = 0$) et le pourcentage pur ($c = 1$) ; c'est souvent un paramètre mal identifié, à ne sortir que si la gamme est large et les données abondantes.

lnorm(s) donne $y = f \cdot e^{s\varepsilon}$ : support strictement positif, distribution asymétrique à droite. Comme $e^{s\varepsilon} \approx 1 + s\varepsilon$ pour $s$ petit, le modèle se comporte comme un proportionnel de CV $\approx s$, mais sans jamais rendre de valeur négative. Le paramètre s vit sur l'échelle log : il est sans dimension. Le lire comme des mg/L produit un nombre plausible et faux. Corollaire : lnorm() produit déjà une erreur proportionnelle à la prédiction — lui empiler un prop() compte deux fois le même effet.

Le bloc ini() : deux opérateurs, deux rôles.

mod <- function() {
  ini({
    tka <- log(1.1)
    tcl <- log(2.8)
    tv  <- log(32)
    eta.cl ~ 0.09              # effet aleatoire : variance
    eta.v  ~ 0.04
    add.err  <- c(0, 0.03)     # ecart-type, mg/L, borne inferieure a 0
    prop.err <- c(0, 0.12)     # ecart-type relatif, fraction
  })
  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 ~ add(add.err) + prop(prop.err) + combined2()
  })
}
fit <- nlmixr2(mod, data, est = "saem")

Le chapitre précédent a posé la règle : dans ini(), le ~ déclare un effet aléatoire par sa variance, le <- déclare un effet fixe. Les paramètres résiduels n'ont pas de troisième voie — ils prennent le <-, exactement comme les theta. Écrire add.err ~ 0.03 ne produit donc pas un terme résiduel : cela crée un ETA nommé add.err dont plus personne ne se sert. C'est la faute de débutant classique, et elle ne fait pas planter grand-chose.

Ce qui mène au point de nommage. nlmixr2 n'a pas de bloc réservé aux paramètres d'erreur. add.err et prop.err sont des noms de convention, pas des mots-clés : appelez-les sigma.floor et cv.assay, rien ne change. Un scalaire déclaré dans ini() devient un paramètre résiduel parce qu'il est consommé par une fonction d'erreur dans la ligne à tilde, et pour aucune autre raison. Le rôle vient de l'usage, pas de la déclaration.

Deux écritures utiles dans ini() : les bornes, add.err <- c(0, 0.03) = borne inférieure 0 et valeur initiale 0,03 (la forme à trois éléments ajoute une borne supérieure) ; et le gel, add.err <- fixed(0.05), pour imposer un plancher de bruit connu du laboratoire au lieu de l'estimer — utile quand les données basses sont trop rares pour l'informer.

Le combiné, et où vit le choix. Deux façons d'assembler un plancher et un pourcentage :

combined2 est la forme qui découle de deux sources de bruit indépendantes : indépendance, donc les variances s'additionnent. C'est celle qui a une justification statistique, c'est l'équivalent exact du Y = F + F*EPS(1) + EPS(2) de NONMEM, et c'est le défaut de nlmixr2. combined1 additionne les écarts-types : aucune paire de sources indépendantes ne produit cela, c'est une paramétrisation commode, rien de plus. Les deux ne diffèrent notablement qu'autour de $f = a/b$, où les régimes s'équilibrent — au plus d'un facteur $\sqrt{2}$.

Le point nlmixr2, lui, n'est pas mathématique. Il est que ce choix se règle aussi par l'option addProp de l'objet de contrôle, hors du modèle :

nlmixr2(mod, data, est = "focei",
        control = foceiControl(addProp = "combined1"))

Écrire + combined2() explicitement sur la ligne à tilde, comme dans le code plus haut, ancre le choix dans le modèle et le rend insensible au contrôle. Ce n'est pas une décoration : voyez le piège fréquent.

Un inhibiteur de kinase oral, 42 patients, 336 concentrations, de 12 mg/L au pic jusqu'à une LOQ à 0,10 mg/L — plus de deux ordres de grandeur. Le SAEM tourne sur le modèle ci-dessus et rend le tableau des paramètres de population :

paramètre	Est.	%RSE	rétro-transformé
tka	0,10	14	$k_a$ = 1,11 h⁻¹
tcl	1,03	5,2	$Cl$ = 2,80 L/h
tv	3,47	3,1	$V$ = 32,1 L
add.err	0,031	9,4	0,031 mg/L
prop.err	0,118	4,6	11,8 %

Trois choses se lisent directement, sans calcul intermédiaire.

Les unités. add.err = 0,031 mg/L est un plancher de bruit, prop.err = 0,118 un CV de 11,8 %. Ce sont des écarts-types : la même erreur exprimée en NONMEM afficherait un $SIGMA de 0,000961 et 0,0139, et il faudrait deux racines carrées pour retrouver ces chiffres.

La cohérence. Le repère coûte deux secondes : add.err doit atterrir dans le voisinage du bruit du dosage, donc du même ordre que la LOQ. Ici 0,031 pour une LOQ à 0,10 — cohérent. Si le SAEM vous rend add.err = 0,9 mg/L avec la même LOQ, le paramètre ne mesure plus le dosage : il éponge un défaut structural. Regardez la courbe avant d'accepter le chiffre.

Le point de bascule. Il est en $f = a/b = 0{,}031/0{,}118 = 0{,}26$ mg/L. Au-dessus, le pourcentage commande ; au-dessous, le plancher. Un creux prédit à 0,15 mg/L reçoit $g = \sqrt{0{,}031^2 + 0{,}0177^2} = 0{,}036$ mg/L. Sous prop(prop.err) seul, il aurait reçu $g = 0{,}118 \times 0{,}15 = 0{,}018$ mg/L : le modèle affirmerait résoudre cette concentration à 18 ng/mL près, six fois mieux que sa propre LOQ. C'est ce que le terme additif empêche.

Le BLQ. Le laboratoire renvoie 47 points sur 336 (14 %) sous la LOQ, presque tous des creux tardifs. Une observation « < LOQ » n'est ni une valeur manquante ni un nombre : c'est une inégalité, et souvent la seule information dont vous disposiez sur la phase terminale.

nlmixr2 la prend au mot, dans le jeu de données, pas dans le modèle. Deux colonnes suffisent : CENS marque la ligne (0 = observation ordinaire, 1 = censurée à gauche, -1 = censurée à droite) et LIMIT fournit l'autre borne de l'intervalle. Sur une ligne censurée, la colonne DV porte la LOQ elle-même.

ID  TIME   DV     AMT  EVID  CMT    CENS  LIMIT
7    0.0   .      200   1    depot   0     .
7    1.0   3.42   .     0    cp      0     .
7    8.0   0.61   .     0    cp      0     .
7   24.0   0.10   .     0    cp      1     .      <- BQL : DV porte la LOQ

La ligne censurée n'entre plus par une densité mais par une probabilité :

C'est la méthode M3 de Beal, obtenue en remplissant une colonne. Ajouter LIMIT = 0 sur ces mêmes lignes borne l'intervalle à $(0;\ 0{,}10)$ au lieu de $(-\infty;\ 0{,}10)$ : c'est M4. Une colonne sépare les deux. Là où NONMEM réclame un F_FLAG, un PHI codé à la main et un LAPLACIAN, la censure est ici un attribut de la donnée — et le coût d'entrée de M3 étant nul, il n'y a plus d'excuse pour un LOQ/2.

Ce que ça change sur ce jeu, en refaisant tourner les deux :

	$Cl$ (L/h)	add.err (mg/L)	%RSE de add.err
M1 — les 47 points écartés	2,55	0,009	71
M3 — CENS = 1	2,80	0,031	9,4

Neuf pour cent d'écart sur la clairance, et ce n'est pas un hasard : les BQL ne manquent pas au hasard. À un temps tardif donné, seuls survivent les patients dont la concentration est au-dessus de la moyenne. On ne supprime pas du bruit, on tronque la queue par le bas — la pente terminale paraît plus plate et $Cl$ sort sous-estimée.

Le tell est dans la dernière colonne. Sans les BQL, add.err n'a plus aucune donnée basse pour l'informer : il s'effondre à 0,009 mg/L et son %RSE passe à 71. nlmixr2 vous imprime gratuitement le signe que le paramètre n'est plus identifié. Un %RSE à trois chiffres sur un terme résiduel additif est presque toujours cette histoire-là.

Les diagnostics. L'IWRES est le résidu divisé par l'écart-type que le modèle d'erreur revendique, donc le seul graphique qui juge $g$ :

Si $g$ est juste, le nuage est centré sur 0, d'écart-type 1, et de même largeur partout. D'où la règle de lecture : IWRES contre prédictions — pas contre le temps — et on lit la largeur, pas le centre. Nuage qui s'ouvre vers la droite : $g$ trop petit en haut de gamme, il manque le pourcentage. Nuage qui se resserre vers la droite : $g$ trop petit en bas de gamme, il manque le plancher. Nuage courbé contre le temps : cela ne parle plus du modèle d'erreur mais du modèle structural.

Le fit est un objet R, donc les deux écosystèmes de diagnostic s'y branchent directement :

fit <- addCwres(fit)    # CWRES/CPRED : non calcules par defaut apres un SAEM
fit <- addNpde(fit)     # NPDE, par simulation depuis le modele de population

library(xpose.nlmixr2)
xpdb <- xpose_data_nlmixr2(fit)
absval_res_vs_pred(xpdb, res = "IWRES")   # la largeur, pas le centre
dv_vs_ipred(xpdb)

library(ggPMX)
ctr <- pmx_nlmixr(fit)
pmx_plot_iwres_ipred(ctr)
pmx_plot_npde_time(ctr)
pmx_report(ctr, name = "diag", save_dir = ".", format = "html")

vpcPlot(fit, n = 500)
sd(fit$IWRES)           # a lire AVANT la forme du nuage

Reprenez la ligne d'erreur la plus banale de nlmixr2, écrite sans modificateur :

cp ~ add(add.err) + prop(prop.err)

Cette ligne ne dit pas comment les deux termes se combinent. Elle nomme deux paramètres et s'arrête là. Le reste — somme des variances ou somme des écarts-types — se décide dans l'objet de contrôle, par l'option addProp, dont le défaut est "combined2" :

fit1 <- nlmixr2(mod, data, est = "focei")
fit2 <- nlmixr2(mod, data, est = "focei",
                control = foceiControl(addProp = "combined1"))

Même fonction mod, même jeu de données, deux modèles d'erreur différents. Sur le jeu de la section précédente, fit1 rend add.err = 0,031 et prop.err = 0,118 ; fit2 rend 0,021 et 0,109 — les estimations se déplacent pour compenser la forme plus large de combined1. Et l'OFV bouge de deux ou trois points, c'est-à-dire de rien : les deux formes ne diffèrent qu'au voisinage de $f = a/b$, donc elles ajustent presque pareil. Rien dans la sortie ne crie.

Le principe déborde largement le cas d'addProp : en nlmixr2, tout ce qui vit dans le contrôle échappe au modèle. C'est vrai aussi de l'OFV.

• La ligne à tilde de model() ne calcule rien : elle déclare comment la prédiction est observée. Sans elle, pas de vraisemblance — juste un simulateur.
• $g(f)$ est le poids de chaque point dans la vraisemblance, pas une description du bruit. Il arbitre quelle partie du profil le SAEM a le droit de rater.
• Les formes : add(a) $= a$ ; prop(b) $= bf$ ; pow(b,c) $= bf^{c}$ ; lnorm(s) $\Rightarrow y = f e^{s\varepsilon}$, où s est un CV sur l'échelle log, pas des mg/L.
• nlmixr2 estime des écarts-types dans les unités de la mesure : prop.err = 0,118 se lit « CV de 11,8 % », sans racine carrée. Pas de SIGMA, pas de 1 FIX, pas de W — le confort que NONMEM demande de construire à la main.
• ini() : ~ pour un effet aléatoire, <- pour tout le reste. Un scalaire devient résiduel parce qu'une fonction d'erreur le consomme, pas parce qu'il est déclaré quelque part de spécial. Les noms sont libres.
• add() + prop() seul est ambigu : la combinaison vient de addProp dans le contrôle (défaut combined2, la somme des variances, équivalente au combiné de NONMEM). Écrivez + combined2() dans la ligne à tilde pour que le modèle se lise seul.
• BLQ : colonnes CENS (1 = censuré à gauche) et LIMIT, la LOQ dans DV. Le point contribue par $\Phi((LOQ-f)/g)$ — le M3 de Beal sans une ligne de code, M4 en ajoutant LIMIT = 0. DV = 0 sur une ligne censurée est une affirmation impossible qui fait dériver le run.
• Jugez $g$ sur les IWRES contre prédictions, en lisant la largeur : qui s'ouvre = pourcentage manquant ; qui se resserre = plancher manquant ; courbé = structural. Mais lisez sd(fit$IWRES) d'abord — sous $\varepsilon$-shrinkage, un beau nuage ne prouve rien.
• Un %RSE énorme sur add.err signale un plancher que plus aucune donnée basse n'informe — typiquement des BQL jetés.