R 警告消息中的调查包

Question

我正在尝试使用调查包从随机调查中找到分类变量的均值。我在使用 svyby() 时遇到了一个问题，它输出： “警告信息： 在矩阵(1:(ns * reps), ncol = reps, nrow = ns, byrow = TRUE) 中： 数据长度[12]不是行数[5]的约数或倍数"

我的问题有两个：是什么导致这个问题以无法解释的表格格式重复出现结果？还有一个更理论的问题，如果我们只是获得手段，在创建设计元素之前是否存在子集问题（据我所知，它只会影响 st. 错误）。

这是我正在使用的代码，尝试使用完整的调查结果和子集结果：

#read-in zip file
library(survey)
library(rio)
td <- tempdir()
tf <- tempfile(tmpdir=td, fileext=".zip")
download.file("https://www.federalreserve.gov/consumerscommunities/files/SHED_public_use_data_2020_(CSV).zip", tf)
file_names <- unzip(tf, list=TRUE)
unzip(tf, exdir=td, overwrite=TRUE)
data <- import(file.path(td, file_names$Name[1]))

#remove weight NAs
data <- data[!is.na(data$weight_pop),]

#create subset
data.subset <- data[data$BK1 == "Yes" & data$afs == "Yes",]

#create svy designs
design <- svydesign(ids = ~CaseID
                    , weights = ~weight_pop
                    , na.rm = TRUE
                    , data = data) #full survey
design2 <- svydesign(ids = ~CaseID
                    , weights = ~weight_pop
                    , na.rm = TRUE
                    , data = data.subset) #subset of survey

svyby(~BK2_a,~race_5cat,design,svymean)
svyby(~BK2_a,~race_5cat,design2,svymean)

svyby(~BK2_b,~race_5cat,design,svymean)
svyby(~BK2_b,~race_5cat,design2,svymean)

svyby(~BK2_c,~race_5cat,design,svymean)
svyby(~BK2_c,~race_5cat,design2,svymean)

#verify observations
table(data$BK2_a,data$race_5cat)
table(data.subset$BK2_a,data.subset$race_5cat)

table(data$BK2_b,data$race_5cat)
table(data.subset$BK2_b,data.subset$race_5cat)

table(data$BK2_c,data$race_5cat)
table(data.subset$BK2_c,data.subset$race_5cat)

有些结果如我所愿，例如 (~BK2_b,~race_5cat,design2,svymean) 或 (~BK2_c,~race_5cat,design,svymean)，但其他结果会抛出警告消息并有无法解释的表格有重复的数字。

使用表函数查看观察结果似乎表明问题可能在于“拒绝”响应，但是当通过“数据 % mutate(BK2_a = str_replace(BK2_a, "Refused", replacement = NA_character_))" 然后在 svyby 函数中 na.rm.all=TRUE 它引入了 NAs 和 NaNs。我也尝试将列更改为因子，但它什么都不做。我是调查包的新手，非常感谢任何帮助。

Answer 1

首先，当您想限制到该子组时，我会使用子集（设计，BK1 ==“是”& afs ==“是”），而不是在子集数据上创建设计2。但这仅在您关心获得正确的标准错误时才重要，正如您所提到的。

其次，要解决您的问题，请尝试将响应变量 BK2_a、BK2_b 和 BK2_c 编码为因子。或者，如果您想要每个组的比例，您可以直接使用svytable()。见下文：

data <- data[!is.na(data$weight_pop),]
data$BK2_a_f <- factor(data$BK2_a)
data$BK2_b_f <- factor(data$BK2_b)
data$BK2_c_f <- factor(data$BK2_c)

#create subset
data.subset <- data[data$BK1 == "Yes" & data$afs == "Yes",]

#create svy designs
design <- svydesign(ids = ~CaseID, weights = ~weight_pop, 
                    na.rm = TRUE, data = data) #full survey
design2 <- svydesign(ids = ~CaseID, weights = ~weight_pop, 
                     na.rm = TRUE, data = data.subset) #subset of survey
design2_better <- subset(design, BK1 == "Yes" & afs == "Yes")


svyby(~BK2_a_f, ~race_5cat, design, svymean)
svyby(~BK2_a_f, ~race_5cat, design2, svymean)
svyby(~BK2_a_f, ~race_5cat, design2_better, svymean)

svyby(~BK2_b_f, ~race_5cat, design, svymean)
svyby(~BK2_b_f, ~race_5cat, design2, svymean)
svyby(~BK2_b_f, ~race_5cat, design2_better, svymean)

svyby(~BK2_c_f, ~race_5cat, design, svymean)
svyby(~BK2_c_f, ~race_5cat, design2, svymean)
svyby(~BK2_c_f, ~race_5cat, design2_better, svymean)

# Alternative approach to getting proportions if you don't care about SEs
prop.table(svytable(~race_5cat + BK2_a, design), 1)
prop.table(svytable(~race_5cat + BK2_b, design), 1)
prop.table(svytable(~race_5cat + BK2_c, design), 1)