Regresión Cuantílica
2024-04-16
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.
2.- EJEMPLO Introductorio
# EJEMPLO 1 simulado
library(ggpubr)## Loading required package: ggplot2library(ggplot2)
# Variable X
x <- seq(0, 100, length.out = 100)
#Varianza heterocedástica
varianza <- 0.1 + 0.75*x
#Parámetros de la ecuación
b_0 <- 1
b_1 <- 2.5
# Error normal estocástico
set.seed(12345)
error <- rnorm(100, mean = 0, sd = varianza)
# respuesta
y <- b_0 + b_1*x + error
datos <- data.frame(x, y)
#Plot
ggplot(data=datos,aes(x=x, y=y)) + geom_point()+theme_classic2()
# Regresión OLS
ggplot(data = datos, aes(x=x, y=y)) +
geom_point() +
geom_smooth(method = "lm", colour = "red", se = TRUE) +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
X > 50 –> La dispersión es más amplia. El intervalo de confianza calculado por OLS es más o menos constante, pero debería ser más amplio cuando hay más dispersión en cuando X es mayor a 30; por lo tanto, el IC no está totalmente bien calculado.
# Regresión de cuartiles
#-------------------------------------------------
library(quantreg)## Warning: package 'quantreg' was built under R version 4.3.3## Loading required package: SparseM##
## Attaching package: 'SparseM'## The following object is masked from 'package:base':
##
## backsolvemodelo_q1 <- rq(formula = y ~x, tau = c(0.25, 0.50, 0.75),
data = datos)
summary(object = modelo_q1, se = "boot")##
## Call: rq(formula = y ~ x, tau = c(0.25, 0.5, 0.75), data = datos)
##
## tau: [1] 0.25
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 1.05855 2.76532 0.38280 0.70270
## x 2.04895 0.15254 13.43247 0.00000
##
## Call: rq(formula = y ~ x, tau = c(0.25, 0.5, 0.75), data = datos)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -0.35595 2.86543 -0.12422 0.90140
## x 2.90813 0.18799 15.46997 0.00000
##
## Call: rq(formula = y ~ x, tau = c(0.25, 0.5, 0.75), data = datos)
##
## tau: [1] 0.75
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -0.10366 2.39117 -0.04335 0.96551
## x 3.24257 0.16537 19.60798 0.00000library(ggplot2)
ggplot(data = datos, aes(x = x, y=y)) +
geom_point()+
geom_quantile(quantiles = c(0.25, 0.50, 0.75), color = "red")+
theme_bw()## Smoothing formula not specified. Using: y ~ x
plot(summary(modelo_q1), 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
Vemos que existen diferencias significativas entre usar un modelo de
cuantiles y OLS, ya que los efectos de la inclinación (parámetrob1) son
diferentes al 95% de confianza.
# Regresión de cuartiles
#-------------------------------------------------
library(quantreg)
modelo_q2 <- 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)## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonuniquesummary(object = modelo_q2, 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) -12.07734 8.90067 -1.35690 0.17793
## x 1.89136 0.18873 10.02141 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.2
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 1.05855 2.44811 0.43240 0.66640
## x 1.90960 0.11450 16.67733 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) 1.58470 2.86543 0.55304 0.58149
## x 2.11467 0.22950 9.21411 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) -0.25413 3.33203 -0.07627 0.93936
## x 2.55878 0.29176 8.77026 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) -0.35595 3.35278 -0.10616 0.91567
## x 2.90813 0.19473 14.93413 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) -0.22484 2.04632 -0.10987 0.91273
## x 3.01891 0.11233 26.87489 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) 0.22956 2.50301 0.09171 0.92711
## x 3.16010 0.10635 29.71513 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) 0.73865 4.87541 0.15151 0.87989
## x 3.36107 0.28741 11.69450 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) 0.16782 4.48556 0.03741 0.97023
## x 3.92620 0.19642 19.98859 0.00000library(ggplot2)
ggplot(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 = "red")+
theme_bw()## Smoothing formula not specified. Using: y ~ x## Warning in rq.fit.br(wx, wy, tau = tau, ...): Solution may be nonunique
plot(summary(modelo_q2), 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
#EJEMPLO 2 simulado
library(ggpubr)
library(ggplot2)
# Variable X
x <- seq(100, 0, lenght.out = 100)## Warning: In seq.default(100, 0, lenght.out = 100) :
## extra argument 'lenght.out' will be disregarded#Varianza heterocedástica
varianza <- 0.1 + 0.8*x
#Parámetros de la ecuación
b_0 <- 5
b_1 <- 2.5
#Error normal estocástico
set.seed(12345)
error <- rnorm(100, mean = 0, sd = varianza)
#respuesta
y <- b_0 - b_1*x + error## Warning in b_0 - b_1 * x + error: longitud de objeto mayor no es múltiplo de la
## longitud de uno menordatos <- data.frame(x, y)
ggplot(data = datos, aes(x = x, y = y)) +
geom_point() +
theme_bw()
Es heterocedástica la varianza.
# Regresion OLS
ggplot(data = datos, aes(x=x, y=y)) +
geom_point()+
geom_smooth(method = "lm", colour = "red", se = TRUE) +
geom_smooth(method = "lm", colour = "red", se = TRUE) +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
# Regresión de deciles
#-----------------------------------
library(quantreg)
modelo_q3 <- 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_q3, 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) 6.96147 7.30674 0.95275 0.34304
## x -3.62776 0.30239 -11.99699 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.2
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 4.07252 3.17579 1.28236 0.20271
## x -3.09051 0.22054 -14.01357 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) 4.59184 2.52954 1.81528 0.07251
## x -2.76074 0.15741 -17.53890 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) 8.43088 4.62566 1.82263 0.07138
## x -2.62256 0.20420 -12.84307 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) 7.52970 4.98562 1.51029 0.13416
## x -2.23182 0.22963 -9.71919 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.26312 3.84001 1.63101 0.10606
## x -2.04362 0.13351 -15.30712 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) 11.34007 3.65457 3.10298 0.00250
## x -1.99575 0.10283 -19.40738 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) 12.27541 8.80110 1.39476 0.16621
## x -1.87628 0.24648 -7.61239 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) 20.66961 14.32706 1.44270 0.15226
## x -1.42548 0.36111 -3.94749 0.00015ggplot(data = datos, aes(x = x, y = y)) +
geom_point()+
geom_smooth(method = "lm", se = TRUE, 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
No es lo mismo hacer una regresión cuantílica que hacer una estimación
de media con ‘mínimos cuadrados ordinarios’. Vemos que va disminuyendo
la renta realmente a medida que aumenta la x.
# Contraste de significación conjunta para saber si hay que hacer regresión cuantílica
#-----------------------------------------------------------------------------------------
modelo_q01 <- rq(formula = y ~x, tau = 0.1, data = datos)
modelo_q02 <- rq(formula = y ~x, tau = 0.2, data = datos)
modelo_q03 <- rq(formula = y ~x, tau = 0.3, data = datos)
modelo_q04 <- rq(formula = y ~x, tau = 0.4, data = datos)
modelo_q05 <- rq(formula = y ~x, tau = 0.5, data = datos)
modelo_q06 <- rq(formula = y ~x, tau = 0.6, data = datos)
modelo_q07 <- rq(formula = y ~x, tau = 0.7, data = datos)
modelo_q08 <- rq(formula = y ~x, tau = 0.8, data = datos)
modelo_q09 <- rq(formula = y ~x, tau = 0.9, data = datos)anova(modelo_q01, modelo_q02, modelo_q03, modelo_q04, modelo_q05, modelo_q06, modelo_q07, modelo_q08, modelo_q09)## 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.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 8 901 8.1955 9.556e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rechazo H0; por tanto, es preferible la regresión de cuantiles.
#Ejemplo3
# Predictor
x <- seq(0, 100, length.out = 100)
# Intercept y pendiente
b_0 <- 6
b_1 <- 0.5
# 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.49821 0.01238 40.233 < 2e-16 ***
## ---
## 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.9429, Adjusted R-squared: 0.9423
## F-statistic: 1619 on 1 and 98 DF, p-value: < 2.2e-16#Regresión OLS
ggplot(data = datos, aes(x = x, y =y)) +
geom_point()+
geom_smooth(method = "lm", colour = "red", se = TRUE)## `geom_smooth()` using formula = 'y ~ x'
theme_bw()## List of 136
## $ line :List of 6
## ..$ colour : chr "black"
## ..$ linewidth : num 0.5
## ..$ linetype : num 1
## ..$ lineend : chr "butt"
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ rect :List of 5
## ..$ fill : chr "white"
## ..$ colour : chr "black"
## ..$ linewidth : num 0.5
## ..$ linetype : num 1
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ text :List of 11
## ..$ family : chr ""
## ..$ face : chr "plain"
## ..$ colour : chr "black"
## ..$ size : num 11
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : num 0
## ..$ lineheight : num 0.9
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ title : NULL
## $ aspect.ratio : NULL
## $ axis.title : NULL
## $ axis.title.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.75points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 2.75points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.bottom : NULL
## $ axis.title.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.75points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.y.left : NULL
## $ axis.title.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 2.75points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey30"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.2points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 2.2points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.bottom : NULL
## $ axis.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 1
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.2points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.y.left : NULL
## $ axis.text.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 2.2points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.theta : NULL
## $ axis.text.r :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0.5
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.2points 0points 2.2points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.ticks :List of 6
## ..$ colour : chr "grey20"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ axis.ticks.x : NULL
## $ axis.ticks.x.top : NULL
## $ axis.ticks.x.bottom : NULL
## $ axis.ticks.y : NULL
## $ axis.ticks.y.left : NULL
## $ axis.ticks.y.right : NULL
## $ axis.ticks.theta : NULL
## $ axis.ticks.r : NULL
## $ axis.minor.ticks.x.top : NULL
## $ axis.minor.ticks.x.bottom : NULL
## $ axis.minor.ticks.y.left : NULL
## $ axis.minor.ticks.y.right : NULL
## $ axis.minor.ticks.theta : NULL
## $ axis.minor.ticks.r : NULL
## $ axis.ticks.length : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ axis.ticks.length.x : NULL
## $ axis.ticks.length.x.top : NULL
## $ axis.ticks.length.x.bottom : NULL
## $ axis.ticks.length.y : NULL
## $ axis.ticks.length.y.left : NULL
## $ axis.ticks.length.y.right : NULL
## $ axis.ticks.length.theta : NULL
## $ axis.ticks.length.r : NULL
## $ axis.minor.ticks.length : 'rel' num 0.75
## $ axis.minor.ticks.length.x : NULL
## $ axis.minor.ticks.length.x.top : NULL
## $ axis.minor.ticks.length.x.bottom: NULL
## $ axis.minor.ticks.length.y : NULL
## $ axis.minor.ticks.length.y.left : NULL
## $ axis.minor.ticks.length.y.right : NULL
## $ axis.minor.ticks.length.theta : NULL
## $ axis.minor.ticks.length.r : NULL
## $ axis.line : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.line.x : NULL
## $ axis.line.x.top : NULL
## $ axis.line.x.bottom : NULL
## $ axis.line.y : NULL
## $ axis.line.y.left : NULL
## $ axis.line.y.right : NULL
## $ axis.line.theta : NULL
## $ axis.line.r : NULL
## $ legend.background :List of 5
## ..$ fill : NULL
## ..$ colour : logi NA
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.margin : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
## ..- attr(*, "unit")= int 8
## $ legend.spacing : 'simpleUnit' num 11points
## ..- attr(*, "unit")= int 8
## $ legend.spacing.x : NULL
## $ legend.spacing.y : NULL
## $ legend.key : NULL
## $ legend.key.size : 'simpleUnit' num 1.2lines
## ..- attr(*, "unit")= int 3
## $ legend.key.height : NULL
## $ legend.key.width : NULL
## $ legend.key.spacing : 'simpleUnit' num 5.5points
## ..- attr(*, "unit")= int 8
## $ legend.key.spacing.x : NULL
## $ legend.key.spacing.y : NULL
## $ legend.frame : NULL
## $ legend.ticks : NULL
## $ legend.ticks.length : 'rel' num 0.2
## $ legend.axis.line : NULL
## $ legend.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.text.position : NULL
## $ legend.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.title.position : NULL
## $ legend.position : chr "right"
## $ legend.position.inside : NULL
## $ legend.direction : NULL
## $ legend.byrow : NULL
## $ legend.justification : chr "center"
## $ legend.justification.top : NULL
## $ legend.justification.bottom : NULL
## $ legend.justification.left : NULL
## $ legend.justification.right : NULL
## $ legend.justification.inside : NULL
## $ legend.location : NULL
## $ legend.box : NULL
## $ legend.box.just : NULL
## $ legend.box.margin : 'margin' num [1:4] 0cm 0cm 0cm 0cm
## ..- attr(*, "unit")= int 1
## $ legend.box.background : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ legend.box.spacing : 'simpleUnit' num 11points
## ..- attr(*, "unit")= int 8
## [list output truncated]
## - attr(*, "class")= chr [1:2] "theme" "gg"
## - attr(*, "complete")= logi TRUE
## - attr(*, "validate")= logi TRUE# Regresión de deciles
#--------------------------------
library(quantreg)
modelo_q4 <- 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_q4, 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) 2.68990 2.05876 1.30656 0.19442
## x 0.48395 0.03006 16.09761 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.2
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 3.49418 1.13471 3.07937 0.00269
## x 0.50172 0.01922 26.09970 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) 5.92341 1.10758 5.34805 0.00000
## x 0.47738 0.02017 23.66365 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.12755 0.85652 7.15398 0.00000
## x 0.48810 0.01804 27.05332 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.74025 1.09245 6.16983 0.00000
## x 0.49438 0.01859 26.58875 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) 7.68284 0.82533 9.30886 0.00000
## x 0.49591 0.01549 32.01511 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) 8.41102 0.72995 11.52276 0.00000
## x 0.49815 0.01375 36.22220 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) 9.17249 0.90674 10.11587 0.00000
## x 0.49793 0.01961 25.39726 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) 10.43941 1.27425 8.19257 0.00000
## x 0.50427 0.02262 22.29470 0.00000ggplot(data = datos, aes(x = x, y = y)) +
geom_point() +
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 inclinación de todos los cuantiles es similar a la de OLS.
# Contraste de significación conjunta para saber si hay que hacer regresión cuantílica
#-----------------------------------------------------------------------------------------
modelo_q01 <- rq(formula = y ~x, tau = 0.1, data = datos)
modelo_q02 <- rq(formula = y ~x, tau = 0.2, data = datos)
modelo_q03 <- rq(formula = y ~x, tau = 0.3, data = datos)
modelo_q04 <- rq(formula = y ~x, tau = 0.4, data = datos)
modelo_q05 <- rq(formula = y ~x, tau = 0.5, data = datos)
modelo_q06 <- rq(formula = y ~x, tau = 0.6, data = datos)
modelo_q07 <- rq(formula = y ~x, tau = 0.7, data = datos)
modelo_q08 <- rq(formula = y ~x, tau = 0.8, data = datos)
modelo_q09 <- rq(formula = y ~x, tau = 0.9, data = datos)anova(modelo_q01, modelo_q02)## Quantile Regression Analysis of Deviance Table
##
## Model: y ~ x
## Joint Test of Equality of Slopes: tau in { 0.1 0.2 }
##
## Df Resid Df F value Pr(>F)
## 1 1 199 0.5501 0.4592Al ser Pr(>F) mayor a 0.05, la inclinación de las ecuaciones del cuantil del 10% y 20% son similares.
anova(modelo_q01, modelo_q02, modelo_q03, modelo_q04, modelo_q05, modelo_q06, modelo_q07, modelo_q08, modelo_q09)## Quantile Regression Analysis of Deviance Table
##
## Model: y ~ x
## Joint Test of Equality of Slopes: tau in { 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 }
##
## Df Resid Df F value Pr(>F)
## 1 8 892 0.676 0.713Las inclinaciones de todos los cuantiles (10%, 20%,…, 90%) son similares.
Da igual hacer regresión de cuantiles que OLS, porque las pendientes son similares.
plot(summary(rq(y~x, tau = c(1:9)/10, data = datos)), parm = 1)## 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
Las constantes son diferentes para cada cuantil. Obviamente cuantiles
inferiores empiezan en valor de Y inferiores. (Esto no es suficiente
para hacer regresiones de cuantiles).
3.- BIBLIOGRAFIA
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