mattfidler
commented
1 year ago I believe this is actually an issue with linCmt() if you change this to an ODE then it seems to work fine.
library(nlmixr2)
#> Loading required package: nlmixr2data
m1 <- function() {
ini({
tka <- .5
tcl <- -3.2
tv <- -1
eta.ka ~ 1
eta.cl ~ 2
eta.v ~ 1
add.err <- 0.1
})
model({
ka <- exp(tka + eta.ka)
cl <- exp(tcl + eta.cl)
v <- exp(tv + eta.v)
kel <- cl/v
d/dt(depot) <- -depot*ka
d/dt(central) <- depot*ka - kel*central
Cc <- central/v
Cc ~ add(add.err)
})
}
theo_md_1 <- theo_md
theo_md_1$EVID[theo_md_1$EVID == 101] <- 1
# make the 24 hour dose a steady-state
theo_md_1$II <- 0
theo_md_1$SS <- 0
# make the 24 hour dose be steady state
theo_md_1$II[theo_md_1$TIME == 24] <- 12
theo_md_1$SS[theo_md_1$TIME == 24] <- 1
theo_md_2 <- theo_md_1
theo_md_2$EVID[theo_md_2$TIME == 24] <- 4
theo_md_2$TIME[theo_md_2$TIME >= 24] <- theo_md_2$TIME[theo_md_2$TIME >= 24] - 24
f <- nlmixr2(m1, theo_md_2, "foce", control=foceiControl(print=0))
#> rxode2 2.0.13 using 4 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`
#> Theta reset (ETA drift)
#> calculating covariance matrix
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 20184
#> → compress parHist in nlmixr2 object, save 3392
print(f)
#> ── nlmixr² FOCE (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCE 546.6585 1045.858 1070.89 -515.929 13.27906 4.723096
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.029 7.544 7.544 0.06 0.01 80.793
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka 0.0528 0.131 247 1.05 (0.816, 1.36) 33.7 11.5%
#> tcl 1.19 0.068 5.7 3.3 (2.89, 3.77) 23.3 4.29%
#> tv 3.25 0.0524 1.61 25.8 (23.3, 28.6) 11.8 33.1%
#> add.err 1.56 1.56
#>
#> 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
#> • last objective function was not at minimum, possible problems in optimization
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> • mu-referenced Thetas were reset during optimization; (Can control by foceiControl(resetThetaP=.,resetThetaCheckPer=.,resetThetaFinalP=.))
#> Censoring ($censInformation): No censoring
#> Minimization message ($message):
#> false convergence (8)
#> In an ODE system, false convergence may mean "useless" evaluations were performed.
#> See https://tinyurl.com/yyrrwkce
#> It could also mean the convergence is poor, check results before accepting fit
#> You may also try a good derivative free optimization:
#> nlmixr2(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 264 × 24
#> ID RESETNO TIME DV PRED RES WRES IPRED IRES IWRES CPRED CRES
#> <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1 0 0.74 0 0.74 0.476 0 0.74 0.476 0 0.74
#> 2 1 1 0.25 2.84 2.82 0.0187 0.0105 3.36 -0.519 -0.333 2.78 0.0627
#> 3 1 1 0.57 6.57 5.37 1.20 0.568 6.31 0.258 0.166 5.31 1.26
#> # ℹ 261 more rows
#> # ℹ 12 more variables: CWRES <dbl>, eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>,
#> # depot <dbl>, central <dbl>, ka <dbl>, cl <dbl>, v <dbl>, kel <dbl>,
#> # tad <dbl>, dosenum <dbl>Created on 2023-08-02 with reprex v2.0.2
I am having an issue with EVID and II/SS, reproduced below using
theo_md. If EVID set to 4 and II/SS are also used to denote a new steady-state dose, nlmixr fails to run.To reproduce I have made 3 datasets:
for github
library(nlmixr2)
dataset 1 we put EVID = 1 and assume that the 24 hour dose was given at steady state
theo_md_1 <- theo_md theo_md_1$EVID[theo_md_1$EVID == 101] <- 1
make the 24 hour dose a steady-state
theo_md_1$II <- 0 theo_md_1$SS <- 0
make the 24 hour dose be steady state
theo_md_1$II[theo_md_1$TIME == 24] <- 12 theo_md_1$SS[theo_md_1$TIME == 24] <- 1
now reset the clock at 24 hours change later doses to reset time at 24 hours
theo_md_2 <- theo_md_1 theo_md_2$EVID[theo_md_2$TIME == 24] <- 4 theo_md_2$TIME[theo_md_2$TIME >= 24] <- theo_md_2$TIME[theo_md_2$TIME >= 24] - 24
drop the II and SS to see effect of EVID = 4
theo_md_3 <- theo_md_2 theo_md_3$SS <- NULL theo_md_3$II <- NULL
m1 <- function() { ini({ tka <- .5 tcl <- -3.2 tv <- -1 eta.ka ~ 1 eta.cl ~ 2 eta.v ~ 1 add.err <- 0.1 }) model({ ka <- exp(tka + eta.ka) cl <- exp(tcl + eta.cl) v <- exp(tv + eta.v) linCmt() ~ add(add.err) }) }
fit1 <- nlmixr2(m1, theo_md_1, est = "foce", table = tableControl(cwres=TRUE, npde=TRUE)) fit2 <- nlmixr2(m1, theo_md_2, est = "foce", table = tableControl(cwres=TRUE, npde=TRUE))
Does not run:
Error in .foceiFitInternal(.ret) :
On initial gradient evaluation, one or more parameters have a zero gradient
Change model, try different initial estimates or use outerOpt="bobyqa")
Error : Could not fit data
On initial gradient evaluation, one or more parameters have a zero gradient
Change model, try different initial estimates or use outerOpt="bobyqa")
Error: Could not fit data
On initial gradient evaluation, one or more parameters have a zero gradient
Change model, try different initial estimates or use outerOpt="bobyqa")
fit3 <- nlmixr2(m1, theo_md_3, est = "foce", table = tableControl(cwres=TRUE, npde=TRUE))
different parameter estimates (expected)
fit1$parFixedDf fit3$parFixedDf