studio - superponer graficas en r ggplot

Regresión por partes con R: trazando los segmentos

Question

Tengo 54 puntos. Representan oferta y demanda de productos. Me gustaría mostrar que hay un punto de quiebre en la oferta.

Primero, ordeno el eje x (oferta) y elimino los valores que aparecen dos veces. Tengo 47 valores, pero elimino el primero y el último (no tiene sentido considerarlos como puntos de ruptura). La rotura es de longitud 45:

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

Luego, para cada uno de estos puntos de ruptura potenciales, estimo un modelo y mantengo en "d" el error estándar residual (sexto elemento en el objeto de resumen del modelo).

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

Al trazar d, observo que mi error estándar residual más pequeño es 34, que corresponde a "Break [34]": 22.4. Así que escribo mi modelo con mi punto de ruptura final:

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

Por último, estoy feliz con mi nuevo modelo. Es significativamente mejor que el lineal simple. Y quiero dibujarlo:

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

Pero tengo un mensaje de advertencia:

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

Y más importante, las líneas son realmente extrañas en mi trama.

Aquí están mis datos:

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)

Answer 1

Aquí es un enfoque más fácil utilizando ggplot2 .

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

EDITAR. También recomendaría echar un vistazo a este paquete segmented que admite la detección automática y la estimación de modelos de regresión segmentada.

ACTUALIZAR:

Aquí hay un ejemplo que hace uso del paquete R segmented para detectar automáticamente los cortes.

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

Answer 2

Las líneas extrañas se deben simplemente al orden en que se trazan los puntos. Lo siguiente debería verse mejor:

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

La advertencia proviene del hecho de que el carácter * es interpretado por 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 

Además, (offer<22.4)*offer es una función discontinua: de aquí proviene la discontinuidad.

Lo siguiente debe estar más cerca de lo que quieres.

model <- lm(
demand ~ ifelse(offer<22.4,offer-22.4,0)
+ ifelse(offer>=22.4,offer-22.4,0) )

Answer 3

Vincent te tiene en el camino correcto. Lo único "extraño" acerca de las líneas en su gráfica resultante es que las lines trazan una línea entre cada punto sucesivo, lo que significa que "salto" se ve si simplemente se conectan los dos extremos de cada línea.

Si no desea ese conector, debe dividir las lines en dos partes separadas.

Además, siento que puedes simplificar un poco tu regresión. Esto es lo que hice:

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

que produce esta trama: