[R] Non-linear regression problem: R vs JMP (long)

About this list Date view Thread view Subject view Author view Attachment view

From: Andrew Robinson (andrewr@uidaho.edu)
Date: Thu 11 Mar 2004 - 03:19:25 EST


Message-id: <200403100819.25760.andrewr@uidaho.edu>

Dear R friends,

I know that this topic has been mulled over before, and that there is a
substantial difference between the convergence criteria for JMP and those for
R. I apologize that this is somwehat raking cold coals.

Summary:

A model/data combination achieves convergence in JMP, and survives a
reasonably rigorous examination (sensible parameter estimates, well-behaved
surface, confidence intervals exclude 0). The same combination fails to
achieve convergence in R despite starting from the estimates reported in JMP.
What can I do?

Full story:

I am collaborating on a project which presently requires the fit of a seven
parameter non-linear function to 874 observations. The function is:

freeze.d <- deriv3(~ dbh^b * u^dbh * a1 * smi^a2 * exp(a3*pbal) *
  (1 + a4*exp(a5*ba)) * cr^a6 * dh5^a7,
  c("a1","a2","a3","a4","a5","a6","a7","b","u"),
  function(dbh, smi, pbal, ba, cr, dh5, a1, a2, a3, a4, a5, a6, a7, b, u){})

The dbh, smi, pbal, ba, cr, and dh5, are known.

The data do not require the level of intricacy reflected in the function. It
is debatable whether they support it. However, my collaborator is anxious to
avoid linear models because of their relatively poor extrapolative
properties. So, that is an ongoing discussion.

He has achieved convergence using non-linear least squares in JMP. JMP uses
the union of three criteria: small change in objective function, small change
in parameter values, and small change in gradient. He can achieve
convergence in any of the three by reducing the other two sufficiently. He
arrives at the same point of convergence from a range of different starting
values. He gets garbage solutions for other starting points, but recognizes
these as such and explores elsewhere. At the point of convergence the model
has reasonable behavior, the confidence intervals for all the parameters
exclude 0, the surface seems to fit the data just fine, and the parameter
estimates all make sense. Some of them are highly correlated (up to -0.95)
but by no means all. In short, it seems to me to be a very thorough and
careful effort.

I'm trying to reproduce the model in R using nls. The convergence criterion
is different than JMP: it's a relative-offset convergence criterion. Even if
I start with his converged parameter estimates, nls will consume as many
iterations as I allow it and produce an error, e.g. "number of iterations
exceeded maximum of 50000". I have the trace on, and for the last n-2
iterations the parameter estimates do not change, at least in the digits that
trace provides. Increasing the tolerance does not seem to help.

So, I wonder: is it possible that JMP could produce results that are garbage
but stand up to reasonable scrutiny? Alternatively, is there some further
element of the fit that I should advise him to examine? I've been reading
Bruce McCullough's work and I know that for example S-plus was one of only
two packages to declare "ns" instead of providing estimates with 0 accurate
digits. But, how does one know that is what is happening in any given
situation, if the model looks good? Finally, it's possible that I am using
nls wrongly: here is my call

freeze.nls.d.1 <-
    nls(dd.5 ~ freeze.d(dbh.0, smi, pbal.0, sba.0, cr.0, dh.5,
                        a1, a2, a3, a4, a5, a6, a7, b, u),
        data = trees, control=ctl.obj, trace=T,
        start = list(b = 0.2, u = 1, a1 = 0.4, a2 = -0.1, a3 = -0.01, a4 = 3,
a5 = -0.02, a6 = 0.25, a7 = 0.40))

Any insights will be appreciated!

Andrew

-- 
Andrew Robinson                      Ph: 208 885 7115
Department of Forest Resources       Fa: 208 885 6226
University of Idaho                  E : andrewr@uidaho.edu
PO Box 441133                        W : http://www.uidaho.edu/~andrewr
Moscow ID 83843                      Or: http://www.biometrics.uidaho.edu
No statement above necessarily represents my employer's opinion.

______________________________________________ R-help@stat.math.ethz.ch mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


About this list Date view Thread view Subject view Author view Attachment view

This archive was generated by hypermail 2.1.3 : Fri 23 Apr 2004 - 18:36:35 EST