我使用复杂的调查设计（来自加州健康访谈调查，CHIS）获得了调查数据，并试图获得四个相互排斥、共同详尽的性取向类别（异性恋、男同性恋/女同性恋、双性恋和异性恋）中每一个的比例估计值。其他）适用于使用
survey例如，在阿拉米达县，92% 的居民可能是异性恋； 3.5%，男同性恋或女同性恋； 3.5%，双性恋； 1% 为“其他”。目标是获得所有县的每个比例的置信区间（例如 92% [91.3%, 92.7%]）。
这是一个可重现的示例：
#Packages
require(survey)
require(dplyr)
require(tidyverse)
#Create data frame
# Set the number of trials and the probabilities for each outcome
set.seed(130)
probs <- c(.94, .03, .02, .01)
n <- 2000
sogi <- data.frame(t(rmultinom(n, 1, prob = probs)))
names(sogi) <- c("straight","gay","bisexual","other")
df_sogi <- sogi %>% mutate(sexual_orientation =
as.factor(case_when(
straight %in% 1 ~ 0,
gay %in% 1 ~ 1,
bisexual %in% 1 ~ 2,
other %in% 1 ~ 3))) %>%
mutate(pw = runif(n, 0.2, 3.5)) %>%
mutate(county = as.factor(rep(LETTERS[1:10], each = n/length(LETTERS[1:10])))) %>%
mutate(ind = row_number()) %>%
select(ind, county, sexual_orientation, pw)
#Examine raw frequencies
table(df_sogi$county, df_sogi$sexual_orientation)
#Make survey design and survey design replicate objects
df_svy_z1 <- svydesign(id=~ind, data = df_sogi, weights = ~pw, strata = ~county)
这是我得到的最接近的，它使用
svyglm#w/ svyglm
s <- svyglm(sexual_orientation == 0 ~ county, design = df_svy_z1, family="quasibinomial")
a <- predict(s, newdata=df_sogi, type="link", se.fit=TRUE)
b <- unique(confint(a))
exp(b)/(1+exp(b))
> exp(b)/(1+exp(b))
2.5 % 97.5 %
1 0.8525092 0.9479540
201 0.8641069 0.9584830
401 0.8766525 0.9607483
601 0.8552640 0.9504926
801 0.8743262 0.9592262
1001 0.9143716 0.9814168
1201 0.8587895 0.9514747
1401 0.9316797 0.9862013
1601 0.8660052 0.9532556
1801 0.8633378 0.9532417
然后我会编写一个函数，循环遍历
sexual_orientationsurvey更令人烦恼的是，我根本不确定通过这种方法获得的置信区间是否正确；它们看起来太宽了。
替代策略可以使用
svymeans <- svymean(~interaction(sexual_orientation, county), design=df_svy_z1)
x <- svycontrast(s, quote(log(`interaction(sexual_orientation, county)0.A` /
(1-`interaction(sexual_orientation, county)0.A`))))
y <- confint(x)
z <- 1/(1+exp(-y))
这里的两个大问题是（1）比例是单元格比例，而不是行比例——即频率表示为整个样本的比例，而不是给定县内样本的比例——以及（2）循环命名对象涉及将表达式
quote(log(`interaction(sexual_orientation, county)0.A`/(1-`interaction(sexual_orientation, county)0.A`))))s我的问题与此类似：
有没有办法使用 R 调查包来估计具有置信区间的多项比例？
但是@Thomas Lumley 教授慷慨地提供的答案仅涉及边际概率，即“整个州”而不是每个县的异性恋、男同性恋/女同性恋等个人的比例.
svyby首先获取按县划分的方法：
library(survey)
svmeans <- svyby(~ sexual_orientation,
design = df_svy_z1,
by = ~ county,
FUN = svymean)
# county sexual_orientation0 sexual_orientation1 sexual_orientation2 sexual_orientation3 se.sexual_orientation0 se.sexual_orientation1
# A A 0.9398657 0.02097910 0.023836687 0.015318555 0.01794317 0.01066598
# B B 0.9433422 0.02165069 0.033570977 0.001436098 0.01845545 0.01154185
# C C 0.9327511 0.03279607 0.022840837 0.011612004 0.01966055 0.01504690
# ....
然后，如果您同意使用标准错误的标准 Wald CI，请将其包装在
confint()中
svCIs <- confint(svmeans)
# 2.5 % 97.5 %
# A:sexual_orientation0 9.046977e-01 0.975033618
# B:sexual_orientation0 9.071702e-01 0.979514244
# C:sexual_orientation0 8.942171e-01 0.971285064
# D:sexual_orientation0 8.881002e-01 0.969468334
# ....
返回的对象有点笨拙，因此为了理解这一切，您可以美化对象并合并它们：
library(dplyr)
library(tidyr)
svmeans_long <- pivot_longer(svmeans[,-grep("se\\.", names(svmeans))], -county,
names_to = "sexual_orientation",
values_to = "mean")
svCIs_long <- data.frame(svCIs) %>%
tibble::rownames_to_column() %>%
separate_wider_delim(rowname, delim = ":",
names = c("county", "sexual_orientation"))
final_data <- left_join(svmeans_long, svCIs_long)
# county sexual_orientation mean X2.5.. X97.5..
# <chr> <chr> <dbl> <dbl> <dbl>
# 1 A sexual_orientation0 0.940 0.905 0.975
# 2 A sexual_orientation1 0.0210 0.0000742 0.0419
# 3 A sexual_orientation2 0.0238 0.00390 0.0438
# 4 A sexual_orientation3 0.0153 -0.00590 0.0365
# 5 B sexual_orientation0 0.943 0.907 0.980
# ....