NLS - ошибка сходимости

Question

Для этого набора данных:

dat = structure(list(x = c(5L, 5L, 5L, 5L, 10L, 10L, 10L, 10L, 15L, 
15L, 15L, 15L, 17L, 17L, 17L, 17L, 20L, 20L, 20L, 20L, 20L, 20L, 
20L, 20L, 22L, 22L, 22L, 22L, 24L, 24L, 24L, 24L, 25L, 25L, 25L, 
25L, 27L, 27L, 27L, 27L, 30L, 30L, 30L, 30L, 35L, 35L, 35L, 35L), 
y = c(2.2, 2.2, 1.95, 1.9, 4.1, 3.95, 3.75, 3.4, 5.15, 4.6, 
4.75, 5.15, 3.7, 4.1, 3.9, 3.5, 7, 6.7, 6.7, 6.95, 4.95, 6, 6.45, 
6.4, 7, 4.45, 6.15, 6.4, 7, 6.6, 6.7, 7, 4.5, 4.7, 5.75, 4.35, 
5.4, 5.15, 5.7, 5.7, 0, 0, 0.5, 0, 0, 0, 0, 0)), .Names = c("x", "y"), 
row.names = c(6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 
15L, 16L, 17L, 34L, 35L, 36L, 37L, 18L, 19L, 20L, 21L, 38L, 39L, 
40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 22L, 23L, 24L, 
25L, 50L, 51L, 52L, 53L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L), 
class = "data.frame") 

Где «х» является температура и «у» переменная отклика биологического процесса

Я пытаюсь приспосабливать эту функцию

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

mod <- nls(y ~ beta.reg(x, Yopt,Tmin,Topt,Tmax, b1), data=dat, 
     start=c(Yopt=6, Tmin=0.1, Topt=24, Tmax=30, b1=1), 
     control=nls.control(maxiter=800)) 

Но, у меня эта ошибка сообщение:

Error en numericDeriv(form[[3L]], names(ind), env) : Missing value or an infinity produced when evaluating the model

Я пробовал ту же функцию с другой такой же набор данных и подходит правильно ...

rnorm<-(10) 
y <- c(20,60,70,49,10) 
rnorm<-(10) 
y <- c(20,60,70,49,10) 
dat<-data.frame(x = rep(c(15,20,25,30,35), times=5), 
       rep = as.factor(rep(1:5, each=5)), 
       y = c(y+rnorm(5), y+rnorm(5),y+rnorm(5),y+rnorm(5),y+rnorm(5))) 

Может кто-нибудь помочь мне с этим?

Session Info:

R version 3.1.1 (2014-07-10) 
Platform: x86_64-pc-linux-gnu (64-bit) 

attached base packages: 
[1] stats  graphics grDevices utils  datasets methods base  

other attached packages: 
[1] nlme_3.1-118  latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29  

loaded via a namespace (and not attached): 
[1] grid_3.1.1 tools_3.1.1 

Answer 1

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

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

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

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

В-четвертых, nls(...) использует метод Гаусса Ньютона по умолчанию. Гаусс Ньютон, как известно, нестабилен со сдвигающими параметрами (параметры, которые добавляются или вычитаются из предиктора, поэтому Tmin и Tmax в вашем случае). К счастью, пакет minpak.lm реализует метод Levenberg Marquardt, который в этих условиях намного более стабилен. Функция 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 намного лучше оценен (и очень отличается от оценки с более ранней модельной функцией). Остатки по-прежнему не нормальны, но в этом случае я хотел бы проверить данные на выбросы.

Answer 2

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

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

library("nls2") 

Exludying 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-нормальный сюжет:

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