Анализ точек изменения с помощью R's nls ()

16

Я пытаюсь реализовать анализ "точки изменения" или многофазную регрессию с использованием nls()R.

Вот некоторые фальшивые данные, которые я сделал . Формула, которую я хочу использовать, чтобы соответствовать данным:

Yзнак равноβ0+β1Икс+β2Максимум(0,Икс-δ)

Предполагается, что это должно соответствовать данным до определенной точки с определенным пересечением и наклоном ( \ beta_0β0 и β1 ), а затем, после определенного значения x ( δ ), увеличить наклон на β2 . Вот о чем вся эта вещь. Перед точкой δ она будет равна 0, а β2 будет обнулен.

Итак, вот моя функция, чтобы сделать это:

changePoint <- function(x, b0, slope1, slope2, delta){ 
   b0 + (x*slope1) + (max(0, x-delta) * slope2)
}

И я пытаюсь подобрать модель таким образом

nls(y ~ changePoint(x, b0, slope1, slope2, delta), 
    data = data, 
    start = c(b0 = 50, slope1 = 0, slope2 = 2, delta = 48))

Я выбрал эти начальные параметры, потому что я знаю, что это начальные параметры, потому что я составил данные.

Тем не менее, я получаю эту ошибку:

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

Я только что сделал неудачные данные? Сначала я попытался подстроить это под реальные данные, и получил ту же ошибку, и я просто подумал, что мои начальные начальные параметры не были достаточно хорошими.

JoFrhwld
источник

Ответы:

12

(Сначала я подумал, что это может быть проблемой из-за того, что maxэто не векторизация, но это не так. Из-за этого трудно работать с changePoint, поэтому следует внести следующие изменения:

changePoint <- function(x, b0, slope1, slope2, delta) { 
   b0 + (x*slope1) + (sapply(x-delta, function (t) max(0, t)) * slope2)
}

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

В любом случае вы можете написать свою целевую функцию и минимизировать ее. Следующая функция выдает квадратичную ошибку для точек данных (x, y) и определенного значения параметров (странная структура аргументов функции должна учитывать, как optimработает):

sqerror <- function (par, x, y) {
  sum((y - changePoint(x, par[1], par[2], par[3], par[4]))^2)
}

Тогда мы говорим:

optim(par = c(50, 0, 2, 48), fn = sqerror, x = x, y = data)

И посмотреть:

$par
[1] 54.53436800 -0.09283594  2.07356459 48.00000006

Обратите внимание, что для моих поддельных данных ( x <- 40:60; data <- changePoint(x, 50, 0, 2, 48) + rnorm(21, 0, 0.5)) существует множество локальных максимумов, в зависимости от заданных вами значений начальных параметров. Я полагаю, что если вы хотите отнестись к этому серьезно, вы бы много раз вызывали оптимизатор со случайными начальными параметрами и проверяли распределение результатов.

Аарон
источник
Этот пост Билла Венаблса хорошо объясняет проблемы, связанные с таким анализом.
Аарон
6
Вместо этого (громоздкого) вызова sapply в вашем первом фрагменте кода вы всегда можете просто использовать pmax .
кардинал
0

Просто хотел добавить, что вы можете сделать это со многими другими пакетами. Если вы хотите получить оценку неопределенности относительно точки изменения (чего не может сделать nls), попробуйте mcpпакет.

# Simulate the data
df = data.frame(x = 1:100)
df$y = c(rnorm(20, 50, 5), rnorm(80, 50 + 1.5*(df$x[21:100] - 20), 5))

# Fit the model
model = list(
  y ~ 1,  # Intercept
  ~ 0 + x  # Joined slope
)
library(mcp)
fit = mcp(model, df)

Давайте построим это с интервалом предсказания (зеленая линия). Синяя плотность - это заднее распределение местоположения точки изменения:

# Plot it
plot(fit, q_predict = T)

Вы можете проверить отдельные параметры более подробно, используя plot_pars(fit)и summary(fit).

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

Йонас Линделёв
источник