nlmixr2 / nlmixr2extra

Extra utilities for nlmixr objects
https://nlmixr2.github.io/nlmixr2extra/
GNU General Public License v3.0
3 stars 1 forks source link

No bootstrap result due to matrix issues #59

Closed andreasmeid closed 1 year ago

andreasmeid commented 1 year ago

Dear Matt,

I am trying to obtain SE or confidence intervals for an nlmixr2 object (named test) after a FOCEi fit using bootstrapping. However, something seems to be wrong with the covariance matrix.

fit2 <- suppressMessages(bootstrapFit(test, nboot = 50)) Warning message: In cov2cor(covMatrix) : diag(.) had 0 or NA entries; non-finite result is doubtful

Can you execute the function on enclosed object without error or where may be the problem?

Many thanks and best wishes

Andreas

test.zip

mattfidler commented 1 year ago

Hi @andreasmeid

Thanks for you patience. I will look at this shortly.

Matt

mattfidler commented 1 year ago

It appears to me that it worked, but didn't update the parameter estimates $parFixed

> test$bootCovMatrix
              tVc tCl tVt tCld tf tRTVeffect tRTVeffectgut     prop.err         eta.f        eta.Vc
tVc             0   0   0    0  0          0  0.0000000000  0.000000000  0.0000000000  0.0000000000
tCl             0   0   0    0  0          0  0.0000000000  0.000000000  0.0000000000  0.0000000000
tVt             0   0   0    0  0          0  0.0000000000  0.000000000  0.0000000000  0.0000000000
tCld            0   0   0    0  0          0  0.0000000000  0.000000000  0.0000000000  0.0000000000
tf              0   0   0    0  0          0  0.0000000000  0.000000000  0.0000000000  0.0000000000
tRTVeffect      0   0   0    0  0          0  0.0000000000  0.000000000  0.0000000000  0.0000000000
tRTVeffectgut   0   0   0    0  0          0  0.0037690186 -0.005706746  0.0002575555  0.0007785759
prop.err        0   0   0    0  0          0 -0.0057067460  0.009631444 -0.0004256940 -0.0013402967
eta.f           0   0   0    0  0          0  0.0002575555 -0.000425694  0.0879265260  0.0184142593
eta.Vc          0   0   0    0  0          0  0.0007785759 -0.001340297  0.0184142593  0.3333001147
eta.Cl          0   0   0    0  0          0 -0.0032338040  0.005392710  0.0016038064  0.0105152980
                    eta.Cl
tVc            0.000000000
tCl            0.000000000
tVt            0.000000000
tCld           0.000000000
tf             0.000000000
tRTVeffect     0.000000000
tRTVeffectgut -0.003233804
prop.err       0.005392710
eta.f          0.001603806
eta.Vc         0.010515298
eta.Cl         0.004676072
> test$bootCorMatrix
              tVc tCl tVt tCld tf tRTVeffect tRTVeffectgut    prop.err       eta.f      eta.Vc
tVc             0   0   0    0  0          0    0.00000000  0.00000000  0.00000000  0.00000000
tCl             0   0   0    0  0          0    0.00000000  0.00000000  0.00000000  0.00000000
tVt             0   0   0    0  0          0    0.00000000  0.00000000  0.00000000  0.00000000
tCld            0   0   0    0  0          0    0.00000000  0.00000000  0.00000000  0.00000000
tf              0   0   0    0  0          0    0.00000000  0.00000000  0.00000000  0.00000000
tRTVeffect      0   0   0    0  0          0    0.00000000  0.00000000  0.00000000  0.00000000
tRTVeffectgut   0   0   0    0  0          0    0.06139233 -0.94717168  0.01414805  0.02196692
prop.err        0   0   0    0  0          0   -0.94717168  0.09813992 -0.01462823 -0.02365579
eta.f           0   0   0    0  0          0    0.01414805 -0.01462823  0.29652407  0.10756638
eta.Vc          0   0   0    0  0          0    0.02196692 -0.02365579  0.10756638  0.57732150
eta.Cl          0   0   0    0  0          0   -0.77029838  0.80356459  0.07909543  0.26635649
                   eta.Cl
tVc            0.00000000
tCl            0.00000000
tVt            0.00000000
tCld           0.00000000
tf             0.00000000
tRTVeffect     0.00000000
tRTVeffectgut -0.77029838
prop.err       0.80356459
eta.f          0.07909543
eta.Vc         0.26635649
eta.Cl         0.06838181
> test
── nlmixr² FOCEi (outer: nlminb) ──

          OBJF      AIC      BIC Log-likelihood
FOCEi 35.29408 418.3831 434.9492      -204.1916

── Time (sec test$time): ──

        setup table compress  other covariance
elapsed 0.003  0.05     0.01 32.277     20.048

── Population Parameters (test$parFixed or test$parFixedDf): ──

                 Parameter  Est. Back-transformed BSV(CV%) Shrink(SD)%
tVc                     Vc  7.83         2.51e+03     89.5      3.13% 
tCl                     Cl  6.29              540     79.4      63.0% 
tVt                     Vt   8.3         4.04e+03                     
tCld                   Cld  7.78         2.39e+03                     
tf                       F     0                1     59.9     0.944% 
tRTVeffect       RTVeffect     0                1                     
tRTVeffectgut RTVeffectgut -0.57            0.566                     
prop.err                     0.5              0.5                     

  Covariance Type (test$covMethod): boot50
    other calculated covs (setCov()): boot5, boot7, boot50, boot50, boot50, boot5, boot50
  No correlations in between subject variability (BSV) matrix
  Full BSV covariance (test$omega) or correlation (test$omegaR; diagonals=SDs) 
  Distribution stats (mean/skewness/kurtosis/p-value) available in test$shrink 
  Information about run found (test$runInfo):
   • gradient problems with initial estimate; see $scaleInfo 
   • parameter estimate near boundary; covariance not calculated: "tRTVeffectgut" use 'getVarCov' to calculate anyway 
   • 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=)) 
  Censoring (test$censInformation): No censoring
  Minimization message (test$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 test is a modified tibble): ──
# A tibble: 203 × 32
  ID     TIME    DV  PRED      RES     WRES IPRED   IRES  IWRES CPRED  CRES  CWRES eta.f eta.Vc eta.Cl
  <fct> <dbl> <dbl> <dbl>    <dbl>    <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>  <dbl>
1 110    0.25  0.8  0.725  0.0752   0.104    1.10 -0.303 -0.549 0.663 0.137 0.120   1.14  0.791  0.474
2 110    0.5   2.34 1.29   1.05     0.865    2.07  0.269  0.260 1.16  1.18  0.568   1.14  0.791  0.474
3 110    0.75  1.73 1.74  -0.00923 -0.00585  2.93 -1.20  -0.819 1.55  0.185 0.0646  1.14  0.791  0.474
# ℹ 200 more rows
# ℹ 17 more variables: AaAtor <dbl>, Ac <dbl>, At <dbl>, RTV_adminis <dbl>, Acyp <dbl>, Acypgut <dbl>,
#   Vc <dbl>, Cl <dbl>, Vt <dbl>, Cld <dbl>, f <dbl>, RTVeffect <dbl>, RTVeffectgut <dbl>, Cc <dbl>,
#   Ct <dbl>, tad <dbl>, dosenum <dbl>
# ℹ Use `print(n = ...)` to see more rows
> 

The zero warning comes from some of the components being fixed (in fact most of the thetas in this case).

I believe that it will update as it is written now if the fit originally tries a covariance and then you update it with the bootstrapped covariance.

mattfidler commented 1 year ago

Added an issue https://github.com/nlmixr2/nlmixr2est/issues/400

andreasmeid commented 1 year ago

Dear Matt,

many thanks, I also think that something happens in the background (the Cov and Cor matrices are formed), but no result is generated. Did I misunderstand something, is there already a solution above (or what exactly do I have to do for this) or are you researching further?

Best wishes

Andreas

mattfidler commented 1 year ago

The "no result" is likely that the $parFixed is not updated. This occurs when the standard error is missing from the first model because it wasn't calculated correctly or the calculation was skipped. The linked bug/feature is where I will put updates on how to fix this in the future.

In the mean time, you can use the confidence interval to manually calculate %RSE, CI etc.

andreasmeid commented 1 year ago

Dear Matt,

ah, I see, so I'will wait for the update. However, I did not understand correctly what could be done in the mean time. Do you have specific functions/code lines to run for this purpose?

Many thanks and best wishes

Andreas