Regresión segmentada en dos trozos

Geyser Old Faithful del parque de Yellowstone.

La regresión segmentada o regresión por pedazos es un método de regresión en que la variable predictora se divide en intervalos de forma que se ajusta una línea o curva a los datos en cada intervalo. La regresión segmentada es útil cuando la variable respuesta muestra una reacción abruptamente diferente a la variable predictora en los diferentes segmentos. En este caso el límite entre los segmentos se llama punto de quiebra. La regresión segmentada se puede aplicar también a la regresión con múltiples variables predictoras. En el documento adjunto se trata únicamente el caso de la regresión lineal simple (una única variable explicativa) segmentada en dos trozos.

Archivo PDF

Regresión local simple

Introducción

La regresión local en el caso simple consiste en predecir la respuesta para un valor particular \(x_0\) de la variable explicativa mediante la recta de regresión ponderada construida con un grupo de puntos de la muestra cuyos valores de la variable explicativa son cercanos a \(x_0\).

En el apartado 2.1.1 del documento Appendix-Nonparametric-Regression.pdf se explica este tipo de regresión no paramétrica de forma clara y con un ejemplo. Aquí vamos a reproducir los gráficos de la figura 1 de este documento de J. Fox & S. Weisberg que es un apéndice del libro An R Companion to Applied Regression.

Los datos

Los datos para la regresión se hallan en el data.frame Prestige del paquete carData. Como ejemplo de regresión local se toman las variables income como regresora y prestige como respuesta.

library(carData)
data("Prestige")
str(Prestige)

## 'data.frame':    102 obs. of  6 variables:
##  $ education: num  13.1 12.3 12.8 11.4 14.6 ...
##  $ income   : int  12351 25879 9271 8865 8403 11030 8258 14163 11377 11023 ...
##  $ women    : num  11.16 4.02 15.7 9.11 11.68 ...
##  $ prestige : num  68.8 69.1 63.4 56.8 73.5 77.6 72.6 78.1 73.1 68.8 ...
##  $ census   : int  1113 1130 1171 1175 2111 2113 2133 2141 2143 2153 ...
##  $ type     : Factor w/ 3 levels "bc","prof","wc": 2 2 2 2 2 2 2 2 2 2 ...

attach(Prestige)
plot(income,prestige,
           xlab="Average income (dollars)",ylab="Prestige")

Para entender como funciona la regresión local en este ejemplo vamos a considerar la predicción de un valor concreto. Exactamente el valor de income que ocupa la posición octogésima que designaremos como \(x_{(80)}\).

# case with the 80th largest income value
rownames(Prestige)[rank(income)==80]

## [1] "chemists"

x_80 <- Prestige["chemists","income"]

Gráfico (a)

En este gráfico se muestran las dos rectas verticales que delimitan los 50 vecinos más próximos a \(x_{(80)}\).

La recta vertical de tipo sólido se corresponde con dicho valor \(x_{(80)}\) y las rectas verticales con guiones son los límites para los 50 vecinos de ese valor. Se ha tomado un conjunto de 50 vecinos ya que \(50/102\approx 1/2\), donde \(102\) es el número total de observaciones en la muestra.

dfr <- data.frame(income, prestige, dist=abs(income-x_80))
dfro <- dfr[order(dfr$dist),]
plot(income, prestige, col="grey",
         xlab="Average income (dollars)", ylab="Prestige")
points(dfro$income[2:51], dfro$prestige[2:51],pch=16)
abline(v=x_80)
abline(v=max(dfro$income[2:51]),lty=2)
abline(v=min(dfro$income[2:51]),lty=2)
axis(3, at=x_80,expression(x[(80)]), tck=0)

Gráfico (b)

Ahora se dibuja la función tricube con los pesos de los 50 vecinos más próximos a \(x_{(80)}\).

tricube <- function(z) ifelse(abs(z)<1, (1-abs(z)^3)^3, 0)
h <- (max(dfro$income[2:51])-min(dfro$income[2:51]))/2
curve(tricube((x-x_80)/h),0,max(income),
                    xlab="Average income (dollars)",
                    ylab="Tricube Kernel Weight")
abline(h=0, col="grey")
abline(h=1, col="grey")
abline(v=max(dfro$income[2:51]),lty=2)
abline(v=min(dfro$income[2:51]),lty=2)
points(dfro$income[2:51],tricube((dfro$income[2:51]-x_80)/h), pch=16)

Gráfico (c)

Aquí se dibuja la regresión lineal local ponderada con los vecinos de \(x_{(80)}\). El punto gordo de color magenta es la predicción de \(x_{(80)}\) por esta regresión.

fit <- lm(prestige ~ income, weights = tricube((income-x_80)/h))
library(shape)
plot(income,prestige,col="grey",
     xlab="Average income (dollars)",ylab="Prestige")
points(dfro$income[2:51],dfro$prestige[2:51],pch=16)
abline(v=x_80)
abline(v=max(dfro$income[2:51]),lty=2)
abline(v=min(dfro$income[2:51]),lty=2)
x0 <- min(dfro$income[2:51])
y0 <- predict(fit, newdata=data.frame(income=min(dfro$income[2:51])))
x1 <- max(dfro$income[2:51])
y1 <- predict(fit, newdata=data.frame(income=max(dfro$income[2:51])))
segments(x0,y0,x1,y1,col="magenta",lwd=3)
points(x=income[rownames(Prestige)=="chemists"], # x_80
       y=predict(fit)[rownames(Prestige)=="chemists"],
       col="magenta",pch=20,cex=3)
x0 <- 18000
y0 <- 45
x1 <- x_80
y1 <- predict(fit)[rownames(Prestige)=="chemists"]  
Arrows(x0,y0,x1,y1,arr.type="triangle",arr.adj = 1)
text(x=20000,y=45,labels = expression(hat(y)[(80)]))

Gráfico (d)

Finalmente, este gráfico muestra la regresión local completa conectando los valores ajustados para todos los datos de \(x\).

plot(prestige ~ income, xlab="Average Income", ylab="Prestige")
lines(lowess(income, prestige, f=0.5, iter=0), lwd=2)

Bibliografía

John Fox & Sanford Weisberg (2018), An R Companion to Applied Regression, Third Edition, SAGE Publications.