Posthoc: two questions

[discovererspirit]

Hi everyone !
I'm a beginner and trying to understand the posthoc fitting and apply it to therapeutic drug monitoring.

This is my code (also attached .txt file)
code posthoc.txt

library(nlmixr2)
mod.dos.cmpt <- function() {
  ini({
    # *1* Efectos fijos "theta" + Error
    # A) Parámetros poblacionales y su valor poblacional
    ltcl <- log(5.42 * 0.001) #Ln_Cl L/Día
    ltv1 <- log(52.4 * 0.001) #Ln_v1 L
    ltQ <-  log(2.26 * 0.001) #Ln_Q  L/Día
    ltv2 <- log(19.6 * 0.001) #Ln_v2 L
    # B) Variabilidad en las observaciones
    add.err  <- 0.15 #error_aditivo
    prop.err <- 0.15 #error_proporcional
    # *2* Efectos aleatorios "omega". Variabilidad entre sujetos "omega"
    # Variabilidad de los parámetros poblacionales (varianza)
    eta.cl  ~ 0.3001 * * 2 #eta_cl 
    eta.v1  ~ 0.1255 * * 2 #eta_v1  
    eta.v2  ~ 0.5165 * * 2 #eta_v2 
  })
  model({
    # First parameters are defined in terms of the initial estimates
    # parameter names.
    cl <- exp(ltcl+LWT -0.313*LW65 -0.855*LALB41+ ATI*log(1.292)+ IMM*log(0.863)  + eta.cl)
    v1 <- exp(ltv1+LWT -0.233*LW65 +eta.v1)
    Q  <- exp(ltQ+LWT)
    v2 <- exp(ltv2 +LWT-0.588*LW65  +eta.v2)
    # After the differential equations are defined
    k <- cl/v1
    k12 <- Q/v1
    k21 <- Q/v2
    d/dt(central)  = -(k+k12)*central+k21*peripheral # masa (mg) de IFX en compart. central 
    d/dt(peripheral) = k12*central-k21*peripheral    #Free Drug second compartment amount
    # And the concentration is then calculated
    conc=central/v1
    conc ~ add(add.err) + prop
(prop.err)
  })
}
nlmixr (mod.dos.cmpt)
  # Calculations required by the model of interest
funcion_transf_logaritmica<-function(nombre_dataframe){
  nombre_dataframe$LWT    <-  log(nombre_dataframe$WT)
  nombre_dataframe$LW65   <-  log(nombre_dataframe$WT/65)
  nombre_dataframe$LALB41 <-  log(nombre_dataframe$ALB/4.1) 
  nombre_dataframe <- data.frame(nombre_dataframe) 
}
  # Read Patient gathered data in Excel from csv file 
datos_pac <-read.csv("pac1_v18_COVARS_ETA.csv", header = TRUE, sep = ";")
datos_pac <- funcion_transf_logaritmica(datos_pac)
  # Estimate EBE's of the patient
fit_pac00 <- nlmixr(mod.dos.cmpt(), datos_pac, est="posthoc")
fit_pac00
fit_pac00$theta
fit_pac00$omega
fit_pac00$omegaR
fit_pac00$eta

The csv file of patient's data it's attached
pac1_v18_COVARS_ETA.csv

Patient's concentrations were chosen to be a perfect fit so...
When I run the model...
I expected a fitting with the etas values once fitted with a result equal to Zero -> and so DV will be equal to PRED and equal to IPRED

Here you can see a screenshot of the result

The questions are?
1- In the fitted object table etas seem to be Zero but if you enter the object$etas a non zero value appears ?
2- Why the fitted IPRED is not equal to the DV?

I don´t know what I`m doing wrong? What I'm missing...?

Your help will be very much appreciated
Thank you very much in advance !

[mattfidler]

Hi @discovererspirit

In the fitted object table etas seem to be Zero but if you enter the object$etas a non zero value appears ?

I would think that the etas in the fit$eta are correct, the eta not showing up in the table is likely a bug.

[mattfidler]

It could also be fixed with nlmixr2. I will check with your code.

[discovererspirit]

Hi Matt,
I'm a beginner and keeping that in mind, I agree that as a general rule, in "real life" patients its "very unlikely" that they present the behaviour of the "typical" patient of the model, and so their IPRED would be different from the PRED.
So I tried to simulate, what I thougth was the easiest patient, a patient, so that patient it's plasma concentrations do not have measurement error and do not deviate from the "typical" patient of the model. This is in an attempt to gain trust in my code and the fitting routine, and that's also the reason why I expected DV to be "almost" equal to PRED and "almost" equal to IPRED.
Is this reasoning at least in theory correct?

[discovererspirit]

I would really appreciate if you could check my code !! and show me the flaws! ;)
Thanks !

[mattfidler]

If the PRED and IPRED are equal it is part of the bug.

The PRED and IPRED and DV values should be different. The PRED should be the furthest from DV and the IPRED should be farthest. However, the PRED isn't optimized for in this case so it may be quite off, while the IPRED should be close to DV.

[mattfidler]

Hi @discovererspirit

Thank you for your patience while I got to this issue.

This is a bug in nlmixr2 posthoc. It currently assumes more than one subject in the data-frame.

As a work-around, if you duplicate the subject and add another ID in your data you will get the right results.

The code change I made was:

datos_pac2 <- datos_pac
datos_pac2$ID <- 2

datos_pac <- rbind(datos_pac, datos_pac2)

Here are the results from my end:

setwd("c:/tmp")
library(nlmixr2)
#> Loading required package: nlmixr2data
mod.dos.cmpt <- function() {
  ini({
    # *1* Efectos fijos "theta" + Error
    # A) Parámetros poblacionales y su valor poblacional
    ltcl <- log(5.42*0.001) #Ln_Cl L/Día
    ltv1 <- log(52.4*0.001) #Ln_v1 L
    ltQ <-  log(2.26*0.001) #Ln_Q  L/Día
    ltv2 <- log(19.6*0.001) #Ln_v2 L
    # B) Variabilidad en las observaciones
    add.err  <- 0.15 #error_aditivo
    prop.err <- 0.15 #error_proporcional
    # *2* Efectos aleatorios "omega". Variabilidad entre sujetos "omega"
    # Variabilidad de los parámetros poblacionales (varianza)
    eta.cl  ~ 0.3001**2 #eta_cl 
    eta.v1  ~ 0.1255**2 #eta_v1  
    eta.v2  ~ 0.5165**2 #eta_v2 
  })
  model({
    # First parameters are defined in terms of the initial estimates
    # parameter names.
    cl <- exp(ltcl+LWT -0.313*LW65 -0.855*LALB41+ ATI*log(1.292)+ IMM*log(0.863)  + eta.cl)
    v1 <- exp(ltv1+LWT -0.233*LW65 +eta.v1)
    Q  <- exp(ltQ+LWT)
    v2 <- exp(ltv2 +LWT-0.588*LW65  +eta.v2)
    # After the differential equations are defined
    k <- cl/v1
    k12 <- Q/v1
    k21 <- Q/v2
    d/dt(central)  = -(k+k12)*central+k21*peripheral # masa (mg) de IFX en compart. central 
    d/dt(peripheral) = k12*central-k21*peripheral    #Free Drug second compartment amount
    # And the concentration is then calculated
    conc=central/v1
    conc ~ add(add.err) + prop(prop.err)
  })
}
nlmixr (mod.dos.cmpt)
#> ℹ parameter labels from comments will be replaced by 'label()'
#>  ── rxode2-based free-form 2-cmt ODE model ────────────────────────────────────── 
#>  ── Initalization: ──  
#> Fixed Effects ($theta): 
#>      ltcl      ltv1       ltQ      ltv2   add.err  prop.err 
#> -5.217659 -2.948849 -6.092390 -3.932226  0.150000  0.150000 
#> 
#> Omega ($omega): 
#>            eta.cl     eta.v1    eta.v2
#> eta.cl 0.09006001 0.00000000 0.0000000
#> eta.v1 0.00000000 0.01575025 0.0000000
#> eta.v2 0.00000000 0.00000000 0.2667722
#> 
#> States ($state or $stateDf): 
#>   Compartment Number Compartment Name
#> 1                  1          central
#> 2                  2       peripheral
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     ini({
#>         ltcl <- -5.21765946353058
#>         label("Ln_Cl L/Día")
#>         ltv1 <- -2.94884868765514
#>         label("Ln_v1 L")
#>         ltQ <- -6.09239046569794
#>         label("Ln_Q  L/Día")
#>         ltv2 <- -3.93222571274567
#>         label("Ln_v2 L")
#>         add.err <- c(0, 0.15)
#>         label("error_aditivo")
#>         prop.err <- c(0, 0.15)
#>         label("error_proporcional")
#>         eta.cl ~ 0.09006001
#>         eta.v1 ~ 0.01575025
#>         eta.v2 ~ 0.26677225
#>     })
#>     model({
#>         cl <- exp(ltcl + LWT - 0.313 * LW65 - 0.855 * LALB41 + 
#>             ATI * log(1.292) + IMM * log(0.863) + eta.cl)
#>         v1 <- exp(ltv1 + LWT - 0.233 * LW65 + eta.v1)
#>         Q <- exp(ltQ + LWT)
#>         v2 <- exp(ltv2 + LWT - 0.588 * LW65 + eta.v2)
#>         k <- cl/v1
#>         k12 <- Q/v1
#>         k21 <- Q/v2
#>         d/dt(central) = -(k + k12) * central + k21 * peripheral
#>         d/dt(peripheral) = k12 * central - k21 * peripheral
#>         conc = central/v1
#>         conc ~ add(add.err) + prop(prop.err)
#>     })
#> }
# Calculations required by the model of interest
funcion_transf_logaritmica<-function(nombre_dataframe){
  nombre_dataframe$LWT    <-  log(nombre_dataframe$WT)
  nombre_dataframe$LW65   <-  log(nombre_dataframe$WT/65)
  nombre_dataframe$LALB41 <-  log(nombre_dataframe$ALB/4.1) 
  nombre_dataframe <- data.frame(nombre_dataframe) 
}
# Read Patient gathered data in Excel from csv file 
datos_pac <-read.csv("pac1_v18_COVARS_ETA.csv", header = TRUE, sep = ";")
datos_pac <- funcion_transf_logaritmica(datos_pac)

# Work-around for single subject data:
datos_pac2 <- datos_pac
datos_pac2$ID <- 2

datos_pac <- rbind(datos_pac, datos_pac2)

# Estimate EBE's of the patient
fit_pac00 <- nlmixr(mod.dos.cmpt(), datos_pac, est="posthoc")
#> rxode2 2.0.13 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 2968
fit_pac00
#> ── nlmixr² FOCE (outer: nlminb) ──
#> 
#>           OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 14.38780 39.73931 34.21596      -10.86965              NA              NA
#> FOCE  14.26093 39.61243 34.08908      -10.80622              NA              NA
#> 
#> ── Time (sec fit_pac00$time): ──
#> 
#>         setup optimize covariance table compress
#> elapsed 0.028        0          0  0.03     0.02
#> 
#> ── Population Parameters (fit_pac00$parFixed or fit_pac00$parFixedDf): ──
#> 
#>           Est. Back-transformed BSV(CV%) Shrink(SD)%
#> ltcl     -5.22            -5.22                     
#> ltv1     -2.95            -2.95                     
#> ltQ      -6.09          0.00226                     
#> ltv2     -3.93            -3.93                     
#> add.err   0.15             0.15                     
#> prop.err  0.15             0.15                     
#> eta.cl                              30.7      100.% 
#> eta.v1                              12.6      100.% 
#> eta.v2                              55.3      100.% 
#>  
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance (fit_pac00$omega) 
#>     or correlation (fit_pac00$omegaR; diagonals=SDs)
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Information about run found (fit_pac00$runInfo):
#>    • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=)) 
#>   Censoring (fit_pac00$censInformation): No censoring
#>   Minimization message (fit_pac00$message):  
#>     Posthoc prediction with provided THETAs 
#> 
#> ── Fit Data (object fit_pac00 is a modified tibble): ──
#> # A tibble: 4 × 31
#>   ID     TIME    DV  PRED         RES     WRES IPRED  IRES  IWRES CPRED     CRES
#>   <fct> <dbl> <dbl> <dbl>       <dbl>    <dbl> <dbl> <dbl>  <dbl> <dbl>    <dbl>
#> 1 1        14  21.4  21.4  0.00000342  3.83e-7  21.1 0.295 0.0933  21.4 -7.84e-5
#> 2 1        42  11.5  11.5 -0.00000199 -2.56e-7  11.2 0.284 0.169   11.5  2.42e-3
#> 3 2        14  21.4  21.4  0.00000342  3.83e-7  21.1 0.295 0.0933  21.4 -7.84e-5
#> # ℹ 1 more row
#> # ℹ 20 more variables: CWRES <dbl>, eta.cl <dbl>, eta.v1 <dbl>, eta.v2 <dbl>,
#> #   central <dbl>, peripheral <dbl>, cl <dbl>, v1 <dbl>, Q <dbl>, v2 <dbl>,
#> #   k <dbl>, k12 <dbl>, k21 <dbl>, tad <dbl>, dosenum <dbl>, LALB41 <dbl>,
#> #   LW65 <dbl>, LWT <dbl>, IMM <dbl>, ATI <dbl>

^{Created on 2023-07-27 with reprex v2.0.2}