基于给定相关矩阵生成二元分布二值变量的方法咨询

问题背景 #

需要生成取值仅为0或1的两个变量(x,y)，二者需满足预设的皮尔逊相关系数。直接对服从指定相关的二元正态变量做截断转二值的方法存在偏差：原本预设0.5的相关系数，截断后实际相关仅为0.3左右，不符合要求。

可行解决方案 #

方案1：使用专门处理相关二分类数据的bindata包（推荐） #

该包可以直接指定变量的边际取1概率、目标相关系数，直接生成符合要求的二值样本，使用逻辑和mvrnorm类似，符合你提到的「通过参数指定分布」的使用习惯。
代码示例：

# 首次使用需安装包
# install.packages("bindata")
library(bindata)
library(tidyverse)

# 设定参数：x取1的概率、y取1的概率，可按需调整
marginal_prob <- c(0.5, 0.5)
# 设定目标相关矩阵
target_cor <- matrix(c(1, 0.5, 0.5, 1), nrow = 2)

# 生成100个样本，固定种子可复现结果
set.seed(123)
df_bin <- rmvbin(
  n = 100,
  margprob = marginal_prob,
  bincorr = target_cor
) %>% as.data.frame() %>% rename(x = V1, y = V2)

# 验证相关系数，输出结果会非常接近预设的0.5
cor(df_bin)

方案2：校正截断正态法的相关系数 #

如果需要沿用正态截断的逻辑，需要先将二值变量的目标相关转换为对应的四分相关系数，再用转换后的相关生成正态样本，最后截断即可得到符合要求的结果。
代码示例：

# 首次使用需安装包
# install.packages(c("MASS", "psych"))
library(MASS)
library(psych)
library(tidyverse)

# 目标二值相关
target_bin_cor <- 0.5
# 转换为底层正态分布需要的四分相关，可按需调整p参数修改x、y的取1概率
tetra_cor <- tetrachoric(
  matrix(c(0.5*(1-target_bin_cor), 0.5*target_bin_cor, 0.5*target_bin_cor, 0.5*(1-target_bin_cor)), 2),
  p = c(0.5,0.5)
)$rho

# 生成二元正态样本
df_norm <- mvrnorm(
  100, mu = c(0,0),
  Sigma = matrix(c(1, tetra_cor, tetra_cor, 1), nrow = 2),
  empirical = TRUE
) %>% as.data.frame() %>% rename(x=V1, y=V2)

# 截断为二值变量
df_bin <- df_norm %>% mutate(
  x = ifelse(x < 0, 0, 1),
  y = ifelse(y < 0, 0, 1)
)

# 验证相关系数
cor(df_bin)

注意事项 #

• 二值变量的相关系数存在理论上限，上限由两个变量的边际取1概率决定，如果设置的目标相关超过理论上限，函数会报错提示，此时需要调整边际概率或者降低目标相关。
• 如果需要生成的二值变量边际取1概率不是0.5，只需要修改对应参数即可，两种方案都支持自定义边际概率。

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