Кусочная регрессия с R: построение отрезков

У меня 54 очка. Они представляют предложение и спрос на продукцию. Я хотел бы показать, что в предложении есть перерыв.

Сначала я сортирую ось x (предложение) и удаляю значения, которые появляются дважды. У меня 47 значений, но я убираю первое и последнее (не имеет смысла рассматривать их как точки останова). Перерыв длиной 45:

Break<-(sort(unique(offer))[2:46])

Затем для каждой из этих потенциальных точек разрыва я оцениваю модель и сохраняю в «d» остаточную стандартную ошибку (шестой элемент в объекте сводки модели).

d<-numeric(45)
for (i in 1:45) {
model<-lm(demand~(offer<Break[i])*offer + (offer>=Break[i])*offer)
d[i]<-summary(model)[[6]] }

Рисуя d, я замечаю, что моя меньшая стандартная остаточная ошибка составляет 34, что соответствует «Разрыву [34]»: 22.4. Итак, я пишу свою модель с моей последней точкой останова:

model<-lm(demand~(offer<22.4)*offer + (offer>=22.4)*offer)

Наконец-то я доволен своей новой моделью. Он значительно лучше простого линейного. И я хочу это нарисовать:

plot(demand~offer)
i <- order(offer)
lines(offer[i], predict(model,list(offer))[i])

Но у меня есть предупреждающее сообщение:

Warning message:
In predict.lm(model, list(offer)) :
  prediction from a rank-deficient fit may be misleading

И что еще важнее, на моем сюжете линии действительно странные.

Мой сюжет с предположительно двумя сегментами, но не соединяющимися

Вот мои данные:

demand <- c(1155, 362, 357, 111, 703, 494, 410, 63, 616, 468, 973, 235,
            180, 69, 305, 106, 155, 422, 44, 1008, 225, 321, 1001, 531, 143,
            251, 216, 57, 146, 226, 169, 32, 75, 102, 4, 68, 102, 462, 295,
            196, 50, 739, 287, 226, 706, 127, 85, 234, 153, 4, 373, 54, 81,
            18)
offer <- c(39.3, 23.5, 22.4, 6.1, 35.9, 35.5, 23.2, 9.1, 27.5, 28.6, 41.3,
           16.9, 18.2, 9, 28.6, 12.7, 11.8, 27.9, 21.6, 45.9, 11.4, 16.6,
           40.7, 22.4, 17.4, 14.3, 14.6, 6.6, 10.6, 14.3, 3.4, 5.1, 4.1,
           4.1, 1.7, 7.5, 7.8, 22.6, 8.6, 7.7, 7.8, 34.7, 15.6, 18.5, 35,
           16.5, 11.3, 7.7, 14.8, 2, 12.4, 9.2, 11.8, 3.9)

r plot piecewise
person Antonin    schedule 06.01.2012    source источник
comment
54 балла - не очень большое количество баллов для обнаружения такого перехода. Вы также можете разместить это на stats.stackexchange.com с конкретным вопросом, достаточно ли очков для обнаружения разрыва данных. Просто мои 2 кар.   -  person Paul Hiemstra    schedule 07.01.2012
comment
Я считаю, что это довольно сомнительно со статистической точки зрения. Лучше оценивать точку останова в самой модели (хотя это делает ее нелинейной). Вы не можете доверять p-значениям или стандартным ошибкам текущего неформального процесса оценки.   -  person hadley    schedule 08.01.2012
comment
Я согласен, 54 балла - это не большая сумма, но моя линейная и кусочно-линейная регрессии значительны. Более того, остаточная стандартная ошибка кусочно-линейной модели составляет 103,9 по сравнению с 121,3 для линейной модели с двумя меньшими степенями свободы. Кусочная модель значительно лучше.   -  person Antonin    schedule 09.01.2012

Ответы (3)


arrow_upward
29
arrow_downward

Вот более простой подход с использованием ggplot2.

require(ggplot2)
qplot(offer, demand, group = offer > 22.4, geom = c('point', 'smooth'), 
   method = 'lm', se = F, data = dat)

РЕДАКТИРОВАТЬ. Я также рекомендовал бы взглянуть на этот пакет segmented, который поддерживает автоматическое обнаружение и оценку моделей сегментированной регрессии.

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

ОБНОВИТЬ:

Вот пример, который использует сегментированный пакет R для автоматически обнаруживать перерывы

library(segmented)
set.seed(12)
xx <- 1:100
zz <- runif(100)
yy <- 2 + 1.5*pmax(xx - 35, 0) - 1.5*pmax(xx - 70, 0) + 15*pmax(zz - .5, 0) + 
  rnorm(100,0,2)
dati <- data.frame(x = xx, y = yy, z = zz)
out.lm <- lm(y ~ x, data = dati)
o <- segmented(out.lm, seg.Z = ~x, psi = list(x = c(30,60)),
  control = seg.control(display = FALSE)
)
dat2 = data.frame(x = xx, y = broken.line(o)$fit)

library(ggplot2)
ggplot(dati, aes(x = x, y = y)) +
  geom_point() +
  geom_line(data = dat2, color = 'blue')

сегментированный

person Ramnath    schedule 06.01.2012
comment
Спасибо за идею использования сегментированного пакета. Muggeo, V.M.R. (2003) Оценка регрессионных моделей с неизвестными точками разрыва. Статистика в Medicine 22, 3055–3071 - интересная статья для понимания того, что происходит в пакете. - person Antonin; 10.01.2012
comment
В частности, у него есть преимущество перед кодом, который я использовал: два сегмента соединены! В The R Book автор не упоминает, что его сегменты не связаны, и даже показывает сюжет со связанными сегментами ... - person Antonin; 10.01.2012
comment
Где пример использования ggplot () с segmented ()? Кажется, я нигде не могу найти. - person Adam Erickson; 19.09.2015
comment
Я добавил пример, в котором ggplot используется с segmented. - person Ramnath; 20.09.2015
comment
После qplot(…) я получаю Error: Unknown parameters: method, se - person sam; 24.03.2016
comment
Подходит ли ваше решение ЗДЕСЬ? - person rnorouzian; 22.11.2020

arrow_upward
7
arrow_downward

Винсент ведет вас на верный путь. Единственное, что «странно» в линиях в вашем результирующем графике, это то, что lines рисует линию между каждой последовательной точкой, что означает, что вы видите «прыжок», если он просто соединяет два конца каждой линии.

Если вам не нужен этот соединитель, вам нужно разделить вызов lines на две отдельные части.

Кроме того, я чувствую, что вы можете немного упростить свою регрессию. Вот что я сделал:

#After reading your data into dat
Break <- 22.4
dat$grp <- dat$offer < Break

#Note the addition of the grp variable makes this a bit easier to read
m <- lm(demand~offer*grp,data = dat)
dat$pred <- predict(m)

plot(dat$offer,dat$demand)
dat <- dat[order(dat$offer),]
with(subset(dat,offer < Break),lines(offer,pred))
with(subset(dat,offer >= Break),lines(offer,pred))

который производит этот сюжет:

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

person joran    schedule 06.01.2012

arrow_upward
4
arrow_downward

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

i <- order(offer)
lines(offer[i], predict(model,list(offer))[i])

Предупреждение возникает из-за того, что символ * интерпретируется lm.

> lm(demand~(offer<22.4)*offer + (offer>=22.4)*offer)
Call:
lm(formula = demand ~ (offer < 22.4) * offer + (offer >= 22.4) * offer)
Coefficients:
            (Intercept)         offer < 22.4TRUE                    offer  
                -309.46                   356.08                    29.86  
      offer >= 22.4TRUE   offer < 22.4TRUE:offer  offer:offer >= 22.4TRUE  
                     NA                   -20.79                       NA

Кроме того, (offer<22.4)*offer является прерывистой функцией: отсюда и происходит разрыв.

Следующее должно быть ближе к тому, что вы хотите.

model <- lm(
  demand ~ ifelse(offer<22.4,offer-22.4,0) 
           + ifelse(offer>=22.4,offer-22.4,0) )
person Vincent Zoonekynd    schedule 06.01.2012
comment
Спасибо за ответ! У меня все еще есть предупреждающее сообщение: Предупреждающее сообщение: В predic.lm (model2, list (offer)): прогноз на основе соответствия с недостаточным рангом может вводить в заблуждение. Результат немного лучше: у меня только один кусочный сегмент, но я не понимаю, почему есть 3 сегмента (т.е. мои два сегмента не связаны естественным образом ...) - person Antonin; 06.01.2012
comment
Не могли бы вы объяснить, почему это делает его плавным? Интуитивно подумал, что в первой части функции будет 22,4-предложение. И предложение-22,4 во второй части. Разве первая часть не означает, что если предложение меньше 22,4, то рассчитывается предложение-22,4, что означает отрицательную ковариату? Следовательно, наш коэффициент будет умножен на эти отрицательные значения, когда мы оценим ввод? - person chris; 19.12.2018