[R] nls() newbie convergence problem

[Katharine Mullen]

dydata
   x1 x2   y
1   9 27 248
2   9 22 213
3   4 23 190
4  11 16 183
5  -6 25 144
6  11 14 169
7  -4 13  72
8   2  8  73
9  10 13 156
10  8 30 263
11 12 10 147
12 -7  5  -2
13  0 10  75
14 12  0  77
15  9  8 115
16 12 24 245
17 34 23 370
18 12  1  84
19 10 37 324
20 26 30 371

weight <- function(y,x1,x2,b0,b1,b2)
{
  pred <-  b0+b1*x1 + b2*x2
  parms <- abs(b1*b2)^(1/3)
  (y-pred)/parms
}

gmfit <-  nls(~weight(y,x1,x2,b0,b1,b2),data=dydata,trace=TRUE,
              start=list(b0=1,b1=1,b2=1), control=list(warnOnly=TRUE))

library(minpack.lm)

weight2 <- function(p, y,x1,x2)
{
  pred <-  p[1]+p[2]*x1 + p[3]*x2
  parms <- abs(p[2]*p[3])^(1/3)
  (y-pred)/parms
}

gmfit2 <-  nls.lm(par = c(1,1,1),
                  fn = weight2, y=dydata$y,x=dydata$x1,x2=dydata$x2,
                  control=list(nprint=1))


dydata is from Norman R. Draper, Yonghong (Fred) Yang, Generalization
of the geometric mean functional relationship, Computational
Statistics & Data AnalysisVolume 23, Issue 3, 9 January 1997, Pages
355-372; I don't think the attachment got through.

nls uses an orthogonality convergence criterium.  Bates talks about this
in the post at
https://stat.ethz.ch/pipermail/r-devel/2000-August/021059.html and it's
also described in the book by Bates and Watts (Nonlinear regression
analysis its applications).  Apparently here the residuals are so small
they are effectively 0, and this criteria does not work; there is a
warning about using zero-residual data in the help page for nls.

nls.lm from the package minpack.lm stops if any of a few different
criteria are met; in this case it stops at the time you think is right.

On Wed, 4 Jun 2008, Bernard Leemon wrote:

> I'm sure this must be a nls() newbie question, but I'm stumped.
> I'm trying to do the example from Draper
> and Yang (1997).  They give this snippet of S-Plus code:
>
> Specify the weight function:
> weight < - function(y,x1,x2,b0,b1,b2)
> {
> pred <-  b0+b1*x1 + b2*x2
> parms <- abs(b1*b2)^(1/3)
> (y-pred)/parms
> }
> Fit the model
> gmfit < -nls(~weight(y,x1,x2,b0,b1,b2), observe,list("starting value"))
>
> in converting this to R, I left the weight function alone and replaced the
> nls() with
>
> gmfit <-
>  nls(~weight(y,x1,x2,b0,b1,b2),data=dydata,trace=TRUE,start=list(b0=1,b1=1,b2=1))
>
> where dydata is the appropriate data.frame.
>
> nls() quickly (6 iterations) finds the exact values from Draper & Yang for
> b0, b1, and b2 but
> despite reporting a discrepancy of only 3.760596e-29 by the 7th iteration,
> it merrily goes on
> to 50 iterations and thinks it never converged.  how do I tell nls() that
> I'm actually quite
> happy with 3.760596e-29 and it need not work further?
>
> thanks for any suggestions.
>
> gary mcclelland (aka bernie)
> univ of colorado
>
> 	[[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.
>