2
$\begingroup$

Using the very nice survey package in R, I am having some trouble with the svyquantile function. From the function documentation it appeared to me that I could obtain confidence intervals for estimated quantiles in addition to SE. That seems to be the case if the survey object is constructed using svydesign(), but not if it is constructed using svrepdesign(). I used a bootstrap svrepdesign object and only get SE returned on my estimated quantiles. Maybe I have a conceptual problem?

From the example in the package documentation:

BRR method

data(scd)

repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1), + c(0,1,0,1,1,0))

scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)

Warning message: In svrepdesign.default(data = scd, type = "BRR", repweights = repweights) : No sampling weights provided: equal probability assumed

svyquantile(~arrests+alive, design=scdrep, quantile=0.5, interval.type="quantile")

Statistic: arrests alive q0.5 185 30 SE: arrests alive q0.5 15.02706 4.756574

confint(svyquantile(~arrests+alive, design=scdrep, quantile=0.5, interval.type="quantile"))

 2.5 % 97.5 %

>

confint(svyquantile(~arrests+alive, design=scdrep, quantile=0.5, confint=TRUE, interval.type="quantile"))

 2.5 % 97.5 %

Thanks in advance for any guidance.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
$\begingroup$

that does look like a bug. here's how to work around it

library(survey)
data(scd)
repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1), + c(0,1,0,1,1,0))
scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)

w <- svyquantile(~arrests, design=scdrep, quantile=0.5, interval.type="quantile")

plusorminus <- SE( w ) * qt( c( 0.025 , 0.975 ) , df = Inf )

# confidence interval for a single statistic
coef( w ) + plusorminus
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks, that certainly is a simple solution...wasn't too sure if basing on a t-distribution with the estimated se was what was supposed to happen in the code, or if it was the best approach. But should at least give give a valid CI estimate, I suppose. Thanks for your insight; I'll use it! $\endgroup$
    – pholck
    Commented Mar 30, 2016 at 2:31
  • $\begingroup$ look inside survey:::tconfint and also type out survey:::confint.sv and then hit tab twice quickly to read what is going on inside the code. good luck! $\endgroup$
    – Anthony Damico
    Commented Mar 30, 2016 at 12:13
1
$\begingroup$

Stumbled across a perhaps better work around to this apparent bug. Turns out that svyby objects are "correctly" handled by confint function, even if svyquantile objects are not. So just use a dummy factor variable in a svyby call with svyquuantile as the passed function:

library(survey) data(scd) repweights<-2*cbind(c(1,0,1,0,1,0), c(1,0,0,1,0,1), c(0,1,1,0,0,1), + c(0,1,0,1,1,0)) scdrep<-svrepdesign(data=scd, type="BRR", repweights=repweights)

quants <- svyby(~arrests,~factor(I(1>0)),svyquantile, quantile=0.5, interval.type="quantile", design=scdrep)

confint(quants)

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.