Open billdenney opened 3 years ago
mattfidler
commented
3 years ago You could use concentrtaion=time.
I understand your point, but it is low priority. I think this is likely a RxODE specific issue, though it could go both ways. This may be a large rewrite for a minor work-around.
billdenney
commented
3 years ago While trying this with my real model, I'm getting an error:
Compiling model...done
unhandled error message: EE:[lsoda] 70000 steps taken before reaching tout
@(lsoda.c:748)e.rds
Error: could not solve the systemI'm trying to generate a reprex, but it's not a simple modification of the above code for some reason. Any suggestions on how to work around this lsoda error?
Edit: to be clear, the real model does use d/dt(EFFECT) = 0 as the only differential equation. So, I'm surprised that any lsoda error occurs since a derivative that always evaluates to zero should be trivial to solve.
billdenney
commented
3 years ago And, I have a thought on the extra code comment above. How is linCmt() handled? Would it be possible to setup a similar escape-from-lsoda hatch for models that are not differential equations?
mattfidler
commented
3 years ago There is already an escape hatch; But it needs to be linked against time, just like linCmt()
mattfidler
commented
3 years ago billdenney
commented
3 years ago My reading of https://nlmixrdevelopment.github.io/RxODE/articles/RxODE-model-types.html#prediction-only-models suggests that I would just need to use ipre and that would shortcut all the compartment selection, etc. A quick try of that (with a few days-old version of nlmixr/RxODE) didn't work due to something that appears to be identical to #514. Is there something that I'm missing?
mattfidler
commented
3 years ago You are right ipre is not needed at all. From that page:
Prediction only models are simple to create. You use the RxODE syntax without any ODE systems in them.
What is missing is:
Currently the independent variable must be time.
RxODE has distinct steps in solving a problem:
When there isn't any ODEs present in the system, step 2 and 3 are skipped. It is a simple test of the number of states that are in the problem. If there are no additional parameters to be calculated (which is rare), then step 4 is skipped. In nlmixr parameter estimation step 5 is skipped until the very end of the problem and step 1 is only done once at the very beginning of the problem.
You probably reproduced a system where the independent variable again is not time.
diabloyg
commented
3 years ago Hi, your discussion seemed difficult for me, but according to Matt's solution a slight modification would be sufficient, though there are still issues:
library(nlmixr)
d_sim <-
data.frame(
concentration=1:1000,
CMT="EFFECT"
)
# Simulate
e0 <- 100
emax <- 100
ec50 <- 200
d_sim$effect_pred <-
e0 + emax*d_sim$concentration/(ec50 + d_sim$concentration)
# Add some exponential and additive error
d_sim$dv <-
d_sim$effect_pred * exp(rnorm(n=nrow(d_sim), sd=0.2)) + rnorm(n=nrow(d_sim))
d_sim_time <- d_sim
d_sim_time$time <- d_sim_time$concentration
test_model <- function() {
ini({
e0 <- 100
emax <- 100
ec50 <- log(200)
prop_err <- 0.2
add_err <- 1
})
model({
d/dt(EFFECT) = 1
effect <- e0 + emax*EFFECT/(exp(ec50) + EFFECT)
effect ~ prop(prop_err) + add(add_err)
})
}Created on 2021-05-07 by the reprex package (v2.0.0)
And the results are as follows:
calculating covariance matrix [====|====|====|====|====|====|====|====|====|====] 0:00:00 done Calculating residuals/tables done Warning message: In (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE, : gradient problems with initial estimate and covariance; see $scaleInfo
fit ...... Minimization message (fit$message):
relative convergence (4)
billdenney
commented
3 years ago @diabloyg , Thanks for updating the example. While there are convergence issues, I'm less concerned about those (since the data would not be good for fitting, I would expect convergence problems). FYI, wrapping your example code with the reprex function in the reprex package can make it easier to understand duplicate what you're showing. For the example below, I put everything inside of the reprex function, and after running, it was on my clipboard so I could paste it into this comment
reprex::reprex({
# All my code
})The challenge is that I want to omit the ODE steps 2 and 3 from Matt's list above (https://github.com/nlmixrdevelopment/nlmixr/issues/517#issuecomment-831605619). I believe that the line d/dt(EFFECT) = 1, requires the ODE system to be setup and solved (even though it is not used), and that is what I want to avoid. I tried using cmt(EFFECT) to overcome that, but it doesn't seem to work.
library(nlmixr)
d_sim <-
data.frame(
concentration=1:1000,
CMT="EFFECT"
)
e0 <- 100
emax <- 100
ec50 <- 200
d_sim$effect_pred <-
e0 + emax*d_sim$concentration/(ec50 + d_sim$concentration)
# Add some exponential and additive error
d_sim$dv <-
d_sim$effect_pred * exp(rnorm(n=nrow(d_sim), sd=0.2)) + rnorm(n=nrow(d_sim))
d_sim_time <- d_sim
d_sim_time$time <- d_sim_time$concentration
test_model <- function() {
ini({
e0 <- 100
emax <- 100
ec50 <- log(200)
prop_err <- 0.2
add_err <- 1
})
model({
d/dt(EFFECT) = 1
effect <- e0 + emax*EFFECT/(exp(ec50) + EFFECT)
effect ~ prop(prop_err) + add(add_err)
})
}
test_model_no_ode <- function() {
ini({
e0 <- 100
emax <- 100
ec50 <- log(200)
prop_err <- 0.2
add_err <- 1
})
model({
cmt(EFFECT)
effect <- e0 + emax*EFFECT/(exp(ec50) + EFFECT)
effect ~ prop(prop_err) + add(add_err)
})
}
fit <- nlmixr(object=test_model, data=d_sim_time, est="focei")
#> RxODE 1.0.8.1 using 4 threads (see ?getRxThreads)
#> Key: U: Unscaled Parameters; X: Back-transformed parameters; G: Gill difference gradient approximation
#> F: Forward difference gradient approximation
#> C: Central difference gradient approximation
#> M: Mixed forward and central difference gradient approximation
#> Unscaled parameters for Omegas=chol(solve(omega));
#> Diagonals are transformed, as specified by foceiControl(diagXform=)
#> calculating covariance matrix
#> done
#> Calculating residuals/tables
#> done
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : gradient problems with initial estimate and covariance; see $scaleInfo
fit <- nlmixr(object=test_model_no_ode, data=d_sim_time, est="focei")
#> Error in nlmixrUIModel(.model, ini, fun): 'names' attribute [22] must be the same length as the vector [21]
#> bad parsed model (cached C:\Users\BILLDE~1\AppData\Local\R\cache\R\RxODE/ui-2bdfd43a36269703008f65cef97cd6f4.bad)Created on 2021-05-06 by the reprex package (v2.0.0)
(Note that the "bad parsed model" error was a names error like #514 when run outside of reprex().)
diabloyg
commented
3 years ago @billdenney , Thanks for your kind notice. From my understanding the goal was to find a way avoiding ODE solving process. Then I thought a solved-system approach would be suitable, as stated (I am not sure whether such approach is ODE-free) in https://nlmixrdevelopment.github.io/nlmixr/reference/foceiFit.html
The model and results are as follows:
library(nlmixr)
d_sim <-
data.frame(
concentration=1:1000,
CMT="EFFECT"
)
# Simulate
e0 <- 100
emax <- 100
ec50 <- 200
d_sim$effect_pred <-
e0 + emax*d_sim$concentration/(ec50 + d_sim$concentration)
# Add some exponential error
d_sim$dv <-
d_sim$effect_pred * exp(rnorm(n=nrow(d_sim), sd=0.2))
d_sim_time <- d_sim
d_sim_time$time <- d_sim_time$concentration
d_sim_time$ID <- 1
#Solved system approach
test_model = function ()
{
e0 = theta[1];
emax = theta[2];
ec50 = theta[3];
}
mod <- RxODE({
ipre = e0 + emax*time/(ec50+time)
})
pred <- function() ipre
errProp <- function(){
return(prop(0.2))
}
inits <- list(THTA=c(120,120,220));
fit <- foceiFit(d_sim_time,inits,test_model,mod,pred,errProp)
#> RxODE 1.0.8 using 4 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`
#> Key: U: Unscaled Parameters; X: Back-transformed parameters; G: Gill difference gradient approximation
#> F: Forward difference gradient approximation
#> C: Central difference gradient approximation
#> M: Mixed forward and central difference gradient approximation
#> Unscaled parameters for Omegas=chol(solve(omega));
#> Diagonals are transformed, as specified by foceiControl(diagXform=)
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> | #| Objective Fun | THETA[1] | THETA[2] | THETA[3] | ERR[1] |
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> | 1| 9179.5499 | 0.08906 | 0.08906 | 1.000 | -1.000 |
#> | U| 9179.5499 | 120.0 | 120.0 | 220.0 | 0.4472 |
#> | X| 9179.5499 | 120.0 | 120.0 | 220.0 | 0.4472 |
#> | G| Gill Diff. | 1815. | 1078. | -354.4 | 1477. |
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> | 2| 47894.769 | -0.6088 | -0.3255 | 1.136 | -1.568 |
#> | U| 47894.769 | 36.26 | 70.25 | 250. | 0.1932 |
#> | X| 47894.769 | 36.26 | 70.25 | 250. | 0.1932 |
#> | 3| 8922.6598 | 0.01927 | 0.04760 | 1.014 | -1.057 |
#> | U| 8922.6598 | 111.6 | 115.0 | 223. | 0.4218 |
#> | X| 8922.6598 | 111.6 | 115.0 | 223. | 0.4218 |
#> | F| Forward Diff. | 1651. | 978.6 | -305.4 | 1597. |
#> |-----+---------------+-----------+-----------+-----------+-----------|
#>.................................................................................................................................................
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> | 25| 8092.2060 | -0.07242 | -0.04317 | 0.9149 | -1.533 |
#> | U| 8092.206 | 100.6 | 104.1 | 201.3 | 0.2089 |
#> | X| 8092.206 | 100.6 | 104.1 | 201.3 | 0.2089 |
#> | 26| 8092.2023 | -0.06837 | -0.04417 | 0.9361 | -1.533 |
#> | U| 8092.2023 | 101.1 | 104.0 | 205.9 | 0.2088 |
#> | X| 8092.2023 | 101.1 | 104.0 | 205.9 | 0.2088 |
#> calculating covariance matrix
#> done
#> Calculating residuals/tables
#> done
#> Warning in .collectWarnings(do.call(foceiFit.data.frame0, call, envir =
#> parent.frame(1))): gradient problems with initial estimate and covariance; see
#> $scaleInfo
fit
#>-- nlmixr Population Only (outer: nlminb) fit -----------------------------------------------------------------------
#>
#> OBJF AIC BIC Log-likelihood Condition Number
#>Pop 8069.386 9915.263 9934.894 -4953.631 3177.729
#>
#>-- Time (sec fit$time): ---------------------------------------------------------------------------------------------
#>
#> setup optimize covariance table
#>elapsed 0.003 0.005 0.005 0.06
#>
#>-- (fit$parFixed or fit$parFixedDf): -------------------------------------------------------------------------------
#>
#> Est. SE %RSE Back-transformed(95%CI)
#>THETA[1] 106 6.63 6.23 106 (93.4, 119)
#>THETA[2] 101 6.55 6.48 101 (88.3, 114)
#>THETA[3] 215 70.1 32.5 215 (78, 353)
#>ERR[1] 0.203 0.203
#>
#> Covariance Type (fit$covMethod): r
#> Some strong fixed parameter correlations exist (fit$cor) :
#> cor:THETA[2],THETA[1] cor:THETA[3],THETA[1] cor:THETA[3],THETA[2]
#> -0.310 0.811 0.263
#>
#>
#> Minimization message (fit$message):
#> relative convergence (4) Created on 2021-05-07 by the reprex package (v2.0.0)
skysir
commented
3 years ago Is there possibly a way to fit closed-form models with estimation methods other than foce? Say laplacian?
mattfidler
commented
3 years ago Somehow I missed comment from @skysir ; It is possible to estimate with saem, and nlme using the same interface and work-around. You may also try gnlmm and gnlmm2 if you want to try Gaussian quadrature
I was trying to fit a relatively simple Emax model today which doesn't have any need for time or other parts. (That is where some of the other issues I just reported are coming up...)
When doing so, I needed to have the model have a compartment (that seems reasonable), but it also appears to need
timewhich does not seem so reasonable. I do realize that nlmixr/RxODE has its origins in integrating ODEs, but it seems like removing the requirement of atimecolumn when nothing really needs to be integrated would be helpful. I would assume that time is not required when there are no ODEs defined byd/dtand there are no uses oflinCmt().(Note that for the model below, the only reason I added the
d/dt()was to define the compartment.)Created on 2021-04-27 by the reprex package (v2.0.0)