nlmixr2 / nlmixr2est

nlmixr2 estimation routines
https://nlmixr2.github.io/nlmixr2est/
GNU General Public License v2.0
7 stars 3 forks source link

Make sure and test that `fit$eta` are non-zero #376

Closed mattfidler closed 1 year ago

mattfidler commented 1 year ago

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.

Originally posted by @mattfidler in https://github.com/nlmixrdevelopment/nlmixr/discussions/643#discussioncomment-6431872

mattfidler commented 1 year ago

Can reproduce with current nlmixr2:

``` r setwd("c:/tmp") library(nlmixr2) #> Warning: package 'nlmixr2' was built under R version 4.3.1 #> 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) # Estimate EBE's of the patient fit_pac00 <- nlmixr(mod.dos.cmpt(), datos_pac, est="posthoc") #> → loading into symengine environment... #> → pruning branches (`if`/`else`) of full model... #> ✔ done #> → calculate jacobian #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> → calculate sensitivities #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> → calculate ∂(f)/∂(η) #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> → finding duplicate expressions in inner model... #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> → optimizing duplicate expressions in inner model... #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> → finding duplicate expressions in EBE model... #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> → optimizing duplicate expressions in EBE model... #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> → compiling inner model... #> using C compiler: 'gcc.exe (GCC) 12.2.0' #> ✔ done #> → finding duplicate expressions in FD model... #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> → optimizing duplicate expressions in FD model... #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> → compiling EBE model... #> using C compiler: 'gcc.exe (GCC) 12.2.0' #> ✔ done #> → compiling events FD model... #> using C compiler: 'gcc.exe (GCC) 12.2.0' #> ✔ done #> rxode2 2.0.13 using 4 threads (see ?getRxThreads) #> no cache: create with `rxCreateCache()` #> → Calculating residuals/tables #> ✔ done #> → compress origData in nlmixr2 object, save 2584 #> → loading into symengine environment... #> → pruning branches (`if`/`else`) of full model... #> ✔ done #> → calculate jacobian #> → calculate sensitivities #> → calculate ∂(f)/∂(η) #> → calculate ∂(R²)/∂(η) #> → finding duplicate expressions in inner model... #> → optimizing duplicate expressions in inner model... #> → finding duplicate expressions in EBE model... #> → optimizing duplicate expressions in EBE model... #> → compiling inner model... #> using C compiler: 'gcc.exe (GCC) 12.2.0' #> ✔ done #> → finding duplicate expressions in FD model... #> → optimizing duplicate expressions in FD model... #> → compiling EBE model... #> using C compiler: 'gcc.exe (GCC) 12.2.0' #> ✔ done #> → compiling events FD model... #> using C compiler: 'gcc.exe (GCC) 12.2.0' #> ✔ done fit_pac00 #> ── nlmixr² FOCE (outer: nlminb) ── #> #> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor) #> FOCEi 7.193900 28.86965 17.10798 -5.434827 NA NA #> FOCE 7.130463 28.80622 17.04454 -5.403109 NA NA #> #> ── Time (sec fit_pac00$time): ── #> #> setup optimize covariance table compress #> elapsed 0.031 0 0 0.03 0 #> #> ── 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 #> eta.v1 12.6 #> eta.v2 55.3 #> #> 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=)) #> • package 'rxode2' was built under R version 4.3.1 #> 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: 2 × 31 #> ID TIME DV PRED RES WRES IPRED IRES IWRES CPRED CRES #> #> 1 1 14 21.4 21.4 0.00000342 3.83e-7 21.1 0.295 0.0933 21.1 0.295 #> 2 1 42 11.5 11.5 -0.00000199 -2.56e-7 11.2 0.284 0.169 11.2 0.284 #> # ℹ 20 more variables: CWRES , eta.cl , eta.v1 , eta.v2 , #> # central , peripheral , cl , v1 , Q , v2 , #> # k , k12 , k21 , tad , dosenum , LALB41 , #> # LW65 , LWT , IMM , ATI fit_pac00$theta #> ltcl ltv1 ltQ ltv2 add.err prop.err #> -5.217659 -2.948849 -6.092390 -3.932226 0.150000 0.150000 fit_pac00$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 fit_pac00$omegaR #> eta.cl eta.v1 eta.v2 #> eta.cl 0.3001 0.0000 0.0000 #> eta.v1 0.0000 0.1255 0.0000 #> eta.v2 0.0000 0.0000 0.5165 fit_pac00$eta #> ID eta.cl eta.v1 eta.v2 #> 1 1 0.01 -0.003878581 -0.0002834194 print(summary(fit_pac00)) #> ID TIME DV PRED RES #> 1:2 Min. :14 Min. :11.46 Min. :11.46 Min. :-1.985e-06 #> 1st Qu.:21 1st Qu.:13.94 1st Qu.:13.94 1st Qu.:-6.344e-07 #> Median :28 Median :16.43 Median :16.43 Median : 7.166e-07 #> Mean :28 Mean :16.43 Mean :16.43 Mean : 7.166e-07 #> 3rd Qu.:35 3rd Qu.:18.92 3rd Qu.:18.92 3rd Qu.: 2.068e-06 #> Max. :42 Max. :21.41 Max. :21.41 Max. : 3.418e-06 #> WRES IPRED IRES IWRES #> Min. :-2.563e-07 Min. :11.17 Min. :0.2836 Min. :0.09334 #> 1st Qu.:-9.643e-08 1st Qu.:13.66 1st Qu.:0.2864 1st Qu.:0.11231 #> Median : 6.346e-08 Median :16.14 Median :0.2892 Median :0.13128 #> Mean : 6.346e-08 Mean :16.14 Mean :0.2892 Mean :0.13128 #> 3rd Qu.: 2.233e-07 3rd Qu.:18.63 3rd Qu.:0.2920 3rd Qu.:0.15024 #> Max. : 3.832e-07 Max. :21.12 Max. :0.2948 Max. :0.16921 #> CPRED CRES CWRES eta.cl eta.v1 #> Min. :11.17 Min. :0.2836 Min. :0.03318 Min. :0 Min. :0 #> 1st Qu.:13.66 1st Qu.:0.2864 1st Qu.:0.03420 1st Qu.:0 1st Qu.:0 #> Median :16.14 Median :0.2892 Median :0.03522 Median :0 Median :0 #> Mean :16.14 Mean :0.2892 Mean :0.03522 Mean :0 Mean :0 #> 3rd Qu.:18.63 3rd Qu.:0.2920 3rd Qu.:0.03624 3rd Qu.:0 3rd Qu.:0 #> Max. :21.12 Max. :0.2948 Max. :0.03726 Max. :0 Max. :0 #> eta.v2 central peripheral cl v1 #> Min. :0 Min. :37.90 Min. :28.29 Min. :0.3558 Min. :3.393 #> 1st Qu.:0 1st Qu.:46.34 1st Qu.:31.14 1st Qu.:0.3558 1st Qu.:3.393 #> Median :0 Median :54.77 Median :33.98 Median :0.3558 Median :3.393 #> Mean :0 Mean :54.77 Mean :33.98 Mean :0.3558 Mean :3.393 #> 3rd Qu.:0 3rd Qu.:63.21 3rd Qu.:36.83 3rd Qu.:0.3558 3rd Qu.:3.393 #> Max. :0 Max. :71.64 Max. :39.67 Max. :0.3558 Max. :3.393 #> Q v2 k k12 #> Min. :0.1469 Min. :1.274 Min. :0.1049 Min. :0.0433 #> 1st Qu.:0.1469 1st Qu.:1.274 1st Qu.:0.1049 1st Qu.:0.0433 #> Median :0.1469 Median :1.274 Median :0.1049 Median :0.0433 #> Mean :0.1469 Mean :1.274 Mean :0.1049 Mean :0.0433 #> 3rd Qu.:0.1469 3rd Qu.:1.274 3rd Qu.:0.1049 3rd Qu.:0.0433 #> Max. :0.1469 Max. :1.274 Max. :0.1049 Max. :0.0433 #> k21 tad dosenum LALB41 LW65 #> Min. :0.1153 Min. : 0 Min. :2 Min. :0 Min. :0 #> 1st Qu.:0.1153 1st Qu.: 7 1st Qu.:2 1st Qu.:0 1st Qu.:0 #> Median :0.1153 Median :14 Median :2 Median :0 Median :0 #> Mean :0.1153 Mean :14 Mean :2 Mean :0 Mean :0 #> 3rd Qu.:0.1153 3rd Qu.:21 3rd Qu.:2 3rd Qu.:0 3rd Qu.:0 #> Max. :0.1153 Max. :28 Max. :2 Max. :0 Max. :0 #> LWT IMM ATI #> Min. :4.174 Min. :0 Min. :0 #> 1st Qu.:4.174 1st Qu.:0 1st Qu.:0 #> Median :4.174 Median :0 Median :0 #> Mean :4.174 Mean :0 Mean :0 #> 3rd Qu.:4.174 3rd Qu.:0 3rd Qu.:0 #> Max. :4.174 Max. :0 Max. :0 ``` Created on 2023-07-12 with [reprex v2.0.2](https://reprex.tidyverse.org)