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

[Andrew Robinson]

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 at 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.