第 91 章 Regression Methods with continuous outcomes 結果變量爲連續型變量
91.1 用於對連續型結果變量做因果推斷的被估計量
邊際潛在結果的差 marginal potential outcomes:
\[ E\{ Y(1) - Y(0) \} \]
或者是條件下潛在結果的差 conditional potential outcomes:
\[ E\{ Y(1) - Y(2) | \mathbf{V = v} \} \]
邊際潛在結果的差，有專門的名字: the Average Causal Effect (ACE) 平均因果效應 或者叫 Average Treatment Effect 平均治療效應。這裏的 “treatment” 其實不是特指治療，而是泛指所有我們想要比較的暴露。
91.2 鑑定 identification - revision
91.2.1 條件因果均差 conditional causal mean difference
\[ \begin{aligned} E\{ Y(1) - Y(0) | \mathbf{C = c} \} & = E\{ Y(1) | \mathbf{C=c} \} - E\{ Y(0) | \mathbf{C=c} \} \\ \text{(By} & \text{ conditional exchangeability given } \mathbf{C}:) \\ &= E\{ Y(1) | X = 1, \mathbf{C=c} \} - E\{ Y(0) | X = 0, \mathbf{C=c} \} \\ \text{(By} & \text{ consistency:)} \\ & = E\{ Y | X = 1, \mathbf{C=c} \} - E\{ Y | X = 0, \mathbf{C=c} \} \\ \end{aligned} \]
91.2.2 簡單分類型條件變量 \(C\) 的 ACE
\[ \begin{aligned} E\{ Y(1) - Y(0)\} & = \sum_cE\{ Y(1) | C=c \}\text{Pr}(C = c) - \sum_c E\{ Y(0) | C=c \}\text{Pr}(C=c) \\ \text{(By} & \text{ the law of total probability }\uparrow) \\ & = \sum_cE\{ Y(1) | X = 1, \mathbf{C=c} \}\text{Pr}(C=c) \\ & \;\;\;\;\;\;\;\;\;- \sum_cE\{ Y(0) | X = 0, \mathbf{C=c} \} \text{Pr}(C = c) \\ \text{(By} & \text{ conditional exchangeability }\uparrow) \\ & = \sum_cE\{ Y | X = 1, \mathbf{C=c} \}\text{Pr}(C=c) \\ & \;\;\;\;\;\;\;\;\;- \sum_cE\{ Y | X = 0, \mathbf{C=c} \} \text{Pr}(C = c) \\ \text{(By} & \text{ consistency }\uparrow) \\ & = \sum_c\{ E(Y|X = 1, C=c) -E(Y|X=0, C=c) \}\text{Pr}(C=c) \end{aligned} \tag{91.1} \]
91.2.3 簡單連續型條件變量 \(C\) 的ACE
\[ \begin{aligned} E\{ Y(1) - Y(0)\} & = \int_cE\{ Y(1) | C=c \}f_C(c)\text{d}c - \int_c E\{ Y(0) | C=c \}f_C(c)\text{d}c \\ & = \int_cE\{ Y(1) | X = 1, \mathbf{C=c} \}f_C(c)\text{d}c \\ & \;\;\;\;\;\;\;\;\;- \int_cE\{ Y(0) | X = 0, \mathbf{C=c} \} f_C(c)\text{d}c \\ & = \int_cE\{ Y | X = 1, \mathbf{C=c} \}f_C(c)\text{d}c \\ & \;\;\;\;\;\;\;\;\;- \int_cE\{ Y | X = 0, \mathbf{C=c} \} f_C(c)\text{d}c \\ & = \int_c\{ E(Y|X = 1, C=c) -E(Y|X=0, C=c) \}f_C(c)\text{d}c \end{aligned} \tag{91.2} \]
91.3 通過線性回歸模型來估計 ACE
91.3.1 條件因果均值差
假設\(Y,X\)分別表示結果變量和暴露變量，有三個變量需要調整 (做爲條件變量): \(C_1\) 連續型，\(C_2\) 二分類型(0/1)，\(C_3\) 分類型(0/1/2/3)，然後我們擬合的線性回歸模型如下:
\[ \begin{aligned} E(Y|X = x, C_1 = c_1, C_2 = c_2, C_3 = c_3) & = \alpha + \beta_Xx + \gamma_{C_1}c_1 + \gamma_{C_2}c_2 \\ \;\;\; +\gamma_{C_{31}}I(c_3 =1)+&\gamma_{C_{32}}I(c_3 =2)+\gamma_{C_{33}}I(c_3 =3) \end{aligned} \tag{91.3} \]
如果，無相互幹擾 no interference，一致性 consistency，條件可置換性 conditional exchangeability三個最主要的前提條件得到滿足，加上，公式 (91.3) 中三個條件變量得到了正確的模型敘述 (specification of the model is correct)。那麼，這個模型估計的回歸系數 \(\beta_Xx\) 就可以被賦予因果關系的解讀:
\[ E\{ Y(1) -Y(0) |\mathbf{C=c}\} \]
bweight，暴露變量是孕期母親是否吸煙 mbsmoke。這裏先只考慮3個條件變量: 懷孕時的年齡 mage，嬰兒是否是該母親的第一個孩子 fbaby，三個懷孕階段中，該母親第一次訪問婦產科醫生的時間段 prenatal，那麼我們可以擬合的最簡單模型其實是這樣的:
cattaneo2 <- read_dta("backupfiles/cattaneo2.dta")
Cat_mod <- lm(bweight ~ as.factor(mbsmoke) + mage + as.factor(fbaby) + as.factor(prenatal), data = cattaneo2)
summary(Cat_mod)##
## Call:
## lm(formula = bweight ~ as.factor(mbsmoke) + mage + as.factor(fbaby) +
## as.factor(prenatal), data = cattaneo2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3062.28 -308.27 28.87 354.08 2000.92
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2735.8442 78.2260 34.9736 < 2.2e-16 ***
## as.factor(mbsmoke)1 -252.2599 21.5677 -11.6962 < 2.2e-16 ***
## mage 5.2681 1.6146 3.2627 0.0011114 **
## as.factor(fbaby)1 -59.9184 17.7004 -3.3851 0.0007173 ***
## as.factor(prenatal)1 578.8464 68.5077 8.4494 < 2.2e-16 ***
## as.factor(prenatal)2 534.2280 70.6032 7.5666 4.592e-14 ***
## as.factor(prenatal)3 458.5222 80.9992 5.6608 1.597e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 562.03 on 4635 degrees of freedom
## Multiple R-squared: 0.0584, Adjusted R-squared: 0.057181
## F-statistic: 47.912 on 6 and 4635 DF, p-value: < 2.22e-16在無相互幹擾 no interference，一致性 consistency，條件可置換性 conditional exchangeability，和該模型是正確模型的前提下，線性回歸的結果 -252.26 可以被賦予因果推斷的解釋: 在懷孕年齡，嬰兒是否是第一胎，第一次訪問婦產科醫生的孕期時期都相同的條件下，如果比較一個懷孕母親全部都在吸煙，和另一個懷孕母親全部都沒有在吸煙的兩個潛在世界，孕期吸煙的世界的母親生的嬰兒平均出生體重比另一個全部都不吸煙的母親生的嬰兒的出生體重輕 252.3 克。且在我們擬合的模型中，認爲這個新生兒體重的差在其他條件變量取任何值時都保持不變。
模型是正確的這個前提其實是可以放寬的，因爲你可能會擬合這樣一個線性回歸模型:
Cat_mod2 <- lm(bweight ~ as.factor(mbsmoke) + mage + I(mage^2) + as.factor(fbaby)*as.factor(prenatal), data = cattaneo2)
summary(Cat_mod2)##
## Call:
## lm(formula = bweight ~ as.factor(mbsmoke) + mage + I(mage^2) +
## as.factor(fbaby) * as.factor(prenatal), data = cattaneo2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3081.576 -307.339 31.472 350.836 2022.096
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2265.09914 174.11509 13.0092 < 2.2e-16 ***
## as.factor(mbsmoke)1 -251.54047 21.58407 -11.6540 < 2.2e-16 ***
## mage 32.66383 11.70667 2.7902 0.0052893 **
## I(mage^2) -0.50575 0.21294 -2.3751 0.0175863 *
## as.factor(fbaby)1 409.55446 153.91082 2.6610 0.0078181 **
## as.factor(prenatal)1 689.11361 79.41616 8.6772 < 2.2e-16 ***
## as.factor(prenatal)2 684.61449 82.81672 8.2666 < 2.2e-16 ***
## as.factor(prenatal)3 555.06907 95.65524 5.8028 6.956e-09 ***
## as.factor(fbaby)1:as.factor(prenatal)1 -460.27080 154.79291 -2.9735 0.0029598 **
## as.factor(fbaby)1:as.factor(prenatal)2 -538.85385 159.42416 -3.3800 0.0007308 ***
## as.factor(fbaby)1:as.factor(prenatal)3 -393.84997 180.50617 -2.1819 0.0291656 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 561.14 on 4631 degrees of freedom
## Multiple R-squared: 0.062187, Adjusted R-squared: 0.060162
## F-statistic: 30.709 on 10 and 4631 DF, p-value: < 2.22e-16這個模型裏，我們給懷孕時年齡擬合了二次項，又允許 fbaby 和 prenatal 之間有交互作用，但是，這並不妨礙我們對我們最關心的因果關系 mbsmke 的回歸系數的解讀，因爲這兩個模型的結果基本沒有差別。
還有別人可能給出的模型是這樣的:
Cat_mod3 <- lm(bweight ~ as.factor(mbsmoke)*as.factor(fbaby) + mage + I(mage^2) + as.factor(prenatal), data = cattaneo2)
summary(Cat_mod3)##
## Call:
## lm(formula = bweight ~ as.factor(mbsmoke) * as.factor(fbaby) +
## mage + I(mage^2) + as.factor(prenatal), data = cattaneo2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3067.913 -309.329 23.246 349.241 2010.918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2370.45016 167.13049 14.1832 < 2.2e-16 ***
## as.factor(mbsmoke)1 -304.70729 27.49473 -11.0824 < 2.2e-16 ***
## as.factor(fbaby)1 -77.54831 19.35568 -4.0065 6.260e-05 ***
## mage 35.28829 11.62801 3.0348 0.002421 **
## I(mage^2) -0.55192 0.21169 -2.6073 0.009156 **
## as.factor(prenatal)1 559.64429 68.58918 8.1594 4.298e-16 ***
## as.factor(prenatal)2 525.31989 70.55368 7.4457 1.144e-13 ***
## as.factor(prenatal)3 455.26790 80.89765 5.6277 1.934e-08 ***
## as.factor(mbsmoke)1:as.factor(fbaby)1 132.78365 43.73770 3.0359 0.002411 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 561.23 on 4633 degrees of freedom
## Multiple R-squared: 0.061474, Adjusted R-squared: 0.059853
## F-statistic: 37.933 on 8 and 4633 DF, p-value: < 2.22e-16模型 Cat_mod3 中，mbsmoke 和新生兒體重之間的因果關系的解釋發生了變化，因爲我們對 mbsmoke 和 fbaby 之間的交互作用進行了檢驗，是有意義的 p = 0.0024*。這時候，在無相互幹擾 no interference，一致性 consistency，條件可置換性 conditional exchangeability，和該模型是正確模型的前提下，且模型 Cat_mod3 是正確的話，數據中的因果關系解釋及就不止一個了: -304.71 的條件因果關系均差 (conditional causal mean difference) 是針對那些已經有過孩子媽媽來說的; -304.71 + 132.784 = -171.9 這一條件因果關系均差 (conditional causal mean difference) 是對那些第一次懷孕當媽媽的人來說的。吸煙這個本來應該十分有害的行爲，對新生兒體重的影響因果關系似乎在第一次當媽媽的人當中影響較小 (這個因果關系陳述以相同懷孕年齡，和有相同的第一次訪問產科醫生時期爲前提)。
91.3.2 效應修正 effect modification 和 交互作用 interaction
在上文中模型 Cat_mod3 中，如果模型是正確的，且無互相幹擾，一致性，條件可置換性前提都得到滿足時，嬰兒是否是第一胎 fbaby 這一變量，對於我們研究的暴露變量 mbsmoke (孕期吸煙) 和結果變量 bweight (新生兒體重) 之間的關系起到了效應修正作用 (effect modification)。因爲我們看到該模型的結果是孕期吸煙對新生兒體重的影響因爲嬰兒是否是第一胎而發生了很大的變化。流行病學中把這個稱爲交互作用 (interaction)。但是，在因果推斷的研究領域中，傾向於把效應修正和交互作用加以區分。效應修正指對我們關心的關系造成效應修正的變量本身，並沒有因果關系的解釋 (effect modification is not causal with respect to the second variable)，對因果關系造成了效應修正的變量本身，沒有“無互相幹擾，一致性，條件可置換性”前提的要求。它只是衆多的條件變量之一。
相反，因果推斷的研究中，把交互作用的專有名詞保留給兩個暴露變量之間，也就是發生了交互作用的兩個變量，都是要研究的暴露變量，都有和結果變量之間因果關系的討論，所以兩個發生了交互作用的暴露變量，都需要滿足“無互相幹擾，一致性，條件可置換性”前提。
假如不光研究孕期吸煙，研究者還想一起研究孕期飲酒習慣 \((X_2)\)，和吸煙習慣 \((X_1)\) 共同對新生兒體重的因果關系影響:
\[ \{ X_1, X_2 \} \perp\perp Y(x_1, x_2) | \mathbf{C}, \forall x_1,x_2 \]
所以，只有當暴露變量有兩個時 (因爲要同時對飲酒習慣和吸煙習慣兩個暴露變量做潛在結果分析 potential outcome)，才會用到交互作用 (interaction)。
91.3.3 分類型條件變量的平均因果效應 (ACE)
Average Causal Effect (ACE) 平均因果效應:
\[ E\{ Y(1) - Y(0) \} \]
在只有一個分類型條件變量的情況下，我們推導過其 ACE (See equations: (91.1)):
\[ \sum_c\{ E(Y|X=1, C=c) - E(Y|X=0, C=c) \}\text{Pr}(C=c) \]
假設分類條件變量 \(C\) 有四個水平 \(0/1/2/3\)，那麼我們可以針對 \(C\) 的每一層水平擬合線性回歸模型:
\[ \begin{aligned} E(Y|X=x, C=c) & = \alpha + \beta_0 x + \gamma_1 I(c=1) + \gamma_2 I(c=2) + \gamma_3 I(c=3) \\ & \;\;\; + \beta_1 xI(c = 1) + \beta_2 x I(c=2) + \beta_3 x I (c=3) \end{aligned} \tag{91.4} \]
模型 (91.4) 是一個飽和模型，因爲 X 和 C 之間一共只有四種分組組合，我們又擬合了一個含有 8 個參數的模型。也就是說，這個模型允許這 8 種 X 和 C 之間的分組，每組都有不同的結果。
\[ \begin{aligned} \beta_0 & = E(Y|X=1,C=0) - E(Y|X=0, C=0) \\ \beta_0 + \beta_1 & = E(Y|X=1,C=1) - E(Y|X=0, C=1) \\ \beta_0 + \beta_2 & = E(Y|X=1,C=2) - E(Y|X=0, C=2) \\ \beta_0 + \beta_3 & = E(Y|X=1,C=3) - E(Y|X=0, C=3) \\ \end{aligned} \]
爲了簡便起見，給他們分別命名:
\[ \begin{aligned} \beta_0 & = \eta_0 \\ \beta_0 + \beta_1 & = \eta_1 \\ \beta_0 + \beta_2 & = \eta_2 \\ \beta_0 + \beta_3 & = \eta_3 \end{aligned} \]
在只有一個分類型條件變量時，當無相互幹擾，一致性，和條件可置換性的前提被滿足，我們可以把公式 (91.1) 中的 \(E(Y|X=1, C=c) - E(Y|X=0, c=c)\) 全部替換成爲 \(\eta_c\):
\[ \begin{aligned} E\{ Y(1) - Y(0) \} & = \sum_c\{ E(Y|X=1, C=c) - E(Y|X=0, c=c)\}\text{Pr}(C = c) \\ & = \sum_c \eta_c \text{Pr}(C=c) \end{aligned} \]
Cat_mod4 <- lm(bweight ~ as.factor(mbsmoke)*as.factor(prenatal), data = cattaneo2)
summary(Cat_mod4)##
## Call:
## lm(formula = bweight ~ as.factor(mbsmoke) * as.factor(prenatal),
## data = cattaneo2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3095.789 -313.854 23.211 360.458 2064.211
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2864.605 85.970 33.3210 < 2.2e-16 ***
## as.factor(mbsmoke)1 -317.160 138.425 -2.2912 0.0220 *
## as.factor(prenatal)1 571.184 86.560 6.5987 4.612e-11 ***
## as.factor(prenatal)2 478.303 89.535 5.3421 9.628e-08 ***
## as.factor(prenatal)3 428.609 102.354 4.1875 2.873e-05 ***
## as.factor(mbsmoke)1:as.factor(prenatal)1 35.913 140.700 0.2552 0.7985
## as.factor(mbsmoke)1:as.factor(prenatal)2 163.470 146.522 1.1157 0.2646
## as.factor(mbsmoke)1:as.factor(prenatal)3 87.177 168.400 0.5177 0.6047
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 563.74 on 4634 degrees of freedom
## Multiple R-squared: 0.052846, Adjusted R-squared: 0.051415
## F-statistic: 36.936 on 7 and 4634 DF, p-value: < 2.22e-16\[ \begin{aligned} \widehat{\eta_0} & = -317.2 \\ \widehat{\eta_1} & = -317.2 + 35.9 = -281.2 \\ \widehat{\eta_2} & = -317.2 + 163.5 = -153.7\\ \widehat{\eta_3} & = -317.2 + 87.2 = -230.0 \\ \end{aligned} \]
爲了估計平均因果效應，我們還需要 prenatal 的分布概率:
with(cattaneo2, tab1(prenatal, graph = F))## prenatal :
## Frequency Percent Cum. percent
## 0 70 1.5 1.5
## 1 3720 80.1 81.6
## 2 697 15.0 96.7
## 3 155 3.3 100.0
## Total 4642 100.0 100.0所以:
\[ \begin{aligned} \widehat{ACE} & = \sum_c \widehat{\eta_c}\widehat{\text{Pr}}(C=c) \\ & = -317.2 \times \frac{70}{4642} -281.2\times\frac{3720}{4642} -153.7\times\frac{697}{4642}-230.0\times\frac{155}{4642} \\ & = -260.9 \end{aligned} \]
91.3.4 Positivity 非零性
當我們用下面的飽和模型的時候，八個可能的分組中，每個格子裏都不能是零，這一前提條件被成爲非零性 (positivity)。
Cat_mod4 <- lm(bweight ~ as.factor(mbsmoke)*as.factor(prenatal), data = cattaneo2)用概率來表達，就是，在所有可能的 \(c\) 層中的對象，其中暴露變量爲 1 的概率必須在 0, 1 之間:
\[ \textbf{Positivity: } \text{if Pr}(C=c) > 0 \text{ then: } 0<\text{Pr}(X=1|C=c) <1 \]
91.3.5 連續型變量的平均因果效應
Cat_mod5 <- lm(bweight ~ factor(mbsmoke) + mage*factor(mbsmoke) + I(mage^2)*factor(mbsmoke), data = cattaneo2)
summary(Cat_mod5)##
## Call:
## lm(formula = bweight ~ factor(mbsmoke) + mage * factor(mbsmoke) +
## I(mage^2) * factor(mbsmoke), data = cattaneo2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3131.920 -309.155 31.845 351.130 2024.364
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2423.99611 164.87921 14.7016 < 2.2e-16 ***
## factor(mbsmoke)1 1121.03497 414.78336 2.7027 0.006903 **
## mage 64.23439 12.37207 5.1919 2.171e-07 ***
## I(mage^2) -0.97679 0.22658 -4.3110 1.659e-05 ***
## factor(mbsmoke)1:mage -92.69003 32.06902 -2.8903 0.003866 **
## factor(mbsmoke)1:I(mage^2) 1.44359 0.60351 2.3920 0.016796 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 564.93 on 4636 degrees of freedom
## Multiple R-squared: 0.04846, Adjusted R-squared: 0.047434
## F-statistic: 47.22 on 5 and 4636 DF, p-value: < 2.22e-16with(cattaneo2, summ(mage, graph = F))## obs. mean median s.d. min. max.
## 4642 26.505 26 5.619 13 45with(cattaneo2, summ(mage^2, graph = F))## obs. mean median s.d. min. max.
## 4642 734.056 676 305.224 169 2025Y <- cattaneo2$bweight
# X <- with(cattaneo2, cbind(fbaby, mmarried, alcohol, fedu, mage))
X <- with(cattaneo2, cbind(mage, mage^2))
treat <- cattaneo2$mbsmoke
fit1<-ATE(Y,treat,X)
summary(fit1)## Call:
## ATE(Y = Y, Ti = treat, X = X)
##
## Estimate StdErr 95%.Lower 95%.Upper Z.value p.value
## E[Y(1)] 3133.7541 20.7403 3093.1038 3174.4044 151.0948 < 2.22e-16 ***
## E[Y(0)] 3409.4859 9.2843 3391.2890 3427.6827 367.2313 < 2.22e-16 ***
## ATE -275.7317 22.7443 -320.3098 -231.1537 -12.1231 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\[ \begin{aligned} \widehat{\beta_0} & = 1121.035 \\ \widehat{\beta_1} & = -92.690 \\ \widehat{\beta_2} & = 1.444 \\ \Rightarrow \widehat{ACE} & = 1121.035 - 92.690\times26.505 + 1.444\times734.056 \\ & = -275.7 \end{aligned} \]
和 STATA 的 teffects ra 結果做個對比:
. teffects ra (bweight mage mage2) (mbsmoke)
Iteration 0: EE criterion = 9.667e-23
Iteration 1: EE criterion = 7.554e-27
Treatment-effects estimation Number of obs = 4,642
Estimator : regression adjustment
Outcome model : linear
Treatment model: none
----------------------------------------------------------------------------------------
| Robust
bweight | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-----------------------+----------------------------------------------------------------
ATE |
mbsmoke |
(smoker vs nonsmoker) | -275.9901 22.74918 -12.13 0.000 -320.5777 -231.4025
-----------------------+----------------------------------------------------------------
POmean |
mbsmoke |
nonsmoker | 3409.482 9.284654 367.22 0.000 3391.284 3427.679
----------------------------------------------------------------------------------------小數點以後的略差異應該是四舍五入的差異。別的估計包括 Robust Std. Err. 都是十分接近的。
91.4 Practical03 - causal inference
注意: 這裏的練習使用的是STATA 因爲，我在 R 裏找不到像 STATA 的 teffects 這樣靈活且方便的命令，如果你知道，歡迎告訴我: abelardccwang@gmail.com。
數據還是吸煙和新生兒體重的關系的數據:
##
## . use "backupfiles/cattaneo2.dta"
## (Excerpt from Cattaneo (2010) Journal of Econometrics 155: 138-154)
##
## . describe
##
## Contains data from backupfiles/cattaneo2.dta
## obs: 4,642 Excerpt from Cattaneo (2010)
## Journal of Econometrics 155:
## 138-154
## vars: 23 21 Feb 2013 14:43
## size: 143,902
## -------------------------------------------------------------------------------
## storage display value
## variable name type format label variable label
## -------------------------------------------------------------------------------
## bweight int %9.0g infant birth weight (grams)
## mmarried byte %10.0g mmarried 1 if mother married
## mhisp byte %9.0g 1 if mother hispanic
## fhisp byte %9.0g 1 if father hispanic
## foreign byte %9.0g 1 if mother born abroad
## alcohol byte %9.0g 1 if alcohol consumed during
## pregnancy
## deadkids byte %9.0g previous births where newborn
## died
## mage byte %9.0g mother's age
## medu byte %9.0g mother's education attainment
## fage byte %9.0g father's age
## fedu byte %9.0g father's education attainment
## nprenatal byte %9.0g number of prenatal care visits
## monthslb int %9.0g months since last birth
## order byte %9.0g order of birth of the infant
## msmoke byte %27.0g smoke cigarettes smoked during
## pregnancy
## mbsmoke byte %9.0g mbsmoke 1 if mother smoked
## mrace byte %9.0g 1 if mother is white
## frace byte %9.0g 1 if father is white
## prenatal byte %9.0g trimester of first prenatal care
## visit
## birthmonth byte %9.0g month of birth
## lbweight byte %9.0g 1 if low birthweight baby
## fbaby float %9.0g YesNo 1 if first baby
## prenatal1 float %9.0g YesNo 1 if first prenatal visit in 1
## trimester
## -------------------------------------------------------------------------------
## Sorted by:
##
## . tab mbsmoke
##
## 1 if mother |
## smoked | Freq. Percent Cum.
## ------------+-----------------------------------
## nonsmoker | 3,778 81.39 81.39
## smoker | 864 18.61 100.00
## ------------+-----------------------------------
## Total | 4,642 100.00
##
## .
## . summ bweight, detail
##
## infant birth weight (grams)
## -------------------------------------------------------------
## Percentiles Smallest
## 1% 1474 340
## 5% 2438 340
## 10% 2693 397 Obs 4,642
## 25% 3033 454 Sum of Wgt. 4,642
##
## 50% 3390 Mean 3361.68
## Largest Std. Dev. 578.8196
## 75% 3726 5188
## 90% 4026 5216 Variance 335032.2
## 95% 4224 5387 Skewness -.784952
## 99% 4621 5500 Kurtosis 5.788678
##
## .
## . *1. 用簡單線性回顧分析一下 `mbsmoke` 和 `bweight` 之間的關系:
## . *a)
## . regress bweight i.mbsmoke
##
## Source | SS df MS Number of obs = 4,642
## -------------+---------------------------------- F(1, 4640) = 164.62
## Model | 53275939.9 1 53275939.9 Prob > F = 0.0000
## Residual | 1.5016e+09 4,640 323622.478 R-squared = 0.0343
## -------------+---------------------------------- Adj R-squared = 0.0341
## Total | 1.5549e+09 4,641 335032.156 Root MSE = 568.88
##
## ------------------------------------------------------------------------------
## bweight | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## mbsmoke |
## smoker | -275.2519 21.4528 -12.83 0.000 -317.3096 -233.1942
## _cons | 3412.912 9.255254 368.75 0.000 3394.767 3431.056
## ------------------------------------------------------------------------------
##
## .
## .
## . *b)
## .
## . regress bweight i.mbsmoke i.fbaby
##
## Source | SS df MS Number of obs = 4,642
## -------------+---------------------------------- F(2, 4639) = 91.56
## Model | 59045489.8 2 29522744.9 Prob > F = 0.0000
## Residual | 1.4958e+09 4,639 322448.533 R-squared = 0.0380
## -------------+---------------------------------- Adj R-squared = 0.0376
## Total | 1.5549e+09 4,641 335032.156 Root MSE = 567.85
##
## ------------------------------------------------------------------------------
## bweight | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## mbsmoke |
## smoker | -281.0638 21.45789 -13.10 0.000 -323.1314 -238.9961
## |
## fbaby |
## Yes | -71.20491 16.8333 -4.23 0.000 -104.2062 -38.20364
## _cons | 3445.178 11.98063 287.56 0.000 3421.69 3468.666
## ------------------------------------------------------------------------------
##
## .
## . *2. 調整了 `fbaby` 之後，暴露和結果之間的關系發生了怎樣的變化？
## . //說明 `fbaby` 是什麼類型的混雜因子？
## .
## . *看兩個結果的報告，吸煙的線性回歸系數從調整 `fbaby` 前的 `-275.25`，
## . // 絕對值變大爲 `-281.06`，這是一種負方向混雜 (negative confounding)。
## . // 這種混雜可以分析 `fbaby` 和 `mbsmoke` 以及 `bweight`
## . // 各自的關系看出，懷第一胎的母親比較少吸煙，且第一胎嬰兒的出生體重
## . // 均值比不是第一胎嬰兒的出生體重要低:
## .
## . tab fbaby mbsmoke, row
##
## +----------------+
## | Key |
## |----------------|
## | frequency |
## | row percentage |
## +----------------+
##
## 1 if first | 1 if mother smoked
## baby | nonsmoker smoker | Total
## -----------+----------------------+----------
## No | 2,066 543 | 2,609
## | 79.19 20.81 | 100.00
## -----------+----------------------+----------
## Yes | 1,712 321 | 2,033
## | 84.21 15.79 | 100.00
## -----------+----------------------+----------
## Total | 3,778 864 | 4,642
## | 81.39 18.61 | 100..
## . tabstat bweight, by(fbaby)
##
## Summary for variables: bweight
## by categories of: fbaby (1 if first baby)
##
## fbaby | mean
## -------+----------
## No | 3386.681
## Yes | 3329.595
## -------+----------
## Total | 3361.68
## ------------------
##
## .
## . tabstat bweight, by(mbsmoke)
##
## Summary for variables: bweight
## by categories of: mbsmoke (1 if mother smoked)
##
## mbsmoke | mean
## ----------+----------
## nonsmoker | 3412.912
## smoker | 3137.66
## ----------+----------
## Total | 3361.68
## ---------------------
##
## .
## . // 這裏需要重新強調的是，通過比較調整新變量前後的回歸系數的變化，
## . // 能且僅僅只能在線性回歸模型 (可壓縮模型) 時使用，邏輯回歸中不適用。
## .
## .
## . *3 在怎樣的假設條件下，這裏的線性回歸模型的回歸系數 `mbsmoke`
## . // 可以被賦予因果關系？
## .
## . // 1. 無相互幹擾 no interference: 一個懷孕母親吸煙與否，和另一個母親
## . // 生下的嬰兒的出生體重之間沒有關系。
## . // 2. 一致性 consistency: 實際觀察到的孕期吸煙母親的嬰兒出生體重，和
## . // 潛在條件下 (當一個懷孕母親被強制吸煙時) 的嬰兒出生體重 (潛在結果)
## . // 是相同的。同樣地，在另一種潛在條件下 (懷孕母親被禁止吸煙時) 的
## . // 嬰兒出生體重 (潛在結果)，和實際觀察到的不吸煙的母親生下的嬰兒體重
## . // 是相同的。
## . // 3. 條件可置換性 conditional exchangeability: 在 `fbaby` 的各個組別中，
## . // 兩種潛在暴露造成的潛在結果，調整了其它共變量之後，和她們真實的暴露情況
## . // (母親是否吸煙)之間是相互獨立的。在這個模型裏，我們只調整了 `fbaby`
## . // 一個共變量，所以如果要給它的回歸系數加上因果關系結論，還必須假設 (雖然
## . // 很可能不合理) 控制 `fbaby` 這個單一的變量，就完全調整了了母親孕期吸煙和
## . // 新生兒體重之間關系的全部混雜因素。
## . // 4. 模型被正確擬合 correct specification of the regression model: 這是指，
## . // 模型中加入的變量與變量之間的關系，被正確地擬合了，因爲目前只有兩個
## . // 分類型變量在模型中，且該模型沒有加入交互作用項，那麼這條前提假設的含義
## . // 就是，我們認爲 `fbaby` 對孕期吸煙和新生兒體重之間的關系沒有交互作用。
## .
## . *4 在前面解釋過的因果關系的前提條件下，要給 `mbsmoke` 一個因果關系的解釋
## . // 的話，(b) 模型的回歸系數該怎麼解釋呢？用潛在結果的概念解釋。
## .
## . // 在 3. 的前提條件下， `mbsmoke` 的回歸系數的因果關系解讀可以是:
## . // 當條件變量 `fbaby` 嬰兒是否是第一胎的變量保持不變時，281.0638 是暴露
## . // (孕期吸煙) 導致的新生兒體重下降的量，其95%信賴區間是 (238.9961, 323.131
## > 4)。
## . // 這是一個潛在結果的差，所以假如所有的媽媽孕期都吸煙，和所有的媽媽孕期都
## . // 不吸煙相比(潛在暴露)，嬰兒的出生平均體重要輕 281.0638 克:
## . // E{Y(1) | C = c} - E{Y(0) | C = c} = 281.0638
## .
## . *5 用 STATA 的 `teffects ra` 命令擬合相同的模型:
## .
## . teffects ra (bweight fbaby) (mbsmoke)
##
## Iteration 0: EE criterion = 2.764e-25
## Iteration 1: EE criterion = 3.079e-26
##
## Treatment-effects estimation Number of obs = 4,642
## Estimator : regression adjustment
## Outcome model : linear
## Treatment model: none
## ------------------------------------------------------------------------------
## | Robust
## bweight | Coef. Std. Err. z P>|z| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## ATE |
## mbsmoke |
## (smoker |
## vs |
## nonsmoker) | -273.1552 20.93503 -13.05 0.000 -314.1872 -232.1233
## -------------+----------------------------------------------------------------
## POmean |
## mbsmoke |
## nonsmoker | 3414.403 9.279263 367.96 0.000 3396.216 3432.59
## ------------------------------------------------------------------------------
##
## .
## . // 因果均差 (ACE) 的估計在 STATA 被叫做 `ATE`，但是估計的結果略低於
## . // 模型 (b) 的結果: -273.1552 vs. -281.0638。
## .
## . *6 在線性回歸模型中加入 `i.mbsmoke##i.fbaby` 的交互作用項，試着計算
## . // `fbaby` 爲 0/1 時各自的 `mbsmoke` 回歸系數:
## .
## . regress bweight i.mbsmoke##i.fbaby
##
## Source | SS df MS Number of obs = 4,642
## -------------+---------------------------------- F(3, 4638) = 65.17
## Model | 62890495.3 3 20963498.4 Prob > F = 0.0000
## Residual | 1.4920e+09 4,638 321689.034 R-squared = 0.0404
## -------------+---------------------------------- Adj R-squared = 0.0398
## Total | 1.5549e+09 4,641 335032.156 Root MSE = 567.18
##
## ------------------------------------------------------------------------------
## bweight | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## mbsmoke |
## smoker | -339.8145 27.35206 -12.42 0.000 -393.4375 -286.1914
## |
## fbaby |
## Yes | -98.18832 18.53668 -5.30 0.000 -134.529 -61.84761
## |
## mbsmoke#|
## fbaby |
## smoker#Yes | 152.2046 44.02482 3.46 0.001 65.89506 238.5142
## |
## _cons | 3457.406 12.47823 277.08 0.000 3432.942 3481.869
## ------------------------------------------------------------------------------
##
## . est store a
##
## .
## . *to get the stratum specific effects:
## .
## . lincom 1.mbsmoke // when baby is not first born
##
## ( 1) 1.mbsmoke = 0
##
## ------------------------------------------------------------------------------
## bweight | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## (1) | -339.8145 27.35206 -12.42 0.000 -393.4375 -286.1914
## ------------------------------------------------------------------------------
##
## .
## . lincom 1.mbsmoke + 1.mbsmoke#1.fbaby // when baby is first born
##
## ( 1) 1.mbsmoke + 1.mbsmoke#1.fbaby = 0
##
## ------------------------------------------------------------------------------
## bweight | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## (1) | -187.6098 34.49709 -5.44 0.000 -255.2405 -119.9791
## ------------------------------------------------------------------------------
##
## .
## . *7 計算 `fbaby` 各組所佔的百分比: (爲了計算孕期吸煙導致的新生兒體重下降
## . // 的邊際效應 marginal effect)
## .
## . tab fbaby
##
## 1 if first |
## baby | Freq. Percent Cum.
## ------------+-----------------------------------
## No | 2,609 56.20 56.20
## Yes | 2,033 43.80 100.00
## ------------+-----------------------------------
## Total | 4,642 100.00
##
## .
## . *8 用 6, 7 的結果，手動計算一下 ACE 的估計量:
## .
## .
## . *now restore model estimates
## . est restore a
## (results a are active now)
##
## .
## . lincom 0.562*1.mbsmoke + 0.438*(1.mbsmoke + 1.mbsmoke#1.fbaby)
##
## ( 1) 1.mbsmoke + .438*1.mbsmoke#1.fbaby = 0
##
## ------------------------------------------------------------------------------
## bweight | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## (1) | -273.1488 21.55453 -12.67 0.000 -315.406 -230.8917
## ------------------------------------------------------------------------------
##
## .
## . margins, dydx(mbsmoke) // 你也可以用這個 margins 的命令，很方便，但是它估計的
## > 標準誤不使用穩健統計學方法，所以略有不同。
##
## Average marginal effects Number of obs = 4,642
## Model VCE : OLS
##
## Expression : Linear prediction, predict()
## dy/dx w.r.t. : 1.mbsmoke
##
## ------------------------------------------------------------------------------
## | Delta-method
## | dy/dx Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## mbsmoke |
## smoker | -273.1552 21.55433 -12.67 0.000 -315.412 -230.8985
## ------------------------------------------------------------------------------
## Note: dy/dx for factor levels is the discrete change from the base level.
##
## .
## . *9 爲什麼沒有加交互作用項的模型 (b) 給出的回歸系數估計和 `teffects ra`
## . // 的結果相差很大？
## .
## . // 這是因爲如果給模型 (b) 的回歸系數賦予因果關系的解釋的話，第四個前提假設
## . // -- 模型選擇正確且變量在模型中的形式也是正確 -- 太過樂觀了。這個前提假設
## . // 認爲沒有交互作用，但是，如果你看加交互作用項的第三個模型中，交互作用項
## . // 的回顧系數其實是有意義的 (有證據顯示交互作用存在):
## .
## . // mbsmoke#|
## . // fbaby |
## . // smoker#Yes | 152.2046 44.02482 3.46 0.001 65.89506 238.5
## > 142
## .
## . *10 現在給模型中加入更多的共變量，用兩種命令分別擬合，比較其結果:
## .
## . regress bweight mbsmoke fbaby mmarried alcohol fedu mage
##
## Source | SS df MS Number of obs = 4,642
## -------------+---------------------------------- F(6, 4635) = 46.86
## Model | 88933839.4 6 14822306.6 Prob > F = 0.0000
## Residual | 1.4660e+09 4,635 316278.403 R-squared = 0.0572
## -------------+---------------------------------- Adj R-squared = 0.0560
## Total | 1.5549e+09 4,641 335032.156 Root MSE = 562.39
##
## ------------------------------------------------------------------------------
## bweight | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## mbsmoke | -221.2625 22.28894 -9.93 0.000 -264.9594 -177.5656
## fbaby | -49.45095 17.62184 -2.81 0.005 -83.99814 -14.90375
## mmarried | 152.3058 21.86593 6.97 0.000 109.4382 195.1734
## alcohol | -64.26606 47.48732 -1.35 0.176 -157.3638 28.83168
## fedu | 4.515101 2.573722 1.75 0.079 -.5306195 9.560821
## mage | 1.282969 1.750586 0.73 0.464 -2.149013 4.714952
## _cons | 3230.456 49.40366 65.39 0.000 3133.601 3327.311
## ------------------------------------------------------------------------------
##
## .
## .
## . teffects ra (bweight fbaby mmarried alcohol fedu mage) (mbsmoke)
##
## Iteration 0: EE criterion = 4.396e-24
## Iteration 1: EE criterion = 7.758e-26
##
## Treatment-effects estimation Number of obs = 4,642
## Estimator : regression adjustment
## Outcome model : linear
## Treatment model: none
## ------------------------------------------------------------------------------
## | Robust
## bweight | Coef. Std. Err. z P>|z| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## ATE |
## mbsmoke |
## (smoker |
## vs |
## nonsmoker) | -231.1408 24.50303 -9.43 0.000 -279.1659 -183.1158
## -------------+----------------------------------------------------------------
## POmean |
## mbsmoke |
## nonsmoker | 3403.054 9.532216 357.01 0.000 3384.372 3421.737
## ------------------------------------------------------------------------------
##
## .
## .
## . // 此時我們發現，簡單現行回歸估計的因果均值差(ATE)總是和考慮了更復雜關系的
## > 模型相比
## . // 相差較多。
## .
## . *11 你可以用下面的非 teffects 代碼還原上面的計算:
## .
## . qui regress bweight mbsmoke fbaby mmarried alcohol fedu mage if mbsmoke==0
##
## . predict Y0
## (option xb assumed; fitted values)
##
## .
## . qui regress bweight mbsmoke fbaby mmarried alcohol fedu mage if mbsmoke==1
##
## . predict Y1
## (option xb assumed; fitted values)
##
## .
## . sum Y0
##
## Variable | Obs Mean Std. Dev. Min Max
## -------------+---------------------------------------------------------
## Y0 | 4,642 3403.054 101.6744 3163.591 3549.308
##
## . gen E0=r(mean)
##
## .
## . sum Y1
##
## Variable | Obs Mean Std. Dev. Min Max
## -------------+---------------------------------------------------------
## Y1 | 4,642 3171.914 67.17853 2839.728 3306.694
##
## . gen E1=r(mean)
##
## .
## . gen ACE = E1-E0
##
## . sum ACE
##
## Variable | Obs Mean Std. Dev. Min Max
## -------------+---------------------------------------------------------
## ACE | 4,642 -231.1409 0 -231.1409 -231.1409
##
## .
## . // 或者使用方便的 margins
## .
## . qui regress bweight i.mbsmoke fbaby mmarried alcohol fedu mage ///
## > i.mbsmoke#i.fbaby i.mbsmoke#i.mmarried i.mbsmoke#i.alcohol i.
## > mbsmoke#c.fedu ///
## > i.mbsmoke#c.mage
##
## .
## . margins, dydx(mbsmoke)
##
## Average marginal effects Number of obs = 4,642
## Model VCE : OLS
##
## Expression : Linear prediction, predict()
## dy/dx w.r.t. : 1.mbsmoke
##
## ------------------------------------------------------------------------------
## | Delta-method
## | dy/dx Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## mbsmoke |
## smoker | -231.1408 24.05129 -9.61 0.000 -278.2928 -183.9888
## ------------------------------------------------------------------------------
## Note: dy/dx for factor levels is the discrete change from the base level.
##
## .
## . // 值得注意的是，即使是使用 `teffects ra` 我們可能對模型形式的指定還是過於簡
## . // 單，例如上面的模型中加入了許多變量，但是 STATA 其實並沒有考慮有三個或者三
## . // 個以上變量之間發生交互作用的情況，而且， `fedu, mage` 被認爲和結果變量
## . // `bweight` 呈簡單一次的線性關系。References
Cattaneo, Matias D. 2010. “Efficient Semiparametric Estimation of Multi-Valued Treatment Effects Under Ignorability.” Journal of Econometrics 155 (2). Elsevier: 138–54.