社會科學統計方法
隨機實驗與因果關係
隨機實驗與因果關係
蔡佳泓
3/1-3/8/2024
1 統計的目的
- 測量：研究者收集資料，經過加總、求平均值、計算離散程度等等，例如：
- 統計一個國家民眾對於政治的興趣、
- 一個地區的空氣品質
- 擲無數次銅板得到正面的比例。有一篇文章探討擲硬幣的實驗結果。
- 預測：
- 收集資料，建立一個模型，然後預測下一個事件發生的機率，如果發生機率接近實際發生的事件，則可以繼續使用該模型。
- 最近流行的人工智慧(AI)，屬於預測的一種，不過生成式AI是產生新的資料而不是預測，MIT有一篇文章Explained: Generative AI討論這些概念。
- 因果關係：我們可以透過隨機實驗，排除所有可能的干擾因素，估計出兩個變數的因果關係。例如對一群人隨機抽樣，有人抽到籤去服役，有人不需要去服役，然後在若干年後詢問這些人對於國家的認同感，這些人之間唯一的差別是有沒有去服役，就可以估計出服役與否對國家認同感的作用。
- 過去我們比較重視預測，但是近年來越來越重視因果關係。
2 因果關係的實例
Figure 2.1: 學前教育
Figure 2.2: 餐廳員工數與最低時薪
- 大家可以想一個可行的實驗設計嗎？
3 實驗設計的原理
3.1 Potential Outcome Model
\(D_{i}\): Indicator of treatment for unit \(i\) \[ \mathrm{D}_{i} = \begin{cases} 1, & \text{if unit} \hspace{3pt} i \hspace{3pt} \text{receives treatment} \\ % & is your "\tab"-like command (it's a tab alignment character) 0, & \text{otherwise} \end{cases} \]
- \(Y_{1i}\)代表我們給予刺激之後的表現，\(Y_{0i}\)代表我們沒有給予刺激之後的表現，例如\(Y_{1i}\)代表我們給盆栽澆水、\(Y_{0i}\)代表沒有澆水的成長情形。
\(Y_{i}\): Observed outcome of interest for unit \(i\)
\[ \tag{1} \mathrm{Y}_{di}=\begin {cases} Y_{0i}, & \text{Potential outcome for unit} \hspace{3pt}i\hspace{3pt} \text{without treatment} \\ Y_{1i}, & \text{Potential outcome for unit } \hspace{3pt} i \hspace{3pt} \text{with treatment} \end{cases} \]
- 為了明確表示刺激有無起見，\(Y_{1i}\)的平均數可以寫成\(E[Y_{1}|D=1]\)，\(Y_{0i}\)的平均數可以寫成\(E[Y_{0}|D=0]\)。
- 要注意到我們無法觀察到\(E[Y_{1}|D=0]\)，也就是我們沒有給予刺激，但是這些個體卻從別的地方接受到刺激。例如沒被澆水的盆栽從其他地方得到水分。同樣的，我們也觀察不到\(E[Y_{0}|D=1]\)，也就是明明澆水，但是盆栽好像被蓋子蓋起來，完全沒接收到水分。換句話說，我們只能觀察到有給予刺激的結果以及沒給予刺激的結果，但是觀察不到有給予刺激卻好像沒給予的結果，也觀察不到沒給予刺激卻好像有得到刺激的結果。
\(Y_{i}=D_{i}\cdot Y_{1i}+(1-D_{i})\cdot Y_{0i}\)
\[\begin{equation} \tag{2} \mathrm{Y}_{i}=\begin {cases} Y_{0i}, & \text{if} \hspace{3pt}D_{i}=0 \\ Y_{1i}, & \text{if} \hspace{3pt} D_{i}=1 \end{cases} \end{equation}\]
\(\alpha_{i} = Y_{1i}-Y_{0i}\)
\[ Y_{i}(D_{1}, D_{2},\ldots,D_{n})=Y_{i}(D'_{1}, D'_{2},\ldots,D'_{n})\quad \text{if}\quad D_{i}=D'_{i} \]
\[\begin{align} \tag{3} ATE & = E[Y_{1}-Y_{0}]\\ & =E[Y_{1}]-E[Y_{0}]\\ & =E[Y_{1}\mid D=1]-E[Y_{0}\mid D=0] \\ & =E[Y_{1}-Y_{0}]\\ & =E[Y_{1}]-E[Y_{0}] \end{align}\]
| \(i\) | \(Y_{1i}\) | \(Y_{0i}\) | \(D_{i}\) | \(\alpha_{i}\) | \(Y_{i}\) |
|---|---|---|---|---|---|
| 1 | 3 | 0 | 1 | 3 | 3 |
| 2 | 1 | 1 | 1 | 0 | 1 |
| 3 | 1 | 0 | 1 | 0 | |
| 4 | 1 | 1 | 0 | 0 | |
| \(E[Y_{1}]\) | 1.5 | ||||
| \(E[Y_{0}]\) | 0.5 | ||||
| \(E[Y_{1}-Y_{0}]\) |
\[ E[Y_{1}] = \frac{1}{N}\Sigma Y_{1i}=1.5 \] \[ E[Y_{0}] = \frac{1}{N}\Sigma Y_{0i}=0.5 \] \[ E[Y_{1}]-E[Y_{0}] = 1 \] \[ \alpha_{ATE}=E[Y_{1}-Y_{0}]=\frac{1}{4}\cdot(3+0+1+0)=1 \]
| \(i\) | \(Y_{1i}\) | \(Y_{0i}\) | \(D_{i}\) | \(\alpha_{i}\) | \(Y_{i}\) |
|---|---|---|---|---|---|
| 1 | 3 | 0 | 3 | 3 | |
| 2 | 1 | 1 | 0 | 1 | |
| \(E[Y_{1}]\) | 2 | ||||
| \(E[Y_{0}]\) | 0.5 | ||||
| \(E[Y_{1}-Y_{0}]\) |
\[ \alpha_{ATT}=E[Y_{1}-Y_{0}|D=1]=\frac{1}{2}\cdot(3+1-0-1)=1.5 \]
3.2 ATE與ATT
- 需要注意的是，ATE可以改寫如下：
\[\begin{align*} \text{ATE} & =E[Y_{1}|D=1]-E[Y_{0}|D=0]\\ & = E[Y_{1}|D=1]-E[Y_{0}|D=1] +E[Y_{0}|D=1]-E[Y_{0}|D=0] \end{align*}\]
\(E[Y_{0}|D=1] - E[Y_{0}|D=1]\)其實是0，但是我們刻意加進去。
\(E[Y_{1}|D=1]-E[Y_{0}|D=1]\)被稱為ATT或者ATET， expected treatment effect given the treatment。而 \(E[Y_{0}|D=1]-E[Y_{1}|D=0]\)代表某種誤差。也就是說：
\[ \text{ATE}=\text{ATT}+Bias \]
如果\(E[Y_{0}|D=1]-E[Y_{0}|D=0]= 0\)，那麼ATE = ATT，代表我們不見得要把受試者分成實驗組跟控制組，只要觀察實驗組就可以。
\(E[Y_{0}|D=1]-E[Y_{0}|D=0]= 0\)表示實際上收到刺激但是「假設」並沒有收到的表現狀況，跟實際上沒收到刺激而且真的沒有收到刺激的表現狀況相等。例如：給植物澆水，但是沒注意到水都澆到盆栽外面去，跟沒澆水的表現狀況相同。或者給頭痛患者服藥，但是因為頭太痛而忘了吃藥，藥效跟沒吃一樣。如果這些情形為真，表示就算我們無法觀察到控制組的受試者「真的」得到刺激的表現，我們可以假設控制組的受試者即使「真的」得到刺激，他們的表現跟沒有得到刺激的受試者相同。
\[ \alpha_{ATT}=E[Y_{1i}]=\frac{1}{2}\cdot(3+0)=1.5 \]
- 理論上\(E[Y_{0}|D=1]\)是觀察不到的，因為它代表在得到刺激的情況下，「假如」他們沒有得到的反事實(counterfactual)狀況。所以我們也觀察不到ATT=\(E[Y_{1}-Y_{0}|D=1]\)。
假設有兩群人，一群人吃頭痛藥以減緩頭痛，另一群人不吃頭痛藥。ATT代表吃頭痛藥的這一群人之後的平均頭痛狀況，ATE則代表吃頭痛藥與不吃頭痛藥之後的平均頭痛狀況差異。如果頭痛藥真的對每個人都有效，ATE會大於ATT，這是因為沒吃藥的人的平均頭痛狀況不會減輕，吃藥的人才會減輕，ATE代表的差異就會很明顯。
- 當\(E[Y_{0}|D=1]>E[Y_{0}|D=0]\)，Bias>0，ATE>ATT。
- 如果控制組「真的」得到刺激，例如控制組的植物真的被澆水了，理論上成長狀況應該比沒澆到水來得好。換句話說，有沒有澆水對植物的影響應該是很顯著的，這時候我們就必須分組觀察有無澆水，而不能只觀察澆過水的植物。
- 對於沒來上課的同學，「假如」要求他們來上這堂課，不一定會比沒上課的同學對於實驗設計了解更多，因為來上這門課可能因為昨天晚上沒睡好、老師講得不清楚等等原因，反而比在家自習的學習效果差，所以Bias < 0，ATE < ATT，代表我們如果只觀察來上課的同學對於實驗設計的了解，有可能會比較好。
- 如果隨機找一群人來參加職業訓練、另一群人請他們自己學習，然後觀察他們的收入。沒來參加訓練的人，因為在同一個時間自己學習或者努力賺錢，反而可能比剛剛參加完職業訓練的人的收入來得高，\(E[Y_{0}|D=1]<E[Y_{0}|D=0]\)。
4 實驗設計與迴歸模型
- 延續之前的符號，只有一個實驗變數D的迴歸模型可以表示如下：
\[ Y=\beta_{0} +\beta_{1}D \]
\[\begin{align} D= \begin{cases}1\\ 0 \end{cases} \end{align}\]
- 如果刺激\(D\)是０與1，那麼迴歸係數\(\beta_{1}\)等於平均值差異的實驗效果。
\[ \hat{\beta}_{OLS}=\frac{\sum_{i=1}^n(Y-\bar{Y})(D-\bar{D})}{\sum_{i=1}^n(D-\bar{D})^2}=\hat{\tau} \]
- 因為樣本迴歸係數是母體的最佳估計，所以我們可以用最小平方法迴歸模型估計母體的實驗效果。
- 但是要應用sandwich這個套件中的\(\texttt{vcovHC}\)函數可以得到相近的變異數估計，而不是直接使用最小平方法估計的標準誤。
- 加上自變數X在模型中，X代表可能的干擾因素，例如D是是否實施新的數學教學方法，X是每個同學之前的數學成績，Y是實施新教學方法後的成績。D的係數\(\beta_{1}\)可以解讀為當X相同時，Y的平均變動程度。
\[ Y=\beta_{0} +\beta_{1}D+\beta_{2}X \]
5 實際操作：小班制與成績
star <- read.csv("~/Dropbox/EastAsia2024/data/DSS/STAR.csv")
head(star)## classtype reading math graduated
## 1 small 578 610 1
## 2 regular 612 612 1
## 3 regular 583 606 1
## 4 small 661 648 1
## 5 small 614 636 1
## 6 regular 610 603 0star <- star %>% mutate(D = as.factor(classtype)) %>%
mutate(graduated = recode_factor(graduated,
'1'='Yes', '0'='No'))
class(star$D)## [1] "factor"levels(star$D)## [1] "regular" "small"mean(star$reading[star$D=='small'])## [1] 632.7mean(star$reading[star$D=='regular'])## [1] 625.5- 畫圖表示實驗組與控制組的平均值：
library(dplyr)
star.sta <- star|> group_by (D) |>
summarise(avg.reading = mean(reading)) |>
mutate(avg.reading =round(avg.reading, 3))
p1 <- ggplot2::ggplot(data=star.sta, aes(x=D, y = avg.reading, fill=D)) +
geom_bar(stat = 'identity') +
geom_text(aes(label = avg.reading) )
p1
library(gtsummary)
star %>% tbl_summary(
by = D,
statistic = list(
all_continuous() ~ "{mean} ({sd})",
graduated ~ "{n} / {N} ({p}%)"
),
digits = all_continuous() ~ 2) %>%
add_overall() %>%
add_n() %>%
modify_header(label ~ "**Variable**") | Variable | N | Overall, N = 1,2741 | regular, N = 6891 | small, N = 5851 |
|---|---|---|---|---|
| classtype | 1,274 | |||
| regular | 689 (54%) | 689 (100%) | 0 (0%) | |
| small | 585 (46%) | 0 (0%) | 585 (100%) | |
| reading | 1,274 | 628.80 (36.73) | 625.49 (35.88) | 632.70 (37.37) |
| math | 1,274 | 631.59 (38.84) | 628.84 (37.94) | 634.83 (39.66) |
| graduated | 1,274 | 1,108 / 1,274 (87%) | 597 / 689 (87%) | 511 / 585 (87%) |
| 1 n (%); Mean (SD); n / N (%) | ||||
mean(star$reading[star$D=='small']) - mean(star$reading[star$D=='regular'])## [1] 7.211- 用廻歸模型估計實驗效果：
m1 <- lm(reading ~ D, data = star)
# Summarize the model
stargazer::stargazer(m1, title = 'OLS and Causal Inference', type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))| Dependent variable: | |
| reading | |
| Dsmall | 7.211*** |
| (2.056) | |
| Constant | 625.500*** |
| (1.393) | |
| Observations | 1,274 |
| R2 | 0.010 |
| Adjusted R2 | 0.009 |
| Residual Std. Error | 36.570 (df = 1272) |
| F Statistic | 12.300*** (df = 1; 1272) |
| Note: | p<0.1; p<0.05; p<0.01 |
6 次團體(subgroup): CATE
- 我們可以估計不同條件或者不同團體的實驗效果，也就是：
\[ \alpha_{ATE(\text{x})}=E[Y_{1}-Y_{0}|X=\text{x}] \]
- 如果是ATT，則寫成：
\[ \alpha_{ATT(\text{x})}=E[Y_{1}-Y_{0}|D=1, X=\text{x}] \]
- 我們可以觀察實驗組中的某一個次團體跟控制組相對應的次團體之間的平均效果差異，這個差異稱為conditional average treatment effect (CATE)
- 例如實驗組服用頭痛藥，控制組多喝水，比較兩組年紀輕的受試者，可能有比較明顯的藥效差異，年紀大的受訪者可能沒有。
- 例如開出大樂透之後的彩券行因為會大大宣傳，所以比沒開過的彩券行應該有比較好的收入，但是在都會區的彩券行可能比較多，附近民眾可能很快就忘記哪家彩券行有開過大樂透。
6.1 估計CATE
- 我們用迴歸模型的交互作用項估計CATE。交互作用項的最小平方法迴歸模型可寫成：
\[ Y = \beta_{0}+\beta_{1}\cdot D+\beta_{2}\cdot X+\beta_{3}\cdot DX \]
最小平方法迴歸模型如果有交互作用，可以估計CATE。這是因為自變數X與刺激變數D之間的交互作用，代表自變數X透過D的作用。當D=1時，代表得到刺激，D=0時，代表屬於控制組。如果D=1且X=x時，D的係數大於當D=0的係數，表示D的作用會隨著X=x而不同。
例如以下資料有X變數與刺激D，X的值是1, 2, 3，D則是0, 1。
# Simulate some data
set.seed(02138)
n <- 250
X <- sample(c(1:3), n, replace = T) # Covariate
D <- rbinom(n, 1, 0.55) # Treatment indicator
Y <- 1 + D * 2 + X * 3 + D * X * 2 + rnorm(n) # Outcome
DF = data.frame(Y=Y, Treatment=D, X=X)- 當受訪者在X=1這一組，實驗效果可計算為：
\[\begin{align} E[Y_{1}-Y_{0}|X=1] =E[Y_{1}|X=1]-E[Y_{0}|X=1]\\ & = (\beta_{0}+\beta_{1}+\beta_{2}+\beta_{3})-(\beta_{0}+\beta_{2})\\ & = \beta_{1}+\beta_{3} \end{align}\]
- 可以看出，\(\beta_{1}\)與\(\beta_{3}\)都與D相關，所以當X變動時，D與Y的關係也會變動。
- 其他組依此類推。用
R估計有交互作用的迴歸模型如下：
# Fit a linear model with interaction between treatment and covariate
res <- lm(Y ~ D * X, data = DF)
# Summarize the model
stargazer::stargazer(res, title = 'Interaction Model', type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))| Dependent variable: | |
| Y | |
| D | 1.967*** |
| (0.344) | |
| X | 2.900*** |
| (0.118) | |
| D:X | 2.028*** |
| (0.156) | |
| Constant | 1.163*** |
| (0.260) | |
| Observations | 250 |
| R2 | 0.955 |
| Adjusted R2 | 0.954 |
| Residual Std. Error | 1.001 (df = 246) |
| F Statistic | 1,721.000*** (df = 3; 246) |
| Note: | p<0.1; p<0.05; p<0.01 |
- 用\(\texttt{interplot}\)函數畫圖了解X變動時，D對於Y的影響，圖6.1顯示，當X從1上升到2，D的係數也從4上升到6，以此類推：
library(ggplot2)
interplot::interplot(res, var1 = "D", var2 = "X", val2 = 1:3,
, ercolor = "#EE00CC", esize = 1.5)+
labs(x = "Covariate (X)") +
geom_point(size = 2, color = "#2211CC") +
theme_bw()
Figure 6.1: Interacton Plot
7 Local Average Treatment Effect (LATE)
在實驗組中每一個個案應該都接收到刺激，有些人可能不確實接受刺激，但是研究者無法強迫。同樣的，在控制組中，每一個個案應該不會接收到刺激，但是有些人可能從其他地方接受到刺激。
願意接受刺激的個案可能跟不願意的個案之間有一些差異，但是我們還沒有觀察到。
換句話說，刺激的效果是因為某些外在的因素而有差異，也就D會被外在的因素影響，造成接收到刺激的差異。我們必須先排除這些外在因素的影響，才能正確估計刺激的效果。因為分成兩個步驟，所以叫做兩階段(two-stage)迴歸模型。
第一階段： \[ D=\pi_{0}+\pi_{1}Z+\nu_{1} \]
第二階段：
\[ Y=\alpha_{0}+\alpha_{1}D+\nu_{2} \]
- 如果兩個假設都成立：\(Cov[\nu_{1}, Z]=0\), 而且\(\pi_{1}\neq 0\), \(Cov[\nu_{3}, Z]=0\)
\[ Y=\gamma_{0}+\gamma_{1} Z +\nu_{3} \]
- 可以導出：
\[ \alpha_{1}=\frac{\gamma_{1}}{\pi_{1}}=\frac{\texttt{Cov}[Y,Z]}{\texttt{Cov}[D,Z]} \]
實驗組中接收到刺激以及控制組中沒接收到刺激的個案，稱為complier。這兩群人之間的平均效果差異，稱為local average treatment effect (LATE)，也稱為complier average causal effect(CACE)。
例如有學生錄取名校的留學顧問中心應該比沒有的中心有比較好的招生收入，但是有些學生要求顧問中心不要宣傳，所以沒有人知道這家中心的學生錄取名校。反而是學生表現普通的留學顧問中心有可能借別人的名字宣傳學生錄取名校，提高招生業績。
LATE的假設是除了接受與不接受刺激之外，沒有其他的干擾因素。例如借學生名字宣傳的留學顧問可能有更多宣傳的財力，這樣就會造成干擾。
7.1 估計LATE
加入工具變數的迴歸模型可以估計LATE。也就是工具變數(instrumental variable, IV)解釋刺激D以及某一個自變數，但是與最後的表現Y無關。工具變數對於自變數的影響代表對於刺激的接受程度。
經由\(\texttt{ivreg}\)函數，可以估計刺激這一個變數的係數，代表LATE。
# Load necessary packages
#install.packages("AER") # Install AER package if not already installed
library(AER)
# Generate example data
set.seed(123)
n <- 1000 # Sample size
Z <- rnorm(n) # Instrumental variable
X <- rnorm(n) # Exogenous variable
D <- rbinom(n, 1, plogis(0.5 + Z + X)) # Endogenous variable (binary treatment)
Y <- 1 + D * 2 + X * 3 + rnorm(n) # Outcome variable
# Create a data frame
data <- data.frame(Y = Y, D = D, Z = Z, X = X)
# Estimate LATE using IV regression (2SLS)
iv_model <- ivreg(Y ~ D | Z + X, data = data)- 在這個模型中，我們假設D被Z影響，我們比較在相同X條件下，D的差異造成Y的差異。模型的估計語法如下：
# Summarize the model
stargazer::stargazer(iv_model, title = 'Instrumental Variable Model', type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))| Dependent variable: | |
| Y | |
| D | 11.500*** |
| (0.581) | |
| Constant | -4.745*** |
| (0.387) | |
| Observations | 1,000 |
| R2 | -0.507 |
| Adjusted R2 | -0.509 |
| Residual Std. Error | 4.550 (df = 998) |
| Note: | p<0.1; p<0.05; p<0.01 |
- 也可以直接輸出LATE:
# Extract LATE coefficient
LATE <- coef(iv_model)["D"]
# Print the LATE estimate
print(paste("Estimated LATE:", round(LATE, 3)))[1] “Estimated LATE: 11.503”
8 因果關係需要注意的議題
8.1 選樣偏誤(selection bias)
- 如果是調查實驗，隨機分配到實驗組與控制組的觀察對象最好可以代表母體。假設有一組有9成男性、1成女性，那麼實驗的結果無法推論到母體。
- 我們也要小心看不見的偏誤，例如小班制的老師跟一般班級的老師表現出來對學生的態度應該沒有太大差異，以免影響實驗結果。
8.2 比對兩組樣本的差異(Balance Check)
#sample gender and interest in study
star<-star %>% mutate(gender=sample(0: 1, 1274, prob = c(0.535, 0.465), replace = T)) %>%
mutate(interest=sample(1:5, 1274, prob = c(0.21,0.17,0.11,0.21,0.3),replace = T))
treat<-subset(star, star$D == 'small')
cont<-subset(star, star$D == 'regular')
roundfnc <- function(x, na.rm = F) round(x, 3)
treat.d <- star %>% select(gender, interest) %>%
dplyr::mutate_if(is.factor, as.numeric) %>%
summarise(across(everything(), list(M=mean, S=sd, max=max, min=min), na.rm=TRUE))
cont.d <- star %>% select(gender, interest) %>%
dplyr::mutate_if(is.factor, as.numeric) %>%
summarise(across(everything(), list(M=mean, S=sd, max=max, min=min), na.rm=TRUE))
treat.d<-roundfnc(treat.d); cont.d<-roundfnc(cont.d)
#png('~/C/Fig/balancestar.png', width=1500, height = 1200, res=200)
par(mfrow=c(1,2))
plot(density(treat$gender, na.rm=TRUE), xlim=c(-0.5,1.5), ylim=c(0,1.6), xlab="Gender", main='')
lines(density(cont$gender, na.rm=TRUE), col="RED")
plot(density(treat$interest, na.rm=TRUE), xlim=c(1,5), ylim=c(0,0.7), xlab="Interest", main='')
lines(density(cont$interest, na.rm=TRUE), col="RED")
8.3 平均值差異的標準誤
估計式除了無偏估計這一個標準外，估計值之間的離散程度應該越小越好。
-
樣本平均值的差異用\(\alpha=\bar{Y_{1}}-\bar{Y_{0}}\)。該估計值的變異數為\(\text{Var}(\hat{\alpha})=\frac{N}{N-1}\Bigl(\frac{\sigma^2_{Y1}}{n1}+\frac{\sigma^2_{Y0}}{n0}\Bigr)\)
Huber-White estimator如下，但是這不是無偏估計: \[ \hat{V}_{hw}=\frac{\frac{1}{n_{1}}\sum_{i}(Y_{i}-\bar{Y_{1}})^2}{n_{1}}+\frac{\frac{1}{n_{0}}\sum_{i}(Y_{i}-\bar{Y_{0}})^2}{n_{0}} \]
無偏估計應該是下面這個估計式：
\[ \hat{V}_{HC2}=\frac{\frac{1}{n_{1}-1}\sum_{i}(Y_{i}-\bar{Y_{1}})^2}{n_{1}}+\frac{\frac{1}{n_{0}-1}\sum_{i}(Y_{i}-\bar{Y_{0}})^2}{n_{0}} \]
- 其中，
\[\frac{1}{n_{1}-1}\sum_{i}(Y_{i}-\bar{Y_{1}})^2=\sigma_{1}^2 \]
\[\frac{1}{n_{0}-1}\sum_{i}(Y_{i}-\bar{Y_{0}})^2=\sigma_{0}^2\]
- 我們可以用
sandwich這個套件裡面的\(\texttt{vcovHC}\)函數，計算ATE的變異數。
# Load necessary library
library(sandwich)
# Estimate the variance of ATE using vcovHC
cov_matrix <- vcovHC(lm(reading ~ D, data = star), type = "HC3")
var_ATE <- cov_matrix[2, 2] # Variance of the coefficient for 'group'
# Print the results
#print(paste("Average Treatment Effect (ATE):", round(ATE, 2)))
print(paste("Variance of ATE:", round(var_ATE, 2)))[1] “Variance of ATE: 4.26”
- 根據以上的估計式，我們用聯合的標準誤檢定是否大到在統計上發生的機率非常小，如果真的非常小而我們又觀察到這樣的差異，那麼我們可以說真的存在這樣的差異。聯合的標準誤公式如下： \[ \sqrt{var(\hat{Y}_{treated}+\hat{Y}_{control})} =\sqrt{var(\hat{Y}_{treated}) + var(\hat{Y}_{control})} =\sqrt{\frac{\sigma^2_{T}}{N_{T}} + \frac{\sigma^2_{C}}{N_{C}}} \]
v.small <- var(star$reading[star$D == 'small']); v.small## [1] 1396v.regular <- var(star$reading[star$D == 'regular']); v.regular## [1] 1287n.group<- star %>% group_by (D) %>%
summarise(n = n())
n.small<-as.numeric(n.group[2,2]); n.regular<-as.numeric(n.group[1,2])
PoolSE<-sqrt((v.small/n.small)+(v.regular/n.regular))
print(paste("標準誤:", round(PoolSE, 2)))## [1] "標準誤: 2.06"- 另外一個估計ATE的變異數的公式為：
\[ \text{SE}(\widehat{ATE})=\sqrt{\frac{1}{N-1} \Bigl\{ \frac{m\times\text{Var_Y0}}{N-m}+\frac{(N-m)\text{Var_Y1}}{m}+2\text{Cov(Y0,Y1)} \Bigr\}} \]
- 其中，\(m\)為實驗組的個案數。如果兩個組的個案數不相等，那麼計算共變量時只能用相對應的個案，或者用插補的方式假定某些不存在的觀察值。
- 上述的公式以
R語法寫成如下。
N <-nrow(star)
m <- n.small #N of treat group
avg.small <- mean(star$reading[star$D=='small'])
avg.regular <- mean(star$reading[star$D=='regular'])
x.small <- star$reading[star$D=='small'] #treated
x.regular <- star$reading[star$D=='regular']
var_Y0 <-sum((x.regular - avg.regular)^2)/(N-m);
var_Y1 <- sum((x.small - avg.small)^2)/(m)
S <- (m * var_Y0 / (N - m)) + (((N - m) * var_Y1) / m) + 2 * (sum((x.small - avg.small) * (x.regular - avg.regular)) / N)
hat.ATE.var <- sqrt((1 / (N - 1)) * S)
print(paste("聯合標準誤:", round(hat.ATE.var, 2)))- 假設有一筆實驗設計的資料如下：
library(sandwich)
# Generate example data
set.seed(123) # for reproducibility
n <- 1000 # total sample size
n_treated <- 500 # sample size of treated group
n_control <- n - n_treated # sample size of control group
treated <- rnorm(n_treated, mean = 70, sd = 10) # simulate test scores for treated group
control <- rnorm(n_control, mean = 65, sd = 10) # simulate test scores for control group
group <- c(rep("Treated", n_treated), rep("Control", n_control))
data <- data.frame(group = group, score = c(treated, control))
# Calculate average outcomes for treated and control groups
avg_treated <- mean(data$score[data$group == "Treated"])
avg_control <- mean(data$score[data$group == "Control"])
# Calculate the ATE
ATE <- avg_treated - avg_control- 我們直接用迴歸模型估計ATE，並且應用sandwich這個套件中的\(\texttt{vcovHC}\)函數，估計變異數。
# Estimate the variance of ATE using vcovHC
cov_matrix <- vcovHC(lm(score ~ group, data = data), type = "HC3")
var_ATE <- cov_matrix[2, 2] # Variance of the coefficient for 'group'
# Print the results
print(paste("Average Treatment Effect (ATE):", round(ATE, 2)))## [1] "Average Treatment Effect (ATE): 5.37"print(paste("Variance of ATE:", round(var_ATE, 2)))## [1] "Variance of ATE: 0.39"print(paste("Standard deviation of ATE:", round(sqrt(var_ATE), 2)))## [1] "Standard deviation of ATE: 0.63"- 套用上面共變量的公式求ATE的標準差：
N <- 1000
m <- n_treated #N of treat group, 50
avg.treated <- avg_treated
avg.control <- avg_control
x.treated <- treated #treated
x.control <- control
var_Y0 <- sum((x.control - avg.control)^2)/(N-m);
var_Y1 <- sum((x.treated - avg.treated)^2)/m;
S <- (m * var_Y0 / (N - m)) + (((N - m) * var_Y1) / m) + 2 * (sum((x.control - avg.control) * (x.treated - avg.treated)) / N)
hat.ATE.var <- sqrt((1 / (N - 1)) * S)
print(paste("標準差:", round(hat.ATE.var, 2)))## [1] "標準差: 0.44"- 兩個公式算出的標準差很接近。當個案數越小，差別越大。
library(gtsummary)
star %>% tbl_summary(
by = D,
statistic = list(
all_continuous() ~ "{mean} ({sd})",
all_categorical() ~ "{n} / {N} ({p}%)"
),
digits = all_continuous() ~ 2) %>%
add_p(pvalue_fun = ~ style_pvalue(.x, digits = 2)) %>%
add_overall() %>%
add_n() %>%
modify_header(label ~ "**Variable**") | Variable | N | Overall, N = 1,2741 | regular, N = 6891 | small, N = 5851 | p-value2 |
|---|---|---|---|---|---|
| classtype | 1,274 | <0.001 | |||
| regular | 689 / 1,274 (54%) | 689 / 689 (100%) | 0 / 585 (0%) | ||
| small | 585 / 1,274 (46%) | 0 / 689 (0%) | 585 / 585 (100%) | ||
| reading | 1,274 | 628.80 (36.73) | 625.49 (35.88) | 632.70 (37.37) | <0.001 |
| math | 1,274 | 631.59 (38.84) | 628.84 (37.94) | 634.83 (39.66) | 0.013 |
| graduated | 1,274 | 1,108 / 1,274 (87%) | 597 / 689 (87%) | 511 / 585 (87%) | 0.71 |
| gender | 1,274 | 586 / 1,274 (46%) | 307 / 689 (45%) | 279 / 585 (48%) | 0.26 |
| interest | 1,274 | 0.91 | |||
| 1 | 256 / 1,274 (20%) | 134 / 689 (19%) | 122 / 585 (21%) | ||
| 2 | 214 / 1,274 (17%) | 114 / 689 (17%) | 100 / 585 (17%) | ||
| 3 | 116 / 1,274 (9.1%) | 64 / 689 (9.3%) | 52 / 585 (8.9%) | ||
| 4 | 297 / 1,274 (23%) | 167 / 689 (24%) | 130 / 585 (22%) | ||
| 5 | 391 / 1,274 (31%) | 210 / 689 (30%) | 181 / 585 (31%) | ||
| 1 n / N (%); Mean (SD) | |||||
| 2 Pearson’s Chi-squared test; Wilcoxon rank sum test | |||||
\[ SE_{p}=S_{p}\sqrt{\frac{1}{n_{1}}+\frac{1}{n_{2}}} \]
R計算如下：
v.small <- var(star$reading[star$D == 'small']); v.small## [1] 1396v.regular <- var(star$reading[star$D == 'regular']); v.regular## [1] 1287n.group<- star %>% group_by (D) %>%
summarise(n = n())
n.small<-n.group[2,2]; n.regular<-n.group[1,2]
p1 <- (((n.small-1)*v.small)+((n.regular-1)*v.regular))
PoolSD<-sqrt(p1/(n.small+n.regular-2));
PoolSE <- PoolSD*sqrt((1/n.small)+(1/n.regular));
print(paste("聯合標準誤:", round(PoolSE, 2)))## [1] "聯合標準誤: 2.06"- 兩個聯合標準誤的公式的計算結果一致，都是2.0左右。
- 因為2倍的標準誤仍然小於實驗效果(2*2.05<7.2)，可以結論實驗效果並不是隨機發生。換句話說，小班制真的會提高成績。
- 套用以上的標準誤公式，可以確定平均數的差異是否大到不是隨機發生。
9 配對（matching）
- 在實驗中，對對照组和實驗组的個體進行匹配，使它们在一些重要的特徵上相似，以消除混雜變量(confounder)的影響，提高研究的内部有效性。
- 將觀察值按照一些關鍵特徵如年龄、性别、教育程度等進行配對，使得實驗组和對照组在這些特徵上盡可能相似。
- 而在非實驗的觀察研究(observational study)中，因為各種原因，我們無法進行隨機分配的實驗。而且在給予刺激之前，實驗組與控制組的特徵有很大的差異，或者有沒有被觀察到的差異。透過匹配，我們可以讓實驗組與控制組接近一致。
- 匹配有兩個方式，第一個是讓每一對的觀察值在各項特徵盡量接近。第二個作法是在已經配對的樣本當中，各項特徵盡量相同。也許有一對樣本的特徵差別很大，但是要是其他對樣本的差異很小，整體的差異還是小。
9.1 傾向分數(propensity score)
- 我們可以透過傾向分數(propensity score)進行配對。傾向分數的作法是用二元勝算對數迴歸模型(logit model)估計觀察到的資料。
- 假設選舉制度越接近比例代表制，政黨數目越多。以往我們用迴歸模型描述選舉制度與政黨數目的關係，但是無法排除文化、其他制度等因素的影響，政黨數目也可能強化選舉制度。現在我們想估計兩者之間的因果關係，也就是如果採取新的制度而且比較接近比例代表制，會不會增加政黨數目？透過傾向分數找出因果關係的方式如下：
- 收集觀察對象（例如100多個國家）的各項特徵以及行為，例如最近有無採用新的選舉制度。
- 透過勝算對數迴歸模型，用這些特徵估計每一個國家是否改變選舉制度的機率。換句話說，我們觀察到每一個國家實際上有沒有改變選舉制度以及改變選舉制度的傾向(0-100%)。
- 然後我們根據傾向分數匹配有無採用新選舉制度的國家，例如有一個國家有高達90%的機率採用新選舉制度而且真的換了選舉制度，另一個國家雖然有85%的機率採用新選舉制度但是沒有換選舉制度，這兩個國家可以配成一對。
- 以此類推，我們可以得到50對左右的國家，然後我們統計有採用新選舉制度的50個國家的平均政黨數目，以及沒有採用新選舉制度的50個國家的平均政黨數目，兩者相減，就得到新選舉制度的平均效果(ATE)。
- 另一個作法是把每一個國家的傾向分數倒過來，乘以每一個國家的政黨數目，然後分別求出有採用新選舉制度以及沒有採用新選舉制度國家的平均政黨數目，相減得到ATE。
10 實驗設計的假設與類型
10.1 實驗假設
- 實驗設計需要隨機分配受試者到實驗組與控制組，並且確保給予刺激時，控制組不會受到影響，也就是SUTVA(stable unit treatment value assumption, SUTVA)的假設必須成立。
- 把學生分成兩組，一組要求每天運動30分鐘，另一組不做任何要求，1個月後測驗體能。如果控制組的受試者意識到他們體能越來越差也去做運動，那麼運動的效果就不會明顯。
- 建議徵求學生同意，一組每天運動30分鐘，另一組告知不要做任何運動，1個月後測驗體能，然後兩組互換，這樣對每個學生比較公平。
- 把學生分成兩組，一組要求每天運動30分鐘，另一組不做任何要求，1個月後測驗體能。如果控制組的受試者意識到他們體能越來越差也去做運動，那麼運動的效果就不會明顯。
10.2 實驗設計類型
- 實驗設計可分以下幾種：
- 研究者可以隨機分配，例如實驗室的實驗、調查實驗、田野實驗等等，研究者可以操弄刺激以估計因果關係，但是要遵守研究倫理的規範。
- 自然實驗：例如8月出生的比較晚入學，9月出生比同年級的同學大幾個月，可以觀察學業的表現。又或者有人被抽到到金門服役，有人在本島服役，到金門服役有可能因為比較接近敵人，所以比較擔心台灣會遭到攻擊。自然實驗雖然可以確定隨機給予刺激，但是要注意研究問題是不是具有社會科學的重要性。
- 統計控制：用匹配或者迴歸模型的控制變項，估計因果關係。
- 注意避免選樣偏誤，也就是讓受試者自行選擇加入實驗組或者控制組，因為可能有不同的因素影響選擇。例如徵求志願者來打疫苗，然後比較有打疫苗跟沒打疫苗的人得到感冒的比例。志願接種疫苗者可能不是因為打疫苗所以比較健康，而是因為他們平常就比較重視自己的身體狀況。
11 作業
請打開哈佛大學的民調資料(HXC23014 Harvard Poll Data.sav)，這筆資料有一個民調實驗，也就是把受訪者分成兩組，其中一組的受訪者先被提示中國在AI的專利申請數領先美國，然後問他們贊不贊成政府花費更多預算在AI上面(Q39)，另一個則是提示受訪者美國在AI專利申請數領先中國，然後問他們贊不贊成政府花費更多預算在AI上面(Q40)。資料中分組的變數是Q39_Q40Split。
- 請以圖形表示受試組與實驗組的次數分佈。
- 請問這兩組受訪者的贊成分數的平均數分別是多少？
- 請問這兩組受訪者的贊成分數的平均數差異為何？
- 請問聯合標準誤為何？
12 更新講義時間
最後更新時間: 2024-03-08 22:08:51