imputation 数据框列表：使用调查包或 R 中的 MIcombine 计算分位数？

Question

我正在使用 R 并遵循 Anthony Damico(@AnthonyDamico) 的代码 (http://asdfree.com/survey-of-consumer-finances-scf.html) 来计算来自消费者金融调查的分位数。代码如下：

## Load necessary libraries (Note: should be installed in the OS)
library(mitools)    # allows analysis of multiply-imputed survey data
library(survey)     # load survey package (analyzes complex design surveys)
library(downloader) # downloads and then runs the source() function on scripts from github
library(foreign)    # load foreign package (converts data files into R)
library(Hmisc)      # load Hmisc package (loads a simple wtd.quantile function)
library(convey)
library(lodown)
library(devtools)
library(srvyr)

## Read SCF data: wave 2016, path.expand("~") is default working directory  
scf_imp <- readRDS("scf 2016.rds" )
scf_rw <- readRDS("scf 2016 rw.rds" )

# define which variables from the five imputed iterations to keep
vars.to.keep <- c( 'y1' , 'yy1' , 'wgt' , 'one' , 'houses', 'oresre', 'mrthel', 'liq', 'income', 'age' )

# restrict each `scf_imp#` data frame to only those variables
scf_imp[[1]] <- scf_imp[[1]][ , vars.to.keep ]
scf_imp[[2]] <- scf_imp[[2]][ , vars.to.keep ]
scf_imp[[3]] <- scf_imp[[3]][ , vars.to.keep ]
scf_imp[[4]] <- scf_imp[[4]][ , vars.to.keep ]
scf_imp[[5]] <- scf_imp[[5]][ , vars.to.keep ]

# clear up RAM
gc()

# turn off scientific notation in most output
options( scipen = 20 )

scf_design <- 
svrepdesign( 
weights = ~wgt , 
repweights = scf_rw[ , -1 ] , 
data = imputationList( scf_imp ) , 
scale = 1 ,
rscales = rep( 1 / 998 , 999 ) ,
mse = TRUE ,
type = "other" ,
combined.weights = TRUE
)

scf_design <-
update(
scf_design,
wealth = liq + houses + oresre - mrthel,
hwealth = houses + oresre,
htoinc  = (houses + oresre)/income,
ltoinc = mrthel / income,
ltv = mrthel + hwealth,
wtoinc = wealth + income
)

### gini coefficient (whole sample) ###
scf_design$designs <- lapply( scf_design$designs , convey_prep )
giniindex <- scf_MIcombine( with( scf_design , svygini( ~ wealth) ) )

在这段代码之后，我尝试了几种方法，包括1）：

q1.wealth <- scf_MIcombine(with(scf_design,svyquantile(~wealth, 0.01, 
method='constant',interval.type='quantile',se=TRUE)))
q5.wealth <- scf_MIcombine(with(scf_design,svyquantile(~wealth, 0.05, method='constant',interval.type='quantile',se=TRUE)))
q20.wealth <- scf_MIcombine(with(scf_design,svyquantile(~wealth, 0.20, method='constant',interval.type='quantile',se=TRUE)))
q40.wealth <- scf_MIcombine(with(scf_design,svyquantile(~wealth, 0.40, method='constant',interval.type='quantile',se=TRUE)))
q60.wealth <- scf_MIcombine(with(scf_design,svyquantile(~wealth, 0.60, method='constant',interval.type='quantile',se=TRUE)))
q80.wealth <- scf_MIcombine(with(scf_design,svyquantile(~wealth, 0.80, method='constant',interval.type='quantile',se=TRUE)))
q90.wealth <- scf_MIcombine(with(scf_design,svyquantile(~wealth, 0.90, method='constant',interval.type='quantile',se=TRUE)))
q95.wealth <- scf_MIcombine(with(scf_design,svyquantile(~wealth, 0.95, method='constant',interval.type='quantile',se=TRUE)))
q99.wealth <- scf_MIcombine(with(scf_design,svyquantile(~wealth, 0.99, method='constant',interval.type='quantile',se=TRUE)))
maxi.wealth <- scf_MIcombine(with(scf_design,svyquantile(~wealth, 1, method='constant',interval.type='quantile',se=TRUE)))  

和下面的方法2）：

quantile.wealth <- scf_MIcombine(with(scf_design,svyquantile(~networth, c(0.01,0.2,0.4,0.6,0.8,0.9,0.99,1), method='constant',interval.type='quantile',se=TRUE)))  

方法1）给我号码，但不准确。首先，它给出了以下警告信号：

Warning messages:
1: In svyquantile.svyrep.design(~networth, 1, method = "constant",  :
Not all replicate weight designs give valid standard errors for quantiles.
2: In svyquantile.svyrep.design(~networth, 1, method = "constant",  :
Not all replicate weight designs give valid standard errors for quantiles.
3: In svyquantile.svyrep.design(~networth, 1, method = "constant",  :
Not all replicate weight designs give valid standard errors for quantiles.
4: In svyquantile.svyrep.design(~networth, 1, method = "constant",  :
Not all replicate weight designs give valid standard errors for quantiles.
5: In svyquantile.svyrep.design(~networth, 1, method = "constant",  :
Not all replicate weight designs give valid standard errors for quantiles.

这说，标准误差不准确......所以我猜平均值也不值得信赖......（不确定）。我尝试使用方法 2)，这只是方法 1) 的一个小调整，使用调查包参数 c[ , ] 在单个命令中计算分位数。然后我收到以下错误。

Error in (1 + 1/m) * evar/vbar : non-conformable arrays
In addition: Warning messages:
1: In svyquantile.svyrep.design(~networth, c(0.01, 0.2, 0.4, 0.6, 0.8,  :
Not all replicate weight designs give valid standard errors for quantiles.
2: In svyquantile.svyrep.design(~networth, c(0.01, 0.2, 0.4, 0.6, 0.8,  :
Not all replicate weight designs give valid standard errors for quantiles.
3: In svyquantile.svyrep.design(~networth, c(0.01, 0.2, 0.4, 0.6, 0.8,  :
Not all replicate weight designs give valid standard errors for quantiles.
4: In svyquantile.svyrep.design(~networth, c(0.01, 0.2, 0.4, 0.6, 0.8,  :
Not all replicate weight designs give valid standard errors for quantiles.
5: In svyquantile.svyrep.design(~networth, c(0.01, 0.2, 0.4, 0.6, 0.8,  :
Not all replicate weight designs give valid standard errors for quantiles.

(ADDITION) 方法 3) 我使用调查包直接引用 (https://github.com/cran/survey/blob/master/man/svyquantile.Rd) 并尝试计算如下分位数：

q <- svyquantile(~networth,scf_design, c(.25,.5,.75),ci=TRUE)

但我收到以下错误。

Error in UseMethod("svyquantile", design) : 
no applicable method for 'svyquantile' applied to an object of class 
 "svyimputationList"

长话短说，我需要以正确的方式计算分位数和基尼系数，以便我可以将数字与模型进行比较。有没有办法更正代码，或者使用其他方法？

请记住，这是数据帧的插补列表，由 5 个数据帧组成。而且......值得注意的是，数据帧不包含相同数量的样本（行）。

Answer 1

我可以看到 3 个问题。第一个很简单。变量networth 不包含在vars.to.keep 中，因此以后无法使用，会引发错误。

其次，警告实际上不是问题。源代码中警告的位置是here。每当您使用interval.type='quantile' 调用svyquantile 并且使用type = 'other' 或任何使用折刀复制权重的类型构建调查设计时，都会引发此警告。您不应更改调查设计的类型。该警告提醒您使用interval.type='quantile' 产生的标准错误在jacknife 调查设计或异常引导设计中无效。 SCF 数据集uses bootstrap，但尚不清楚标准错误是否有效。或者，interval.type='probability' 不会引发错误，但可能不是您想要的。

q1.wealth <- scf_MIcombine(with(scf_design,
    svyquantile(~wealth,0.01,
                method='constant',
                interval.type='probability',
               se=TRUE, na.rm=T)));

第三，最后一个错误是讨厌的。我指的是这个

> quantile.wealth <- scf_MIcombine(with(
      scf_design,
      svyquantile(~networth,
                  c(0.01,0.2,0.4,0.6,0.8,0.9,0.99,1),
                  method='constant',
                  interval.type='quantile',
                  se=T,na.rm=T)));
# Error in (1 + 1/m) * evar/vbar : non-conformable arrays

似乎scf_MIcombine 和MIcombine 无法将svyquantile 调用的结果与多个分位数结合起来。有一个变通方法，svyby 函数不受此 bug 影响。

scf_MIcombine( with(scf_design,
    svyby(~networth,~as.factor(1),
          svyquantile,
          c(0.01,0.2,0.4,0.6,0.8,0.9,0.99,1),
          se=T, na.rm=T,
          method='constant',
          interval.type='quantile')))

在上面的示例中，函数svyby 用于计算调查数据子集的统计数据。在这种情况下，我使用了一个包含所有样本的虚假子集 ~as.factor(1)。