Regresión Cuantílica (Quantile Regression)
Joaquín Amat Rodrigo
Julio 2017
Versión PDF: Github
Más sobre ciencia de datos: cienciadedatos.net
Introducción
La regresión lineal por mínimos cuadrados, así como muchas de sus adaptaciones, se emplean como método para estimar la media de una variable respuesta condicionada a uno o varios predictores. Es decir, parte de la premisa de que la media de la variable respuesta depende del valor que tomen otras variables. En determinadas situaciones, puede ocurrir que la media no sea el parámetro más informativo y que en su lugar se desee conocer un determinado cuantil, por ejemplo el cuantil 0.5 (la mediana). Es aquí donde la regresión de cuantiles (quantile regression) interviene. Si bien la matemática es diferente a la empleada en regresión por mínimos cuadrados ordinarios, la interpretación final es muy similar. Se obtienen coeficientes de regresión que estiman el efecto que tiene cada predictor sobre un cuantil específico de la variable respuesta.
El poder realizar regresión sobre cualquier parte de la distribución permite conocer la influencia de los predictores desde el mínimo al máximo rango de la variable respuesta. Esto es especialmente útil en modelos de regresión en los que no se cumple la condición de varianza constante, ya que esto significa que no hay un único ratio de cambio (pendiente) que represente bien a toda la variable respuesta a la vez.
Cuantil: El cuantil de orden \(\tau\) \((0 < \tau < 1)\) de una distribución es el valor de la variable \(X\) que marca un corte de modo que una proporción \(\tau\) de valores de la población es menor o igual que dicho valor. Por ejemplo, el cuantil de orden 0,36 deja un 36% de valores por debajo y el cuantil de orden 0,50 se corresponde con la mediana de la distribución.
En la práctica, la proporción de observaciones menor o igual a un determinado plano de regresión del cuantil no es exactamente igual a \(\tau\), debido a la varianza muestral.
Ejemplo
Supóngase las siguientes dos variables simuladas de tal forma que su varianza no es constante.
library(ggplot2)
# Predictor
x <- seq(0, 100, length.out = 100)
# Varianza no constante
varianza <- 0.1 + 0.05*x
# Intercept y pendiente
b_0 <- 6
b_1 <- 0.1
# Error aleatorio normal con varianza no constante
set.seed(1)
error <- rnorm(100, mean = 0, sd = varianza)
# Variable respuesta
y <- b_0 + b_1*x + error
datos <- data.frame(x, y)
ggplot(data = datos, aes(x = x, y = y)) +
geom_point() +
theme_bw()
La representación gráfica muestra una clara falta de homocedasticidad, la varianza de la variable respuesta aumenta con forme aumenta el valor del predictor. Esto supone una violación de las condiciones necesarias para la regresión lineal por mínimos cuadrados ordinales, haciendo que sus resultados tengan una utilidad limitada.
ggplot(data = datos, aes(x = x, y = y)) +
geom_point() +
geom_smooth(method = "lm", colour = "firebrick", se = TRUE) +
theme_bw()
A pesar de que la regresión lineal por mínimos cuadrados ordinarios va a calcular la estimación de la media de forma no sesgada, generando una estimación tan buena como es posible, para valores altos de X (por ejemplo X = 75) los límites de confianza van a ser mucho más optimistas de lo que realmente muestran los datos. Existen diferentes formas de intentar minimizar este problema, como por ejemplo la regresión por Weighted Least Squares. Esta sección se centra en la regresión por cuantiles.
La función rq() del paquete quantreg permite ajustar modelos de regresión para cualquier cuantil.
library(quantreg)
# regresión del quantil 0.9
modelo_q09 <- rq(formula = y ~ x, tau = 0.9, data = datos)
summary(modelo_q09)##
## Call: rq(formula = y ~ x, tau = 0.9, data = datos)
##
## tau: [1] 0.9
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 6.04168 5.93265 6.39328
## x 0.16125 0.15591 0.17826El summary() de un modelo rq devuelve los coeficientes de regresión de cada predictor junto con los correspondientes intervalos de confianza calculados por el método rank (existen 4 métodos adicionales). Acorde al modelo, por cada unidad que se incrementa el predictor X, el cuantil 0.9 de la variable respuesta aumenta en promedio 0.16 unidades. Además, como el intervalo de confianza no contiene el 0, hay evidencias de que el predictor influye de forma significativa.
Para obtener el p-value asociado a cada predictor, se tiene que recurrir a los otros métodos disponibles:
summary(object = modelo_q09, se = "boot")##
## Call: rq(formula = y ~ x, tau = 0.9, data = datos)
##
## tau: [1] 0.9
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 6.04168 0.18565 32.54414 0.00000
## x 0.16125 0.01164 13.85130 0.00000summary(object = modelo_q09, se = "nid")##
## Call: rq(formula = y ~ x, tau = 0.9, data = datos)
##
## tau: [1] 0.9
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 6.04168 0.26581 22.72912 0.00000
## x 0.16125 0.01632 9.88007 0.00000Empleando las estimaciones de intercept y del predictor se puede representar la recta de regresión del cuantil.
ggplot(data = datos, aes(x = x, y = y)) +
geom_point() +
geom_abline(intercept = coef(modelo_q09)[1], slope = coef(modelo_q09)[2],
colour = "firebrick") +
theme_bw()
El paquete gráfico ggplot2 contiene la función geom_quantile() que ajusta y representa directamente la recta de regresión.
ggplot(data = datos, aes(x = x, y = y)) +
geom_point() +
geom_quantile(quantiles = 0.9) +
theme_bw()rq permite ajustar varios cuantiles a la vez indicando la secuencia en el argumento tau.
modelo_q <- rq(formula = y ~ x, tau = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),
data = datos)
summary(object = modelo_q, se = "boot")##
## Call: rq(formula = y ~ x, tau = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
## 0.8, 0.9), data = datos)
##
## tau: [1] 0.1
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 5.95736 0.52287 11.39357 0.00000
## x 0.03796 0.01559 2.43534 0.01669
##
## Call: rq(formula = y ~ x, tau = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
## 0.8, 0.9), data = datos)
##
## tau: [1] 0.2
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 5.93735 0.26311 22.56605 0.00000
## x 0.06943 0.00909 7.63490 0.00000
##
## Call: rq(formula = y ~ x, tau = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
## 0.8, 0.9), data = datos)
##
## tau: [1] 0.3
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 6.19112 0.20374 30.38690 0.00000
## x 0.07395 0.00821 9.00988 0.00000
##
## Call: rq(formula = y ~ x, tau = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
## 0.8, 0.9), data = datos)
##
## tau: [1] 0.4
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 6.08146 0.21617 28.13322 0.00000
## x 0.09208 0.01083 8.49844 0.00000
##
## Call: rq(formula = y ~ x, tau = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
## 0.8, 0.9), data = datos)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 6.02410 0.20302 29.67212 0.00000
## x 0.10351 0.01071 9.66019 0.00000
##
## Call: rq(formula = y ~ x, tau = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
## 0.8, 0.9), data = datos)
##
## tau: [1] 0.6
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 6.16841 0.19185 32.15225 0.00000
## x 0.11507 0.00878 13.10467 0.00000
##
## Call: rq(formula = y ~ x, tau = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
## 0.8, 0.9), data = datos)
##
## tau: [1] 0.7
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 6.09508 0.18344 33.22737 0.00000
## x 0.12784 0.00791 16.16605 0.00000
##
## Call: rq(formula = y ~ x, tau = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
## 0.8, 0.9), data = datos)
##
## tau: [1] 0.8
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 6.10539 0.21215 28.77934 0.00000
## x 0.13783 0.01339 10.29678 0.00000
##
## Call: rq(formula = y ~ x, tau = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
## 0.8, 0.9), data = datos)
##
## tau: [1] 0.9
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 6.04168 0.17196 35.13381 0.00000
## x 0.16125 0.01083 14.88335 0.00000ggplot(data = datos, aes(x = x, y = y)) +
geom_point() +
geom_quantile(quantiles = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),
color = "firebrick") +
theme_bw()
Se puede observar que la pendiente de la recta de cada cuantil es distinta, lo que significa que el predictor X influye de forma distinta a cada cuantil de la variable respuesta. Para estudiar como varia esta influencia y su significancia, se pueden representar la pendiente para cada cuantil.
plot(summary(modelo_q), parm = "x")
Cada punto representa el coeficiente de regresión estimado de un cuantil. La línea continua roja es el coeficiente de regresión estimado para el predictor utilizando mínimos cuadrados ordinarios (lm(y ~ x, data = datos)) y las líneas discontinuas sus límites de confianza del 95%. Para este ejemplo, los cuantiles 0.1, 0.2, 0.3, 0.7, 0.8 y 0.9 se diferencian significativamente del valor obtenido por mínimos cuadrados ordinarios.
Además de conocer si un cuantil específico se ve significativamente influenciado por un predictor (p-value calculado en el summary), puede ser de interés determinar si los coeficientes de regresión de dos cuantiles son distintos. Para ello se ajustan modelos individuales con cada cuantil y se comparan mediante la función anova().
modelo_q05 <- rq(formula = y ~ x, tau = 0.5, data = datos)
modelo_q045 <- rq(formula = y ~ x, tau = 0.45, data = datos)
modelo_q07 <- rq(formula = y ~ x, tau = 0.7, data = datos)
anova(modelo_q05, modelo_q045)## Quantile Regression Analysis of Deviance Table
##
## Model: y ~ x
## Joint Test of Equality of Slopes: tau in { 0.5 0.45 }
##
## Df Resid Df F value Pr(>F)
## 1 1 199 2.3477 0.1271anova(modelo_q05, modelo_q07)## Quantile Regression Analysis of Deviance Table
##
## Model: y ~ x
## Joint Test of Equality of Slopes: tau in { 0.5 0.7 }
##
## Df Resid Df F value Pr(>F)
## 1 1 199 10.03 0.001783 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Los test F confirman que no hay diferencia significativa entre los coeficientes de regresión de los cuantiles 0.45 y 0.5 pero sí entre 0.5 y 0.7.
Regresión de cuantiles para comparar medianas
Cuando se quiere comparar dos poblaciones pero los datos no siguen una distribución normal, por ejemplo por que presentan colas asimétricas, y además el tamaño muestral es pequeño, suele ser adecuado comparar medianas en lugar de medias. Con frecuencia se recurre al test de Mann–Whitney–Wilcoxon para contrastar la hipótesis nula de dos medianas son iguales. Este test solo puede utilizarse con este fin si se cumple la exigente condición de que la única diferencia entre las poblaciones es su localización, el resto de características (asimetría, dispersión…) tienen que ser idénticas. La regresión de cuantiles permite comparar medianas (cuantil 0.5) sin necesidad de que se cumpla esta condición.
Supóngase que se dispone de las siguientes muestras y de que se desea conocer si existe diferencia significativa entre las poblaciones de origen.
grupo_a <- c(5, 5, 5, 5, 5, 5, 7, 8, 9, 10)
grupo_b <- c(1, 2, 3 ,4, 5, 5, 5, 5, 5, 5)
datos <- data.frame(grupo = rep(c("a", "b"), each = 10),
valor = c(grupo_a, grupo_b))
par(mfrow = c(1,2))
hist(grupo_a, col = "cyan4", main = "")
hist(grupo_b, col = "firebrick", main = "")
En vista de que el tamaño muestral es pequeño y de que ambos grupos muestran una clara asimetría, el t-test queda descartado. Una posible alternativa es comparar las medianas. Como la asimetría de los grupos tiene dirección opuesta, el test de Mann–Whitney–Wilcoxon no puede emplearse con esta finalidad. En su lugar se recurre a la regresión de cuantiles.
modelo_q50 <- rq(valor ~ grupo, tau = 0.5, data = datos)
summary(modelo_q50, se = "boot")##
## Call: rq(formula = valor ~ grupo, tau = 0.5, data = datos)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 5.00000 1.12438 4.44691 0.00031
## grupob 0.00000 1.32145 0.00000 1.00000El resultado no muestra evidencias en contra de la hipótesis nula de que las medianas de ambos grupos son iguales. De hecho, la media de ambos grupos es la misma.
median(grupo_a)## [1] 5median(grupo_b)## [1] 5Si se emplea de Mann–Whitney–Wilcoxon, debido a la diferencia en forma de las distribuciones, el resultado es significativo, lo que llevaría a conclusiones erróneas.
wilcox.test(grupo_a, grupo_b, paired = FALSE)##
## Wilcoxon rank sum test with continuity correction
##
## data: grupo_a and grupo_b
## W = 82, p-value = 0.007196
## alternative hypothesis: true location shift is not equal to 0Regresión de cuantiles y regresión OLS en distribuciones normales
Una de las propiedades de la distribución normal es que el valor de la media y mediana son iguales. Esto significa que el resultado de comparar las medias (regresión por mínimos cuadrados ordinarios) y medianas (regresión del cuantil 0.5) tiende a ser el mismo a medida que aumenta el tamaño muestral. Las pequeñas diferencias se deben a variaciones del muestreo.
# Predictor
x <- seq(0, 100, length.out = 100)
# Intercept y pendiente
b_0 <- 6
b_1 <- 0.1
# Error aleatorio normal
set.seed(1)
error <- rnorm(100, mean = 0, sd = 4)
# Variable respuesta
y <- b_0 + b_1*x + error
datos <- data.frame(x, y)
modelo_ols <- lm(y ~ x, data = datos)
summary(modelo_ols)##
## Call:
## lm(formula = y ~ x, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.360 -2.423 0.062 2.341 9.190
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.52486 0.71676 9.103 1.07e-14 ***
## x 0.09821 0.01238 7.931 3.56e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.611 on 98 degrees of freedom
## Multiple R-squared: 0.3909, Adjusted R-squared: 0.3847
## F-statistic: 62.9 on 1 and 98 DF, p-value: 3.558e-12modelo_q50 <- rq(y ~ x, tau = 0.5, data = datos)
summary(modelo_q50)##
## Call: rq(formula = y ~ x, tau = 0.5, data = datos)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 6.74025 5.03377 8.54760
## x 0.09438 0.06034 0.12350De hecho, si la distribución es normal, la pendiente de la regresión de todos los cuantiles es la misma e igual a la obtenida mediante regresión por mínimos cuadrados ordinarios. Cualquier desviación se debe a variaciones de muestreo.
ggplot(data = datos, aes(x = x, y = y)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, color = "firebrick") +
geom_quantile(quantiles = c(1:9)/10) +
theme_bw()
plot(summary(rq(y ~ x, tau = c(1:9)/10, data = datos)), parm = "x")
La superposición de los coeficientes de regresión de cada cuantil con el coeficiente de regresión de la media por OLS y su intervalo del 95% muestra que no hay diferencias significativas entre ellos. De hecho, la comparación entre cualquier par de cuantiles no debería ser significativa si la población es de tipo normal.
modelo_q01 <- rq(formula = y ~ x, tau = 0.1, data = datos)
modelo_q05 <- rq(formula = y ~ x, tau = 0.5, data = datos)
modelo_q09 <- rq(formula = y ~ x, tau = 0.9, data = datos)
anova(modelo_q01, modelo_q05)## Quantile Regression Analysis of Deviance Table
##
## Model: y ~ x
## Joint Test of Equality of Slopes: tau in { 0.1 0.5 }
##
## Df Resid Df F value Pr(>F)
## 1 1 199 0.1164 0.7333anova(modelo_q01, modelo_q09)## Quantile Regression Analysis of Deviance Table
##
## Model: y ~ x
## Joint Test of Equality of Slopes: tau in { 0.1 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 1 199 0.2477 0.6192Bibliografía
A gentle introduction to quantile regression for ecologists, Brian S. Cade and Barry R. Noon
http://data.library.virginia.edu/getting-started-with-quantile-regression/
Quantile regression in r: a vignette, Roger Koenker
This work by Joaquín Amat Rodrigo is licensed under a Creative Commons Attribution 4.0 International License.