-
Notifications
You must be signed in to change notification settings - Fork 17
/
quantile-regression.Rmd
113 lines (80 loc) · 2.18 KB
/
quantile-regression.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
# Quantile Regression
## Data Setup
We'll use the <span class="pack" style = "">quantreg</span> package for comparison, and the classic data set on Belgian household income and food expenditure. Scale income if you want a meaningful 'centercept'.
```{r qr-setup}
library(tidyverse)
library(quantreg)
data(engel)
# engel$income = scale(engel$income)
X = cbind(1, engel$income)
colnames(X) = c('Intercept', 'income')
```
## Function
Loss function. It really is this simple.
```{r qr-func}
qreg <- function(par, X, y, tau) {
lp = X%*%par
res = y - lp
loss = ifelse(res < 0 , -(1 - tau)*res, tau*res)
sum(loss)
}
```
## Estimation
We'll estimate the median to start. Compare `optim` output with <span class="pack" style = "">quantreg</span> package.
```{r qr-est}
optim(
par = c(intercept = 0, income = 0),
fn = qreg,
X = X,
y = engel$foodexp,
tau = .5
)$par
rq(foodexp ~ income, tau = .5, data = engel)
```
### Other quantiles
Now we will add additional quantiles to estimate.
```{r qr-quants}
# quantiles
qs = c(.05, .1, .25, .5, .75, .9, .95)
fit_rq = coef(rq(foodexp ~ income, tau = qs, data = engel))
fit_qreg = map_df(qs, function(tau)
data.frame(t(
optim(
par = c(intercept = 0, income = 0),
fn = qreg,
X = X,
y = engel$foodexp,
tau = tau
)$par
)))
```
## Comparison
Compare results.
```{r qr-compare, echo =FALSE}
rbind(fit_rq, t(fit_qreg)) %>%
as.data.frame() %>%
rownames_to_column('coef') %>%
mutate(` ` = rep(c('fit_rq', 'fit_qreg'), each = 2)) %>%
select(` `, everything()) %>%
kable_df()
```
## Visualize
Let's visualize the results.
```{r qr-vis, echo=FALSE}
engel %>%
qplot(data = .,
income,
foodexp,
color = I(scales::alpha('orange', .25))) +
geom_abline(aes(
intercept = intercept,
slope = income,
color = group
),
data = data.frame(fit_qreg, group = factor(qs))) +
scico::scale_color_scico_d(palette = 'vik')
```
## Python
The above is available as a Python demo in the [supplemental section](#python-qreg).
## Source
Original code available at https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/quantile_regression.Rmd