拓端tecdat|Stata辅导估算观测数据的风险比

原文链接：http://tecdat.cn/?p=6419

 

在分析二元结果时，逻辑回归是分析师对回归建模的默认方法。随机研究中，当然很容易估计比较两个治疗组的风险比。对于观察数据，治疗不是随机分配的，估计治疗效果的风险比有点棘手。

 

理想情况 - 随机治疗分配
理想情况下，我们首先模拟（在Stata中）一个大型数据集，该数据集可能在随机试验中出现：

gen x = rnormal（）
gen z =（runiform（）<0.5）
gen xb = x + z
gen pr = exp（xb）/（1 + exp（xb））
gen y =（runiform（）<pr）

此代码为10,000个人生成数据集。每个都有一个基线变量x的值，它是从标准N（0,1）分布模拟的。接下来，根据随机研究，我们模拟一个二进制变量z，概率0.5为1，概率0.5为0.然后生成二元结果y，我们从逻辑回归模型生成它，对数几率为1等于x + z。因此，对于x，调整z = 1与z = 0的真实优势比是exp（1）= 2.72。

由于处理是随机分配的，我们可以忽略x并使用带有日志链接的GLM命令估计比较z = 1到z = 0的风险比：

Generalized linear models                          No. of obs      =     10000
Optimization     : ML                              Residual df     =      9998
                                                   Scale parameter =         1
Deviance         =  12960.39176                    (1/df) Deviance =  1.296298
Pearson          =        10000                    (1/df) Pearson  =    1.0002

Variance function: V(u) = u*(1-u)                  [Bernoulli]
Link function    : g(u) = ln(u)                    [Log]

                                                   AIC             =  1.296439
Log likelihood   = -6480.195879                    BIC             = -79124.59

------------------------------------------------------------------------------
             |                 OIM
           y | Risk Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           z |   1.429341   .0241683    21.13   0.000     1.382748    1.477503
       _cons |    .495886   .0070829   -49.11   0.000     .4821963    .5099643
------------------------------------------------------------------------------

风险比估计为1.43，因为数据集很大，95％置信区间非常窄。

 

估算观测数据的风险比

现在让我们考虑观测数据的情况。为此，我们模拟了一个新的数据集 ：

gen x=rnormal()
gen tr_xb=x
gen tr_pr=exp(tr_xb)/(1+exp(tr_xb))
gen z=(runiform() < tr_pr)
gen xb=x+z
gen pr=exp(xb)/(1+exp(xb))
gen y=(runiform() < pr)

如果我们为y运行相同的GLM模型，忽略x，我们得到：

. glm y z, family(binomial) link(log) eform

Generalized linear models                          No. of obs      =     10000
Optimization     : ML                              Residual df     =      9998
                                                   Scale parameter =         1
Deviance         =  12014.93919                    (1/df) Deviance =  1.201734
Pearson          =        10000                    (1/df) Pearson  =    1.0002

Variance function: V(u) = u*(1-u)                  [Bernoulli]
Link function    : g(u) = ln(u)                    [Log]

                                                   AIC             =  1.201894
Log likelihood   = -6007.469594                    BIC             = -80070.04

------------------------------------------------------------------------------
             |                 OIM
           y | Risk Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           z |   1.904111   .0353876    34.65   0.000        1.836    1.974748
       _cons |   .4101531   .0069811   -52.36   0.000      .396696    .4240667
------------------------------------------------------------------------------

 

使用对数广义线性模型

最明显的方法是在我们的GLM命令中添加x：

. glm y z x, family(binomial) link(log) eform

Iteration 0:   log likelihood = -9271.4631  
Iteration 1:   log likelihood = -5833.7585  
Iteration 2:   log likelihood = -5733.9167  (not concave)
Iteration 3:   log likelihood = -5733.9167  (not concave)
...

然而，这无法收敛  。

 

通过逻辑模型估计风险比率

一个相对简单的替代方案是使用逻辑模型来估计调整x的治疗风险比。

 

Logistic regression                               Number of obs   =      10000
                                                  LR chi2(2)      =    2797.60
                                                  Prob > chi2     =     0.0000
Log likelihood = -5343.6848                       Pseudo R2       =     0.2075

------------------------------------------------------------------------------
           y | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           z |   2.947912   .1442639    22.09   0.000     2.678297    3.244668
           x |   2.693029   .0811728    32.87   0.000     2.538542    2.856918
       _cons |   .9910342   .0322706    -0.28   0.782      .929761    1.056345
------------------------------------------------------------------------------

然而，我们可以使用该模型来计算每个人的预测概率，首先假设所有个体都未被治疗（z = 0），然后假设所有个体都被治疗（z = 1）：

replace z=0
predict pr_z0, pr
replace z=1
predict pr_z1, pr

此代码首先生成一个新变量zcopy，它保留原始治疗分配变量的副本。然后我们将所有个体z设置为0，并询问y = 1的预测概率。然后我们将所有个体设置为z = 1，并再次计算P（y = 1）。现在估计z = 1与z = 0相比的风险比，我们简单地考虑这两个条件下的边际风险比，即P（y = 1 | z = 1）/ P（y = 1） | Z = 0）：

summ pr_z0 pr_z1

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       pr_z0 |     10000    .4970649    .2060767   .0166056   .9765812
       pr_z1 |     10000    .7094445    .1765126   .0474181   .9919309

di .7094445/.4970649
1.4272673
replace z=zcopy

 给出了我们估计的风险比，比较z = 1到z = 0，为1.43，与我们第一次模拟数据时估计的风险比相同，其中治疗分配是完全随机的（特别是独立于x）。

 

 置信区间

 

我们已经找到了风险比的点估计，但我们当然也喜欢置信区间，以指示估计的精确度。

首先，当z设置为0然后设置为1时，我们使用teffects给出我们对二元结果的边际均值的估计（相当于y = 1的概率）：

 
Iteration 0:   EE criterion =  1.898e-21  
Iteration 1:   EE criterion =  5.519e-34  

Treatment-effects estimation                    Number of obs      =     10000
Estimator      : regression adjustment
Outcome model  : logit
Treatment model: none
------------------------------------------------------------------------------
             |               Robust
           y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
POmeans      |
           z |
          0  |   .4957607   .0072813    68.09   0.000     .4814897    .5100318
          1  |   .7078432   .0071566    98.91   0.000     .6938164    .7218699
------------------------------------------------------------------------------

 

为了计算风险比和置信区间，我们首先使用teffects ra，coeflegend来查找Stata保存估算值的名称：

Treatment-effects estimation                    Number of obs      =     10000
Estimator      : regression adjustment
Outcome model  : logit
Treatment model: none
------------------------------------------------------------------------------
           y |      Coef.  Legend
-------------+----------------------------------------------------------------
POmeans      |
           z |
          0  |   .4957607  _b[POmeans:0bn.z]
          1  |   .7078432  _b[POmeans:1.z]
------------------------------------------------------------------------------

我们现在可以使用nlcom计算风险比及其置信区间。但是，由于这将为我们提供基于Wald的对称置信区间，因此最好找到对数风险比的这个区间，然后将得到的区间反向转换为风险比例：

 

       _nl_1:  log(_b[POmeans:1.z] / _b[POmeans:0bn.z])

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _nl_1 |   .3561291   .0172327    20.67   0.000     .3223536    .3899047
------------------------------------------------------------------------------

. di exp(.3561291)
1.4277919

. di exp(.3223536)
1.3803728

. di exp(.3899047)
1.47684

因此，我们估计风险比为1.43,95％置信区间为1.38至1.48。

 

如果您有任何疑问，请在下面发表评论。