Replies: 3 comments 7 replies
-
|
Hi @chansgary I would suggest using The following example shows the exploration I did with your example: library(nlmixr2est)
#> Loading required package: nlmixr2data
sample.data = read.csv("~/Downloads/sample.data.csv")[,-1]
mod.2CMT <- function() {
ini({
tVc <- c(-2.30, 2.68, 3.91) # log_V1 (L)
tCl <- c(0, 2.21, 4.60) # log_Cl (L/hr)
tVp <- c(-2.30, 1.53, 3.00) # log_V2 (L)
tQ <- c(-2.30, 0.46, 2.30) # log_Q (L/hr)
eta.Vc + eta.Cl ~ c(0.373973,
0.333469, 0.530645) # IIV
Aeff.Cl <- c(0.1,1.178,5) # Effect of A on CL
Beff.Vc <- c(0.1,1.56,2) # Effect of B on Vc
prop.err = c(0,0.133,1)
})
model({
Vc <- exp(tVc + eta.Vc + log(WT70)) * (1 + Beff.Vc * B)
Vp <- exp(tVp + eta.Vc + log(WT70))
Cl <- exp(tCl + eta.Cl + 0.75 * log(WT70) + Aeff.Cl * A)
Q = exp(tQ + 0.75 * log(WT70))
linCmt() ~ prop(prop.err)
})
}
nlmixr2(mod.2CMT)
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ── rxode2-based solved PK 2-compartment model ──────────────────────────────────
#> ── Initalization: ──
#> Fixed Effects ($theta):
#> tVc tCl tVp tQ Aeff.Cl Beff.Vc prop.err
#> 2.680 2.210 1.530 0.460 1.178 1.560 0.133
#>
#> Omega ($omega):
#> eta.Vc eta.Cl
#> eta.Vc 0.373973 0.333469
#> eta.Cl 0.333469 0.530645
#> ── μ-referencing ($muRefTable): ──
#> theta eta level covariates
#> 1 tCl eta.Cl id A*Aeff.Cl
#>
#> ── Model (Normalized Syntax): ──
#> function() {
#> ini({
#> tVc <- c(-2.3, 2.68, 3.91)
#> label("log_V1 (L)")
#> tCl <- c(0, 2.21, 4.6)
#> label("log_Cl (L/hr)")
#> tVp <- c(-2.3, 1.53, 3)
#> label("log_V2 (L)")
#> tQ <- c(-2.3, 0.46, 2.3)
#> label("log_Q (L/hr)")
#> Aeff.Cl <- c(0.1, 1.178, 5)
#> label("Effect of A on CL")
#> Beff.Vc <- c(0.1, 1.56, 2)
#> label("Effect of B on Vc")
#> prop.err <- c(0, 0.133, 1)
#> eta.Vc + eta.Cl ~ c(0.373973, 0.333469, 0.530645)
#> })
#> model({
#> Vc <- exp(tVc + eta.Vc + log(WT70)) * (1 + Beff.Vc *
#> B)
#> Vp <- exp(tVp + eta.Vc + log(WT70))
#> Cl <- exp(tCl + eta.Cl + 0.75 * log(WT70) + Aeff.Cl *
#> A)
#> Q = exp(tQ + 0.75 * log(WT70))
#> linCmt() ~ prop(prop.err)
#> })
#> }
fit.example = try(nlmixr2(mod.2CMT,sample.data,"nlme"))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of nlme model...
#> ✔ done
#> → finding duplicate expressions in nlme model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in nlme model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#>
#> **Iteration 1
#> LME step: Loglik: -79.20682, nlminb iterations: 45
#> reStruct parameters:
#> ID1 ID2 ID3
#> 4.495867 -5.518972 8.059559
#> Beginning PNLS step: .. Error in nlme.formula(model = DV ~ .nlmixrNlmeFun(pars = list(tCl = tCl + :
#> Singularity in backsolve at level 0, block 1
#> Error : Singularity in backsolve at level 0, block 1
mod2 <- mod.2CMT |>
model(Vc <- exp(tVc + log(WT70)) * (1 + Beff.Vc * B)) |>
model(Vp <- exp(tVp + log(WT70))) |>
model(Cl <- exp(tCl + 0.75 * log(WT70) + Aeff.Cl * A))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ! remove between subject variability `eta.Vc`
#> ! remove between subject variability `eta.Cl`
fit.example <- nlmixr2(mod2, sample.data, "nlm", nlmControl(print=0))
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of population log-likelihood model...
#> ✔ done
#> → calculate jacobian
#> → calculate ∂(f)/∂(θ)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04
#> → finding duplicate expressions in nlm llik gradient...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in nlm llik gradient...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in nlm pred-only...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in nlm pred-only...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> rxode2 2.1.2 using 8 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> → calculating covariance
#> ✔ done
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> ✔ done
#> → loading llik model into symengine environment...
#> → pruning branches (`if`/`else`) of llik full model...
#> ✔ done
#> → finding duplicate expressions in Llik EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in Llik EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling Llik EBE model...
#> ✔ done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 4432
#> → compress parHistData in nlmixr2 object, save 9448
mod3 <- mod2 |>
model(Vc <- exp(tVc + eta.Vc + log(WT70)) * (1 + Beff.Vc * B)) |>
model(Vp <- exp(tVp + eta.Vc + log(WT70))) |>
model(Cl <- exp(tCl + eta.Cl + 0.75 * log(WT70) + Aeff.Cl * A)) |>
ini(eta.Vc + eta.Cl ~ c(0.373973,
0.333469, 0.530645))
#> ℹ add between subject variability `eta.Vc` and set estimate to 1
#> ℹ add between subject variability `eta.Cl` and set estimate to 1
#> ℹ change initial estimate of `eta.Vc` to `0.373973`
#> ℹ add covariance between `eta.Cl` and `eta.Vc` with initial estimate `0.333469`
#> ℹ change initial estimate of `eta.Cl` to `0.530645`
fit.example = try(nlmixr2(mod3, sample.data, "nlme"))
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of nlme model...
#> ✔ done
#> → finding duplicate expressions in nlme model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in nlme model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#>
#> **Iteration 1
#> LME step: Loglik: -79.20682, nlminb iterations: 45
#> reStruct parameters:
#> ID1 ID2 ID3
#> 4.495867 -5.518972 8.059559
#> Beginning PNLS step: .. Error in nlme.formula(model = DV ~ .nlmixrNlmeFun(pars = list(tCl = tCl + :
#> Singularity in backsolve at level 0, block 1
#> Error : Singularity in backsolve at level 0, block 1
fit.exampleF <- nlmixr2(mod.2CMT, sample.data, "focei", foceiControl(print=0))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:08
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 4432
#> → compress parHistData in nlmixr2 object, save 7960
print(fit.exampleF)
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 1498.389 1593.742 1610.878 -786.871 5.733959 1.112537
#>
#> ── Time (sec $time): ──
#>
#> setup optimize covariance table compress
#> elapsed 0.000926 8.136331 8.136333 0.031 0.03
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tVc log_V1 (L) 2.7 0.0029 0.107 14.9 (14.8, 15)
#> tCl log_Cl (L/hr) 2.23 0.00304 0.137 9.29 (9.24, 9.35) 83.7
#> tVp log_V2 (L) 1.78 0.00227 0.127 5.92 (5.9, 5.95)
#> tQ log_Q (L/hr) 0.46 0.00127 0.277 1.58 (1.58, 1.59)
#> Aeff.Cl Effect of A on CL 1.18 0.00176 0.15 1.18 (1.17, 1.18)
#> Beff.Vc Effect of B on Vc 1.56 0.00195 0.125 1.56 (1.56, 1.56)
#> prop.err 0.133 0.133
#> eta.Vc 67.3
#> Shrink(SD)%
#> tVc
#> tCl 24.2%
#> tVp
#> tQ
#> Aeff.Cl
#> Beff.Vc
#> prop.err
#> eta.Vc 92.6%
#>
#> Covariance Type ($covMethod): r
#> Correlations in between subject variability (BSV) matrix:
#> cor:eta.Cl,eta.Vc
#> 0.749
#>
#>
#> 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
#> • using R matrix to calculate covariance
#> • S matrix non-positive definite
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> 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: 41 × 23
#> ID TIME DV PRED RES WRES IPRED IRES IWRES CPRED CRES CWRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 101 16.3 0.436 208. -207. -1.25 16.1 -15.7 -7.32 56.6 -56.2 -4.71
#> 2 101 16.7 42.3 225. -183. -1.06 25.9 16.4 4.75 79.3 -37.0 -2.06
#> 3 101 17.3 22.5 220. -197. -1.17 17.1 5.40 2.37 64.3 -41.8 -3.17
#> # ℹ 38 more rows
#> # ℹ 11 more variables: eta.Vc <dbl>, eta.Cl <dbl>, Vc <dbl>, Vp <dbl>,
#> # Cl <dbl>, Q <dbl>, tad <dbl>, dosenum <dbl>, A <dbl>, B <dbl>, WT70 <dbl>Created on 2024-02-04 with reprex v2.1.0 However, I find the assumption curious that you have exactly the same between subject varaibility between central volume and peripheral volume curious. I have not seen this very often is something to think a bit more about if you want to find a solution in If you wish to use library(nlmixr2est)
#> Loading required package: nlmixr2data
sample.data = read.csv("~/Downloads/sample.data.csv")[,-1]
mod.2CMT <- function() {
ini({
tVc <- c(-2.30, 2.68, 3.91) # log_V1 (L)
tCl <- c(0, 2.21, 4.60) # log_Cl (L/hr)
tVp <- c(-2.30, 1.53, 3.00) # log_V2 (L)
tQ <- c(-2.30, 0.46, 2.30) # log_Q (L/hr)
eta.Vc + eta.Cl ~ c(0.373973,
0.333469, 0.530645) # IIV
Aeff.Cl <- c(0.1,1.178,5) # Effect of A on CL
Beff.Vc <- c(0.1,1.56,2) # Effect of B on Vc
prop.err = c(0,0.133,1)
})
model({
Vc <- exp(tVc + eta.Vc + log(WT70)) * (1 + Beff.Vc * B)
Vp <- exp(tVp + log(WT70))
Cl <- exp(tCl + eta.Cl + 0.75 * log(WT70) + Aeff.Cl * A)
Q = exp(tQ + 0.75 * log(WT70))
linCmt() ~ prop(prop.err)
})
}
nlmixr2(mod.2CMT)
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ── rxode2-based solved PK 2-compartment model ──────────────────────────────────
#> ── Initalization: ──
#> Fixed Effects ($theta):
#> tVc tCl tVp tQ Aeff.Cl Beff.Vc prop.err
#> 2.680 2.210 1.530 0.460 1.178 1.560 0.133
#>
#> Omega ($omega):
#> eta.Vc eta.Cl
#> eta.Vc 0.373973 0.333469
#> eta.Cl 0.333469 0.530645
#> ── μ-referencing ($muRefTable): ──
#> theta eta level covariates
#> 1 tVc eta.Vc id
#> 2 tCl eta.Cl id A*Aeff.Cl
#>
#> ── Model (Normalized Syntax): ──
#> function() {
#> ini({
#> tVc <- c(-2.3, 2.68, 3.91)
#> label("log_V1 (L)")
#> tCl <- c(0, 2.21, 4.6)
#> label("log_Cl (L/hr)")
#> tVp <- c(-2.3, 1.53, 3)
#> label("log_V2 (L)")
#> tQ <- c(-2.3, 0.46, 2.3)
#> label("log_Q (L/hr)")
#> Aeff.Cl <- c(0.1, 1.178, 5)
#> label("Effect of A on CL")
#> Beff.Vc <- c(0.1, 1.56, 2)
#> label("Effect of B on Vc")
#> prop.err <- c(0, 0.133, 1)
#> eta.Vc + eta.Cl ~ c(0.373973, 0.333469, 0.530645)
#> })
#> model({
#> Vc <- exp(tVc + eta.Vc + log(WT70)) * (1 + Beff.Vc *
#> B)
#> Vp <- exp(tVp + log(WT70))
#> Cl <- exp(tCl + eta.Cl + 0.75 * log(WT70) + Aeff.Cl *
#> A)
#> Q = exp(tQ + 0.75 * log(WT70))
#> linCmt() ~ prop(prop.err)
#> })
#> }
fit.example = nlmixr2(mod.2CMT,sample.data,"saem", saemControl(print=0))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> rxode2 2.1.2 using 8 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`
#> NaN in prediction; Consider: relax atol & rtol; change initials; change seed; change structural model
#> warning only issued once per problem
#> Calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem predOnly model 0...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in saem predOnly model 1...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem predOnly model 1...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in saem predOnly model 2...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem predOnly model 2...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 4432
#> → compress phiM in nlmixr2 object, save 152256
#> → compress parHistData in nlmixr2 object, save 22064
#> → compress saem0 in nlmixr2 object, save 30560
print(fit.example)
#> ── nlmixr² SAEM OBJF by FOCEi approximation ──
#>
#> Gaussian/Laplacian Likelihoods: AIC() or $objf etc.
#> FOCEi CWRES & Likelihoods: addCwres()
#>
#> ── Time (sec $time): ──
#>
#> setup covariance saem table compress other
#> elapsed 0.00137 0.011007 8.335 0.075 0.026 1.423623
#>
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#>
#> Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%)
#> tVc log_V1 (L) 2.57 0.6 23.4 13 (4.02, 42.3) 8.75
#> tCl log_Cl (L/hr) 2.41 0.134 5.56 11.1 (8.57, 14.5) 13.4
#> tVp log_V2 (L) 1.12 0.515 45.9 3.07 (1.12, 8.42)
#> tQ log_Q (L/hr) 0.875 0.191 21.8 2.4 (1.65, 3.49)
#> Aeff.Cl Effect of A on CL 2.18 0.208 9.56 2.18 (1.77, 2.58)
#> Beff.Vc Effect of B on Vc 0.79 0.836 106 0.79 (-0.849, 2.43)
#> prop.err 0.295 0.295
#> Shrink(SD)%
#> tVc -21.9%
#> tCl -21.8%
#> tVp
#> tQ
#> Aeff.Cl
#> Beff.Vc
#> prop.err
#>
#> Covariance Type ($covMethod): linFim
#> Correlations in between subject variability (BSV) matrix:
#> cor:eta.Cl,eta.Vc
#> -1.00
#>
#>
#> Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs)
#> Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink
#> Information about run found ($runInfo):
#> • log-scale mu referenced time varying covariates (Aeff.Cl) may have better results on no log-transformed scale (https://github.com/nlmixr2/nlmixr2est/issues/348), check results for plausibility
#> Censoring ($censInformation): No censoring
#>
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 41 × 19
#> ID TIME DV PRED RES IPRED IRES IWRES eta.Vc eta.Cl Vc Vp
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 101 16.3 0.436 114. -114. 103. -102. -3.37 -0.0684 0.105 5.12 0.438
#> 2 101 16.7 42.3 144. -102. 134. -91.6 -2.32 -0.0684 0.105 5.12 0.438
#> 3 101 17.3 22.5 129. -106. 116. -93.4 -2.73 -0.0684 0.105 5.12 0.438
#> # ℹ 38 more rows
#> # ℹ 7 more variables: Cl <dbl>, Q <dbl>, tad <dbl>, dosenum <dbl>, A <dbl>,
#> # B <dbl>, WT70 <dbl>Created on 2024-02-04 with reprex v2.1.0 |
Beta Was this translation helpful? Give feedback.
-
|
This is more of a question than a bug so I am recategorizing this. |
Beta Was this translation helpful? Give feedback.
-
|
@chansgary I am sure that the other folks are busy; If you reach out to my email address maybe I can setup a MS Teams call to help you with your setup issue. (matt dot fidler at novartis dot com) |
Beta Was this translation helpful? Give feedback.
Uh oh!
There was an error while loading. Please reload this page.
-
Hi Dr. Fidler and nlmixr2 users,
I have trouble fitting the following model with the data attached here. Previously, I had no trouble fitting using the same model and the same data. I came up with this error recently with even the simplest model(without any covariate).
This is crucial in one of our research questions in building PK Models. I wonder if you could help me look into this and your help will be greatly appreciated.
`library(dplyr)
sample.data.csv
sample.data = read.csv("sample.data.csv")[,-1]
mod.2CMT <- function() {
ini({
tVc <- c(-2.30, 2.68, 3.91) # log_V1 (L)
tCl <- c(0, 2.21, 4.60) # log_Cl (L/hr)
tVp <- c(-2.30, 1.53, 3.00) # log_V2 (L)
tQ <- c(-2.30, 0.46, 2.30) # log_Q (L/hr)
eta.Vc + eta.Cl ~ c(0.373973,
0.333469, 0.530645) # IIV
Aeff.Cl <- c(0.1,1.178,5) # Effect of A on CL
Beff.Vc <- c(0.1,1.56,2) # Effect of B on Vc
prop.err = c(0,0.133,1)
})
model({
Vc <- exp(tVc + eta.Vc + log(WT70)) * (1 + Beff.Vc * B)
Vp <- exp(tVp + eta.Vc + log(WT70))
Cl <- exp(tCl + eta.Cl + 0.75 * log(WT70) + Aeff.Cl * A)
Q = exp(tQ + 0.75 * log(WT70))
linCmt() ~ prop(prop.err)
})
}
nlmixr2(mod.2CMT)
fit.example = nlmixr2(mod.2CMT,sample.data,"nlme")`
The following is the error it prompts:
`ℹ parameter labels from comments will be replaced by 'label()'
→ loading into symengine environment...
→ pruning branches (
if/else) of nlme model...✔ done
→ finding duplicate expressions in nlme model...
→ optimizing duplicate expressions in nlme model...
✔ done
ld: warning: -single_module is obsolete
ld: warning: -multiply_defined is obsolete
Error in MEEM(object, conLin, control$niterEM) :
Singularity in backsolve at level 0, block 1
Error: Singularity in backsolve at level 0, block 1`
Beta Was this translation helpful? Give feedback.
All reactions