Comparing nls() to nls2() - what am I doing wrong

Question

I am trying to emulate an nls() fit with nls2() via brute-force, when nls() works, so that I can look to a second option when it doesn't.

What am I doing wrong in how I have specified nls2() below? Is it my grid values?

nls() output gives:

Formula: y ~ C0 * exp(-lambda1 * x) + C1 * exp(-lambda2 * x)

Parameters:
                            Estimate Std. Error t value Pr(>|t|)   
C0.intercept_initial       -0.108025   0.012703  -8.504  0.00105 **
lambda1.log_slope_initial   0.438553   0.072185   6.075  0.00371 **
C1.intercept_terminal       0.107742   0.012531   8.598  0.00101 **
lambda2.log_slope_terminal  0.057911   0.008045   7.199  0.00197 **

with plot:

whereas nls2() gives:

Formula: y ~ C0 * exp(-lambda1 * x) + C1 * exp(-lambda2 * x)

Parameters:
          Estimate Std. Error t value Pr(>|t|)
C0       1.000e+01  2.043e+14       0        1
C1      -1.000e+01  5.285e+11       0        1
lambda1  0.000e+00  2.043e+14       0        1
lambda2  0.000e+00  5.285e+11       0        1

Code:

# Data
dat <-  data.frame(x = c(0,2,4,8,12,24,48,72), y = c(0.000,0.05,0.0671,0.068,0.05,0.0250,0.0103,0.0043))

# Get initial parameters from SSbiexp
params <-  getInitial(y ~ SSbiexp(x, intercept_initial, log_slope_initial, intercept_terminal, log_slope_terminal), data = dat)
# Fit nls model
nls_mod_bi <-  nls(y ~ C0*exp(-lambda1*x) + C1*exp(-lambda2*x), start = list(C0 = params[1], lambda1 = exp(params[2]), C1 = params[3], lambda2 = exp(params[4])), data = dat)
summary(nls_mod_bi)
xNew <- seq(0, 240, length.out = 100) # new grid of times
yNew <- predict(nls_mod_bi, list(x = xNew)) # predict activity
dfNew <-  data.frame(x = xNew, y = yNew)
dfNew <-  dfNew |> slice(-1) # remove first row which contains -ve prediction for time 0
ggplot(dfNew, aes(x = x, y = y)) +
  geom_line() +
  geom_point(data = dat, aes(x = x, y = y), size = 3) +
  xlab("x") + ylab("y") +
  theme_bw(base_size = 20)

# Now fit nls2 model by brute-force
fo <- y ~ C0*exp(-lambda1*x) + C1*exp(-lambda2*x)
grd <- data.frame(C0      = c(-10, 10),
                  C1      = c(-10, 10),
                  lambda1 = c(0, 1),
                  lambda2 = c(0, 1))
nls_mod_bi2 <- nls2(fo,
                    data = dat,
                    start = grd,
                    algorithm = "brute-force")
summary(nls_mod_bi2)

Answer 1

The default number of points is 50 and with 4 parameters and given that 50^(1/4) = 2.65 you are essentially using a grid of less than 3x3x3x3 which is not enough.

Given that this function is partially linear it would make more sense to use the "plinear-brute-force" algorithm. It does not require starting values for the linear parameters (C0, C1) so it reduces the problem size down to 2 parameters (lambda1, lambda2). To use that algorithm supply a formula whose right hand side is a matrix whose i-th column multiplies the implicit i-th linear parameter.

fo3 <- y ~ cbind(C0 = exp(-lambda1 * x), C1 = exp(-lambda2 * x))
nls_mod_bi3.0 <- nls2(fo3, data = dat, start = grd[3:4], 
  algorithm = "plinear-brute-force")
coef(nls_mod_bi3)
##    lambda1    lambda2    .lin.C0    .lin.C1 
##  0.2857143  0.1428571 -0.3030034  0.2995545 
 
nls(fo3, data = dat, start = coef(nls_mod_bi3.0)[1:2], alg = "plinear")
## Nonlinear regression model
##   model: y ~ cbind(C0 = exp(-lambda1 * x), C1 = exp(-lambda2 * x))
##    data: dat
##  lambda1  lambda2  .lin.C0  .lin.C1 
##  0.43855  0.05791 -0.10802  0.10774 
##  residual sum-of-squares: 4.679e-05
##
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 7.392e-06