第 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}\} \]

Example 91.1 孕期吸煙和嬰兒出生體重的關系: 數據來自(Cattaneo 2010) 結果變量是出生體重 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

這個模型裏,我們給懷孕時年齡擬合了二次項,又允許 fbabyprenatal 之間有交互作用,但是,這並不妨礙我們對我們最關心的因果關系 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 和新生兒體重之間的因果關系的解釋發生了變化,因爲我們對 mbsmokefbaby 之間的交互作用進行了檢驗,是有意義的 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-16
with(cattaneo2, summ(mage, graph = F))
##  obs. mean   median  s.d.   min.   max.  
##  4642 26.505 26      5.619  13     45
with(cattaneo2, summ(mage^2, graph = F))
##  obs. mean    median  s.d.    min.   max.  
##  4642 734.056 676     305.224 169    2025
Y <- 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.