SAEM fit has very large parameter jump at first parameter to a "wrong" value and I don't understand why.

[iamstein]

The model fit below, using the SAEM method, gives the output further below. The initial guess for resp_add_sd is 0.005 and it jumps at iteration 1 to 0.996, and then stays around 1 the whole rest of the estimation routine. The other model parameters all stay right around their initial values. Based on monolix fits of the same model and the nlmixr2 focei fit, 0.005 is around the "right" answer. I'm attaching here both the saem and focei fits.

The SAEM fit has other strange properties as well. For instance, the SE and %RSE of lec50 is absurdly small, and this also does not reflect fits from other software.

Do you have suggestions for how to debug this issue?

Input:

  fit = qs::qread("Task09_saem_v2.qs")
  fit$parHist[1,]
  fit$iniUi$theta

Output:

>   fit$iniUi$theta
       lkout          lr0        lec50        lemax  resp_add_sd resp_prop_sd 
  -5.0359531    0.6151856    4.7706846    6.0867747    0.0050000    0.3000000 
>   fit$parHist[1,]
  iter     lkout       lr0    lec50    lemax V(omega2.r0) resp_add_sd resp_prop_sd
1    1 -5.252822 0.4647881 4.748225 5.817413      0.54872   0.9963458    0.6096235

Small value of SE for lec50:

fit$parFixed["lec50",]
             Parameter Est.       SE    %RSE Back-transformed(95%CI)  BSV(CV%) Shrink(SD)%
lec50 EC50 for B cells 4.54 0.000203 0.00446       93.6 (93.6, 93.6) fix(10.0)      96.6%>

Let me know too if questions like this belong in Discussions rather than Issues. I wasn't entirely sure.

Andy

Task09_saem_v2.qs.zip

Task09_focei.qs.zip

[mattfidler]

The initial estimate for residual error in additive and proprotional errors are ignored in nlmixr2 (and nlmixr) for saem. Wenping was clever and calculated the standard deviations for each of these components from the predicted residuals. This is likely why there are differences between other software and this software. This changes the path of the saem algorithm.

I would have to read the Monolix/NONMEM manuals to see if this is something that they do (but I don't think they do, nor is it always in the manual).

A (implied) feature request is to allow the mcmc process to optimize regardless of if it can actually calculate the residual errors. This also means thinking of a good one-dimensional optimizer, though Monolix may still use Nelder-Meade...

[iamstein]

This could be a case where nlmixr2:focei+Monolix give similar "good" answers while nlmixr2:SAEM gives a "bad" answer, much farther from the global optimum. I'll follow up with you on this 1-1, but I was just wondering, would it be helpful to you for me to create a little case study of this issue?

[mattfidler]

It could be; In general, Wenping's approach may be better, though there are situations where optimization of the residual components may be better.

[iamstein]

Independent of which approach is better, if you cannot specify an initial estimate for the error term with the saem method, then I think nlmixr should at least give a prominent warning that this initial estimate is being ignored and set to a calculated value.

[mattfidler]

It should be in the fit information (when you print it).

[mattfidler]

I'm unsure what "prominant" is too.

Also saem ignores boundaries of parameters (in general), which is also in the warnings (I believe)

[mattfidler]

Hm. On p. 6 and 7 of http://lixoft.com/wp-content/uploads/2016/03/monolixMethodology.pdf it says that the variance is not estimated but calculated in Monolix unless it is a add+prop model (just like nlmixr2/nlmixr).

Perhaps the initial parameter is used as the first estimate of the parameter (which may not be done in saem, or perhaps the calculation method is different). I believe they also say power models are calculated with sufficient statistics (which I don't believe is true in nlmixr2 currently)