Здесь так много проблем, что я сомневаюсь, что они могут быть адекватно освещены в сообщении SO, но это должно помочь вам начать.

Во-первых, похоже, что вы хотите Tmax < max(dat$x), например ‹35. Это вызывает проблему, потому что тогда Tmax - x < 0 для некоторых значений x и когда вы пытаетесь возвести отрицательное число в степень (во втором члене вашей формулы), вы получаете NA's. Это причина сообщения об ошибке.

Во-вторых, сходимость нелинейной модели зависит от формулы модели, а также от данных, поэтому тот факт, что процесс сходится с одним набором данных, но не с другим, совершенно не имеет значения.

В-третьих, нелинейное моделирование итеративно минимизирует остаточную сумму квадратов как функцию параметров. Если поверхность RSS имеет локальные минимумы, а ваш start близок к единице, алгоритмы его найдут. Но только глобальный минимум является верным решением. У вашей проблемы много-много локальных минимумов.

В-четвертых, nls(...) по умолчанию использует метод Гаусса-Ньютона. Гаусс Ньютон, как известно, нестабилен со смещением параметров (параметры, которые добавляются или вычитаются из предиктора, поэтому Tmin и Tmax в вашем случае). К счастью, в пакете minpak.lm реализован метод Левенберга-Марквардта, который в этих условиях гораздо более устойчив. Функция nlsLM(...) в этом пакете использует ту же последовательность вызовов, что и nls(...), и возвращает объект типа nls, поэтому все методы для этого класса объекта также работают. Используйте это.

В-пятых, фундаментальное предположение нелинейной регрессии (фактически, все регрессии методом наименьших квадратов) состоит в том, что остатки распределены нормально. Таким образом, вы должны проверить любое решение, используя график Q-Q.

В-шестых, ваша модель обладает извращенным набором характеристик. Когда Tmin -> -Inf первый член модели приближается к 1. Оказывается, это дает более низкий RSS, чем любое другое значение Tmin меньше min(dat$x), поэтому все алгоритмы имеют тенденцию приводить Tmin к большим отрицательным значениям. В этом легко убедиться:

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))

Это выглядит довольно неплохо, но это не так: график Q-Q показывает, что остатки далеко не нормальны. Тот факт, что и Tmin, и b1 очень плохо оценены, а значение Tmin не имеет физического смысла, это проблемы с данными, а не с соответствием.

В-седьмых, оказывается, что приведенное выше соответствие на самом деле является локальным минимумом. Мы можем убедиться в этом, выполнив поиск по сетке по Tmin, Tmax и b1 (исключая Yopt и Topt, чтобы сэкономить время, и потому, что эти параметры хорошо оцениваются независимо от начальной точки).

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))

С математической точки зрения это явно превосходное соответствие: RSS ниже, а остатки распределены гораздо ближе к нормальному. Опять же, тот факт, что параметры плохо оценены и не имеют физического смысла, является проблемой с данными (и, возможно, с формулой модели), а не с процессом подбора.

Все вышесказанное говорит о том, что с вашей моделью что-то не так. С математической точки зрения одна проблема заключается в том, что функция не определена для x за пределами (Tmin,Tmax). Поскольку у вас есть данные до x=35, алгоритм подгонки никогда не даст Tmax < 35 (если он сходится). Подход, чтобы справиться с этим, немного изменяет функцию вашей модели, чтобы обрезать до 0 за пределами этого диапазона. (Хотя я понятия не имею, законно ли это, исходя из физики вашей проблемы ...).

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
  ))
}

Выполнение приведенного выше кода с этой функцией дает:

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))

Фактически поиск по сетке дает точно такой же результат независимо от начальной точки. Обратите внимание, что RSS ниже, чем любой из результатов с более ранней моделью, и что b1 гораздо лучше оценен (и очень отличается от оценки с более ранней функцией модели). Остатки по-прежнему ненормальные, но в этом случае я хотел бы проверить данные на наличие выбросов.

Добавление еще одного возможного решения к @jlhoward ...

Я нашел этот nls2 пакет:

library("nls2")

Исключение x~17,35 из исходного набора данных:

newdat <- subset(dat, x!=17 & x!=35 )

Применение функции к сокращенному набору данных:

beta.reg<-with(newdat,  
           y ~ Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / Tmax-Topt))^b1
           )

Создание набора стартеров:

st1 <- expand.grid(Yopt = seq(4, 8, len = 4),
                   Tmin = seq(0, 4, len = 4), 
                   Topt = seq(15, 25, len = 4),
                   Tmax= seq(28, 38, len = 4),
                   b1 = seq(0, 4, len = 4))

Примерка модели:

mod <- nls2(beta.reg, start = st1, algorithm = "brute-force")

Коэффициенты извлечения:

round(coef(summary(mod)),3)

#     Estimate Std. Error t value Pr(>|t|)
# Yopt    6.667      0.394  16.925    0.000
# Tmin    0.000     12.023   0.000    1.000
# Topt   21.667      0.746  29.032    0.000
# Tmax   31.333      1.924  16.289    0.000
# b1      1.333      1.010   1.320    0.197

Диагностика:

sum(residuals(mod)^2)

# [1] 50.18246

И, наконец, скорректированная функция и график QQ-normal:

par(mfrow=c(1,2))
with(newdat,plot(y~x,xlim=c(0,35))) 
points(fitted(mod)~I(newdat$x), pch=19)
with(as.list(coef(mod)),
 curve(
  Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1,
   add=TRUE, col="red"))

qqnorm(residuals(mod))
qqline(residuals(mod))