原文連結：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。