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

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

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

gen x = rnormal（）gen z =（runiform（）<0.5）gen xb = x + zgen 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      =     10000Optimization     : ML                              Residual df     =      9998                                                   Scale parameter =         1Deviance         =  12960.39176                    (1/df) Deviance =  1.296298Pearson          =        10000                    (1/df) Pearson  =    1.0002Variance function: V(u) = u*(1-u)                  [Bernoulli]Link function    : g(u) = ln(u)                    [Log]                                                   AIC             =  1.296439Log 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=xgen tr_pr=exp(tr_xb)/(1+exp(tr_xb))gen z=(runiform() < tr_pr)gen xb=x+zgen pr=exp(xb)/(1+exp(xb))gen y=(runiform() < pr)

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

. glm y z, family(binomial) link(log) eformGeneralized linear models                          No. of obs      =     10000Optimization     : ML                              Residual df     =      9998                                                   Scale parameter =         1Deviance         =  12014.93919                    (1/df) Deviance =  1.201734Pearson          =        10000                    (1/df) Pearson  =    1.0002Variance function: V(u) = u*(1-u)                  [Bernoulli]Link function    : g(u) = ln(u)                    [Log]                                                   AIC             =  1.201894Log 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) eformIteration 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.0000Log 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=0predict pr_z0, prreplace z=1predict 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   .9919309di .7094445/.49706491.4272673replace 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      =     10000Estimator      : regression adjustmentOutcome model  : logitTreatment 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      =     10000Estimator      : regression adjustmentOutcome model  : logitTreatment 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。

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。