nlmixr2 — ini, model et la variabilité

La grammaire de ini({}) tient en trois signes. Ce ne sont pas des détails de style : chacun déclare une nature de paramètre, et c'est le seul moyen que nlmixr2 a de savoir à quoi il a affaire.

• <- (ou =) déclare un effet fixe : tcl <- log(4.5). Un theta, une valeur estimée pour toute la population.
• ~ déclare un effet aléatoire : eta.cl ~ 0.1. Un oméga, une dispersion autour du typique.
• | déclare un niveau de variabilité : iov.cl ~ 0.03 | occ. Un tirage par occasion et non plus par patient.

Le nom, lui, ne déclare rien. tcl n'est pas un mot réservé : le t est une convention humaine pour « theta de cl », lisible par vous, ignorée par la machine. Vous pourriez l'appeler bidule. Ce qui fait de tcl un paramètre log-transformé, ce n'est pas son nom ni sa valeur initiale log(4.5) — c'est l'exp() que vous écrivez autour de lui dans model({}).

Le modèle de base

mod <- function() {
  ini({
    tcl <- log(4.5)        # clairance typique, sur l echelle log
    tv  <- log(32)
    tka <- log(1.1)
    eta.cl ~ 0.1           # VARIANCE de l eta, pas un ecart-type
    eta.v  ~ 0.1
    eta.ka ~ 0.2
    prop.err <- 0.15
  })
  model({
    cl <- exp(tcl + eta.cl)
    v  <- exp(tv  + eta.v)
    ka <- exp(tka + eta.ka)
    d/dt(depot) = -ka*depot
    d/dt(centr) =  ka*depot - (cl/v)*centr
    cp = centr/v
    cp ~ prop(prop.err)
  })
}
fit <- nlmixr2(mod, dat, est = "saem")

soit, pour la clairance :

La convention du t, et pourquoi elle existe

Comparez à l'écriture NONMEM du même modèle, CL = THETA(1)*EXP(ETA(1)). Les deux décrivent la même log-normale, mais ils n'estiment pas le même nombre :

Chez NONMEM, l'EXP n'enveloppe que l'eta : THETA(1) est une clairance, en L/h, sur l'échelle naturelle. Chez nlmixr2, l'exp() enveloppe theta et eta ensemble : tcl n'est pas une clairance, c'est un logarithme de clairance, sans unité. tcl = 1.52 ne veut rien dire tant qu'on ne l'a pas exponentié en 4,57 L/h.

Cette convention n'est pas obligatoire — cl <- tcl * exp(eta.cl) avec tcl <- 4.5 fonctionne aussi. De nombreux exemples placent néanmoins les theta sur l'échelle log : un paramètre exponentié n'a pas besoin d'une borne uniquement pour garantir sa positivité. La forme du modèle assure cette propriété et laisse l'optimiseur travailler sur une échelle non bornée.

La valeur après le tilde est une variance

eta.cl ~ 0.1 déclare $\omega_{cl}^2 = 0{,}1$. C'est la même échelle que l'$OMEGA de NONMEM, et l'inverse de la sd de Monolix. Pour une log-normale, le CV exact du paramètre vaut :

ini ($\omega^2$)	$\omega$	$CV$ approché	$CV$ exact
0,04	0,20	20 %	20,2 %
0,09	0,30	30 %	30,7 %
0,10	0,32	32 %	32,4 %
0,16	0,40	40 %	41,7 %
0,50	0,71	71 %	80,5 %

Le ~ 0.1 qui traîne dans tous les tutoriels — y compris dans le chapitre précédent — n'est donc pas un chiffre magique : c'est 32 % de CV, une supposition de départ délibérément raisonnable pour un paramètre PK. Vous avez maintenant de quoi la remplacer par la vôtre.

Le bloc de covariance

Par défaut, les etas sont indépendants. On les corrèle en les additionnant à gauche du tilde :

ini({
  tcl <- log(4.5)
  tv  <- log(32)
  eta.cl + eta.v ~ c(0.1,
                     0.05, 0.1)
  eta.ka ~ 0.2
})

Le c() donne le triangle inférieur de la matrice de covariance, ligne par ligne — exactement comme le $OMEGA BLOCK de NONMEM :

L'ordre est donc var(eta.cl), cov(eta.cl, eta.v), var(eta.v) : la covariance est au milieu, pas à la fin. Si vous préférez déclarer ce que vous lisez plutôt que ce qui est estimé, cor() accepte les écarts-types sur la diagonale et la corrélation hors diagonale :

  eta.cl + eta.v ~ cor(0.32,
                       0.5, 0.32)

Les deux blocs ci-dessus décrivent la même matrice. Seule change la façon dont vous écrivez la valeur initiale ; le paramètre estimé reste la covariance.

L'IOV se déclare comme un niveau

nlmixr2 ne code pas l'IOV avec des etas répétés. Il la déclare comme un niveau de variabilité, avec une barre verticale :

ini({
  tcl <- log(4.5)
  eta.cl ~ 0.1
  iov.cl ~ 0.03 | occ      # entre occasions, pas entre patients
})
model({
  cl <- exp(tcl + eta.cl + iov.cl)
  ...
})

| occ se lit « variabilité entre occ », par symétrie avec l'IIV qui est la variabilité entre id. La colonne occ est une colonne ordinaire du jeu de données, que vous construisez vous-même (par exemple dat$occ <- ifelse(dat$TIME < 336, 1, 2)).

Les covariables sont du R ordinaire

Ici nlmixr2 est plus direct que ses deux concurrents : il n'y a pas de bloc de covariables. Toute colonne du jeu de données est directement utilisable dans model({}), et la transformation s'écrit sur place.

ini({
  tcl <- log(4.5)
  beta.cl.crcl <- 0.7
  beta.cl.sex  <- 0.1
  eta.cl ~ 0.1
})
model({
  cl <- exp(tcl + beta.cl.crcl*log(CRCL/90) + beta.cl.sex*SEX + eta.cl)
})

avec SEX codée 0/1 dans les données. Regardez ce que l'échelle log fait gratuitement : la covariable entre additivement sur le log, donc multiplicativement sur le paramètre.

Une covariable vise à expliquer une partie de la variabilité systématique portée jusque-là par l'eta. exp(tcl) devient la clairance d'un patient de référence — ici 90 mL/min de CRCL et SEX = 0 — et eta.cl porte le reste. Une baisse de l'oméga peut soutenir son utilité, mais doit être interprétée avec sa plausibilité, sa précision, les diagnostics et les performances prédictives.

44 patients, administration orale, un compartiment avec absorption d'ordre 1, erreur proportionnelle. Deux occasions de prélèvement par patient : un profil riche au jour 1, deux points résiduels au jour 21. On construit le modèle statistique par étapes, en SAEM.

Run	Modèle de variabilité	Paramètres	OFV	$\Delta$
1	trois etas indépendants	7	1682,4	—
2	+ bloc eta.cl + eta.v ~ c(...)	8	1670,1	−12,3
3	+ IOV sur cl au niveau occ	9	1651,8	−18,3
4	+ CRCL sur cl	10	1633,6	−18,2

Run 1 → 2. La corrélation sort à $r = 0{,}45$. Un paramètre de plus, l'OFV baisse de 12,3, très au-dessus du seuil de 3,84 à 1 degré de liberté : le bloc est retenu.

Run 2 → 3. iov.cl s'estime à 19,2 % de CV, et le BSV de la clairance descend de 41,2 % à 35,1 %. C'est le résultat le plus instructif du tableau : une part de ce qu'on attribuait à « ce patient élimine fort » était en réalité « cette visite-là était différente ». Sans niveau d'occasion, l'IIV absorbe l'IOV et se retrouve surestimée — et un modèle qui surestime l'IIV surestime aussi la dispersion des simulations dont vous tirerez une recommandation de dose.

Run 3 → 4. L'exposant sur la clairance de la créatinine sort à 0,71, et le BSV de la clairance passe de 35,1 % à 29,9 %. La fonction rénale explique donc environ 5 points de CV. C'est cette phrase-là, et non le $\Delta$ de 18,2, qui a un sens clinique et qui ira dans le rapport.

La table des paramètres

C'est la vraie force ergonomique de nlmixr2 : tout est dans une seule table, déjà rétro-transformée.

── nlmixr² SAEM(ODE); OBJF by FOCEi approximation ──

  Parameter      Est.     SE  %RSE  Back-transformed(95%CI)  BSV(CV%)  Shrink(SD)%
  tcl            1.52  0.043   2.8       4.57 (4.20, 4.97)       29.9        7.6%
  tv             3.46  0.048   1.4       31.8 (28.9, 35.0)       39.1       12.8%
  tka           0.131  0.104    79       1.14 (0.93, 1.40)       52.4       45.3%
  beta.cl.crcl   0.71  0.128    18       0.71 (0.46, 0.96)
  prop.err      0.116  0.009   7.8                    0.116
  iov.cl                                                          19.2

Trois colonnes méritent qu'on s'y arrête.

Est. est sur l'échelle estimée, Back-transformed sur l'échelle utile. tcl = 1,52 n'est pas reportable ; $e^{1{,}52} = 4{,}57$ L/h l'est. Cette colonne n'a rien de magique : nlmixr2 l'a produite en lisant l'exp() de votre model({}). La ligne beta.cl.crcl le montre bien — comme aucune exponentielle ne l'entoure, la valeur rétro-transformée est identique à l'estimation, seul l'intervalle de confiance est ajouté.

BSV(CV%) vous épargne la racine carrée. 29,9 % est le $\sqrt{e^{\omega^2}-1}$ de l'oméga de eta.cl. Notez que iov.cl apparaît sur sa propre ligne, avec sa dispersion mais sans estimation ni erreur standard : c'est un niveau de variabilité, pas un effet fixe.

Les objets à interroger

La table imprimée résume ; l'objet fit contient tout, sous forme de matrices R directement réutilisables.

fit$omega     # covariance complete des effets aleatoires
fit$omegaR    # meme matrice : SD sur la diagonale, correlations hors diagonale
fit$shrink    # shrinkage et statistiques de distribution des etas
fit$eta       # les etas individuels, un par patient
fit$iov       # les kappa individuels, par patient et par occasion

> fit$omega
       eta.cl  eta.v eta.ka
eta.cl 0.0856 0.0497 0.0000
eta.v  0.0497 0.1422 0.0000
eta.ka 0.0000 0.0000 0.2426

> fit$omegaR
       eta.cl  eta.v eta.ka
eta.cl 0.2926 0.4500 0.0000
eta.v  0.4500 0.3771 0.0000
eta.ka 0.0000 0.0000 0.4925

fit$omega et fit$omegaR décrivent la même chose dans deux langues. La première est celle que le modèle estime ; la seconde est celle que vous savez lire — on y voit d'un coup d'œil $\omega_{cl} = 0{,}29$, $\omega_{v} = 0{,}38$ et $r_{cl,v} = 0{,}45$, sans diviser quoi que ce soit. Les zéros de la colonne eta.ka ne sont pas des estimations proches de zéro : ce sont des cases jamais estimées, puisque eta.ka a été déclaré hors du bloc.

Et le shrinkage, déjà présent dans la colonne Shrink(SD)% :

$\eta_{cl}$ est à 7,6 % : informés par la courbe entière, ses diagnostics individuels sont ici plus informatifs, sans que ce chiffre garantisse leur exactitude. $\eta_{ka}$ est à 45,3 %, et la raison se lit dans le protocole — seul le jour 1 porte des prélèvements précoces. Chez les patients concernés, l'estimation individuelle retombe vers la population.

Vous traduisez en nlmixr2 un modèle NONMEM qui tourne depuis des années. Le control stream contient :

$THETA (0, 1.2)          ; KA
$PK
  KA = THETA(3)*EXP(ETA(3))

Le (0, ...) est ici parfaitement correct et parfaitement idiomatique : THETA(3) est une constante d'absorption en h⁻¹, elle doit rester positive, on le lui dit. Vous transposez fidèlement :

ini({
  tka <- c(0, log(1.2))     # « comme dans NONMEM : borne basse a 0 »
  eta.ka ~ 0.2
})
model({
  ka <- exp(tka + eta.ka)
})

Vous lancez en est = "focei". Le run se termine. tka sort à 0,0000, ka rétro-transformée à 1,00 h⁻¹ — pile. La SE est aberrante, le %RSE part en vrille. Vous relancez le même fichier en est = "saem" : cette fois tka = -0,478, soit 0,62 h⁻¹, et les diagnostics de la phase d'absorption redeviennent propres.

Deux méthodes, un seul fichier, deux réponses. La conclusion tentante — « SAEM et FOCEI ne sont pas d'accord, c'est un artefact d'algorithme » — est fausse, et elle vous ferait passer des jours à comparer des moteurs alors que le bug est dans une parenthèse.

Reste l'asymétrie, qui est le vrai danger. SAEM ne gère pas les bornes et les ignore — nlmixr2 vous en avertit d'ailleurs dans la console. FOCEI, lui, les respecte : son optimiseur externe accepte des contraintes de boîte, et il ne dit rien puisqu'il fait exactement ce qu'on lui a demandé. Le moteur qui vous prévient est donc celui qui n'avait pas de problème, et le moteur qui vous donne le mauvais chiffre est celui qui se tait. Un avertissement noyé dans un mur de sortie R, et l'affaire est jouée.

Le correctif tient en trois caractères :

tka <- log(1.2)           # aucune borne, et c est volontaire

• La grammaire de ini({}) tient en trois signes : <- déclare un effet fixe, ~ un effet aléatoire, | un niveau de variabilité. Le nom tcl ne déclare rien — le t est une convention humaine.
• Ce qui rend tcl log-transformé, c'est l'exp() de model({}), pas sa valeur initiale. nlmixr2 relit votre model({}) pour rétro-transformer ses sorties : il constate ce que vous avez écrit, il ne corrige jamais un désaccord entre les deux blocs.
• La valeur après le ~ est une variance, comme l'$OMEGA de NONMEM et à l'inverse de la sd de Monolix. $CV = \sqrt{e^{\omega^2}-1}$ ; le ~ 0.1 des tutoriels vaut 32 % de CV ; une IIV usuelle vit entre 0,04 et 0,25.
• eta.cl + eta.v ~ c(0.1, 0.05, 0.1) donne le triangle inférieur de la covariance : la covariance est au milieu. cor() permet de saisir SD et corrélation à la place.
• L'IOV est un niveau (iov.cl ~ 0.03 | occ) : une seule variance quel que soit le nombre d'occasions, donc pas de SAME à écrire. Sans elle, l'IIV absorbe l'IOV et se retrouve surestimée.
• Les covariables sont du R ordinaire dans model({}), sans bloc dédié. Sur l'échelle log, beta*log(CRCL/90) est le modèle puissance. Une baisse d'oméga peut soutenir l'effet, avec les autres critères de sélection et de validation.
• Sur un theta log-transformé, le %RSE est trompeur près de zéro : jugez sur l'intervalle rétro-transformé. fit$omegaR se lit sans calcul ; l'information des graphiques eta-covariable décroît avec le shrinkage, 20–30 % restant une heuristique d'alerte.
• Ne mettez aucune borne sur un paramètre enveloppé d'exp() : elle est inutile et son échelle est trompeuse. SAEM les ignore en vous prévenant, FOCEI les applique en silence — le moteur dangereux est celui qui se tait.