Синтаксис трехэлементной сегментированной регрессии с использованием NLS в R при вогнутой

Question

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

Я использовал пакет segmented для оценки местоположений точек останова (используемых в качестве начальных значений в NLS), но хотел бы сохранить свои модели в формате NLS, в частности, nlsLM {minipack.lm}, потому что я использую другие типы кривые к моим данным с использованием NLS, хочу позволить NLS оптимизировать значения узлов, я иногда использую переменные веса и мне нужно иметь возможность легко вычислять доверительные интервалы Монте-Карло из propagate. Хотя я очень близок к тому, чтобы иметь правильный синтаксис для формулы, я не получаю ожидаемого / необходимого поведения около точек останова. Сегменты ДОЛЖНЫ встречаться непосредственно в точках останова (без каких-либо скачков), но, по крайней мере, на этих данных я получаю странный локальный минимум в точке останова (см. Графики ниже).

Ниже приведен пример моего данные и общий процесс. Я считаю, что моя проблема заключается в формуле NLS.

library(minpack.lm)
library(segmented)

y <- c(-3.99448113, -3.82447011, -3.65447803, -3.48447030, -3.31447855, -3.14448753, -2.97447972, -2.80448401, -2.63448380, -2.46448069, -2.29448796, -2.12448912, -1.95448783, -1.78448797, -1.61448563, -1.44448719, -1.27448469, -1.10448651, -0.93448525, -0.76448637, -0.59448626, -0.42448586, -0.25448588, -0.08448548,  0.08551417,  0.25551393,  0.42551411,  0.59551395,  0.76551389,  0.93551398)

x <- c(61586.1711, 60330.5550, 54219.9925, 50927.5381, 48402.8700, 45661.9175, 37375.6023, 33249.1248, 30808.6131, 28378.6508, 22533.3782, 13901.0882, 11716.5669, 11004.7305, 10340.3429,  9587.7994,  8736.3200,  8372.1482,  8074.3709,  7788.1847,  7499.6721,  7204.3168,  6870.8192,  6413.0828,  5523.8097,  3961.6114,  3460.0913,  2907.8614, 2016.1158,   452.8841)


df<- data.frame(x,y)


#Use Segmented to get estimates for parameters with 2 breakpoints
my.seg2 <- segmented(lm(y ~ x, data = df), seg.Z = ~ x, npsi = 2)


#extract knot, intercept, and coefficient values to use as NLS start points
my.knot1 <- my.seg2$psi[1,2]
my.knot2 <- my.seg2$psi[2,2]
my.m_2 <- slope(my.seg2)$x[1,1]
my.b1 <- my.seg2$coefficients[[1]]
my.b2 <- my.seg2$coefficients[[2]]
my.b3 <- my.seg2$coefficients[[3]]

#Fit a NLS model to ~replicate segmented model. Presumably my model formula is where the problem lies
my.model <- nlsLM(y~m*x+b+(b2*(ifelse(x>=knot1&x<=knot2,1,0)*(x-knot1))+(b3*ifelse(x>knot2,1,0)*(x-knot2-knot1))),data=df, start = c(m = my.m_2, b = my.b1, b2 = my.b2, b3 = my.b3, knot1 = my.knot1, knot2 = my.knot2))

Как это должно выглядеть

plot(my.seg2)

How it does look

plot(x, y)
lines(x=x, y=predict(my.model), col='black', lty = 1, lwd = 1)

I was pretty sure I had it "right", but when the 95% confidence intervals are plotted with the line and prediction resolution (e.g., the density of x points) is increased, things seem dramatically incorrect.

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

Всем спасибо за помощь.

Answer 1

Определите g как вектор группировки, имеющий ту же длину, что и x, который принимает значения 1, 2, 3 для 3 частей оси X, и создайте из них модель nls. Полученный график выглядит нормально.

my.knots <- c(my.knot1, my.knot2)
g <- cut(x, c(-Inf, my.knots, Inf), label = FALSE)
fm <- nls(y ~ a[g] + b[g] * x, df, start = list(a = c(1, 1, 1), b = c(1, 1, 1)))

plot(y ~ x, df)
lines(fitted(fm) ~ x, df, col = "red")

(продолжение после графика)

Ограничения

Хотя приведенное выше выглядит нормально и может быть Достаточно, это не гарантирует, что отрезки пересекаются в узлах. Для этого мы должны наложить ограничения, что обе стороны равны в узлах:

a[2] + b[2] * my.knots[1] = a[1] + b[1] * my.knots[1]
a[3] + b[3] * my.knots[2] = a[2] + b[2] * my.knots[2]

поэтому

a[2] = a[1] + (b[1] - b[2]) * my.knots[1]
a[3] = a[2] + (b[2] - b[3]) * my.knots[2]
     = a[1] + (b[1] - b[2]) * my.knots[1] + (b[2] - b[3]) * my.knots[2]

давая:

# returns a vector of the three a values
avals <- function(a1, b) unname(cumsum(c(a1, -diff(b) * my.knots)))

fm2 <- nls(y ~ avals(a1, b)[g] + b[g] * x, df, start = list(a1 = 1, b = c(1, 1, 1)))

Чтобы получить мы можем использовать три значения a:

co <- coef(fm2)
avals(co[1], co[-1])

Чтобы получить остаточную сумму квадратов:

deviance(fm2)
## [1] 0.193077

Многочлен

Хотя он включает в себя большое количество параметров, полиномиальная аппроксимация может использоваться вместо сегментированной линейной регрессии. Полином 12-й степени включает 13 параметров, но имеет меньшую остаточную сумму квадратов, чем сегментированная линейная регрессия. Более низкая степень может использоваться с соответствующим увеличением остаточной суммы квадратов. Полином 7-й степени включает 8 параметров и визуально выглядит неплохо, хотя имеет более высокую остаточную сумму квадратов.

fm12 <- nls(y ~ cbind(1, poly(x, 12)) %*% b, df, start = list(b = rep(1, 13)))

deviance(fm12)
## [1] 0.1899218

Answer 2

Это может частично отражать ограничение в segmented. segmented возвращает значение одной точки изменения без количественной оценки связанной неопределенности. Повторяя анализ с использованием mcp, который возвращает байесовские апостериорные значения, мы видим, что вторая точка изменения распределяется бимодально:

library(mcp)
model = list(
  y ~ 1 + x,  # Intercept + slope in first segment
  ~ 0 + x,  # Only slope changes in the next segments
  ~ 0 + x
)

# Fit it with a large number of samples and plot the change point posteriors
fit = mcp(model, data = data.frame(x, y), iter = 50000, adapt = 10000)
plot_pars(fit, regex_pars = "^cp*", type = "dens_overlay")

FYI, mcp также может отображать достоверные интервалы (красные пунктирные линии):

plot(fit, q_fit = TRUE)