如何生成符合指定DGP和组内相关系数的聚类结构二元结果变量

我们可以通过多层logistic模型的框架来同时满足指定DGP和组内相关系数要求，具体实现逻辑和代码如下：

核心原理 #

二元结果的组内相关系数ICC定义为聚类层面的方差占总方差的比例：
ICC = σ²_u / (σ²_u + σ²_e)
其中对于logistic链接，个体层面残差服从固定方差为π²/3≈3.29的logistic分布，我们可以通过调整聚类（groupid）层面的随机效应方差σ²_u来达到目标ICC。

实现步骤 #

步骤1：参数设置 #

首先确定自定义参数：

• 目标ICC值（示例取0.05）
• 常数项a：即attrib1取参照组（attrib1_A=0、attrib1_B=0）时y=1的概率对应的logit转换值，比如你要参照组概率为20%，则a = qlogis(0.2)
• 给定的变量系数：attrib1_A=0.3、attrib1_B=-0.5
• 如果你的DGP里包含attrib2变量，也可以同步设置对应系数

步骤2：计算聚类随机效应方差 #

根据目标ICC反解得到groupid层面的随机效应标准差：
σ_u = sqrt( (ICC * π²/3) / (1 - ICC) )

步骤3：生成随机效应与结果变量 #

为每个groupid生成一个独立的随机效应，合并到原数据集后计算线性预测值，转换为概率后生成二元结果y。

完整R代码 #

# 加载依赖包
library(dplyr)

# ----------------------
# 自定义参数调整区
# ----------------------
target_icc <- 0.05          # 目标组内相关系数
ref_prob <- 0.2             # 参照组（所有虚拟变量取0）y=1的概率
coef_attrib1_A <- 0.3       # attrib1_A的系数
coef_attrib1_B <- -0.5      # attrib1_B的系数
# 若有attrib2变量可在此处添加对应系数，例如 coef_attrib2_A <- 0.2

# ----------------------
# 计算随机效应参数
# ----------------------
a <- qlogis(ref_prob)                       # 常数项logit转换
resid_var_logit <- pi^2 / 3                 # logistic分布固定个体残差方差
sigma_u_sq <- (target_icc * resid_var_logit) / (1 - target_icc)
sigma_u <- sqrt(sigma_u_sq)                 # group层面随机效应标准差

# ----------------------
# 生成结果变量y
# ----------------------
# 为每个groupid生成随机效应
group_random_eff <- d %>%
  distinct(groupid) %>%
  mutate(group_u = rnorm(n(), mean = 0, sd = sigma_u))

# 合并回原数据集生成y
d <- d %>%
  left_join(group_random_eff, by = "groupid") %>%
  mutate(
    # 若有attrib2变量，直接在xb计算中添加对应项即可
    xb = a + coef_attrib1_A * attrib1_A + coef_attrib1_B * attrib1_B + group_u,
    y_prob = plogis(xb),                     # 将线性预测转换为概率
    y = rbinom(n(), size = 1, prob = y_prob) # 生成二元结果
  )

效果验证 #

你可以通过拟合空多层模型验证生成的ICC是否符合预期：

library(lme4)
null_model <- glmer(y ~ 1 + (1|groupid), family = binomial, data = d)
var_group <- as.numeric(VarCorr(null_model)$groupid)
estimated_icc <- var_group / (var_group + pi^2/3)
cat("生成数据的ICC估计值为：", round(estimated_icc, 3))

注意事项 #

• 如果你后续分析用probit模型，只需把上述代码中的plogis替换为pnorm，个体残差方差替换为1即可
• 由于随机抽样的存在，单次生成的ICC可能和目标值有微小偏差，多次生成取符合要求的样本即可

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