如何在R中结合复杂抽样设计获取加权多列相关矩阵？

目前R中没有直接支持survey对象一键生成完整加权Polychoric相关矩阵的现成函数，不过可以通过手动循环或结合bootstrap的方式实现，以下是两种可行方案：

方案1：手动循环计算变量对相关并构建矩阵 #

这种方法直接从survey对象中提取变量和抽样权重，循环计算每对有序变量的Polychoric相关，最后拼接成完整矩阵。

步骤示例： #

1. 加载所需包并提取survey对象中的数据与权重

library(survey)
library(wCorr) # 或psych包，二选一即可

# 假设你的复杂抽样设计对象为svy_design
target_vars <- c("ordered_var1", "ordered_var2", "ordered_var3") # 替换为你的有序变量名
survey_data <- svy_design$variables[, target_vars, drop = FALSE]
sampling_weights <- svy_design$weights

2. 初始化空矩阵并循环计算

var_count <- length(target_vars)
poly_corr_mat <- matrix(NA, nrow = var_count, ncol = var_count)
rownames(poly_corr_mat) <- colnames(poly_corr_mat) <- target_vars

# 循环计算每对变量的Polychoric相关
for (i in 1:var_count) {
  for (j in i:var_count) {
    # 使用wCorr包的weightedCorr函数
    corr_val <- weightedCorr(survey_data[[i]], survey_data[[j]], 
                            method = "polychoric", weights = sampling_weights)
    # 若用psych包，替换为：
    # corr_val <- psych::polychoric(cbind(survey_data[[i]], survey_data[[j]]), weights = sampling_weights)$rho
    
    poly_corr_mat[i, j] <- corr_val
    poly_corr_mat[j, i] <- corr_val # 利用对称矩阵特性填充
  }
}

# 查看结果
print(poly_corr_mat)

方案2：结合Bootstrap考虑复杂抽样的设计效应 #

如果需要考虑复杂抽样中的分层、聚类等设计效应（仅用权重无法完全体现），可以通过survey包的svybootstrap函数对抽样设计进行bootstrap重复抽样，再计算每轮样本的Polychoric相关矩阵，最终取均值作为结果。

步骤示例： #

library(survey)
library(wCorr)

# 定义计算Polychoric相关矩阵的函数
calc_poly_corr <- function(design, vars) {
  data_subset <- design$variables[, vars, drop = FALSE]
  wts <- design$weights
  n <- length(vars)
  corr_mat <- matrix(NA, n, n)
  rownames(corr_mat) <- colnames(corr_mat) <- vars
  
  for (i in 1:n) {
    for (j in i:n) {
      corr_mat[i,j] <- corr_mat[j,i] <- weightedCorr(data_subset[[i]], data_subset[[j]], 
                                                    method = "polychoric", weights = wts)
    }
  }
  return(corr_mat)
}

# 执行bootstrap（replicates可根据需求调整，一般至少1000次）
boot_replicates <- svybootstrap(svy_design, calc_poly_corr, vars = target_vars, replicates = 1000)

# 提取平均相关矩阵与标准差矩阵
mean_poly_corr <- apply(boot_replicates$replicates, c(1,2), mean)
sd_poly_corr <- apply(boot_replicates$replicates, c(1,2), sd)

# 查看结果
print(mean_poly_corr)
print(sd_poly_corr)

注意事项： #

• 确保用于计算的变量是有序分类变量，Polychoric相关仅适用于这类变量。
• 若使用psych包的polychoric函数，需注意该函数对缺失值的处理默认是删除含缺失的观测，可通过na.rm参数调整。

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