如何通过svyVGAM的svy_vglm获取分类响应的加权百分比及置信区间

对于你用svy_vglm拟合的无协变量多分类模型，要提取各分类的加权百分比和置信区间，可通过以下两种方式实现：

方法一：从拟合模型中提取计算 #

1. 提取模型参数与方差矩阵 #

从cat_svy_fit中提取对数几率系数和对应的方差协方差矩阵：

# 提取B、C相对于参考类别A的对数几率系数
coefs <- coef(cat_svy_fit)
# 提取系数的方差协方差矩阵
vcov_mat <- vcov(cat_svy_fit)

2. 将对数几率转换为百分比 #

利用多分类模型的概率公式，把对数几率转换为总体百分比：

# 计算各分类的相对几率（A的几率固定为1）
odds <- c(1, exp(coefs))
# 转换为概率后转为百分比
pct <- (odds / sum(odds)) * 100
names(pct) <- c("A", "B", "C")

3. 用delta方法计算置信区间 #

由于概率是系数的非线性函数，需用delta方法推导标准误，再构建95%置信区间：

# 定义概率对系数的偏导数函数
delta_fun <- function(beta) {
  denom <- 1 + exp(beta[1]) + exp(beta[2])
  dA <- -(exp(beta[1]) + exp(beta[2])) / denom^2
  dB <- exp(beta[1]) * (1 + exp(beta[2])) / denom^2
  dC <- exp(beta[2]) * (1 + exp(beta[1])) / denom^2
  rbind(dA, dB, dC)
}

# 计算雅可比矩阵与概率的方差
jacobian <- delta_fun(coefs)
prob_vcov <- jacobian %*% vcov_mat %*% t(jacobian)
pct_se <- sqrt(diag(prob_vcov)) * 100

# 生成置信区间结果
pct_ci <- data.frame(
  category = names(pct),
  pct_estimate = round(pct, 2),
  lower_ci = round(pct - qnorm(0.975) * pct_se, 2),
  upper_ci = round(pct + qnorm(0.975) * pct_se, 2)
)

# 查看结果
print(pct_ci)

方法二：直接用svymean计算（更简便） #

当模型仅含截距项时，等价于直接计算加权比例，用svymean效率更高：

# 将分类变量转换为哑变量矩阵
dummy_vars <- model.matrix(~ outcome - 1, data = cat_data)
# 计算加权均值（即百分比）及置信区间
svy_pct <- svymean(dummy_vars, design = cat_svy_design)

# 转换为直观的数据框格式
result_df <- as.data.frame(svy_pct) %>% 
  mutate(
    category = c("A", "B", "C"),
    pct_estimate = round(mean * 100, 2),
    lower_ci = round((mean - 1.96 * SE) * 100, 2),
    upper_ci = round((mean + 1.96 * SE) * 100, 2)
  ) %>% 
  select(category, pct_estimate, lower_ci, upper_ci)

# 查看结果
print(result_df)

两种方法的结果完全一致，后者无需依赖模型拟合结果，操作更简洁。

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