# Clear work space: 
rm(list = ls()) 

# Load quantreg for Quantile Regression: 

library(quantreg)

# Use engel data set provided by Koenker and Bassett (1982): 

data(engel)

# Conduct QR: 

my_qr <- rq(foodexp ~ income, data = engel, tau = seq(0.25, 0.75, 0.25))

# Show results: 
summary(my_qr) 

## 
## Call: rq(formula = foodexp ~ income, tau = seq(0.25, 0.75, 0.25), data = engel)
## 
## tau: [1] 0.25
## 
## Coefficients:
##             coefficients lower bd  upper bd 
## (Intercept)  95.48354     73.78608 120.09847
## income        0.47410      0.42033   0.49433
## 
## Call: rq(formula = foodexp ~ income, tau = seq(0.25, 0.75, 0.25), data = engel)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             coefficients lower bd  upper bd 
## (Intercept)  81.48225     53.25915 114.01156
## income        0.56018      0.48702   0.60199
## 
## Call: rq(formula = foodexp ~ income, tau = seq(0.25, 0.75, 0.25), data = engel)
## 
## tau: [1] 0.75
## 
## Coefficients:
##             coefficients lower bd  upper bd 
## (Intercept)  62.39659     32.74488 107.31362
## income        0.64401      0.58016   0.69041

# Alternatives: 

q25 <- rq(foodexp ~ income, data = engel, tau = 0.25)
q50 <- rq(foodexp ~ income, data = engel, tau = 0.50)
q75 <- rq(foodexp ~ income, data = engel, tau = 0.75)


# Compare results: 

library(stargazer)
stargazer(q25, q50, q75, title = "Quantile Regression Results", type = "text")

## 
## Quantile Regression Results
## ==========================================
##                   Dependent variable:     
##              -----------------------------
##                         foodexp           
##                 (1)       (2)       (3)   
## ------------------------------------------
## income       0.474***  0.560***  0.644*** 
##               (0.029)   (0.028)   (0.023) 
##                                           
## Constant     95.484*** 81.482*** 62.397***
##              (21.392)  (19.251)  (16.305) 
##                                           
## ------------------------------------------
## Observations    235       235       235   
## ==========================================
## Note:          *p<0.1; **p<0.05; ***p<0.01

# Visualization for Quantile Regression with some tau values: 
library(tidyverse)

intercept_slope <- my_qr %>% 
  coef() %>% 
  t() %>% 
  data.frame() %>% 
  rename(intercept = X.Intercept., slope = income) %>% 
  mutate(quantile = row.names(.))


ggplot() + 
  geom_point(data = engel, aes(income, foodexp), alpha = 0.5) + 
  geom_abline(data = intercept_slope, aes(intercept = intercept, slope = slope, color = quantile)) + 
  theme_minimal() + 
  labs(x = "Income", y = "Food Expenditure", 
       title = "Quantile Regression with tau = 0.25, 0.50 and 0.75", 
       caption = "Data Source: Koenker and Bassett (1982)")

LS0tDQp0aXRsZTogIlF1YW50aWxlIFJlZ3Jlc3Npb24gVXNpbmcgUiINCmF1dGhvcjogIk5ndXllbiBDaGkgRHVuZyINCnN1YnRpdGxlOiAiRGFpbHkgR3JhcGggU2VyaWVzIg0Kb3V0cHV0Og0KICBodG1sX2RvY3VtZW50Og0KICAgIGNvZGVfZG93bmxvYWQ6IHllcw0KICAgICMgY29kZV9mb2xkaW5nOiBoaWRlDQogICAgaGlnaGxpZ2h0OiB6ZW5idXJuDQogICAgdGhlbWU6IGZsYXRseQ0KICAgIHRvYzogeWVzDQogICAgdG9jX2Zsb2F0OiB5ZXMNCiAgd29yZF9kb2N1bWVudDoNCiAgICB0b2M6IHllcw0KLS0tDQoNCmBgYHtyIHNldHVwLGluY2x1ZGU9RkFMU0V9DQprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUsIHdhcm5pbmcgPSBGQUxTRSwgbWVzc2FnZSA9IEZBTFNFLCBmaWcucmV0aW5hPTIpDQpgYGANCg0KIyBSZWZlcmVuY2VzDQoNCmh0dHA6Ly93d3cuZGVwb2Nlbi5vcmcvdm4vb3VyLXB1YmxpY2F0aW9uL2FydGljbGVzLzE4OS1jYWMteWV1LXRvLWFuaC1odW9uZy1kZW4ta2V0LXF1YS1ob2MtdGFwLW8taG9hLWt5Lmh0bWwNCg0KIyBSIENvZGVzIGZvciBDb25kdWN0aW5nIFF1YW50aWxlIFJlZ3Jlc3Npb24NCg0KYGBge3J9DQoNCiMgQ2xlYXIgd29yayBzcGFjZTogDQpybShsaXN0ID0gbHMoKSkgDQoNCiMgTG9hZCBxdWFudHJlZyBmb3IgUXVhbnRpbGUgUmVncmVzc2lvbjogDQoNCmxpYnJhcnkocXVhbnRyZWcpDQoNCiMgVXNlIGVuZ2VsIGRhdGEgc2V0IHByb3ZpZGVkIGJ5IEtvZW5rZXIgYW5kIEJhc3NldHQgKDE5ODIpOiANCg0KZGF0YShlbmdlbCkNCg0KIyBDb25kdWN0IFFSOiANCg0KbXlfcXIgPC0gcnEoZm9vZGV4cCB+IGluY29tZSwgZGF0YSA9IGVuZ2VsLCB0YXUgPSBzZXEoMC4yNSwgMC43NSwgMC4yNSkpDQoNCiMgU2hvdyByZXN1bHRzOiANCnN1bW1hcnkobXlfcXIpIA0KDQojIEFsdGVybmF0aXZlczogDQoNCnEyNSA8LSBycShmb29kZXhwIH4gaW5jb21lLCBkYXRhID0gZW5nZWwsIHRhdSA9IDAuMjUpDQpxNTAgPC0gcnEoZm9vZGV4cCB+IGluY29tZSwgZGF0YSA9IGVuZ2VsLCB0YXUgPSAwLjUwKQ0KcTc1IDwtIHJxKGZvb2RleHAgfiBpbmNvbWUsIGRhdGEgPSBlbmdlbCwgdGF1ID0gMC43NSkNCg0KDQojIENvbXBhcmUgcmVzdWx0czogDQoNCmxpYnJhcnkoc3RhcmdhemVyKQ0Kc3RhcmdhemVyKHEyNSwgcTUwLCBxNzUsIHRpdGxlID0gIlF1YW50aWxlIFJlZ3Jlc3Npb24gUmVzdWx0cyIsIHR5cGUgPSAidGV4dCIpDQoNCg0KIyBWaXN1YWxpemF0aW9uIGZvciBRdWFudGlsZSBSZWdyZXNzaW9uIHdpdGggc29tZSB0YXUgdmFsdWVzOiANCmxpYnJhcnkodGlkeXZlcnNlKQ0KDQppbnRlcmNlcHRfc2xvcGUgPC0gbXlfcXIgJT4lIA0KICBjb2VmKCkgJT4lIA0KICB0KCkgJT4lIA0KICBkYXRhLmZyYW1lKCkgJT4lIA0KICByZW5hbWUoaW50ZXJjZXB0ID0gWC5JbnRlcmNlcHQuLCBzbG9wZSA9IGluY29tZSkgJT4lIA0KICBtdXRhdGUocXVhbnRpbGUgPSByb3cubmFtZXMoLikpDQoNCg0KZ2dwbG90KCkgKyANCiAgZ2VvbV9wb2ludChkYXRhID0gZW5nZWwsIGFlcyhpbmNvbWUsIGZvb2RleHApLCBhbHBoYSA9IDAuNSkgKyANCiAgZ2VvbV9hYmxpbmUoZGF0YSA9IGludGVyY2VwdF9zbG9wZSwgYWVzKGludGVyY2VwdCA9IGludGVyY2VwdCwgc2xvcGUgPSBzbG9wZSwgY29sb3IgPSBxdWFudGlsZSkpICsgDQogIHRoZW1lX21pbmltYWwoKSArIA0KICBsYWJzKHggPSAiSW5jb21lIiwgeSA9ICJGb29kIEV4cGVuZGl0dXJlIiwgDQogICAgICAgdGl0bGUgPSAiUXVhbnRpbGUgUmVncmVzc2lvbiB3aXRoIHRhdSA9IDAuMjUsIDAuNTAgYW5kIDAuNzUiLCANCiAgICAgICBjYXB0aW9uID0gIkRhdGEgU291cmNlOiBLb2Vua2VyIGFuZCBCYXNzZXR0ICgxOTgyKSIpDQoNCg0KDQoNCg0KYGBgDQoNCg==