There are so many issues here that I doubt it can be adequately addressed in the SO post, but that should get you started.

First, it looks like you want Tmax < max(dat$x)

for example <35. This causes a problem because then Tmax - x < 0

for some values x

, and when you try to raise a negative number to cardinality (in the second term of the formula), you get NA

. This is the cause of the error message.

Second, the convergence of a nonlinear model depends on the model formula as well as on the data, so the fact that the process converges on one dataset but not another is completely irrelevant.

Third, nonlinear modeling iteratively minimizes the residual sum of squares depending on the parameters. If there are local minima on the RSS surface , and yours start

is close to one, the algorithms will find it. But the only true solution is the global minimum . Your problem has many, many local lows.

Fourth, it nls(...)

uses the default Gaussian Newton method. Gauss Newton, known to be unstable with shearing parameters (parameters that are added to or subtracted from the predictor, therefore Tmin

, and Tmax

in your case). Fortunately, the package minpak.lm

implements the Levenberg Marquardt method, which is much more stable under these conditions. The function nlsLM(...)

in this package uses the same call sequence nls(...)

as and returns an object of type nls

, so all methods for this object class also work. Use this.

Fifth, the fundamental assumption in nonlinear regression (in fact, least squares regression) is that the residuals are usually distributed. Therefore, you need to test some solution using the QQ graph.

Sixth, your model has a perverse set of characteristics. When Tmin -> -Inf

, the first term in the model approaches 1

. It turns out that this gives a lower RSS than any other value Tmin

less than min(dat$x)

, so algorithms tend to lead Tmin

to large negative values. You can easily see it like this:

library(minpack.lm)
mod <- nlsLM(y ~ beta.reg(x, Yopt,Tmin,Topt,Tmax, b1), data=dat,
             start=c(Yopt=6,Tmin=0,Topt=24,Tmax=50, b1=1),
             control=nls.lm.control(maxiter=1024,maxfev=1024))
coef(summary(mod))
#         Estimate   Std. Error     t value     Pr(>|t|)
# Yopt    6.347019    0.2919686 21.73870235 8.055342e-25
# Tmin -155.530098 2204.0011003 -0.07056716 9.440694e-01
# Topt   21.157545    0.6702713 31.56564484 2.240134e-31
# Tmax   35.000000   11.4838614  3.04775537 3.933164e-03
# b1      3.321326    9.1844548  0.36162468 7.194035e-01
sum(residuals(mod)^2)
# [1] 50.24696

par(mfrow=c(1,2))
plot(y~x,dat)
with(as.list(coef(mod)),curve(beta.reg(x, Yopt,Tmin,Topt,Tmax, b1),add=TRUE))
qqnorm(residuals(mod))

      

        
        
        
      

    

It looks like a pretty decent fit , but it is not : the QQ graph shows that the residuals are not remotely normal. The fact that both Tmin

is b1

very poorly evaluated and the value for is Tmin

physically meaningless is data issues, not compliance.

Seventh, it turns out that the fit above is actually a local minimum . We can see this by searching for the grid on Tmin

, Tmax

and b1

(leaving Yopt

and Topt

to save time and because these parameters are well estimated regardless of the starting point).

init <- c(Yopt=6, Topt=24)
grid <- expand.grid(Tmin= seq(0,4,len=100),
                    Tmax= seq(35,100,len=10),
                    b1  = seq(1,10,len=10))
mod.lst <- apply(grid,1,function(gr){
  nlsLM(y ~ beta.reg(x, Yopt,Tmin,Topt,Tmax, b1), data=dat,
        start=c(init,gr),control=nls.control(maxiter=800)) })
rss <- sapply(mod.lst,function(m)sum(residuals(m)^2))
mod <- mod.lst[[which.min(rss)]]   # fit with lowest RSS
coef(summary(mod))
#        Estimate   Std. Error      t value     Pr(>|t|)
# Yopt   6.389238    0.2534551 25.208557840 2.177168e-27
# Topt  22.636505    0.5605621 40.381798589 7.918438e-36
# Tmin  35.000002  104.6221159  0.334537316 7.396005e-01
# Tmax  36.234602  133.4987344  0.271422809 7.873647e-01
# b1   -41.512912 7552.0298633 -0.005496921 9.956395e-01
sum(residuals(mod)^2)
# [1] 34.24019

plot(y~x,dat)
with(as.list(coef(mod)),curve(beta.reg(x, Yopt,Tmin,Topt,Tmax, b1),add=TRUE))
qqnorm(residuals(mod))

      

        
        
        
      

    

Mathematically, this is clearly a better fit: the RSS is lower and the residuals are much more normally distributed. Again, the fact that the parameters are poorly estimated and not physically significant is a problem with the data (and possibly the model formula), not a fitting process.

All of the above assumes that something is wrong with your model. One problem from a mathematical point of view is that the function is undefined for the x

outside (Tmin,Tmax)

. Since you have the data before x=35

, the fitting algorithm will never yield Tmax < 35

(if it converges). The approach to solving this problem slightly changes your modeling function for a clip to 0 outside this range. (I have no idea if this is legitimately based on the physics of your problem though ...).

beta.reg<-function(x, Yopt,Tmin,Topt,Tmax, b1) {
  ifelse(x>Tmax,0,
    ifelse(x<Tmin,0,
      Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1
  ))
}

      

        
        
        
      

    

Executing the above code with this function gives:

coef(summary(mod))
#         Estimate   Std. Error     t value     Pr(>|t|)
# Yopt   6.1470413   0.21976766   27.970636 3.202940e-29
# Tmin -52.8172658 184.16899439   -0.286787 7.756528e-01
# Topt  23.0777898   0.63750721   36.200045 7.638121e-34
# Tmax  30.0039413   0.02529877 1185.984187 1.038918e-98
# b1     0.5966129   0.32439982    1.839128 7.280793e-02

sum(residuals(mod)^2)
# [1] 28.10144

par(mfrow=c(1,2))
plot(y~x,dat)
with(as.list(coef(mod)),curve(beta.reg(x, Yopt,Tmin,Topt,Tmax, b1),add=TRUE))
qqnorm(residuals(mod))
qqline(residuals(mod))

      

        
        
        
      

    

In fact, searching for a grid gives exactly the same result regardless of the starting point. Note that the RSS is lower than any of the results with the earlier model, and that b1

it scores much better (and is very different from the score with the earlier model function). The residuals are still not normal, but in this case I would like to check the data for outliers.