Svyby confidence intervals

Question

I was wondering if anyone came up with a function that can create confidence intervals from an svyby object for proportions (in my case, crosstabs for the binary element in the survey package). I often compare proportions between groups, and it would be very convenient to have a function that can extract confidence intervals (with the svyciprop polling function, not with the config). The example below shows what I would like to achieve.

Data loading

library(survey) library(weights) data(api) apiclus1$both<-dummify(apiclus1$both)[,1]#Create dummy variable dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc) 

Create an svyby object that compares the proportion of the variable "both" by type

 b<-svyby(~both, ~stype, dclus1, svymean) confint(b)#This works, but svyciprop is best in other cases, especially when proportion is close to 0 or 1 svyciprop(b)#This requires that you specify each level and a design object 

Is it possible to create a function (for example, using CI (b, method = "likelihood") that achieves the same as confint (b), but using svyciprop? Basically this should go through each level of the svyby object and create a trust interval My attempts have been unsuccessful so far.

There may be another way, but I like to use svyby () as it is fast and intuitive.

Thanks in advance.

Answer 1

svyby () has a vartype = argument to indicate how you want to determine the sampling uncertainty. Use vartype = "ci" to get confidence intervals, for example

 svyby(~I(ell>0),~stype,design=dclus1, svyciprop,vartype="ci",method="beta") 

It is easy to verify that this gives the same thing as each level manually, for example,

 confint(svyciprop(~I(ell>0), design=subset(dclus1,stype=="E"),method="beta")) 

Answer 2

interesting .. these two commands should not give the same result .. the first should probably throw an error or warning:

 svyby( ~both , ~stype , dclus1 , svyciprop , method = 'likelihood' ) svyby( ~both , ~stype , dclus1 , svymean ) 

you can warn Dr. Lumley on this issue - code close to line 80 of surveyby.R may be slightly modified to get svyciprop working inside svyby too .. but I can ignore something (and he may have noticed this somewhere in the documentation), so be sure to carefully read everything before contacting him about it

anyway, here is a workaround that might solve your problem

 # create a svyby-like function specific for svyciprop svyciby <- function( formula , by , design , method = 'likelihood' , df = degf( design ) ){ # steal a bunch of code from the survey package source # stored in surveyby.R.. byfactors <- model.frame( by , model.frame( design ) , na.action = na.pass ) byfactor <- do.call( "interaction" , byfactors ) uniquelevels <- sort( unique( byfactor ) ) uniques <- match( uniquelevels , byfactor ) # note: this may not work for all types.. # i only tested it out on your example. # run the svyciprop() function on every unique combo all.cis <- lapply( uniques , function( i ){ svyciprop( formula , design[ byfactor %in% byfactor[i] ] , method = method , df = df ) } ) # transpose the svyciprop confidence intervals t.cis <- t( sapply( all.cis , attr , "ci" ) ) # tack on the names dimnames( t.cis )[[1]] <- as.character( sort( unique( byfactor ) ) ) # return the results t.cis } # test out the results svyciby( ~both , ~stype , dclus1 , method = 'likelihood' ) # pretty close to your b, but not exact (as expected) confint(b) # and this one does match (as it should) svyciby( ~both , ~stype , dclus1 , method = 'mean' , df = Inf )