从频率学派Logistic回归转贝叶斯Logistic回归的代码问题求助

问题背景 #

使用proc surveylogistic执行频率学派Logistic回归时，因变量部分类别样本量过低出现完全分离现象，计划切换到贝叶斯Logistic回归以稳定模型，同时保留原调查数据的抽样权重（主权重+BRR重复权重），目标输出ParameterEstimates和OddsRatios表，但编写proc mcmc代码时出现报错。

原频率学派代码 #

proc surveylogistic data=final_cross varmethod=BRR (fay=0.3);     
ods select  ParameterEstimates   OddsRatios;     
weight A_PWGT;     
repweights A_PWGT1 - A_PWGT100;     
class ecignicotinelevel3_7 (ref="3") age2_7 (ref="1") gender_7 (ref="2") race_7 (ref="1") education_7 (ref="1") income_7 (ref="1") ecigdaily_7 (ref="1") cigformer_7 (ref="1")   othertob_7 (ref="1")  ecig_device2_7 (ref="1")/ param = ref;      
model ecignicotinelevel3_7=age2_7 gender_7 race_7 education_7 income_7 cigformer_7 othertob_7 ecigdaily_7 ecig_device2_7 ecig_time_7 ecig_puff_7/link=glogit;     
where dual_7=3 and (ecignicotinelevel3_7=3 | ecignicotinelevel3_7=4);     
run; 

尝试的贝叶斯代码（报错版本） #

proc mcmc data=final_cross outpost=posterior seed=1234 nbi=1000 nmc=2000; 
ods select ParameterEstimates OddsRatios; 
weight A_PWGT; 
repweights A_PWGT1 - A_PWGT100; 
model ecignicotinelevel3_7(event='3') = age2_7 gender_7 race_7 education_7 income_7cigformer_7 othertob_7 ecigdaily_7 ecig_device2_7ecig_time_7 ecig_puff_7 / link=logit; 
prior all ~ normal(0, var=1000);  
where dual_7=3 and (ecignicotinelevel3_7=3 | ecignicotinelevel3_7=4);
run; 

报错信息（翻译后） #

559     weight A_PWGT;
        ------
        180
560     repweights A_PWGT1 - A_PWGT100;
        ----------
        180
561     model ecignicotinelevel3_7(event='3') = age2_7 gender_7 race_7 education_7 income_7
                                  -
                                  22
                                  76
561! cigformer_7
ERROR 180-322: 语句无效或使用顺序错误。
ERROR 22-322: 语法错误，预期以下内容之一：名称、-、:、_ALL_、_CHARACTER_、_CHAR_、_NUMERIC_、~。
ERROR 76-322: 语法错误，语句将被忽略。

错误原因分析 #

1. 语法错误：模型语句中变量名连写，income_7cigformer_7应为income_7 cigformer_7，ecig_device2_7ecig_time_7应为ecig_device2_7 ecig_time_7，导致22/76号语法报错。
2. 过程不支持的语句：proc mcmc不支持weight和repweights语句，这是180号报错的核心原因——贝叶斯MCMC过程处理复杂抽样权重的逻辑与频率学派调查过程完全不同，不能直接照搬surveylogistic的权重声明方式。

修正方案 #

步骤1：修复语法错误 #

先修正模型语句中的变量分隔问题，修正后的模型行如下：

model ecignicotinelevel3_7(event='3') = age2_7 gender_7 race_7 education_7 income_7 cigformer_7 othertob_7 ecigdaily_7 ecig_device2_7 ecig_time_7 ecig_puff_7 / link=logit;

步骤2：处理抽样权重与重复权重 #

针对贝叶斯框架下的复杂调查数据，提供两种可行思路：

思路1：加权似然法（整合主权重）

手动将主权重A_PWGT整合到Logistic似然函数中，通过ll语句自定义似然。同时需要声明分类变量（与原代码一致）：

proc mcmc data=final_cross outpost=posterior seed=1234 nbi=1000 nmc=2000;
  ods select ParameterEstimates OddsRatios;
  /* 声明分类变量，与原代码保持一致的参考水平 */
  class ecignicotinelevel3_7 (ref="3") age2_7 (ref="1") gender_7 (ref="2") race_7 (ref="1") 
        education_7 (ref="1") income_7 (ref="1") ecigdaily_7 (ref="1") cigformer_7 (ref="1") 
        othertob_7 (ref="1") ecig_device2_7 (ref="1")/ param = ref;
  /* 定义弱信息先验，缓解完全分离问题 */
  prior all ~ normal(0, var=1000);
  /* 构造线性预测器 */
  eta = intercept + age2_7 + gender_7 + race_7 + education_7 + income_7 + cigformer_7 + 
        othertob_7 + ecigdaily_7 + ecig_device2_7 + ecig_time_7 + ecig_puff_7;
  /* 计算加权Logistic似然 */
  ll = sum(A_PWGT * ( (ecignicotinelevel3_7='3') * log(invlogit(eta)) + 
                      (ecignicotinelevel3_7='4') * log(1 - invlogit(eta)) ));
  /* 指定自定义似然 */
  model general(ll);
  where dual_7=3 and (ecignicotinelevel3_7=3 | ecignicotinelevel3_7=4);
run;

注意：该方法仅整合了主抽样权重，未处理BRR重复权重。若需基于重复权重估计方差，需额外编写循环对每个重复权重数据集运行MCMC，再合并结果。

思路2：使用贝叶斯调查专用过程

SAS的proc bayes（部分版本支持）提供了调查数据权重的处理能力，可以尝试替换proc mcmc，语法更接近surveylogistic：

proc bayes data=final_cross seed=1234 nbi=1000 nmc=2000;
  ods select ParameterEstimates OddsRatios;
  class ecignicotinelevel3_7 (ref="3") age2_7 (ref="1") gender_7 (ref="2") race_7 (ref="1") 
        education_7 (ref="1") income_7 (ref="1") ecigdaily_7 (ref="1") cigformer_7 (ref="1") 
        othertob_7 (ref="1") ecig_device2_7 (ref="1")/ param = ref;
  model ecignicotinelevel3_7(event='3') = age2_7 gender_7 race_7 education_7 income_7 cigformer_7 
        othertob_7 ecigdaily_7 ecig_device2_7 ecig_time_7 ecig_puff_7 / link=logit;
  weight A_PWGT;
  /* 弱信息先验缓解完全分离 */
  prior all ~ normal(0, var=1000);
  where dual_7=3 and (ecignicotinelevel3_7=3 | ecignicotinelevel3_7=4);
run;

补充说明 #

• 因变量已筛选为二分类（仅3和4），原surveylogistic的glogit链接可替换为贝叶斯中的logit链接，逻辑一致。
• 弱信息先验normal(0, var=1000)可有效缓解完全分离问题，为参数提供正则化约束，避免频率学派中出现的无穷大估计值。

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