4th parameter of Boltzmann sigmoid must be greater than .9 in R

Question

I'm trying to fit a 4 parameter boltzmann sigmoid and get an error: "Error in nls(y ~ a0 + (a1 - a0)/(1 + exp((a2 - x)/a3)), start = list(a0 = max(y), : singular gradient"

I have figured out that the code runs if the a3 parameter is set to .9 or greater, or -.9 or less. Does anyone have the reason this is? I want to provide a starting parameter for a3 as the slope according to the description on this website: http://www.originlab.com/doc/Origin-Help/Boltzmann-FitFunc . That is why I have the linear fit coefficient a3.s, but the result is < .9 and I get the error. Is there a way to estimate a3.s prior to use as starting parameter for nls? I am simply using a linear fit of the midpoint of the sigmoid +/- 10 x units - is that the correct interpretation of the a3 parameter?

Here is my code:

#fit boltzman sigmoid
a0.s=max(y); a1.s=min(y); a2.i=which.min( abs(((a0.s+a1.s)/2) - y) ); a2.s=x[a2.i]
lin.x.i=x<a2.s+10 & x>a2.s-10
a3.s=unname(coef(lm(y[lin.x.i]~so[lin.x.i]))[2])
fit <- nls(y ~ a0 + (a1-a0)/(1+exp((a2-x)/a3)), 
start=list(a0=max(y), a1=min(y), a2=a2.s,a3=.9) , trace=TRUE)
params=coef(fit)
curve(params[1]+(params[2]-params[1])/(1+exp((params[3]-x)/params[4])), 1,100,col='black',add=T,type='l')

Here is the data:

x=c( 75,  40,  90,  55, 15, 100,  10,  70,  90,  50,  15,   5,   5,  70, 100,  20,  60,  65,  20,  50,  30,  85,  60,  80,  55,  40,  45,  95,  10,  55, 60,  10,  35,  80,  75,  25,  30,   5,  35,  50, 100,  40,  30,  80,  20,  45,  25,  25,  95,  95,  65,  35,  90,  85,  70,  15,  75,  45,  85,  65);

y=c(4.673686, 0.034781, 5.014355, 0.843847, 0.013337, 4.214557, 0.015299, 5.017280, 4.327815, 0.041139, 0.008704, 0.007437, 0.005125, 4.725786, 3.869776, 0.018725, 4.514051, 3.232932, 0.012979, 0.257651, 0.028170, 4.723512, 2.676991, 5.018232, 0.633399, 0.040133, 0.051864, 5.019395, 0.006505, 0.642376, 2.752317, 0.010827, 0.029303, 4.050711, 3.698887, 0.018385, 0.029491, 0.013894, 0.032034, 0.053761, 5.029349, 0.038272, 0.032619, 5.030450, 0.022356, 0.053421, 0.025370, 0.024763, 4.948973, 3.254528, 1.149153, 0.038530, 4.612227, 4.048692, 4.809153, 0.016246, 5.014711, 0.062841, 5.026961, 2.951881)

Related to this question: the formula on the linked to website has a slight variation in the equation, where the the a2 parameter is used in the form "x-a2", while the equation I provided, and got from my data acquisition software's curve fitting function is the one I provided in the code with "a2-x". Which form of the Boltzmann is correct? Does the difference matter?

Answer 1

Why are you reinventing the wheel? Use the native function SSfpl for your model.

fitr <- nls(y ~ SSfpl(x, a1, a0, a2, ma3))

The parameter ma3 is -a3 in your notation, but otherwise the parametrization is identical, and you get slightly better convergence.

You should probably be using weighted least squares, since ordinary least squares assumes the variability of $y-E(y)$ does not depend on the value of $x$; which is clearly violated in your data.

Answer 2

I guess this would also give an idea:

df <- data.frame(x,y)

# using drm() in "aomisc" package should give a solution
model.df <- drm(y ~ x, fct = L.4(), data = df)
summary(model.df)
plot(model.df, log="", main = "df")