1.- 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 τ (0<τ<1) de una distribución es el valor de la variable X que marca un corte de modo que una proporción τ 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.

#** 2.- REGRESIÓN DE CUANTILES**

library(ggplot2) 
## Warning: package 'ggplot2' was built under R version 4.2.3
# 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(size=1, col="blue") +
theme_bw()

La representación gráfica muestra una clara falta de homocedasticidad, la varianza de la variable respuesta aumenta conforme 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(size = 1, col="blue") + 
  geom_smooth(method = "lm", colour = "red", se = TRUE) + 
  theme_bw() 
## `geom_smooth()` using formula = 'y ~ x'

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.

library(quantreg)
## Warning: package 'quantreg' was built under R version 4.2.3
## Loading required package: SparseM
## 
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
## 
##     backsolve
# regresión del cuantil 0.9
modelo_q09 <- rq(formula = y ~x, tau = 0.9, data = datos)
summary(modelo_q09)
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique
## 
## 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.17826

El 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.00000
summary(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.00000

Empleando 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(size=1, col = "blue") +
  geom_abline(intercept = coef(modelo_q09)[1], slope = coef(modelo_q09)[2],
              colour = "red") +
  theme_bw()

ggplot(data=datos, aes(x = x, y = y)) +
  geom_point(size = 1, col = "blue") +
  geom_quantile(quantiles = 0.9, col= "red") +
  theme_bw()
## Smoothing formula not specified. Using: y ~ x

Para varios cuantiles de la distribución

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.00000
ggplot(data = datos, aes(x = x, y=y)) +
  geom_point(size=1, col="blue") +
  geom_quantile(quantiles = c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),
                color = "red") +
  theme_bw()
## Smoothing formula not specified. Using: y ~ x

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 cómo varía esta influencia y su significancia, se puede representar la pendiente para cada cuantil.

plot(summary(modelo_q), parm="x")
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

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)
## Warning in summary.rq(x, se = se, R = R, covariance = TRUE): 2 non-positive fis
## 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.1271

Los test F confirman que no hay diferencia significativa entre los coeficientes de regresión de re los cuantiles 0.45 y 0.5, pero sí entre 0.5 y 0.7.

#3.- 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 = "cyan", main = "")
hist(grupo_b, col = "red", 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)
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
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.00000

El 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] 5
median(grupo_b)
## [1] 5

Si 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)
## Warning in wilcox.test.default(grupo_a, grupo_b, paired = FALSE): cannot
## compute exact p-value with ties
## 
##  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 0

#** 4.- REGRESIÓN DE CUANTILES Y REGRESIÓN OLS EN DIST. 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-12
modelo_q50 <- rq(y ~ x, tau = 0.5, data = datos)
summary(modelo_q50)
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique
## 
## 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.12350
ggplot(data = datos, aes(x = x, y = y)) +
  geom_point(size = 1, col="blue") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  geom_quantile(quantiles = c(1:9)/10) +
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
## Smoothing formula not specified. Using: y ~ x

plot(summary(rq(y ~ x, tau = c(1:9)/10, data = datos)), parm = "x")
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique

La superposición (están dentor de las líneas discontinuas) 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.7333
anova(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.6192