{"word_count": 482, "character_count": 4994}
[an updated version is now online here]After my series of post on classification algorithms, it's time to get back to R codes, this time for quantile regression. Yes, I still want to get a better understanding of optimization routines, in R. Before looking at the quantile regression, let us compute the median, or the quantile, from a sample.MedianConsider a sample [latex]\{y_1,\cdots,y_n\}[/latex]. To compute the median, solve[latex display="true"]\min_\mu \left\lbrace\sum_{i=1}^n|y_i-\mu|\right\rbrace[/latex]which can be solved using linear programming techniques. More precisely, this problem is equivalent to[latex display="true"]\min_{\mu,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^na_i+b_i\right\rbrace[/latex]with [latex]a_i,b_i\geq 0[/latex] and [latex]y_i-\mu=a_i-b_i[/latex], [latex]\forall i=1,\cdots,n[/latex].To illustrate, consider a sample from a lognormal distribution,n = 101 set.seed(1)y = rlnorm(n)median(y)[1] 1.077415For the optimization problem, use the matrix form, with [latex]3n[/latex] constraints, and [latex]2n+1[/latex] parameters,library(lpSolve)A1 = cbind(diag(2*n),0) A2 = cbind(diag(n), -diag(n), 1)r = lp("min", c(rep(1,2*n),0),rbind(A1, A2),c(rep("&#62;=", 2*n), rep("=", n)), c(rep(0,2*n), y))tail(r$solution,1) [1] 1.077415It looks like it's working well...QuantileOf course, we can adapt our previous code for quantilestau = .3quantile(x,tau)      30% 0.6741586The linear program is now[latex display="true"]\min_{\mu,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^n\tau a_i+(1-\tau)b_i\right\rbrace[/latex]with [latex]a_i,b_i\geq 0[/latex] and [latex]y_i-\mu=a_i-b_i[/latex], [latex]\forall i=1,\cdots,n[/latex]. The R code is nowA1 = cbind(diag(2*n),0) A2 = cbind(diag(n), -diag(n), 1)r = lp("min", c(rep(tau,n),rep(1-tau,n),0),rbind(A1, A2),c(rep("&#62;=", 2*n), rep("=", n)), c(rep(0,2*n), y))tail(r$solution,1) [1] 0.6741586So far so good...Quantile Regression (simple)Consider the following dataset, with rents of flat, in a major German city, as function of the surface, the year of construction, etc.base=read.table("http://freakonometrics.free.fr/rent98_00.txt",header=TRUE)The linear program for the quantile regression is now[latex display="true"]\min_{\mu,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^n\tau a_i+(1-\tau)b_i\right\rbrace[/latex]with [latex]a_i,b_i\geq 0[/latex] and [latex display="true"]y_i-[\beta_0^\tau+\beta_1^\tau x_i]=a_i-b_i[/latex][latex]\forall i=1,\cdots,n[/latex]. So use hererequire(lpSolve) tau = .3n=nrow(base)X = cbind( 1, base$area)y = base$rent_euroA1 = cbind(diag(2*n), 0,0) A2 = cbind(diag(n), -diag(n), X) r = lp("min",       c(rep(tau,n), rep(1-tau,n),0,0), rbind(A1, A2),       c(rep("&#62;=", 2*n), rep("=", n)), c(rep(0,2*n), y)) tail(r$solution,2)[1] 148.946864   3.289674Of course, we can use R function to fit that modellibrary(quantreg)rq(rent_euro~area, tau=tau, data=base)Coefficients:(Intercept)        area  148.946864    3.289674Here again, it seems to work quite well. We can use a different probability level, of course, and get a plotplot(base$area,base$rent_euro,xlab=expression(paste("surface (",m^2,")")),     ylab="rent (euros/month)",col=rgb(0,0,1,.4),cex=.5)sf=0:250yr=r$solution[2*n+1]+r$solution[2*n+2]*sflines(sf,yr,lwd=2,col="blue")tau = .9r = lp("min",       c(rep(tau,n), rep(1-tau,n),0,0), rbind(A1, A2),       c(rep("&#62;=", 2*n), rep("=", n)), c(rep(0,2*n), y)) tail(r$solution,2)[1] 121.815505   7.865536yr=r$solution[2*n+1]+r$solution[2*n+2]*sflines(sf,yr,lwd=2,col="blue")Quantile Regression (multiple)Now that we understand how to run the optimization program with one covariate, why not try with two ? For instance, let us see if we can explain the rent of a flat as a (linear) function of the surface and the age of the building.require(lpSolve) tau = .3n=nrow(base)X = cbind( 1, base$area, base$yearc )y = base$rent_euroA1 = cbind(diag(2*n), 0,0,0) A2 = cbind(diag(n), -diag(n), X) r = lp("min",       c(rep(tau,n), rep(1-tau,n),0,0,0), rbind(A1, A2),       c(rep("&#62;=", 2*n), rep("=", n)), c(rep(0,2*n), y)) tail(r$solution,3)[1] 0.000000 3.257562 0.077501Unfortunately, this time, it is not working well...library(quantreg)rq(rent_euro~area+yearc, tau=tau, data=base)Coefficients: (Intercept)         area        yearc -5542.503252     3.978135     2.887234Results are quite different. And actually, another technique can confirm the later (IRLS - Iteratively Reweighted Least Squares)eps = residuals(lm(rent_euro~area+yearc, data=base))for(s in 1:500){  reg = lm(rent_euro~area+yearc, data=base, weights=(tau*(eps&#62;0)+(1-tau)*(eps&#60;0))/abs(eps))  eps = residuals(reg)}reg$coefficients (Intercept)         area        yearc -5484.443043     3.955134     2.857943I could not figure out what went wrong with the linear program. Not only coefficients are very different, but also predictions...yr = r$solution[2*n+1]+r$solution[2*n+2]*base$area+r$solution[2*n+3]*base$yearcplot(predict(reg),yr)abline(a=0,b=1,lty=2,col="red")It's now time to investigate....