R: Не удалось подобрать 4-параметрическую кривую хоккейной клюшки с помощью nls

Мой набор данных:

mydata<-structure(list(t = c(0.208333333, 0.208333333, 0.208333333, 0.208333333, 
1, 1, 1, 1, 2, 2, 2, 2, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 
16, 16, 0.208333333, 0.208333333, 0.208333333, 0.208333333, 1, 
1, 1, 1, 2, 2, 2, 2), parent = c(1.2, 1.4, 0.53, 1.2, 1, 0.72, 
0.93, 1.1, 0.88, 0.38, 0.45, 0.27, 0.057, 0.031, 0.025, 0.051, 
0.027, 0.015, 0.034, 0.019, 0.017, 0.025, 0.024, 0.023, 0.29, 
0.22, 0.34, 0.19, 0.12, 0.092, 0.41, 0.28, 0.064, 0.05, 0.058, 
0.043)), .Names = c("t", "Ct"), row.names = c(325L, 326L, 
327L, 328L, 341L, 342L, 343L, 344L, 357L, 358L, 359L, 360L, 373L, 
374L, 375L, 376L, 389L, 390L, 391L, 392L, 401L, 402L, 403L, 404L, 
805L, 806L, 807L, 808L, 821L, 822L, 823L, 824L, 837L, 838L, 839L, 
840L), class = "data.frame")

Подлежащая подгонке функция - кривая клюшки; т.е. после точки изгиба tb она расплющивается:

hockeystick<-function (t, C0, k1, k2, tb) 
{
  Ct = ifelse(t <= tb, C0 -k1 * t, C0 -k1*tb -k2*t)
}

Монтаж с использованием nls:

start.hockey<-c(C0=3,k1=1,k2=0.1,tb=3)
nls(log(Ct)~hockeystick(t,C0,k1,k2,tb),start=start.hockey,data=mydata)

Независимо от того, какие начальные значения я использую, я всегда получаю эту ошибку:

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

Я пробовал как port, так и стандартные методы nls. Я пробовал как линеаризованное (показано здесь), так и нормальное состояние модели, но, похоже, ни один из них не работает.

Изменить: в соответствии с предложением Карла я попытался подогнать модель к набору данных, где я сначала усреднил значения Ct на значение t и все еще получаю ошибку.

edit: Немного изменена модель, поэтому значение k2 положительное, а не отрицательное. Отрицательное значение кинетически не имеет смысла.


r curve-fitting nls
person Pinemangoes    schedule 11.09.2014    source источник
comment
Попробуйте построить график hockeystick(mydata$t,C0,k1,k2,tb) против mydata$t. Это не хоккейная клюшка. Кроме того, ваши повторяющиеся значения t почти гарантированно приведут к сбою регрессии.   -  person Carl Witthoft    schedule 11.09.2014
comment
В общем, мне было бы лучше подогнать регрессию к усредненным значениям на значение t? Хоккейная клюшка имеет линеаризованную форму, поэтому ось Y находится в логарифмических единицах.   -  person Pinemangoes    schedule 11.09.2014
comment
tb - «точка изгиба» модели хоккейной клюшки, точка, в которой кривая изменяет «скорость спуска».   -  person Pinemangoes    schedule 11.09.2014

Ответы (1)


arrow_upward
2
arrow_downward

Я еще не решил проблему nls(), но у меня есть несколько предложений.

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

hockeystick<-function (t, C0, k1, k2, tb) 
{
   Ct <- ifelse(t <= tb, C0 -k1 * t, C0 -k1*t -k2*(t-tb))
}

Глазное яблоко:

par(las=1,bty="l") ## cosmetic
plot(log(Ct)~t,data=mydata)
curve(hockeystick(x,C0=0,k1=0.8,k2=-0.7, tb=3),add=TRUE)

введите описание изображения здесь

Я сделал k2 отрицательным, поэтому наклон на втором этапе меньше, чем на первом.

start.hockey <- c(C0=0,k1=0.8,k2=-0.7, tb=3)
nls(log(Ct)~hockeystick(t,C0,k1,k2,tb),
                        start=start.hockey,data=mydata)

Модели с точками останова часто не дифференцируются по параметрам, но я не совсем понимаю, в чем здесь проблема ...

Это действительно работает:

library(bbmle)
m1 <- mle2(log(Ct)~dnorm(hockeystick(t,C0,k1,k2,tb),
                  sd=exp(logsd)),
          start=c(as.list(start.hockey),list(logsd=0)),
          data=mydata)

Параметры приемлемые (и отличаются от начальных значений):

coef(summary(m1))
##         Estimate Std. Error   z value        Pr(z)
## C0    -0.4170749  0.2892128 -1.442104 1.492731e-01
## k1     0.6720120  0.2236111  3.005271 2.653439e-03
## k2    -0.5285974  0.2400605 -2.201934 2.766994e-02
## tb     2.0007688  0.1714292 11.671108 1.790751e-31
## logsd -0.2218745  0.1178580 -1.882558 5.976033e-02

Сюжетные прогнозы:

pframe <- data.frame(t=seq(0,15,length=51))
pframe$pred <- predict(m1,newdata=pframe)
with(pframe,lines(t,pred,col=2))

введите описание изображения здесь

person Ben Bolker    schedule 11.09.2014
comment
Спасибо за исчерпывающий ответ. Да, ваша модификация с tb на t-tb - большое улучшение. Я проверю оставшуюся часть вашего ответа, когда вернусь домой сегодня вечером. - person Pinemangoes; 11.09.2014
comment
Не могли бы вы объяснить, почему вы оборачиваете функцию dnorm вокруг hockeystick() для оценки? Кроме того, для этого конкретного соединения мне, вероятно, потребуется получить данные в диапазоне 2-10 для _3 _... - person Pinemangoes; 11.09.2014
comment
Кроме того, похоже, что параметр tb вообще не изменился во время подгонки. Любое значение, которое я использую в hockey.start, сохраняется в последнем наборе параметров. Хотя я обычно могу сделать обоснованное предположение на основе данных, я не могу сделать это в бумаге. - person Pinemangoes; 11.09.2014
comment
параметр tb для тестового набора изменился (см. правки выше). Если вы используете разумные начальные значения, я не знаю, почему оптимизатор зависнет. Я использую dnorm(), потому что mle2 не выполняет метод наименьших квадратов, он выполняет общую оценку максимального правдоподобия, а dnorm() - это способ привести его к эквиваленту аппроксимации методом наименьших квадратов. - person Ben Bolker; 11.09.2014
comment
было бы более эффективно и надежно подогнать кусочно-линейную модель с использованием lm для указанной точки останова, а затем выполнить одномерную нелинейную оптимизацию (optimize) для выбора точки останова. Или используйте пакет strucchange. - person Ben Bolker; 11.09.2014