[R] Nonlinear regression aborting due to error

[Dennis Fisher]

Colleagues,

I am fitting an Emax model using nls.  The code is:
	START			<- list(EMAX=INITEMAX, EFFECT=INITEFFECT, C50=INITC50)
	CONTROL		<- list(maxiter=1000, warnOnly=T)
	#FORMULA		<- as.formula(YVAR ~ EMAX - EFFECT * XVAR^GAMMA / (XVAR^GAMMA + C50^GAMMA))	## alternate version of formula
	FORMULA		<- as.formula(YVAR ~ EMAX - EFFECT / (1 + (C50/XVAR)^GAMMA))
	FIT				<- nls(FORMULA, start=START, control=CONTROL, trace=T)

If GAMMA equals 10-80, nls converges successfully and the fit tracks the fit from a smoother (Supersmoother).  However, if I attempt to estimate GAMMA using:
	START			<- list(EMAX=INITEMAX, EFFECT=INITEFFECT, C50=INITC50, GAMMA=INITGAMMA)
GAMMA increases rapidly to > 500 and nls terminates with:
	Error in chol2inv(object$m$Rmat()) : 
	  element (4, 4) is zero, so the inverse cannot be computed
	In addition: Warning message:
	In nls(FORMULA, start = START, control = CONTROL, trace = T) :
	  singular gradient

I also tried fixing GAMMA to > 1000 and I get a similar error message:
	Error in chol2inv(object$m$Rmat())
	 : element (2, 2) is zero, so the inverse cannot be computed 
	In addition: Warning message: 
	In nls(FORMULA, start = START, control = CONTROL, trace = T)
	 : singular gradient 

The data do not suggest a very large value for GAMMA so I am surprised that the estimate is increasing so rapidly.  I attempted to use the port algorithm with an upper bound on GAMMA but the upper bound is reached rapidly, suggesting that the data support a large  value for GAMMA.
 
A subset of the data (with added noise) is shown below.  A GAMMA value of 1280 triggers the error with this subset

XVAR	<- c(26, 31.3, 20.9, 24.8, 22.9, 4.79, 19.6, 18, 19.6, 9.69, 21.7, 26.6, 27.8, 9.12, 10.5, 20.1, 16.7, 14.1, 10.2, 19.2, 24.7, 34.6, 
26.6, 25.1, 5.98, 13.4, 15.7, 9.59, 7.39, 21.5, 15.7, 12.4, 19.2, 17.8, 19.7, 27.1, 25.6, 36.4, 22.9, 8.68, 27, 25.9, 33.3, 24.2, 
21.4, 31, 19.1, 18.7, 23.5, 19.4, 10.3, 12.8, 13.9, 18.5, 21, 15.2, 18.9, 9.12, 16.9, 12.9, 29.5, 15.5, 7.34, 8.97, 8.04, 23.7, 
16.3, 37.6, 35.2, 13.7, 28.1, 29.5, 15.1, 26, 6.52)


YVAR	<- c(-34.2, -84.2, -71.1, -91.9, -104.1, -23.2, -27.2, -13.4, -143.2,  24.7, -72.1, -38, 25.2, -8, -34.1, -15.1, -112.6, -93.5, -130.9, 
-127.8, -118.7, -53.5, -29.8, 98, 0, -37.6, -99.4, 57.9, 0.2, -62.2, -27.3, 8.3, -51.6, -111.6, -25.6, -51.7, -106.4, -85.1, 
-63.1, -60.8, -27.7, -20.7, 22.9, -49.4, -85.7, -90.9, -107,  -20.6, -36.3, -40.2, 39.8, -55, -54.5, -103.9, -53.1, -2.3, -72.3, 
-65.6, -57.8, -64.4, -129.1, 10.4, -9.9, -29.6, -40.8, 52, -94, 8.8, -98.8, 28, -16.3, -99.2, -48.5, -111.9, -15.4)

I suspect that I am making a conceptual error in the use of nls.  Any help would be appreciated.  If a different function to fit nonlinear regression would work better, please direct me.

Dennis

Dennis Fisher MD
P < (The "P Less Than" Company)
Phone: 1-866-PLessThan (1-866-753-7784)
Fax: 1-866-PLessThan (1-866-753-7784)
www.PLessThan.com