Error with nlmixr 2.0.4 for model that ran OK with earlier versions

[nholford]

I have a warfarin pd model that I have used to test previous versions of nlmixr (latest was 1.1.1.-3). When I run this with nlmixr 2.0.4 I get the following errors.

pd.fit.CEF <- nlmixr(pd.ce.emax, data.pd, est="focei", control=foceiControl(print=100000))

Any suggestions on what to do?

Thanks Nick

unhandled error message: EE:lsoda -- at t = 6.03775e-009 and step size _C(h) = 2.30322e-014, the corrector convergence failed repeatedly or with fabs(_C(h)) = hmin @(lsoda.c:887) Theta reset (ETA drift) Theta reset (zero gradient values); Switch to bobyqa Error in .fitFun(.ret) : Evaluation error: Starting values violate bounds. Restart 1 Error in .fitFun(.ret) :

Evaluation error: Starting values violate bounds. Restart 2 Error in .fitFun(.ret) : Evaluation error: Starting values violate bounds. Restart 3 Error in .fitFun(.ret) : Evaluation error: Starting values violate bounds. Error in (function (data, inits, PKpars, model = NULL, pred = NULL, err = NULL, : Could not fit data. Calls: source ... nlmixr.function -> do.call -> -> .collectWarnings Execution halted

[mattfidler]

Hi Nick,

Thank you for trying the new version on nlmixr. In this version of nlmixr we use symngine instead of sympy. Because of this, the symbolic expressions are in different a different order. This can affect the rounding errors that accumulate in the optimization which causes a different solution for the same problem.

In general, you can try the following to have the same problem run:

1. Use bobyqa instead of the default optimization. In this case the outer problem is optimized by a gradient free approach. This will likely avoid this specific problem. This can be done with pd.fit.CEF <- nlmixr(pd.ce.emax, data.pd, est="focei", control=foceiControl(print=100000, outerOpt="bobyqa"))
2. Another option is to change the initial estimates, this can help in some of the problems that can be seen.

We appreciate your support. I know that the joint PKPD model works still since it is run before every release, its results are here:

https://nlmixrdevelopment.github.io/nlmixr/articles/multiple-endpoints.html#focei-fits-1

I know that @RikSchoemaker has also tested a group of models on the warfarin dataset so he may have some insights on this particular problem.

[nholford]

Hi Matt,

Thanks for the advice. I tried using the magic bobyqa control and the model ran. It did so with these warnings:

Warning messages: 1: In (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE, : initial ETAs were nudged; (can control by fnythoceiControl(etaNudge=., etaNudge2=)) 2: In (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE, : ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.)) 3: In (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE, : last objective function was not at minimum, possible problems in optimization 4: In (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE, : parameter estimate near boundary; covariance not calculated use 'getVarCov' to calculate anyway 5: In (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE, : Undefined 'dvid' integer values in data: 1, 2

Most of them look like the kind of thing I get from NONMEM and largely ignore. Warning number 5 doesn’t mean anything to me. ‘dvid’ is a covariate in the data set. Why is saying that values that are defined ‘1, 2’ are undefined?

I am glad to hear you are using some warfarin data and models in your testing. I was not aware that you were using this warfarin example to show how to use multiple end points. Nice 😊

Perhaps you might add the effect compartment model that caused the default controls to fail as an example to show how bobyqa might help?

pd.ce.emax <- function() { ini({ popemax <- c(0,0.9,0.999) tc50 <- log(0.5) tkout <- log(0.05) te0 <- log(100)

eta.emax ~ .5
eta.c50  ~ .5
eta.kout ~ .5
eta.e0 ~ .5

eps.pdadd <- 100

}) model({ poplogit = log(popemax/(1-popemax)) logit=poplogit+eta.emax emax = 1/(1+exp(-logit)) c50 = exp(tc50 + eta.c50) kout = exp(tkout + eta.kout) e0 = exp(te0 + eta.e0)

# PK parameters from input dataset
ktr = IKTR
ka = IKA
cl = ICL
v = IV

cp = center/v

d/dt(depot) = -ktr * depot
d/dt(gut) =  ktr * depot -ka * gut
d/dt(center) =  ka * gut - cl / v * center

d/dt(ce) = kout*(cp - ce)

effect=e0*(1-emax*ce/(c50+ce))

effect ~ add(eps.pdadd)

}) }

-- Nick Holford, Professor Clinical Pharmacology Dept Pharmacology & Clinical Pharmacology, Bldg 503 Room 302A University of Auckland,85 Park Rd,Private Bag 92019,Auckland,New Zealand office:+64(9)923-6730 mobile:NZ+64(21)46 23 53 FR+33(6)62 32 46 72 email: @.*** http://holford.fmhs.auckland.ac.nz/ http://orcid.org/0000-0002-4031-2514 Read the question, answer the question, attempt all questions

From: Matthew Fidler @.> Sent: Monday, 26 April 2021 5:15 PM To: nlmixrdevelopment/nlmixr @.> Cc: Nick Holford @.>; Author @.> Subject: Re: [nlmixrdevelopment/nlmixr] Error with nlmixr 2.0.4 for model that ran OK with earlier versions (#512)

Hi Nick,

Thank you for trying the new version on nlmixr. In this version of nlmixr we use symngine instead of sympy. Because of this, the symbolic expressions are in different a different order. This can affect the rounding errors that accumulate in the optimization which causes a different solution for the same problem.

In general, you can try the following to have the same problem run:

1. Use bobyqa instead of the default optimization. In this case the outer problem is optimized by a gradient free approach. This will likely avoid this specific problem. This can be done with pd.fit.CEF <- nlmixr(pd.ce.emax, data.pd, est="focei", control=foceiControl(print=100000, outerOpt="bobyqa"))
2. Another option is to change the initial estimates, this can help in some of the problems that can be seen.

We appreciate your support. I know that the joint PKPD model works still since it is run before every release, its results are here:

https://nlmixrdevelopment.github.io/nlmixr/articles/multiple-endpoints.html#focei-fits-1https://nlmixrdevelopment.github.io/nlmixr/articles/multiple-endpoints.html#focei-fits-1

I know that @RikSchoemakerhttps://github.com/RikSchoemaker has also tested a group of models on the warfarin dataset so he may have some insights on this particular problem.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/nlmixrdevelopment/nlmixr/issues/512#issuecomment-826920889, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AGI42QEOHVYKKLJW6YRGVFTTKV7OPANCNFSM43P7J5VA.

[mattfidler]

You are right; It may be a good idea to add it as well to see what is going on here. Maybe more of an example of a full model development. I'm currently at the upper limit of the time allowed for automatic building and testing of the models. I need to figure out how to run these longer models in a more automatic way to make sure something we add doesn't break nlmixr or change things too much.

[mattfidler]

Most of them look like the kind of thing I get from NONMEM and largely ignore. Warning number 5 doesn’t mean anything to me. ‘dvid’ is a covariate in the data set. Why is saying that values that are defined ‘1, 2’ are undefined?

In this case you have a dvid in your dataset that is ignored, since multiple endpoints are not used. That is what this warning means. As you state, in this case it is meangless since it is a single-endpoint model.

[nholford]

Hi Matt, Thanks for this explanation. Whether this is a single or multiple endpoint model should be detectable by the software and thus avoid this warning which serves no purpose except to confuse the user. Best wishes, Nick

From: Matthew Fidler @.> Sent: Wednesday, 28 April 2021 2:23 AM To: nlmixrdevelopment/nlmixr @.> Cc: Nick Holford @.>; Author @.> Subject: Re: [nlmixrdevelopment/nlmixr] Error with nlmixr 2.0.4 for model that ran OK with earlier versions (#512)

Most of them look like the kind of thing I get from NONMEM and largely ignore. Warning number 5 doesn’t mean anything to me. ‘dvid’ is a covariate in the data set. Why is saying that values that are defined ‘1, 2’ are undefined?

In this case you have a dvid in your dataset that is ignored, since multiple endpoints are not used. That is what this warning means. As you state, in this case it is meangless since it is a single-endpoint model.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/nlmixrdevelopment/nlmixr/issues/512#issuecomment-828048107, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AGI42QFYXSEWJP4CYQZVBZTTK5IOVANCNFSM43P7J5VA.

[nholford]

Hi Matt, As you probably know this warfarin example was developed to show 3 models for warfarin PKPD. The data is simulated so the true model is known. Comparing the 3 models (turnover (true), effect compartment and immediate effect) can be done using nlmixr diagnostics. The lack of usefulness of the numerous residual plots can be easily appreciated. The VPC shows which is the correct model. So I agree extending the example to a more complete picture of the typical modelling process would be useful. I am happy to help. Best wishes, Nick

From: Matthew Fidler @.> Sent: Wednesday, 28 April 2021 2:18 AM To: nlmixrdevelopment/nlmixr @.> Cc: Nick Holford @.>; Author @.> Subject: Re: [nlmixrdevelopment/nlmixr] Error with nlmixr 2.0.4 for model that ran OK with earlier versions (#512)

You are right; It may be a good idea to add it as well to see what is going on here. Maybe more of an example of a full model development. I'm currently at the upper limit of the time allowed for automatic building and testing of the models. I need to figure out how to run these longer models in a more automatic way to make sure something we add doesn't break nlmixr

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/nlmixrdevelopment/nlmixr/issues/512#issuecomment-828046358, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AGI42QFLATKXNWZOTIEL45DTK5H4LANCNFSM43P7J5VA.

[MissLesley2]

Hi Matt,

I have the same issues using nlmxir2 package. The error is in below. Do you have any suggestions on what I can do? Could nlmxir2 be able to deal with effect compartment or transit compartments? And also, could lagtime will be an issue for sequential PKPD modeling?

Error in .foceiFitInternal(.ret) : Starting values violate bounds
Restart 1
Theta reset (ETA drift)
Error : Could not fit data"

myPDFit <- nlmixr(pk.turnover.emax, df, est="focei", control=foceiControl(print=100000, outerOpt="bobyqa"))

pk.turnover.emax <- function() {
  ini({

    tke <- log(1)
    temax <- logit(0.5)
    tec50 <- log(100)
    tkout <- log(0.05)
    tkin <- log(0.05*70000)
    tkd <-log(2)

    ##
    eta.ke ~ 1
    eta.emax ~ 1
    eta.ec50  ~ 1
    eta.kout ~ 1
    eta.kin ~ 1
    eta.kd ~ 1

    ## 
    prop.err <- 0.2
  })
  model({

    ke = exp(tke + eta.ke)
    emax=expit(temax+eta.emax)
    ec50 =  exp(tec50 + eta.ec50)
    kout = exp(tkout + eta.kout)
    kin = exp(tkin + eta.kin)
    kd = exp(tkd + eta.kd)

    # PK parameters from input dataset
    Ka = IKA
    CL = ICL
    Vp = IVP
    Vc = IVC
    F  = IF
    Q = IQ
    e0 = IE0
    lagtime = ILAG

    #PK
    d/dt(A1) = -Ka* A1
    d/dt(A2) = Ka*A1 - CL / Vc * A2 + Q/Vp * A3 - Q/Vc * A2
    d/dt(A3) = Q/Vc * A2 - Q/Vp * A3
    cp = A2 / Vc

    lag(A1) = lagtime
    f(A1)= F

    ##PD
    effect(0) = e0 
    d/dt(ce) = ke*(cp-ce)
    d/dt(cd)= -kd*cd
    d/dt(effect) = kin*(1+emax*ce/(ec50+ce)) -kout*effect +kd*cd

    effect ~ prop(prop.err)
  })
}

[mattfidler]

Do you have any suggestions on what I can do?

I always suggest doing your analysis in a reproducible manner. If you need the same results, you should use the same version of R, same version of compilers, same platform. I hope to add some information to a nlmixr2 fit so this can be discovered after the time of the fit. I believe you should always at least print the sessionInfo() at the end so you can manually setup the R environment as it was at the time of your fit. You can use tools like renv or hedgehog.

Could nlmxir2 be able to deal with effect compartment or transit compartments?

Yes.

And also, could lagtime will be an issue for sequential PKPD modeling?

Possibly, but it depends on the model. Often it isn't. In fact, the treatment of lag times changed in nlmixr2 and is more likely to give a reasonable estimate of the parameters.

This is discussed here:

https://blog.nlmixr2.org/blog/2022-11-10-nlmixr2-nonmem/