12

是否有一个现有的函数可以从一个svyby对象的比例创建置信区间（在我的例子中是包中二进制项目的交叉表survey）。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)svycipropsvyby

可能有另一种解决方法，但我喜欢使用svyby()它，因为它既快速又直观。

4

2 回答 2

12

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"))
于 2013-01-02T21:19:14.807 回答
2

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

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

您可能想提醒 Lumley 博士注意这个问题 - 第 80 行附近的代码surveyby.R可能会稍作修改以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 )
于 2013-01-01T18:51:35.340 回答