SAS rand与R rgamma/rpois生成负二项数据的差异排查

我采用Molenbergh、Kalema和Iddi的第二种组合模型，通过Gamma随机效应引入相关性与离散度，结合Poisson分布构建Poisson-Gamma混合模型模拟相关负二项数据，再使用SAS PROC GENMOD进行经验功效计算。但发现R生成的数据计算出的95%置信区间更宽（1500次重复中R均值为1.01，SAS为0.9），导致功效更低（目标样本量下R为81%，SAS为89%）。已确认数据生成公式一致，Gamma随机效应及负二项数据的均值、方差收敛结果相近，且功效计算采用同一SAS代码，推测差异源于数据生成环节，寻求排查方法以确定准确性并缩小结果差异。

%let shapex = %sysevalf(&a1 - (&min_a * &rho1));
%let shapey = %sysevalf(&a2 - (&min_a * &rho1));
%let shapez = %sysevalf(&min_a * &rho1);
* Generate independent gammas;
data gamma_gen;
    call streaminit(12345); 
    do i = 1 to &n;
        x = rand('Gamma', &shapex, 1);
        y = rand('Gamma', &shapey, 1);
        z = rand('Gamma', &shapez, 1);
    output;
end;
run;
* Obtain correlated gammas;
data corr_gamma;
    set gamma_gen;
    G1 = &b1 * (x + z);
    G2 = &b2 * (y + z);
run;

* Correlated NB;
data corr_nb;
    set corr_gamma;
    lambda1 = &m1 * G1;
    lambda2 = &m2 * G2;
    call streaminit(12345);
    V1 = rand('Poisson', lambda1);
    V2 = rand('Poisson', lambda2);
    drop i x y z G1 G2 lambda1 lambda2;
run;

rG = function(N, a, b, rho) {
  rho1 = rho * sqrt(max(a)/min(a))
  if (rho1 > 1) return(NA)
  else {
    X = rgamma(n = N, shape = a[1] - min(a) * rho1, scale = 1)
    Y = rgamma(n = N, shape = a[2] - min(a) * rho1, scale = 1)
    Z = rgamma(n = N, shape = min(a) * rho1, scale = 1)
    cbind(b[1] * (X + Z), b[2] * (Y + Z))
  }
}
rNB = function(N, k, m, r) {
  rho = getrho(k,m,r) # this value is the exact same between both programs
  gamma_sample = rG(N, 1/k, k, rho)
  if (!any(is.na(gamma_sample))){
    return(cbind(rpois(n = N, lambda = m[1] * gamma_sample[,1]),
                 rpois(n = N, lambda = m[2] * gamma_sample[,2])))
  }
}

• 强制随机数序列完全一致：SAS中两次调用call streaminit(12345)分别生成Gamma和Poisson数据；R需在开头显式设置set.seed(12345)，且确保Gamma生成与Poisson生成的顺序和SAS严格对应——先一次性生成所有Gamma样本，再按相同顺序生成Poisson数据。可仅生成10个样本，对比两边V1、V2的具体数值，这是验证数据一致性最直接的方法。
• 核对Gamma分布的所有参数：将SAS宏中的&a1、&a2、&min_a、&rho1与RrG函数的a、rho1逐一对比。比如用SAS的%put输出&shapex、&shapey、&shapez，在R中打印a[1]-min(a)*rho1等值，确保每一个参数完全匹配。
• 对比Gamma变量的分布特征：各生成10000个样本，计算SAS的G1、G2与R生成的对应变量的均值、方差、协方差（重点是G1和G2的相关性），确认两者的分布特征几乎无差异（浮点精度差异除外）。
• 验证Poisson生成环节的一致性：把SAS生成的G1、G2导出为CSV，导入R后用相同种子生成Poisson数据，看结果是否与SAS的V1、V2一致；反之，将R生成的Gamma样本导入SAS生成Poisson，对比结果，定位差异是否出在Poisson生成步骤。
• 排查浮点精度累积差异：SAS与R的浮点运算精度可能存在细微差别，可对比两边生成的G1、G2的小数点后多位数值，看是否因精度差异导致后续Poisson结果不同。
• 确认参数传递的准确性：检查SAS宏中的&b1、&b2、&m1、&m2是否与R的b、m参数完全对应，比如&b1是否等于R中的k[1]，&m1是否等于m[1]，打印参数值逐一验证。

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