nlmixr2 — le modèle structural en R

La vraie question, quand on découvre un outil, est de savoir où il coupe. nlmixr2 ne coupe pas là où Monolix coupe.

• ini({...}) ne contient que des nombres : où part la recherche, et quel paramètre est aléatoire. Aucune équation n'y a sa place.
• model({...}) contient toutes les équations : la valeur typique, l'effet aléatoire, les EDO, la concentration, et la ligne d'observation.

La conséquence est immédiate et déroute tous ceux qui arrivent de Monolix : cl <- exp(tcl + eta.cl) s'écrit dans model, à trois lignes des d/dt(). La variabilité inter-individuelle — de la statistique pure — vit dans le même bloc que la structure. Il n'y a pas de [INDIVIDUAL] à chercher : la loi lognormale s'écrit à la main. C'est le découpage de NONMEM ($THETA/$OMEGA d'un côté, $PK/$DES de l'autre), transposé en R.

Le 1 compartiment oral tient en deux quantités — le dépôt $A_d$, le central $A_c$ — et une division :

Version EDO. Le système câblé à la main :

oral1cpt <- function() {
  ini({
    tka <- log(0.9)      # <- : parametre de POPULATION. ka = 0.9 /h, pose sur
    tcl <- log(6.0)      #      l'echelle log pour garantir une valeur positive
    tv  <- log(45)
    eta.cl ~ 0.09        # ~  : effet ALEATOIRE ; 0.09 est une VARIANCE
    eta.v  ~ 0.04
    prop.sd <- 0.15      # <- : erreur residuelle, sur l'echelle ECART-TYPE
  })
  model({
    ka <- exp(tka)             # pas d'eta ici : ka est suppose sans IIV
    cl <- exp(tcl + eta.cl)    # la loi lognormale s'ecrit A LA MAIN
    v  <- exp(tv  + eta.v)

    d/dt(depot)   = -ka * depot                     # declarer d/dt CREE le compartiment
    d/dt(central) =  ka * depot - (cl/v) * central  # les EDO portent des QUANTITES

    cp = central / v           # la division par V est A VOTRE charge
    cp ~ prop(prop.sd)         # ligne d'OBSERVATION : cp est confronte a DV
  })
}

Ligne à ligne :

• Dans ini, <- déclare un paramètre de population, ~ déclare un effet aléatoire et sa variance. Il n'y a pas de matrice OMEGA à écrire : elle se déduit des lignes en ~. Pour des etas corrélés, on donne le triangle inférieur : eta.cl + eta.v ~ c(0.09, 0.03, 0.04).
• Toujours dans ini, tcl <- c(-Inf, log(6.0), Inf) fixe borne basse, valeur initiale et borne haute ; label("Clairance (L/h)") nomme la ligne dans la table de sortie ; fix() gèle le paramètre.
• Dans model, écrire d/dt(depot) suffit à créer le compartiment depot — il n'y a rien à déclarer ailleurs. Les compartiments sont numérotés dans leur ordre d'apparition, ce qui rend un cmt numérique fragile : réordonnez deux lignes et le jeu de données ne dose plus au bon endroit. Utilisez les noms dans la colonne cmt.
• cp n'est pas un mot-clé : c'est le nom que vous donnez à la prédiction.

Version linCmt(). Le bloc ini ne change pas ; seul model maigrit :

  model({
    ka <- exp(tka)
    cl <- exp(tcl + eta.cl)
    v  <- exp(tv  + eta.v)
    cp = linCmt()            # ka + cl + v en portee -> 1 cpt oral, rend une CONCENTRATION
    cp ~ prop(prop.sd)
  })

Le modèle écrit, le reste du script est du R banal :

library(nlmixr2)

dat <- read.csv("pk_oral.csv")     # ID TIME AMT DV EVID CMT ; dose 200 mg dans depot

fit <- nlmixr2(oral1cpt, dat, est = "saem",     # oral1cpt SANS parentheses
               control = saemControl(nBurn = 300, nEm = 400, seed = 42))

fit$parFixed    # table des parametres de population, RSE, CV des etas
head(fit)       # data.frame : ID TIME DV PRED IPRED CWRES eta.cl eta.v cl v cp ...

Avant de lire quoi que ce soit, vérifiez que la machine fait ce que vous croyez. Pour un sujet typique — dose orale 200 mg, $k_a = 0{,}9$ h⁻¹, $V = 45$ L, $Cl = 6{,}0$ L/h, $F = 1$ — la constante d'élimination vaut $k_e = Cl/V \approx 0{,}133$ h⁻¹, d'où un rapport $k_a/k_e = 6{,}75$ et un pic à

Comme $C_{max} = \frac{F \cdot D}{V}e^{-k_e t_{max}}$ pour ce modèle, on attend $C_{max} \approx 4{,}44 \times e^{-0{,}332} \approx 3{,}19$ mg/L, pour une exposition $AUC = F \cdot D / Cl \approx 33{,}3$ mg·h/L. Trois nombres calculés en trente secondes, qui disqualifient un modèle mal câblé avant qu'il ne vous coûte une demi-journée de diagnostics.

fit n'est pas un fichier de sortie à ouvrir dans un éditeur : c'est un data.frame augmenté, une ligne par observation, avec PRED, IPRED, CWRES et les etas déjà en colonnes, plus des tables (fit$parFixed, fit$omega) qui partent directement dans ggplot2. C'est là que se paie le choix du modèle-objet.

Les fautes de syntaxe, dans nlmixr2, sont bruyantes : un paramètre déclaré dans ini mais absent de model fait échouer l'analyse avant même le premier run. Le piège coûteux est ailleurs, et il est muet.

Le symptôme, lui, ne ressemble pas du tout à la cause. Ce que vous verrez, si vous voyez quelque chose, c'est une étape de covariance qui échoue ou deux paramètres de population aux RSE aberrantes — la signature de leur confusion. Vous partirez chercher un problème de données ou d'identifiabilité. L'erreur est un tilde.

Le même bloc ini cache une seconde asymétrie, purement de lecture celle-là. eta.cl ~ 0.09 donne une variance : l'écart-type vaut 0,3 et le CV lognormal environ 31 %, pas 9 %. Deux lignes plus bas, prop.sd <- 0.15 donne un écart-type : 15 % d'erreur proportionnelle, là où NONMEM vous ferait écrire 0,0225 dans $SIGMA. Même bloc, même syntaxe, deux échelles. D'où le double réflexe : convertir avant de recopier une valeur initiale trouvée dans un run NONMEM, et convertir avant d'annoncer une IIV dans un rapport.

• Le modèle est une fonction R à deux blocs : ini({...}) ne contient que des nombres, model({...}) que des équations. nlmixr2 lit le corps de la fonction, il ne l'exécute pas — d'où nlmixr2(oral1cpt, ...), sans parenthèses.
• Le découpage est celui de NONMEM, pas celui de Monolix : cl <- exp(tcl + eta.cl) s'écrit à la main, dans le même bloc que les EDO. Aucun [INDIVIDUAL] à chercher.
• Dans ini, <- = paramètre de population, ~ = effet aléatoire et sa variance ; l'erreur résiduelle, elle, se déclare en écart-type.
• d/dt(depot) crée le compartiment par sa seule déclaration, les numéros suivent l'ordre d'apparition — dosez par nom — et la division par $V$ est à votre charge.
• linCmt() ne prend aucun argument : il choisit le modèle d'après les noms des paramètres en portée (ka, cl, v = 1 cpt oral), route la dose et rend une concentration.
• cp ~ prop(prop.sd) est la ligne d'observation ; le tilde change de sens selon son bloc.
• Un <- mis à la place d'un ~ dans ini supprime une IIV sans rien casser — et le symptôme désigne un autre coupable.