svyby比例的置信区间

Question

是否有一个现有的函数可以从survey包中的svyby对象中创建比例的置信区间（在我的情况下，这是二元项的交叉表）。我经常比较不同组之间的比例，所以拥有一个可以提取置信区间的函数会非常方便（使用调查函数svyciprop而不是confint函数）。下面的示例显示了我想要实现的内容。

加载数据

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)

创建一个svyby对象，比较stype中变量"both"的比例。

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

是否有可能创建一个函数（例如byCI(b，method="likelihood")），它能够像使用confint(b)一样使用svyciprop来实现相同的功能？它基本上需要遍历svyby对象的每个级别并创建置信区间。到目前为止，我的尝试都没有成功。

也许还有另一种方法解决这个问题，但我喜欢使用svyby()，因为它快速且直观。

- maycobra

Answer 1

svyby()有一个vartype=参数，用于指定您希望如何指定抽样不确定性。使用vartype="ci"来获取置信区间，例如

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

很容易检查这个方法与手动逐级操作得到的结果相同，例如：

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

- Thomas Lumley

---

@maycobra 我之前没有意识到这是可能的 —— 这显然更加合理，你应该修改接受的答案 :) - Anthony Damico

Answer 2

有趣的是，这两个命令不应该返回相同的结果。第一个命令可能会抛出错误或警告：

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

你可能需要向卢姆利博士报告这个问题 - 在surveyby.R的第80行附近的代码可能需要进行轻微修改，以使svyciprop也能在svyby中工作。 但我可能忽略了某些细节（他可能在文档中有说明），所以在联系他之前请仔细阅读所有内容。 无论如何，这里有一个临时解决方案，可能可以解决你的问题。

# 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's 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 )

- Anthony Damico