[R] Fwd: Making rq and bootcov play nice

[Frank E Harrell Jr]

roger koenker wrote:
> 
> John,
> 
> You can make a local version of bootcov which either:
> 
>     deletes these arguments from the call to fitter, or
>     modify the switch statement to include rq.fit,
> 
> the latter would need to also modify rq() to return a fitFunction
> component, so the first option is simpler.   One of these days I'll
> incorporate clustered se's into summary.rq, but meanwhile
> this seems to be a good alternative.
> 
> Roger

With some luck I'll have a Design package interface to quantreg for the 
case where one value of tau is requested, along with a corresponding 
addition to bootcov.  With some luck I'll have this by tomorrow.  This 
would give you nomograms, effect plots, LaTeX expression of model fits, etc.

As a side issue, as I'm just getting into quantile regression, I've done 
some simulations to better understand the operating characteristics. 
I've seen that the efficiency of rq estimates at fixed covariate values, 
when you have a normal distribution, relative to least squares 
estimates, is about the same as the efficiency of the sample median 
compared to the sample mean, or a bit worse.  Code is below.  As is 
always the case, if you are willing to assume a lot (e.g., normality 
with constant variance) you can do better.  On the other hand, quantile 
regression is more invariant (but not perfectly so) to transformations 
of Y as studied using the 2nd bit of code below.

On another issue, I had trouble installing quantreg on one Ubuntu 
machine, as I've forgotten which LAPACK/ATLAS ubuntu/debian package 
needs to be installed.  I got an error related to the absence of 
llapack.  If anyone has a pointer I'd appreciate getting it.

This gives me an opportunity to thank Roger for his wonderful work in 
this area, and his excellent vignette that goes with quantreg.

Frank Harrell


n <- 50
m <- 1000
a <- b <- aq <- bq <- numeric(m)

for(i in 1:m)
   {
     cat(i,'\r')
     x1 <- rnorm(n)
     y <- x1 + rnorm(n)/4
     f <- lm(y ~ x1)
     k <- coef(f)
     a[i] <- k[1]; b[i] <- k[2]

     f <- rq(y ~ x1)
     k <- coef(f)
     aq[i] <- k[1]; bq[i] <- k[2]
   }
cat('\n')
z <- c(mean(a^2), mean(aq^2))
zq <- c(mean((b-1)^2), mean((bq-1)^2))
c(z, z[1]/z[2])
c(zq, zq[1]/zq[2])

-------------------

n <- 50
par(mfrow=c(5,4), mar=c(2,1,1,1))
for(i in 1:20)
   {
     x1 <- rnorm(n)
     y <- exp(x1 + rnorm(n)/4)
     range(y)

     f <- lm(y ~ exp(x1))
     x1s <- seq(-2.25, 2.25, length=100)
     d <- data.frame(x1=x1s)
     plot(x1s, log(predict(f,d)), type='n', ylab='',
          xlim=c(-2.5,2.5), ylim=c(-4, 3.5))
     lines(x1s, log(predict(f,d)), col='blue')
     k <- coef(f)
     a <- k[1];  b <- k[2]
     c(a^2, 2*a*b, b^2)
     f <- lm(log(y) ~ x1)
     lines(x1s, predict(f,d), col='green')


     f <- rq(y ~ exp(x1))
#    plot(x1s, log(predict(f,d)), type='n', ylab='')
     lines(x1s, log(predict(f,d)), col='pink')
     k <- coef(f)
     a <- k[1];  b <- k[2]
     c(a^2, 2*a*b, b^2)
     f <- rq(log(y) ~ x1)
     lines(x1s, predict(f,d), col='red')
   }





> 
> 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 Jul 23, 2009, at 8:10 PM, John Gardner wrote:
> 
>> I have a quick question, and I apologize in advance if, in asking, I
>> expose my woeful ignorance of R and its packages. I am trying to use
>> the bootcov function to estimate the standard errors for some
>> regression quantiles using a cluster bootstrap. However, it seems that
>> bootcov passes arguments that rq.fit doesn't like, preventing the
>> command from executing. Here is an example:
>>
>> e<-bootcov(rq(y~x),clust,B=10,fitter=rq.fit)
>>
>> (where clust is my clustering variable) results in
>>
>> Error in rq.fit.br(x, y, tau = tau, ...) :
>> unused argument(s) (maxit = 15, penalty.matrix = NULL)
>>
>> In contrast, the lm.fit function seems to just ignore these arguments,
>> resulting in the following warning:
>>
>> 10: In fitter(X[obs, , drop = FALSE], Y[obs, , drop = FALSE], maxit = 
>> maxit,  :
>> extra arguments maxitpenalty.matrix are just disregarded.
>>
>> Is there a way that I can either (a) modify bootcov so that it doesn't
>> pass these arguments or (b) modify rq so that it ignores them?
>>
>> Thanks in advance,
>> John Gardner
>>
>> ______________________________________________
>> 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.
> 
> ______________________________________________
> 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.
> 


-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University