4

In my statistics class we use Stata and since I'm an R user I want to do the same things in R. I've gotten the right results but it seems like a somewhat awkward way of getting something as simple as confidence intervals.

Here's my crude solution:

library(quantreg)
na = round(runif(100, min=127, max=144))
f <- rq(na~1, tau=.5, data=ds)
s <- summary.rq(f, se="boot", R=1000)
coef(s)[1]
coef(s)[1]+ c(-1,1)*1.96*coef(s)[2]

I've also experimented a little at the boot package but I haven't gotten it to work:

library(boot)
b <- boot(na, function(w, i){ 
        rand_bootstrap_sample = w[i]
        f <- rq(rand_bootstrap_sample~1, tau=.5)
        return(coef(f)) 
        }, R=100)
boot.ci(b)

Gives an error:

Error in bca.ci(boot.out, conf, index[1L], L = L, t = t.o, t0 = t0.o, : estimated adjustment 'a' is NA

My questions:

  • What I wan't is to know if there is another better way of getting the confidence interval
  • why is the bootstrap code complaining?
lmo
  • 37,904
  • 9
  • 56
  • 69
Max Gordon
  • 5,367
  • 2
  • 44
  • 70

1 Answers1

4

Your example does not give an error message for me (Windows 7/64,R 2.14.2), so it could be a problem of random seeds. So if you post an example using some random method, better add a line set.seed; see example.

Note that the error message refers to the bca type of boot.ci; since this one often complains, deselect it by giving type explicitly.

I do not know exactly why you use the rather complex rq in the bootstrap. If you really wanted to profile rq, forget the simple example below, but please give some more details.

library(boot)
set.seed(4711)
na = round(runif(100, min=127, max=144))

b <- boot(na, function(w, i) median(w[i]), R=1000)
boot.ci(b,type=c("norm","basic","perc"))
Dieter Menne
  • 10,076
  • 44
  • 67
  • Thank you for your answer. It actually runs fine now on my laptop (R 2.15.0) - strange. The reason I'm using rq is simply due to that I'm trying to translate from a Stata excercise that we did. I'm mostly curious if there is a quantile regression that gives a bootstrap confidence interval – Max Gordon Apr 17 '12 at 07:31
  • Could you explain what you expect from rq(something~1)? In my view, this is a rather degenerate "regression". – Dieter Menne Apr 17 '12 at 08:58
  • As in your example I get the median but it's a stat course example that I just want to play around with - I agree that it barely makes sense... – Max Gordon Apr 17 '12 at 09:02
  • 3
    It seems that I get the bootstrap error when the resamplings are low and I use the bca instead of the "percentile", with 1000 resamplings I have no issue even using the bca but with 100 the error seems pretty predictable – Max Gordon Apr 17 '12 at 13:06