[R] Finding starting values for the parameters using nls() or nls2()

[Pinglei Gao]

Hi,

I have some data that i'm trying to fit a double exponential model: data.
Frame (Area=c (521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91, 989.05,
1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),

Retention=c (95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
9.72, 97.92, 71.44, 44.52, 24.44, 15.26) ) and the formula of the double
exponential is: exp (b0*exp (b1*x^th)).



I failed to guess the initial parameter values and then I learned a measure
to find starting values from Nonlinear Regression with R (pp. 25-27):



> cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),

+ Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )

> expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}

> grid.Disperse <- expand.grid(list(b0 = seq(0.01,4, by = 0.01), th =
c(0.02),b1 = seq(0.01, 4, by = 0.01)))

> Disperse.m2a <- nls2(Retention ~expFct(Area, b0, b1,th), data = cl, start
= grid.Disperse, algorithm = "brute-force")

> Disperse.m2a

Nonlinear regression model

model: Retention ~ expFct(Area, b0, th, b1)

data: cl

b0 th b1

3.82 0.02 0.01

residual sum-of-squares: 13596

Number of iterations to convergence: 160000

Achieved convergence tolerance: NA



I got no error then I use the output as starting values to nls2 ():

> nls.m2<- nls2(Retention ~ expFct(Area, b0, b1, th), data = cl, start =
list(b0 = 3.82, b1 = 0.02, th = 0.01))

Error in (function (formula, data = parent.frame(), start, control =
nls.control(), :

Singular gradient



Why? Did I do something wrong or misunderstand something?



Later, I found another measure from Modern Applied Statistics with S (pp.
216-217):



> negexp <- selfStart(model = ~ exp(b0*exp(b1*x^th)),initial =
negexp.SSival, parameters = c("b0", "b1", "th"),

+ template = function(x, b0, b1, th) {})

> Disperse.ss <- nls(Retention ~ negexp(Area, b0, b1, th),data = cl, trace =
T)

b0 b1 th

4.208763 144.205455 1035.324595

Error in qr.default(.swts * attr(rhs, "gradient")) :

NA/NaN/Inf (arg1) can not be called when the external function is called.



Error happened again. How can I fix it? I am desperate.



Best regards,



Pinglei Gao




[[alternative HTML version deleted]]

______________________________________________
R-h...@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Andrew Robinson]

Here are some things to try. Maybe divide Area by 1000 and retention
by 100. Try plotting the data and superimposing the line that
corresponds to the 'fit' from nls2. See if you can correct it with
some careful guesses.

Getting suitable starting parameters for non-linear modeling is one of
the black arts of statistical fitting. Good luck! And don't forget to
check for sensitivity.

Andrew

--
Andrew Robinson
Deputy Director, CEBRA, School of Biosciences
Reader & Associate Professor in Applied Statistics Tel: (+61) 0403 138 955
School of Mathematics and Statistics Fax: +61-3-8344 4599
University of Melbourne, VIC 3010 Australia
Email: a.rob...@ms.unimelb.edu.au
Website: http://www.ms.unimelb.edu.au/~andrewpr

MSME: http://www.crcpress.com/product/isbn/9781439858028
FAwR: http://www.ms.unimelb.edu.au/~andrewpr/FAwR/
SPuR: http://www.ms.unimelb.edu.au/spuRs/

[Bert Gunter]

Well... (inline -- and I hope this isn't homework!)

On Sun, Oct 9, 2016 at 3:05 PM, Andrew Robinson
<A.Rob...@ms.unimelb.edu.au> wrote:
> Here are some things to try. Maybe divide Area by 1000 and retention
> by 100. Try plotting the data and superimposing the line that
> corresponds to the 'fit' from nls2. See if you can correct it with
> some careful guesses.
>
> Getting suitable starting parameters for non-linear modeling is one of

> the black arts of statistical fitting. ...
>
> Andrew

True. But it's usually worthwhile thinking about the math a bit before guessing.

Note that the model can be linearized to:

log(log(Retention)) = b0 + b1*Area^th

So a plot of log(log(Retention)) vs Area may be informative and useful
for finding starting values. e.g., for a grid of th's, do linear
regression fits .

However, when I look at that plot, it seems pretty linear with a
negative slope. This suggests that you may have an overparametrization
problem . i.e. fix th =1 and use the b0 and b1 from the above
regression for starting values.

Do note that this strategy isn't foolproof, as it ignores that the
error term is additive in the above transformed metric, rather than
the original. This can sometimes mislead. But this is just a
heuristic.

Cheers,
Bert

[peter dalgaard]

> On 10 Oct 2016, at 00:40 , Bert Gunter <bgunte...@gmail.com> wrote:
>
> Well... (inline -- and I hope this isn't homework!)
>

Pretty much same as I thought.

Fixing th=0.02 in the grid search looks wrong. Bert's plot is pretty linear, so th=1 is a good guesstimate. There's a slight curvature but to reduce it, you would increase th, not decrease it. Running the regression, as Bert suggests, indicates that b0=5.16 and b1= -0.00024 could work as reasonable starting values. Notice that the grid search had "b1 = seq(0.01, 4, by = 0.01)" which is wrong in both sign and scale.

Andrew's suggestion of dividing Retention by 100 is tempting, since it looks like a percentage, but that would make all Y values less than 1 and the double exponential function as written has values that are always bigger than 1. (It is conceivable that the model itself is wrong, though. E.g. it could be that Retention on a scale from 0 to 1 could be modeled as exp(-something), but we really have no idea of the context here.)

(If this was in fact homework, you should now go and write a proper SelfStart initializer routine for this model. Even if it isn't homework, you do need to study the text again, because you have clearly not understood how self-starting models work.)

-pd

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd....@cbs.dk Priv: PDa...@gmail.com

[ProfJCNash]

I didn't try very hard, but got a solution from .1, 1, .1 with nlxb() from nlmrt. It took a lot
of iterations and looks to be pretty ill-conditioned. Note nlmrt uses analytic derivatives if it
can, and a Marquardt method. It is designed to be a pit bull -- tenacious, not fast.

I'm working on a replacement for this and nls(), but it will be a while. However, I welcome
short scripts like this as tests. My script below

Best, JN

cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),

Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )

expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}

expf2 <- "Retention ~ exp(b0*exp(b1*Area^th))"

# grid.Disperse <- expand.grid(list(b0 = seq(0.01,4, by = 0.01), th =
#c(0.02),b1 = seq(0.01, 4, by = 0.01)))

#> Disperse.m2a <- nls2(Retention ~expFct(Area, b0, b1,th), data = cl, start
#= grid.Disperse, algorithm = "brute-force")

# Disperse.m2a

library(nlmrt)
test <- nlxb(expf2, start= c(b0=.1, b1=1, th=.1), trace=TRUE, data=cl)

[Gabor Grothendieck]

If you are not tied to that model the SSasymp() model in R could be
considered and is easy to fit:

# to plot points in order
o <- order(cl$Area)
cl.o <- cl[o, ]

fm <- nls(Retention ~ SSasymp(Area, Asym, R0, lrc), cl.o)
summary(fm)

plot(Retention ~ Area, cl.o)
lines(fitted(fm) ~ Area, cl.o, col = "red")

[ProfJCNash]

Despite nlmrt "solving" the OP's problem, I think Gabor's suggestion likely gives a more sensible approach
to the underlying modelling problem.

It is, of course, sometimes important to fit a particular model, in which case nls2 and nlmrt are set up to grind away.
And hopefully the follow-up to nlmrt I'm working on will have enough capability in getting analytic derivatives
to work for a wider class of models. Note that functional approaches in nlmrt and minpack.lm allow users to
provide derivatives. Too many users think numerical approximations are a panacea, but my experience is that
most problems benefit from very accurate derivatives, of which analytic expressions are generally the best.

JN

[ProfJCNash]

The key lines are

library(nlmrt)
test <- nlxb(expf2, start= c(b0=.1, b1=1, th=.1), trace=TRUE, data=cl)

Thus I started with .1 1 and .1. The "solution" from nlxb, which is using analytic derivatives
and a very aggressive Marquardt code to keep trying even in bad situations, was
as you included below. Note that the singular values of the Jacobian are given (they are
recorded on the same table as the parameters, but do NOT correspond to the parameters.
The placement was simply a tidy place to put these numbers.)

The ratio of these sv's is 1.735e+16/0.004635 or approx 4E+18, so the condition number
of the traditional Gauss Newton approach is about 1E+37. Not a nice problem!

You probably should reformulate.

JN

On 16-10-10 10:41 AM, Pinglei Gao wrote:
> Thanks very much for your kindness help. I run your script then came out
> lots of outputs and I also studied the solution you posted. Forgive my
> ignorance, I still can't find the suitable starting values. Did I
> misunderstand something?
>
> Best,
>
> Pinglei Gao
>
> -----邮件原件-----
> 发件人: ProfJCNash [mailto:profj...@gmail.com]
> 发送时间: 2016年10月10日 10:41
> 收件人: Gabor Grothendieck; Pinglei Gao
> 主题: Re: [R] Finding starting values for the parameters using nls() or
> nls2()
>
> I forgot to post the "solution" found by nlmrt:
>
> nlmrt class object: x
> residual sumsquares = 1086.8 on 15 observations
> after 5001 Jacobian and 6991 function evaluations
> name coeff SE tstat pval gradient
> JSingval
> b0 5.3274e-14 NA NA NA -6.614e+13
> 1.735e+16
> b1 33.5574 NA NA NA -3.466
> 11518
> th -0.00721203 NA NA NA -740.8
> 0.004635
>
>
> Note the singular values -- this is the worst SV(max)/SV(min) ratio I've
> observed!
>
> JN

[Pinglei Gao]

Thanks very much for taking time on this. Your assistances are very much
appreciated. But, I am afraid that I still have a question to bother you.

I am working on a paper about weed seeds dispersal with harvest machine. I
found three general models for seed dispersal and retention after a review
of relevant literature. All models were optimized using nonlinear least
squares via the nls function in the statistical package R. The model that
best described the data will be determined by comparing Akaike Information
Criterion (AIC) values and the model with the lowest AIC score will be
selected.

The first general model incorporated simple exponential and power
exponential functions, its starting value was easily to be found. But, I am
stuck with model 2 (which was mentioned previously) and model 3 with the
form: Retention = (b0*Area^th+1)^b1. The model 3 is totally different to
others. I tried the measures that you were mentioned. But I still can’t
find suitable starting values because of my limited knowledge. I hope you
can do me the favor again. I can send the draft to you when I finished the
paper, if it is necessary. Maybe you can give me some constructive
suggestion about statistic and model construction and I can name you as a
coauthor for your contributions.

Best,

Pinglei Gao

-----邮件原件-----
发件人: peter dalgaard [mailto:pda...@gmail.com]
发送时间: 2016年10月10日 7:41
收件人: Pinglei Gao
抄送: Andrew Robinson; R help (r-h...@r-project.org); Bert Gunter

主题: Re: [R] Finding starting values for the parameters using nls() or
nls2()

-pd

>>> + 33.04, 23.46,

[Pinglei Gao]

Thanks very much for your kindness help. I run your script then came out
lots of outputs and I also studied the solution you posted. Forgive my
ignorance, I still can't find the suitable starting values. Did I
misunderstand something?

Best,

Pinglei Gao

-----邮件原件-----

发件人: ProfJCNash [mailto:profj...@gmail.com]
发送时间: 2016年10月10日 10:41
收件人: Gabor Grothendieck; Pinglei Gao

主题: Re: [R] Finding starting values for the parameters using nls() or
nls2()

I forgot to post the "solution" found by nlmrt:

nlmrt class object: x
residual sumsquares = 1086.8 on 15 observations
after 5001 Jacobian and 6991 function evaluations
name coeff SE tstat pval gradient
JSingval
b0 5.3274e-14 NA NA NA -6.614e+13
1.735e+16
b1 33.5574 NA NA NA -3.466
11518
th -0.00721203 NA NA NA -740.8
0.004635

Note the singular values -- this is the worst SV(max)/SV(min) ratio I've
observed!

JN

______________________________________________

[dave fournier]

Actually this converges very nicely if you use these starting values
that I obtained with
AD Model Builder

th 9.1180e-01
b0 5.2104e+00
b1 -4.6725e-04

The R result looks like

nls.m2
Nonlinear regression model
model: Retention ~ expFct(Area, b0, b1, th)
data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,

2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,

1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,

18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,

15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")
b0 b1 th
5.2104466 -0.0004672 0.9118029
residual sum-of-squares: 686.8

Number of iterations to convergence: 1

Achieved convergence tolerance: 1.75e-06

[dave fournier]

I've spent quite a bit of time trying to convince people on various
lists that the solution to these kinds of
problems lies in the stable parameterization of the model. I write the
solutions in AD Model Builder because it
is easy. But R people are generally stuck in R (or mired) so the
message somehow always gets lost.

I thought I would bite the bullet and figure out how to do this in R.
It turns out that one can fit this model
with only 6 function evaluations using nls2. I apologize in advance
for my execrable R code. But it does do
the job.

The data were fit using 3 calls to nls2. For the first call only the
parameter th was estimated. This converged
in 4 function evaluations. For the second call all three of the
parameters were estimated
(but the other 2 parameters are different from b0, b1) .This converged
to the solution in four function evaluations.

The final call to nls2 uses the converged values to calculate b0,b1 and
theta and starts from there.
As you can see the model is already converged. This final call to nls2
is used to get the standard error
estimates for the parameters.

> nls.m1
Nonlinear regression model
model: Retention ~ expFct1(Area, y1, yn, th)
data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,

2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,

1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,

18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,

15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")

th
0.9862
residual sum-of-squares: 730

Number of iterations to convergence: 2
Achieved convergence tolerance: 1.616e-06

> nls.m2
Nonlinear regression model
model: Retention ~ expFct2(Area, y1, yn, c, d, th)
data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,

2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,

1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,

18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,

15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")

c d th
0.9716 1.1010 0.9118
residual sum-of-squares: 686.8

Number of iterations to convergence: 4
Achieved convergence tolerance: 2.398e-06

> nls.m
Nonlinear regression model

model: Retention ~ expFct(Area, b0, b1, th)

data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,

2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,

1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,

18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,

15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")
b0 b1 th

5.2104452 -0.0004672 0.9118040
residual sum-of-squares: 686.8

Number of iterations to convergence: 0
Achieved convergence tolerance: 1.384e-06

Parameters:
Estimate Std. Error t value Pr(>|t|)
b0 5.2104452 0.4999594 10.422 2.29e-07 ***
b1 -0.0004672 0.0013527 -0.345 0.7358
th 0.9118040 0.3583575 2.544 0.0257 *

So what is the stable parameterization for this model. To simplify let
x be the independent variable and y be the
dependent variable and write t instead of th So the model is

y = exp(b0*exp(b1*x^t) (1)

Now it would be nice to reorder the x's so that they are monotone
increasing or decreasing, but in this case the
first x is almost the largest and the last x is almost the smallest
(slight reach) so they will do. Call them x1 and xn.
The new parameters of the model are the predicted y's for x1 and xn.
Call them y1 and yn. Note that y1 and yn
are parameters, not observations.


The data were fit using 3 calls to nls2. For the first call only the
parameter th was estimated. this converged
in 4 function evaluestions. for the second call all three of the
parameters were estimated
(but the other 2 parameters are different from b0, b1) .This converged
to the solution in four function evaluations.

The final call to nls2 uses the converged values to calculate b0,b1 and
theta and starts from there.
as you can see the model is already converged. this final call to nls2
is used to get the standard error
estimates for the parameters.

> nls.m1
Nonlinear regression model
model: Retention ~ expFct1(Area, y1, yn, th)
data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,

2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,

1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,

18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,

15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")

th
0.9862
residual sum-of-squares: 730

Number of iterations to convergence: 2
Achieved convergence tolerance: 1.616e-06

> nls.m2
Nonlinear regression model
model: Retention ~ expFct2(Area, y1, yn, c, d, th)
data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,

2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,

1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,

18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,

15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")

c d th
0.9716 1.1010 0.9118
residual sum-of-squares: 686.8

Number of iterations to convergence: 4
Achieved convergence tolerance: 2.398e-06

> nls.m
Nonlinear regression model

model: Retention ~ expFct(Area, b0, b1, th)

data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8,

2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52,

1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36,

18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44,

15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class
= "data.frame")
b0 b1 th

5.2104452 -0.0004672 0.9118040
residual sum-of-squares: 686.8

Number of iterations to convergence: 0
Achieved convergence tolerance: 1.384e-06

Parameters:
Estimate Std. Error t value Pr(>|t|)
b0 5.2104452 0.4999594 10.422 2.29e-07 ***
b1 -0.0004672 0.0013527 -0.345 0.7358
th 0.9118040 0.3583575 2.544 0.0257 *

So what is the stable parameterization for this model. To simplify letx
be the independent variable and y be the
dependent variable and write t insted of th So the model is

y = exp(b0*exp(b1*x^t) (1)

Now it would be nice to reorder the x's so that they are monotone
increasing or decreasing, but in this case the
first x is almost the largest and the last x is almost the smallest
(slight reach) so they will do. Call them x1 and xn.
The new parameters of the model are the predicted y's for x1 and xn.
Call them y1 and yn. Note that y1 and yn
are parameters, not observations.

y1 = exp(b0*exp(b1*x1^t)
(2)
yn = exp(b0*exp(b1*xn^t)

One can solve for b1 and b0 in terms of y1 and yn.

b1=log(log(y1)/log(yn))/(x1^t-xn^t);
(3)
b0=log(y1)*exp(-b1*x1^t);


To use this we do the fitting of the model in two stages. In the first
stage we use the obvious
estimates of y[1] for y1 and y[n] for yn and fix these values and
estimate t using the obvious
starting value of 1 for t.

In the second stage we set

y1=c*y[1]

yn=d*y[n]

and estimate c, d, and t using the obvious starting values of 1 for c
and d.
That's it. this converges as stated in 6 function evaluations. This
technique works for
may curves such as 3,4,5 parameter logistic double exponential
vonBertalanffy
(where Schnute and I first discovered it in the early 1980's before
anyone was born).

One caveat is that usually not all values of y1 and yn are
permissible. For example in this case if b>0
then y1 and yn are both >1. To make this routine bulletproof one needs
to impose these kinds of
constraints. But this was not needed here.

Here is my R code for this model.

library("nls2")
cl=data.frame(Area=c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,

989.05,
1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),
Retention=c (95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )

y1 <<- cl$Retention[1];
yn <<- cl$Retention[length(cl$Retention)]

expFct1 <- function(x, y1, yn,th)
{
x1=x[1]
xn=x[length(x)]
b1=log(log(y1)/log(yn))/(x1^th-xn^th)
b0=log(y1)*exp(-b1*x1^th)
exp(b0*exp(b1*x^th));
}

expFct2 <- function(x, y_1,y_n,c,d,th)
{
x1=x[1]
xn=x[length(x)]
y1=c*y_1
yn=d*y_n
b1=log(log(y1)/log(yn))/(x1^th-xn^th)
b0=log(y1)*exp(-b1*x1^th)
exp(b0*exp(b1*x^th));
}

expFct <- function(x, b0, b1,th)
{
exp(b0*exp(b1*x^th))
}


nls.m1<- nls2(Retention ~ expFct1(Area, y1, yn, th), data = cl,
start = list(th = 1))

th_start=coef(nls.m1)


nls.m2<- nls2(Retention ~ expFct2(Area, y1, yn, c, d, th), data = cl,
start = list(c=1, d=1, th = th_start))

x=cl$Area
x1=x[1]
xn=x[length(x)]
th_start=coef(nls.m2)[3]
cy1=coef(nls.m2)[1]*y1

dyn=coef(nls.m2)[2]*yn
b1_start=log(log(cy1)/log(dyn))/(x1^th_start-xn^th_start)
b0_start=log(cy1)*exp(-b1_start*x1^th_start)

th_start
b1_start
b0_start
nls.m<- nls2(Retention ~ expFct(Area, b0, b1, th), data = cl,
start = list(b0 = b0_start, b1 = b1_start, th = th_start))

nls.m1
nls.m2
nls.m

[Pinglei Gao]

Dear Dave
Thanks for your kindness help. I am sorry, I am on a filed survey these days. I did not check my email for a long time. Please forgive me for the misunderstanding brought to you. Your answer works well for the model: Retention = exp (b0*exp (b1*Area^th)). Could you find the starting values for the another model: Retention = (b0*Area^(th+1))^b with the same data for me. I also asked the same question on R-help list before.

Bests,

PG

-----邮件原件-----
发件人: dave fournier [mailto:da...@otter-rsch.com]
发送时间: 2016年10月20日 5:47
收件人: r-h...@r-project.org; gaopi...@163.com

主题: Re: [R] Finding starting values for the parameters using nls() or nls2()

Actually this converges very nicely if you use these starting values that I obtained with AD Model Builder

th 9.1180e-01
b0 5.2104e+00
b1 -4.6725e-04

The R result looks like

nls.m2
Nonlinear regression model

model: Retention ~ expFct(Area, b0, b1, th)
data: structure(list(Area = c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91, 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2), Retention = c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46, 9.72, 97.92, 71.44, 44.52, 24.44, 15.26)), .Names = c("Area", "Retention"), row.names = c(NA, -15L), class = "data.frame")
b0 b1 th

5.2104466 -0.0004672 0.9118029
residual sum-of-squares: 686.8

Number of iterations to convergence: 1
Achieved convergence tolerance: 1.75e-06

______________________________________________

[dave fournier]

Unfortunately this problem does not appear to be well posed.

Retention = (b0*Area^(th+1))^b

If b0, th, and b are the parameter only the product (th+1)*b is determined.

This comes from noting that powers satisfy



(a^b)^c = a^(b*c)



So your model can be written as



(b0^b) * Area^((th+1)*b)

[David Winsemius]

> On Oct 25, 2016, at 9:29 AM, dave fournier <da...@otter-rsch.com> wrote:
>
>
>
> Unfortunately this problem does not appear to be well posed.
>
> Retention = (b0*Area^(th+1))^b
>
> If b0, th, and b are the parameter only the product (th+1)*b is determined.
>
> This comes from noting that powers satisfy
>
>
>
> (a^b)^c = a^(b*c)
>
>
>
> So your model can be written as
>
>
>
> (b0^b) * Area^((th+1)*b)
>

... which nicely completes the thread since one model had:

th1 = 9.1180e-01
b01= 5.2104e+00
b11 = -4.6725e-04
(th1+1)*b11
[1] -0.0008932885


b0 = 5.2104466 ; b1 = -0.0004672 ; th = 0.9118029
((th+1)*b1)
[1] -0.0008931943

So both the R's nls2 and AD_MOdel_Builder results yield that same predictions for any given data point at least up to four decimal places.

So perhaps there was some other physical interpretation of those parameters and there exists an underlying theory that would support adding some extra constraints?

--
David Winsemius
Alameda, CA, USA

[dave fournier]

>
>> On Oct 25, 2016, at 9:29 AM, dave fournier <davef at otter-rsch.com>

I'm not really sure what your point is. The OP has two models. One is well
posed and the other is not. I was discussing solution of the former model
which is well posed. The form of the model, using a, b, and t and x,y to
simplify the expression is

y = exp( a * exp(b * x^t) )

My point is that there are many models like this where the obvious
parameterization (something like the parameters the user is interested
in) leads to parameter estimates with highly correlated errors.
This does not necessarily mean that
the model is badly determined. The model exists independent of the
particular parameterization. To fix the ideas assume there are n
observations
(x_i,y_i) and simplify by assuming that x_1<=x_2<=...<=x_n
(but not all equal lets hope)

A stable parameterization of the model can often be obtained by
picking as new parameters y1 and yn where y1 and yn are the predicted
values for y_1 and y_n, that is

y1 = exp( a * exp(b * x_1^t) )
yn = exp( a * exp(b * x_n^t) )

Then one solves for some of the original model parameters in terms of y1
and yn.
The particular way this is done depends on the model. Often some has some
linear parameters and the procedure is easy. For this model as I showed
one can solve for a and b in terms of y1 and yn.
Then one can fit the model easily with AD Model Builder or
nls2 using this parameterization. nls2 provides the standard errors for
the parameter estimates. The AD Model Builder solution provides the
estimated variance covariance matrix of the parameter estimates via the
standard maximum likelihood Hessian calculations. It also provides the
covariance for any desired dependent variable. So one can fit the model
using y1,yn, and t and get the covariance matrix for a,b, and t
in one step. In nls2 one needs to fit the model using y1,yn and then
solve for a and b and run the model again from that point. That is no big
deal, and I showed how to do it, but it is one more step for the user.

It is interesting to see the correlation matrix for the parameter estimates
and the dependent variables.
std dev correlation
1 logt -9.233e-02 3.4026e-01 1.0000
2 c 9.7164e-01 3.7006e-02 -0.2690 1.0000
3 d 1.1010e+00 1.7852e-01 -0.7495 0.0325 1.000
4 t 9.1180e-01 3.1025e-01 1.0000 -0.2690 -0.749 1.0000
5 a 5.2104e+00 4.3404e-01 -0.9781 0.4191 0.621 -0.9781 1.000
6 b -4.6725e-04 1.1714e-03 0.9994 -0.2836 -0.726 0.9994 -0.984 1.00

Here the independent variable are c, d, and logt
where y1=c*y_1 y2=d*y2
That is just an additional thing so that one can start with c=1 and d=1
Also logt is used where t=exp(logt) which makes sure t>0.
Notice that the correlation between c and d is 0.0325 which is small.


If a and b were the parameters of the model

4 t 9.1180e-01 3.1025e-01 1.0000
5 a 5.2104e+00 4.3404e-01 -0.9781 1.000
6 b -4.6725e-04 1.1714e-03 0.9994 -0.984 1.00

One can see how close to 1 and -1 the correlations are.

Notice that the parameter b is very badly determined.
So rather than saying the model is no good one sees that
the model is no good if one want to get information about
the parameter b. In contrast the parameter a is fairly well
determined and the parameter t is "kind of" determined.

Because of this parameter confounding nls2 is extremely sensitive
to the starting parameter values using this parameterization.
If I change the parameters by about 15% or less as below

#b0start=5.210445
#b1_start=-0.0004672373
#th_start=0.911804
b0_start=5.7
b1_start=-0.00055
th_start=1.05

I get the dreaded

Error in (function (formula, data = parent.frame(), start, control =
nls.control(), :

singular gradient

So there is nothing wrong with either the data or the model. The model
just needs to be given a stable parameterization.