2

I am trying to fit oral pharmacokinetic function using R nls. The function can be seen here: enter image description here

I tried the following:

dat <- data.frame(
        Time = c(10, 15, 20, 30, 40, 60, 90, 120, 180, 210, 240, 300, 360),
        C_PO = c(0,0.28,0.55,1.2,2,1.95,1.85,1.6,
                 0.86,0.78,0.6,0.21,0.18))

plot(dat$C_PO ~ dat$Time, data = dat, log = "y")

fit <- nls(C_PO~ka*100* (exp(-k*Time) - exp(-ka*Time))/(v*(ka- k)),
    data = dat,
    start = list(ka = 0.09, k = 0.01, v = 50))
#> Error in nls(C_PO ~ ka * 100 * (exp(-k * Time) - exp(-ka * Time))/(v * : singular gradient

The code is trying to fit a non-linear model according to the concentration values C_PO over a time period Time. The initial values were calculated using graph stripping methods. However, the error message reads "singular gradient." How can this be fixed?

UseR10085
  • 7,120
  • 3
  • 24
  • 54
PPenton
  • 25
  • 5
  • You can use `minpack.lm` or `gslnls` for better non-linear model fitting. – UseR10085 Jun 16 '23 at 06:26
  • As it is, one issue is that division by zero is possible. As a first step, you should substitute `deltak = ka - k` and `ka = deltak + k`. Then you should set lower boundaries for the parameters (and `algorithm = "port"` if you use `nls`). However, I was not able to get a successful fit where `deltak` was not at the lower boundary. Even if you get a successfull fit with one of the optimizers suggested by comments and answers, you should be wary. Better repeat the experiment a few times to get more (reliable) data. – Roland Jun 16 '23 at 06:33
  • Are you sure you have the equation right? Often difficulty in fitting suggests that the wrong model is being used. Usually here one uses the biexponential function. See `?SSbiexp` and in PK package see `?biexp` . – G. Grothendieck Jun 17 '23 at 14:50
  • @G.Grothendieck I find these self starting functions interesting. But In this case, the equation I have should be able to descibe it. My question is if those functions can give also the formula? – PPenton Jun 17 '23 at 15:12
  • @UseR10085 Can you elaborate briefly how do they differ? – PPenton Jun 17 '23 at 15:12
  • @Roland Thank you for your suggestion. I will try to avoid zeros by adding constraints from not on. – PPenton Jun 17 '23 at 15:13
  • @PPenton, The biexponential does fit your data better than the model in the question. – G. Grothendieck Jun 17 '23 at 17:05
  • @G.Grothendieck The thing is that in PK I don't always easily understanding the data relationships. So, these derived forms are very helpful. – PPenton Jun 18 '23 at 00:12
  • The problem is that what you believe is the model isn't really how it works. Good models don't depend on the 7th decimal and they fit easily. – G. Grothendieck Jun 18 '23 at 12:30

1 Answers1

1

You can use gslnls package for fitting non-linear model like

library(gslnls)

dat <- data.frame(
  Time = c(10, 15, 20, 30, 40, 60, 90, 120, 180, 210, 240, 300, 360),
  C_PO = c(0,0.28,0.55,1.2,2,1.95,1.85,1.6,
           0.86,0.78,0.6,0.21,0.18))

plot(dat$C_PO ~ dat$Time, data = dat, log = "y")

fit <- gsl_nls(C_PO ~ ka*100* (exp(-k*Time) - exp(-ka*Time))/(v*(ka- k)), 
                    data = dat, 
                    start = c(ka = 0.09, k = 0.01, v = 50), 
                    control = gsl_nls_control(maxiter = 500))

summary(fit)

#> Formula: C_PO ~ ka * 100 * (exp(-k * Time) - exp(-ka * Time))/(v * (ka - 
#>     k))
#> 
#> Parameters:
#>     Estimate Std. Error t value Pr(>|t|)
#> ka    0.0132     4.0027   0.003    0.997
#> k     0.0132     4.0025   0.003    0.997
#> v    21.0181  6372.7704   0.003    0.997
#> 
#> Residual standard error: 0.3363 on 10 degrees of freedom
#> 
#> Number of iterations to convergence: 18 
#> Achieved convergence tolerance: 0
UseR10085
  • 7,120
  • 3
  • 24
  • 54