Kevin Chang
2012-Jun-07 18:49 UTC
[R] Quantile regression: Discrepencies Between optimizer and rq()
Hello Everyone, I'm currently learning about quantile regressions. I've been using an optimizer to compare with the rq() command for quantile regression. When I run the code, the results show that my coefficients are consistent with rq(), but the intercept term can vary by a lot. I don't think my optimizer code is wrong and suspects it has something to do with the starting values. The results seems very sensitive to different starting values and I don't know how to make sense of it. Advice from the community would be greatly appreciated. Sincerely, Kevin Chang ###################### CODE Below ########################### library(quantreg) data(engel) y<-cbind(engel[,2]) x<-cbind(rep(1,length(engel[,1])),engel[,1]) x1<-cbind(engel[,1]) nn<-nrow(engel) nn bhat.ls<-solve(t(x)%*%x)%*%t(x)%*%y #bhat.ls # QUANTILES quant=.25 fr.1=function(bhat.opt) { uu=y-x%*%bhat.opt sample.cond.quantile=quantile(uu,quant) w.less=rep(0,nn) for(ii in 1:nn){if(uu[ii]<sample.cond.quantile) w.less[ii]=1} sum((quant-1)*sum((y-x%*%bhat.opt)*w.less) #negative residuals +quant*sum((y-x%*%bhat.opt)*(1-w.less))) #positive residuals } start<-c(0,0) result=optim(start,fr.1) bhat.cond=result$par #Quantile Command Results fit.temp=rq(y~x1,tau=quant) fit.temp #OPTIMIZER Results bhat.cond #OLS Command Results mean=lm(y~x1) mean [[alternative HTML version deleted]]
Roger Koenker
2012-Jun-07 19:40 UTC
[R] Quantile regression: Discrepencies Between optimizer and rq()
Optim() by default is using Nelder-Mead which is an extremely poor way to do linear programming, despite the fact that ?optim says that: "It will work reasonably well for non-differentiable functions." I didn't check your coding of the objective function fully, but at the very least you should explicitly pass the arguments y, x, and quant. and you need to replace what you call sample.cond.quantile by 0 in the definition of w.less. url: www.econ.uiuc.edu/~roger Roger Koenker email rkoenker at uiuc.edu Department of Economics vox: 217-333-4558 University of Illinois fax: 217-244-6678 Urbana, IL 61801 On Jun 7, 2012, at 1:49 PM, Kevin Chang wrote:> Hello Everyone, > > > > I'm currently learning about quantile regressions. I've been using an > optimizer to compare with the rq() command for quantile regression. > > When I run the code, the results show that my coefficients are consistent > with rq(), but the intercept term can vary by a lot. > > I don't think my optimizer code is wrong and suspects it has something to do > with the starting values. > > The results seems very sensitive to different starting values and I don't > know how to make sense of it. > > > > Advice from the community would be greatly appreciated. > > > > Sincerely, > > > Kevin Chang > > > > ###################### CODE Below ########################### > > > > library(quantreg) > > data(engel) > > y<-cbind(engel[,2]) > > x<-cbind(rep(1,length(engel[,1])),engel[,1]) > > x1<-cbind(engel[,1]) > > nn<-nrow(engel) > > nn > > > > bhat.ls<-solve(t(x)%*%x)%*%t(x)%*%y > > #bhat.ls > > > > # QUANTILES > > quant=.25 > > > > fr.1=function(bhat.opt) > > { > > uu=y-x%*%bhat.opt > > sample.cond.quantile=quantile(uu,quant) > > w.less=rep(0,nn) > > for(ii in 1:nn){if(uu[ii]<sample.cond.quantile) w.less[ii]=1} > > > > sum((quant-1)*sum((y-x%*%bhat.opt)*w.less) #negative residuals > > +quant*sum((y-x%*%bhat.opt)*(1-w.less))) #positive residuals > > } > > start<-c(0,0) > > result=optim(start,fr.1) > > bhat.cond=result$par > > > > #Quantile Command Results > > fit.temp=rq(y~x1,tau=quant) > > fit.temp > > > > #OPTIMIZER Results > > bhat.cond > > > > #OLS Command Results > > mean=lm(y~x1) > > mean > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.
Apparently Analagous Threads
- quantile regression: plotting coefficients on only one variable (rq)
- applying quantile to a list using values of another object as probs
- problems with generation of quantiles under For ()
- package quantreg behaviour in weights in function rq,
- why won't rq draw lines?