R help - Aug 2009 - Help with nls and error messages singular gradient

Hi All,

I'm trying to run nls on the data from the study by Marske (Biochemical
Oxygen Demand Interpretation Using Sum of Squares Surface. M.S. thesis,
University of Wisconsin, Madison, 1967) and was reported in Bates and Watts
(1988).

Data is as follows, (stored as mydata)

  time  bod
1    1 0.47
2    2 0.74
3    3 1.17
4    4 1.42
5    5 1.60
6    7 1.84
7    9 2.19
8   11 2.17

I then run the following;
#Plot initial curve
plot(mydata$time, mydata$bod,xlab="Time (in
days)",ylab="biochemical oxygen
demand (mg/l) ")

model <- nls(bod ~ beta1/(1 - exp(beta2*time)), data mydata, start=list(beta1
= 3, beta2 = -0.1),trace=T)

The start values are recommended, (I have used these values in SPSS without
any problems, SPSS returns values of Beta1 = 2.4979 and Beta2 = -2.02 456)

but return the error message,
Error in nls(bod ~ beta1/(1 - exp(beta2 * time)), data = mydata, start
list(beta1 = 3,  : singular gradient

Can anyone offer any advice?

Thanks in advance

Mike









-- 
Michael Pearmain
Senior Analytics Research Specialist

“Statistics are like women; mirrors of purest virtue and truth, or like
whores to use as one pleases”

Google UK Ltd
Belgrave House
76 Buckingham Palace Road
London SW1W 9TQ
United Kingdom
t +44 (0) 2032191684
mpearmain@google.com

If you received this communication by mistake, please don't forward it to
anyone else (it may contain confidential or privileged information), please
erase all copies of it, including all attachments, and please let the sender
know it went to the wrong person. Thanks.

	[[alternative HTML version deleted]]

Hi, Michael,

I think the SPSS answer is wrong. Your starting values are way off.
Look at this plot for verification:

con <- textConnection("time  bod
1    1 0.47
2    2 0.74
3    3 1.17
4    4 1.42
5    5 1.60
6    7 1.84
7    9 2.19
8   11 2.17")
mydata <- read.table(con, header = TRUE)
close(con)

beta <- c(3, -0.1) # your initial values
beta <- c(2.4979, -2.02456)  # SPSS answer
mycurve <- function(x) {
  beta[1]/(1 - exp(beta[1] * x))
}
curve(mycurve, from = 1, to = 11,
      ylim = range(mydata$bod, mycurve(mydata$time)))
points(mydata$time, mydata$bod)

You might want to look at package nls2 which allows a brute force grid
search to find some starting values. Or rethink the equation you're
trying to fit.

HTH,

--sundar

On Tue, Aug 25, 2009 at 9:10 AM, Michael Pearmain<mpearmain at google.com>
wrote:
> Hi All,
>
> I'm trying to run nls on the data from the study by Marske (Biochemical
> Oxygen Demand Interpretation Using Sum of Squares Surface. M.S. thesis,
> University of Wisconsin, Madison, 1967) and was reported in Bates and Watts
> (1988).
>
> Data is as follows, (stored as mydata)
>
> ?time ?bod
> 1 ? ?1 0.47
> 2 ? ?2 0.74
> 3 ? ?3 1.17
> 4 ? ?4 1.42
> 5 ? ?5 1.60
> 6 ? ?7 1.84
> 7 ? ?9 2.19
> 8 ? 11 2.17
>
> I then run the following;
> #Plot initial curve
> plot(mydata$time, mydata$bod,xlab="Time (in
days)",ylab="biochemical oxygen
> demand (mg/l) ")
>
> model <- nls(bod ~ beta1/(1 - exp(beta2*time)), data > mydata,
start=list(beta1 = 3, beta2 = -0.1),trace=T)
>
> The start values are recommended, (I have used these values in SPSS without
> any problems, SPSS returns values of Beta1 = 2.4979 and Beta2 = -2.02 456)
>
> but return the error message,
> Error in nls(bod ~ beta1/(1 - exp(beta2 * time)), data = mydata, start >
list(beta1 = 3, ?: singular gradient
>
> Can anyone offer any advice?
>
> Thanks in advance
>
> Mike
>
>
>
>
>
>
>
>
>
> --
> Michael Pearmain
> Senior Analytics Research Specialist
>
> ?Statistics are like women; mirrors of purest virtue and truth, or like
> whores to use as one pleases?
>
> Google UK Ltd
> Belgrave House
> 76 Buckingham Palace Road
> London SW1W 9TQ
> United Kingdom
> t +44 (0) 2032191684
> mpearmain at google.com
>
> If you received this communication by mistake, please don't forward it
to
> anyone else (it may contain confidential or privileged information), please
> erase all copies of it, including all attachments, and please let the
sender
> know it went to the wrong person. Thanks.
>
> ? ? ? ?[[alternative HTML version deleted]]
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> 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.
>
>