[R] nls, gnls, starting values, and covariance matrix

[Spencer Graves]

	  Question:  I note that b1 is close to 0.  In other words, do you get 
any significant reduction in the quality of the fit by removing it from 
the regression equation?  If no, doing so resolved the identifiability 
problem mentioned by Tom Blackwell.

	  If b1 needs to be nonzero in the equation, then I suggest another 
parameterization:  Have you considered replacing a0=c00+c01, b0=c00-c01, 
a1=c10+c11, and b1=c10-c11, so c00=(a0+b0)/2, c01=(a0-b0)/2, 
c10=(a1+b1)/2, and c11=(a1-b1)/2?  With this, we get the following:

	  Y ~ c00+c10*x + log(exp(c01+c11*x)+exp(-c01-c11*x)).

Since b1 is close to 0, both c10 and c11 are approximately a1/2.  Since 
a1 is always positive, you could further replace c11 by exp(c111), and 
break the identifiability problem that way.

hth, spencer graves
################
Also, not all algorithms are equally robust.  McCullougn (either 1998 or
1999, "Assessing the reliability of statistical software: Part I or
II'', The American Statistician,  52: 358-366 or 53:  149-159) reported
that he often had trouble getting "nls" in S-Plus to converge.  He
recommended using "nlminb" to get convergence and then using the output
from nlminb as starting values in "nls" to get confidence intervals that
nls provided but nlminb did not.  Today, I use "optim" in the MASS
library in place of nlminb.

hth, spencer graves
#################
	  I'm not familiar with gnls, but Bates and Watts (1988, Nonlinear
Regression Analysis and Its Applications, Wiley, esp. pp. 256-259)
establish that changing the way a problem is parameterized can have a
major impact on numerical stability and on the adequacy of approximate
normal theory for confidence intervals, etc.

	  If you let a0 = c0+c1 and b0 = c0-c1, then your model can be
rewritten as follows:

	  Y ~ c0 + log(exp(c1-a1*X)+exp(-c1-b1*X))

hth, spencer graves

Thomas W Blackwell wrote:
  > Simon  -
  >
  > There's a symmetry in the model you are fitting, and the error
  > message returned sounds to me as though it is referring to that.
  >
  > Could you try a model formula of the form
  >
  > Y ~ a0 + log(exp(-c0-a1*X) + exp(+c0-b1*X))
  >
  > Maybe you will need to make the two slopes identifiable in much
  > the same fashion that I've done for the two intercepts, in order
  > to get it to work.  Or maybe not.  I'll leave that up to you.
  >
  > HTH  -  tom blackwell  -  u michigan medical school  -  ann arbor  -
  >
  > On Sat, 19 Apr 2003 sdfrost at ucsd.edu wrote:
  >
  >
  >>Dear R-Help,
  >>
  >>I'm trying to fit a model of the following form using gnls. I've
fitted it
  >>using nlsList with the following syntax:
  >>
  >>nlsList(Y~log(exp(a0-a1*X)+exp(b0-b1*X))|K,start=list
  >>(a0=6,a1=0.2,b0=4.5,b1=0.001),data=data.frame(Y=y,X=X,K=k)))
  >>
  >>which works just fine:
  >>
  >><snip>
  >>
  >>Coefficients:
  >>         a0        a1       b0            b1
  >> 1 5.459381 0.5006811 5.137458 -0.0040548687
  >> 2 5.761496 0.1716723 6.359151 -0.0022802595
  >> 3 5.683510 0.5436838 5.906742 -0.0007788076
  >> 4 6.225745 0.2807003 5.875803 -0.0008351051
  >> 5 6.558350 0.1388707 5.071080  0.0014594212
  >> 6 5.483639 0.2757080 2.406683 -0.0003282243
  >> 7 5.746064 0.4354105 5.883882 -0.0002577279
  >> 8 5.448679 0.3385350 2.851571  0.0011627360
  >> 9 5.259762 0.5654369 5.498967  0.0015381718
  >>10 6.546022 0.8008781 4.913085  0.0051150166
  >>12 5.602982 1.1538595 5.008253 -0.0006087786
  >>13 6.452605 0.1752357 6.229393  0.0007899073
  >>15 5.937199 0.2214811 4.980386  0.0081102533
  >>16 5.998689 0.2925840 6.077816  0.0062388250
  >>
  >>However, I'd like to be able to fit the model using gnls. The format 
is a
  >>little different, but I get an error when I use the following syntax:
  >>
 
 >>gnls(Y~log(exp(a0-a1*X)+exp(b0-b1*X)),params=a0+a1+b0+b1~K,start=list(rep(c
  >>(6.02,0.2,4.5,0.001),16)),data=data.frame(Y=y,X=x,K=k),control=list
  >>(msVerbose=TRUE,apVar=FALSE,returnObject=TRUE))
  >>
  >>Error in gnls(Y ~ log(exp(a0 - a1 * X) + exp(b0 - b1 * X)), params =
a0 +  :
  >>        Approx. covariance matrix for parameter estimates not of full
rank
  >>
  >>I assume that I'm getting the format of my starting values wrong. Any
  >>suggestions would be greatly appreciated.
  >>
  >>Best wishes
  >>Simon
  >
  >
  > ______________________________________________
  > R-help at stat.math.ethz.ch mailing list
  > https://www.stat.math.ethz.ch/mailman/listinfo/r-help