How to use logit error in rxode2/nlmixr2

[mattfidler]

I have a very simple pd model with

model({

    …

    pd ~ prop(prop.err) + add(add.sd)

  })

Which works okay but gives values just above and below the natural range of pd (0, 10).

Obviously some sort of logit would be helpful. And normally I’d hand-code the solution but it feels like I’d be duplicating the effort put into logitNorm() etc. I can’t find any code examples in the package or online. Can you point me to something?

[mattfidler]

library(nlmixr2)
#> Loading required package: nlmixr2data

f <- function() {
  ini({
    Emax <- 10
    E0 <- 0.001
    tEC50 <- log(10)
    h <- 1
    eta.ec50 ~ 0.1
    eff.sd <- 0.5
  })
  model({
    EC50 <- exp(tEC50+eta.ec50)
    Cp <- time
    eff  <- E0 + (Emax-E0) * Cp^h / (EC50^h + Cp^h)
    eff ~ logitNorm(eff.sd, 0, 10)
  })
}

et <- et(seq(0, 20, by=0.1)) %>%
  et(id=1:10)

s <- rxSolve(f, et)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’

plot(s, eff, sim)

s <- s %>%
  dplyr::rename(dv=sim) %>%
  dplyr::select(id, time, dv)


fit <- nlmixr2(f, s, "focei", control=list(print=0))
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → calculate jacobian
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 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 (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.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 (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> → compiling events FD model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> rxode2 3.0.2.9000 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 26368
#> → compress parHistData in nlmixr2 object, save 5888

print(fit)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#> 
#>           OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 1601.564 5307.697 5341.332      -2647.849        35171619        43.98629
#> 
#> ── Time (sec $time): ──
#> 
#>            setup optimize covariance table compress    other
#> elapsed 0.034237 0.301738   0.301741  0.21    0.018 3.210284
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#> 
#>            Est.      SE %RSE      Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> Emax         10   0.376 3.76              10 (9.27, 10.7)                     
#> E0     0.000974 9.1e-05 9.35 0.000974 (0.000795, 0.00115)                     
#> tEC50      2.47   0.129 5.24            11.8 (9.17, 15.2)     34.3     -5.18% 
#> h          1.02  0.0214 2.11           1.02 (0.973, 1.06)                     
#> eff.sd    0.506                                     0.506                     
#>  
#>   Covariance Type ($covMethod): |r|,|s|
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Information about run found ($runInfo):
#>    • gradient problems with initial estimate and covariance; see $scaleInfo 
#>    • since sandwich matrix is corrected, you may compare to $covR or $covS if you wish 
#>    • S matrix non-positive definite but corrected by S = sqrtm(S%*%S) 
#>    • R matrix non-positive definite but corrected by R = sqrtm(R%*%R) 
#>    • bad solve during optimization 
#>    • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.)) 
#>    • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=)) 
#>    • Hessian reset during optimization; (can control by foceiControl(resetHessianAndEta=.)) 
#>   Censoring ($censInformation): No censoring
#>   Minimization message ($message):  
#>     relative convergence (4) 
#> 
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 2,010 × 15
#>   ID     TIME       DV     PRED       RES   WRES   IPRED     IRES  IWRES   CPRED
#>   <fct> <dbl>    <dbl>    <dbl>     <dbl>  <dbl>   <dbl>    <dbl>  <dbl>   <dbl>
#> 1 1       0   0.000650 0.000974 -0.000324 -0.800 9.74e-4 -3.24e-4 -0.800 9.74e-4
#> 2 1       0.1 0.0681   0.0792   -0.0111   -0.252 9.89e-2 -3.08e-2 -0.745 7.92e-2
#> 3 1       0.2 0.250    0.158     0.0922    0.773 1.97e-1  5.31e-2  0.483 1.58e-1
#> # ℹ 2,007 more rows
#> # ℹ 5 more variables: CRES <dbl>, CWRES <dbl>, eta.ec50 <dbl>, EC50 <dbl>,
#> #   Cp <dbl>

^{Created on 2024-12-03 with reprex v2.1.1}