在R中计算协变量分层内的加权患病率比及置信区间

针对加权调查数据中，需要计算各协变量分层内不同时间段的避孕套使用患病率比及95%置信区间的需求，以下提供两种可行方法：

方法1：分层计算患病率后用delta法推导置信区间 #

先通过svyby按协变量分层，得到各层内不同年份的患病率及标准误，再利用delta法计算患病率比的置信区间。

library(survey)
library(tidyr)
set.seed(123)

# 示例数据集
dat <- data.frame(
  year = factor(sample(c("2019", "2020"), 100, replace = TRUE)),
  sv_weight = sample(1:150, 100, replace = TRUE),
  strat = sample(531011:532149, 100, replace = TRUE),
  condom_use = sample(c(0, 1), 100, replace = TRUE),
  age = factor(sample(c("18-29", "30-39", "40-49"), 100, replace = TRUE)),
  rurality = factor(sample(c("Rural", "Urban"), 100, replace = TRUE)),
  sexual_orientation =  factor(sample(c("Straight", "Not straight"), 100, replace = TRUE)))

# 构建调查设计
options(survey.lonely.psu = "adjust")
design <- svydesign(data = dat,
                    id = ~1,
                    strata = ~strat,
                    weights = ~sv_weight)

# 定义函数：计算患病率及标准误
get_prevalence <- function(d) {
  prop <- svyciprop(~condom_use, d, method = "logit")
  data.frame(prev = coef(prop), se = SE(prop))
}

# 年龄分层的患病率比及CI
age_prev <- svyby(~year, ~age, design, FUN = get_prevalence, keep.var = FALSE)
age_pr <- age_prev %>%
  pivot_wider(names_from = year, values_from = c(prev, se)) %>%
  mutate(
    pr = prev_2020 / prev_2019,
    # delta法计算对数尺度的标准误
    se_log_pr = sqrt((se_2020/prev_2020)^2 + (se_2019/prev_2019)^2),
    lci = exp(log(pr) - 1.96 * se_log_pr),
    uci = exp(log(pr) + 1.96 * se_log_pr)
  ) %>%
  select(age, pr, lci, uci)
print(age_pr)

# 城乡分层的患病率比及CI
rurality_prev <- svyby(~year, ~rurality, design, FUN = get_prevalence, keep.var = FALSE)
rurality_pr <- rurality_prev %>%
  pivot_wider(names_from = year, values_from = c(prev, se)) %>%
  mutate(
    pr = prev_2020 / prev_2019,
    se_log_pr = sqrt((se_2020/prev_2020)^2 + (se_2019/prev_2019)^2),
    lci = exp(log(pr) - 1.96 * se_log_pr),
    uci = exp(log(pr) + 1.96 * se_log_pr)
  ) %>%
  select(rurality, pr, lci, uci)
print(rurality_pr)

# 性取向分层的患病率比及CI
so_prev <- svyby(~year, ~sexual_orientation, design, FUN = get_prevalence, keep.var = FALSE)
so_pr <- so_prev %>%
  pivot_wider(names_from = year, values_from = c(prev, se)) %>%
  mutate(
    pr = prev_2020 / prev_2019,
    se_log_pr = sqrt((se_2020/prev_2020)^2 + (se_2019/prev_2019)^2),
    lci = exp(log(pr) - 1.96 * se_log_pr),
    uci = exp(log(pr) + 1.96 * se_log_pr)
  ) %>%
  select(sexual_orientation, pr, lci, uci)
print(so_pr)

说明：svyciprop用logit方法计算患病率的置信区间，保证比例在0-1范围内；delta法通过对数变换将比值的方差转换为两个患病率对数方差之和，再转换回原始尺度得到置信区间。

方法2：分层log-binomial模型直接输出结果 #

通过拟合分层的log-binomial模型，直接得到各层内2020年相对2019年的患病率比及95%置信区间，步骤更简洁。

library(survey)
set.seed(123)

# 复用之前的数据集和调查设计（无需重复定义）

# 年龄分层的模型结果
age_pr <- svyby(~year, ~age, design, FUN = function(d) {
  # 拟合log-binomial模型
  mod <- svyglm(condom_use ~ year, d, family = quasibinomial(link = "log"))
  coef_year <- coef(mod)["year2020"]
  se_year <- SE(mod)["year2020"]
  # 转换为患病率比及CI
  data.frame(
    pr = exp(coef_year),
    lci = exp(coef_year - 1.96 * se_year),
    uci = exp(coef_year + 1.96 * se_year)
  )
})
print(age_pr)

# 城乡分层的模型结果
rurality_pr <- svyby(~year, ~rurality, design, FUN = function(d) {
  mod <- svyglm(condom_use ~ year, d, family = quasibinomial(link = "log"))
  coef_year <- coef(mod)["year2020"]
  se_year <- SE(mod)["year2020"]
  data.frame(
    pr = exp(coef_year),
    lci = exp(coef_year - 1.96 * se_year),
    uci = exp(coef_year + 1.96 * se_year)
  )
})
print(rurality_pr)

# 性取向分层的模型结果
so_pr <- svyby(~year, ~sexual_orientation, design, FUN = function(d) {
  mod <- svyglm(condom_use ~ year, d, family = quasibinomial(link = "log"))
  coef_year <- coef(mod)["year2020"]
  se_year <- SE(mod)["year2020"]
  data.frame(
    pr = exp(coef_year),
    lci = exp(coef_year - 1.96 * se_year),
    uci = exp(coef_year + 1.96 * se_year)
  )
})
print(so_pr)

说明：log-binomial模型的对数链接函数使得暴露变量（year2020）的系数指数即为患病率比；使用quasibinomial可处理数据可能存在的过度离散问题，若数据无过度离散，替换为binomial即可。

内容的提问来源于stack exchange，提问作者notasfarwest