"Missing value or an infinity produced when evaluating the model"

Question

I'm trying to solve a two-component decay model in R using the nls function, but running into errors. The equation is:

Where t is time, Ctot is C1+C2, and p1 and p2 are known proportions of Ctot.

my data (dd) is:

> head(dd,n=15)
      t    Ctot
1  0.00 6.62
2  0.33 6.45
3  0.50 6.38
4  0.67 6.44
5  0.83 6.38
6  1.00 6.39
7  1.17 6.35
8  1.33 6.33
9  1.50 6.33
10 1.67 6.28
11 1.83 6.17
12 2.00 6.11
13 2.17 6.07
14 2.33 5.89
15 2.50 5.86

Using nls I have tried:

p1 <- 0.3
p2 <- 0.7   
z <- nls(Ctot~(p1*C1*(exp(-k1*t)))+(p2*C2*(exp(-k2*t))), data=dd, start=list(C1=6, C2=0.1, k1=0.01, k2=0.01))

However I am getting:

z <- nls(Ctot~(p1*C1*(exp(-k1*t)))+(p2*C2*(exp(-k2*t))), data=dd, start=list(C1=6, C2=0.1, k1=0.01, k2=0.01))
Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model

I would be grateful if anyone has suggestions!

It may be that the optimum values of `k1` and `k2` are close to each other so that a model with only one term would be better. — G. Grothendieck, Sep 29 '15 at 20:51
Seems the p's and C's might not be separately estimable. If `Ctot == C1+C2`, then wouldn't p1 and p2 need to be 1 or be in a fixed ratio? — IRTFM, Sep 29 '15 at 23:31
I guess that `Ctot==C1+C2` only holds at `t==t0`. Otherwise `C1` and `C2` couldn't be constants, since `Ctot` is not constant over time. — mra68, Sep 29 '15 at 23:43
@mra68: That's what the exp(-k*t) terms are there for .... to drive the two components to zero over time. Ctot[1] == C1[1]+C2[1] puts constraints on what p1 and p2 can be. — IRTFM, Sep 29 '15 at 23:47
What is "t0" in the formula? If it were 0 then `Ctot==p1*C1+p2*C2` at t=0. But `Ctot==C1+C2`. — mra68, Sep 30 '15 at 00:08
@mra68 I assume you are asking this of `squishy`, since I assumed that `t0` meant time=0, so you are now repeating my initial question. — IRTFM, Sep 30 '15 at 00:12

IRTFM · Answer 1 · 2015-09-30T00:29:46.710

The data seems fairly limited and clearly incomplete since it only the head. If we make up some data for testing methods ... and leave out the confusing p1 and p2:

 t=seq(0, 20, by=.3)
 Ctot = 3 * exp( -1 * t) + 4 * exp(-5*t)
 # following hte example on gnm::gnm's help page:

saved.fits <- list(); library(gnm)
for (i in 1:10) {
      saved.fits[[i]] <- suppressWarnings( gnm(Ctot ~ Exp(1 + t, inst = 1) +
                                                      Exp(1 + t, inst = 2),
                    verbose=FALSE))}
plot(Ctot~t)
lines(saved.fits[[3]]$fitted~t)
lines(saved.fits[[3]]$fitted~t,col="red")

I wasn't familiar with the gnm package and so ended up reading the first few sections and then the worked 2 component data fitting example in its vignette: https://cran.r-project.org/web/packages/gnm/vignettes/gnmOverview.pdf . Most of the fits will be as expected, but some will find a local maximum in likelihood that is not a global max:

> saved.fits[[1]]$coefficients
                     (Intercept) Exp(. + t, inst = 1).(Intercept) 
                    1.479909e-12                     1.098612e+00 
          Exp(1 + ., inst = 1).t Exp(. + t, inst = 2).(Intercept) 
                   -1.000000e+00                     1.386294e+00 
          Exp(1 + ., inst = 2).t 
                   -5.000000e+00 
attr(,"eliminated")
[1] 0
> exp( saved.fits[[1]]$coefficients[4] )
Exp(. + t, inst = 2).(Intercept) 
                               4 
> exp( saved.fits[[1]]$coefficients[2] )
Exp(. + t, inst = 1).(Intercept) 
                               3

score 0 · Answer 2 · answered Sep 25 '18 at 18:00

With the data shown in the question it does not seem to work well but if you are open to other parametric models then this 3 parameter model seems reasonable.

fm <- nls(Ctot ~ 1 / (a + b * t^c), dd, st = list(a = 1, b = 1, c = 1))
plot(dd)
lines(fitted(fm) ~ t, dd, col = "red")

"Missing value or an infinity produced when evaluating the model"

2 Answers2