Кусочная регрессия с 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, что соответствует "Break [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

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

My plot with the supposedly two segments, but not joining

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

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)

Ответы

Ответ 1

Вот более простой подход, используя ggplot2.

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

ИЗМЕНИТЬ. Я бы также рекомендовал взглянуть на этот пакет segmented, который поддерживает автоматическое обнаружение и оценку сегментированных регрессионных моделей.

enter image description here

UPDATE:

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

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

segmented

Ответ 2

У Винсента вы на правильном пути. Единственное, что "странно" в строках вашего результирующего графика состоит в том, что 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))

который производит этот график:

enter image description here

Ответ 3

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

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