社會科學統計方法

隨機實驗與因果關係

1 統計的目的

有三種統計的用途。

1. 測量：研究者收集資料，根據不同的分布形態，找出中央趨勢、計算離散程度等等，例如：

• 統計不同國家民眾對於政治的興趣；
• 統計排隊多久才能打到流感疫苗；
• 統計擲銅板得到正面的比例。有一篇文章探討擲硬幣的實驗結果。
• 統計文章出現的關鍵字與主題的機率。

2. 預測：

• 收集資料，建立一個模型，然後預測下一個事件發生的機率，如果發生機率接近實際發生的事件，則可以繼續使用該模型。
• 最近流行的人工智慧(AI)，屬於預測的一種，不過生成式AI是產生新的資料而不是預測，MIT有一篇文章Explained: Generative AI討論這些概念。

3. 因果關係：我們可以透過隨機實驗，排除所有可能的干擾因素，估計出兩個變數的因果關係。

• 例如：對一群人隨機抽樣，有人抽到籤去服役，有人不需要去服役，然後在若干年後詢問這些人對於國家的認同感，這些人之間唯一的差別是有沒有去服役，就可以估計出服役與否對國家認同感的作用(Erikson, 2011)。

• 例如：對叛亂組織轟炸後，如果這些組織存活，轟炸會降低還是提高這些組織的地位？比對有受到與沒受到轟炸的叛亂組織的後續發展，可以估計轟炸的效果(Lyall, 2009)。

• 例如：在1997，英國太陽報從支持保守黨轉而支持工黨，而在2009，太陽報從支持工黨又轉回保守黨。藉由統計報紙的立場轉變前後的讀者態度，可以估計媒體的影響 (Reeves, McKee, Stuckler, 2016)。

• 過去統計比較重視預測，但是近年來越來越重視因果關係，因為因果關係的研究可以告訴我們相關背後的機制，進而延伸到其他領域的研究。但是因果關係需要好的研究設計以及適當的資料，有時候可遇不可求。

• 另一方面，實驗設計雖然可以確認因果關係，但是是不是具有理論意涵？政策意涵？還是只是發現造成差異的條件？需要進一步地研究。

2 因果關係的實例

本週將介紹因果關係與隨機實驗。

我們用因果關係來解釋兩個變數的關係。例如在美國，家庭收入越高，學童的成績越好，這之間的因果關係是兒童越早開始學習，越早啟發腦力，但是美國的學前教育收費昂貴，只有一定收入的家庭才負擔得起，所以越能負擔得起的家庭，學童的表現越好。

Figure 2.1: 學前教育

過去認為如果A發生在前、B發生在後，所以如果A發生而且B也發生，A與B之間可能有因果關係。

但是這樣的關係其實沒有排除可能的干擾因素。例如父母身高跟子女身高的關係看起來是因果關係，但是我們並無法觀察到未生育民眾的子女的身高，只觀察到已生育民眾與其子女的身高，所以身高來自遺傳？還是來自後天營養？如何進行遺傳的實驗並排除環境的影響？

另外一個例子是喜歡吃辣是不是有益健康？吃辣會增加排汗，幫助身體排毒？理論上我們可以請一組人常常吃辣，另外一組人吃清淡食物，然候觀察他們的身心狀況。即使我們可以隨機分配兩組人吃辣以及不吃辣，但是實際生活中喜歡吃辣可能本身就是身體健康，越吃越過癮，不喜歡的人可能自己覺得身體無法負荷，看到辣就全身冒汗。所以吃辣的行為跟身體健康與否互相影響，難以估計因果關係。

在勞動經濟學，有一個著名的研究是提高最低工資是否會造成勞工失業？也就是雇主不堪最低工資的負擔，寧可犧牲服務品質也要減少雇用勞工？或者說，每提高多少最低工資會提高多少失業率？Card and Krueger (1994) 調查兩個相鄰的州–New Jersey, Pennsylvania的餐廳員工人數以及時薪 ，其中賓州的最低時薪一直維持不變，而紐澤西把最低工資從4.25美元提高到5.05美元，假設其他條件都一樣，紐澤西的失業率在調高最低工資後增加了嗎？而賓州的失業率都不變嗎？圖2.2顯示兩州的地圖。

Figure 2.2: 餐廳員工數與最低時薪

圖2.3顯示，兩個州的政策前後的平均雇用人數的差異，以及標準誤。研究發現，紐澤西餐廳的員工人數在提高時薪之後增加了，同一時間賓州的餐廳的員工人數減少了。雖然這個結論聽起來不合乎常理，然而有可能賓州的員工因為最低時薪選擇到紐澤西就業，或者員工更加努力工作使得餐廳的生意彌補了提高時薪的成本，使得餐廳並沒有因此解雇員工。不論如何，這個difference-in-difference的實驗設計，啟發了許多後續研究。David Card與兩位學者因為勞動經濟學的研究，獲得2021年的諾貝爾經濟學獎。

Figure 2.3: 餐廳員工數與最低時薪

• 大家可以想一個可行的實驗設計嗎？

3 實驗設計的原理

從前面的例子可以看出實驗設計的基本原理是我們觀察兩組群體在固定一個差異情況下，比較群體之間的差異。這是Ronald Fisher 提出的實驗設計，透過隨機分組到控制組與實驗組 ，然後用變異數分析或者區塊(blocking)設計，分析實驗組與控制組的差異。

所以我們把一群人隨機分組，然後假設這兩組人除了給予的刺激之外都沒有不同之處。而刺激的效果就是這兩組人表現差異的平均。

如果我們只觀察一個人或是一群人，然後比較有無給予某個刺激或者某個政策的差異，稱為反事實的統計(counterfactual)，但是事實上我們不可能做到，例如請一位病患既服用安慰劑又服用高血壓藥，也不可能找一個班級學生同時增加數學課又減少數學課的上課時間，然後比較教學方法的差異。

Donald Rubin提出潛在反應(potential outcome)的概念，每個人在接受刺激跟沒接受刺激情況下，分別有不同的反應，這兩者之間的差異就是因果關係。

潛在反應的模型特別是用在無法隨機分組的觀察資料，透過配對(matching)、傾向分數(propensity score)、工具變數(instrumental variable)估計實驗效果。

3.1 Potential Outcome Model

我們用Potential Outcome Model定義以上的原理。

首先是實驗刺激(D)：

\(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]\)，也就是明明澆水，但是盆栽好像被蓋子蓋起來，完全沒接收到水分。換句話說，我們只能觀察到有給予刺激的結果以及沒給予刺激的結果，但是觀察不到有給予刺激卻好像沒給予的結果，也觀察不到沒給予刺激卻好像有得到刺激的結果。

根據方程式1，Y有兩個結果，一個是有接收到刺激，一個是沒接收到刺激，這兩個Y可以寫在一個方程式裡面：

\(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}\)是有無刺激的反應相減，得到：

\[\alpha_{i} = Y_{1i}-Y_{0i}\]

等一下我們會顯示\(\alpha_{i}\)是\(Y_{1i}\)的平均數與\(Y_{0i}\)的平均數的差異。

為了確保以上的模型得到正確的估計，我們必須假設\(D_{i}\)對每一個個體都有同樣的影響，也沒有任何的外溢效果。寫成以下的等式：

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

我們定義平均實驗效果(Average Treatment Effect, ATE)為：

\[\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] \end{align}\]

假設我們有4個受試者，其中兩個屬於實驗組，兩個在控制組，在第一個時間點測量結果表示為\(Y_{0i}\)，第2個時間點測量結果為\(Y_{1i}\)，整理如下表：

\(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	0	1	0
4	1	1	0	0	0
\(E[Y_{1}]\)	1.5
\(E[Y_{0}]\)		0.5
\(E[Y_{1}-Y_{0}]\)				1

• 計算第2個時間點或者假設受到刺激的測量結果的平均值\(E[Y_{1}]\)： \[ E[Y_{1}] = \frac{1}{N}\Sigma Y_{1i}=1.5 \]

• 計算第1個時間點或者假設沒得到刺激的測量結果的平均值\(E[Y_{0}]：\)$ E[Y_{0}] = 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 \]

• 因為 \[Y_{i}=Y_{1}\cdot D+(1-D)\cdot Y_{0}\]

• 所以\(\sum Y_{i}=3+1+0+0=4\)。\(E[Y_{i}]=1\)。

ATE也稱為difference-in-means估計式(estimator)，在隨機實驗的條件成立下，可以正確估計因果關係。

3.2 Naive TE

如果我們只看\(Y_{11}\)與\(Y_{12}\)這兩個觀察值，計算\(E[Y_{i}|D=1]\)得到2，計算\(E[Y_{i}|D=0]\)得到0.5，兩者相差為1.5，與\(\bar{\alpha_{i}}\)以及\(E[Y_{1}]-E[Y_{0}] = 1\)不同。我們稱\(E[Y_{i}|D=1]-E[Y_{i}|D=0]\)為Naive treatment effect(NATE)。

在上面表格中，\(Y_{11}=3\)，而\(Y_{12}=1\)。也就是說，\(D_{i}\)有不同的效果。這時候會造成NATE不等於ATE。

如果把\(Y_{11}\)改成1，那麼ATE=NATE。換句話說，NATE有誤差，需要先假設刺激對於每個受試者的效果都相同。

3.3 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_{0}|D=0]\)代表某種誤差。也就是說：

\[ \text{ATE}=\text{ATT}+Bias \]

• 如果\(E[Y_{0}|D=1]-E[Y_{0}|D=0]= 0\)，那麼ATE = ATT，代表我們不見得要把受試者分成實驗組跟控制組，只要觀察實驗組在兩個時間點的變化就可以估計實驗效果。

計算Average Treatment Effects on the Treated (ATT)如下：

\(i\)	\(Y_{1i}\)	\(Y_{0i}\)	\(D_{i}\)	\(\alpha_{i}\)	\(Y_{i}\)
1	3	0	1	3	3
2	1	1	1	0	1
\(E[Y_{1}]\)	2
\(E[Y_{0}]\)		0.5
\(E[Y_{1}-Y_{0}]\)				1.5

\[ \alpha_{ATT}=E[Y_{1}-Y_{0}|D=1]=\frac{1}{2}\cdot(3+1-0-1)=1.5 \]

• \(\alpha_{ATT}\)與\(\alpha_{i}\)相差0.5。

• \(\alpha_{ATT}\)與\(\alpha_{i}\)相等的前提是\(E[Y_{0}|D=1]-E[Y_{0}|D=0]= 0\)。也就是說：實際上收到刺激但是「假設」並沒有收到的表現狀況，跟實際上沒收到刺激而且真的沒有收到刺激的表現狀況相等。例如：給植物澆水，但是沒注意到水都澆到盆栽外面去，跟沒澆水的表現狀況相同。或者給頭痛患者服藥，但是因為頭太痛而忘了吃藥，藥效跟沒吃一樣。如果這些情形為真，表示就算我們無法觀察到控制組的受試者「真的」得到刺激的表現，我們可以假設控制組的受試者即使「真的」得到刺激，他們的表現跟沒有得到刺激的受試者相同。

ATE指的是所有參與者的因果關係的估計，ATT指的是對於給予刺激者的因果關係的估計：

\[ \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。因為受到刺激的個人或者集體對於刺激的反應不是完全相同，所以我們可以用ATT說明刺激的作用。但是ATE代表有無接收到刺激的平均差異，比較具有說服力。

• 假設有兩群人，一群人吃頭痛藥以減緩頭痛，另一群人不吃頭痛藥。ATT代表吃頭痛藥的這一群人之後的平均頭痛狀況，ATE則代表吃頭痛藥與不吃頭痛藥之後的平均頭痛狀況差異。如果頭痛藥真的對每個人都有效，ATE會大於ATT，這是因為沒吃藥的人的平均頭痛狀況不會減輕，吃藥的人才會減輕，ATE代表的差異就會很明顯。

  • 當\(E[Y_{0}|D=1]>E[Y_{0}|D=0]\)，Bias>0，ATE>ATT。
  • 如果控制組「真的」得到刺激，例如控制組的植物真的被澆水了，理論上成長狀況應該比沒澆到水來得好。換句話說，有沒有澆水對植物的影響應該是很顯著的，這時候我們就必須分組觀察有無澆水，而不能只觀察澆過水的植物。
  • 對於沒來上課的同學，「假如」要求他們來上這堂課，不一定會比沒上課的同學對於實驗設計了解更多，因為來上這門課可能因為昨天晚上沒睡好、老師講得不清楚等等原因，反而比在家自習的學習效果差，所以Bias < 0，ATE < ATT，代表我們如果只觀察來上課的同學對於實驗設計的了解，有可能會比較好。
  • 例如政府提供免費的職業訓練，所有人可以免費報名參加，我們只能觀察實際上有參加職訓的人之後的收入，但是觀察不到沒參加的人，要估計ATT只能假設或者觀察過去的資料。沒來參加訓練的人，有可能本身不愁收入，也可能在同一個時間自己學習或者努力賺錢，反而可能比剛剛參加完職業訓練的人的收入來得高，也就是\(E[Y_{0}|D=1]<E[Y_{0}|D=0]\)，Bias > 0，ATE>ATT。

在實際生活中我們不見得可以隨機分配受試者為實驗組、控制組。如果我們可以找到一個隨機的事件，觀察研究對象在事件發生時間前後的變化，同樣可以進行實驗研究。

• 例如每次選舉中都有兩種候選人，一種是現任者，另一種是挑戰者。假設成為現任者是隨機發生的，例如在里長選舉，有些候選人可能差3票當選，有些則是差2票落選。那麼我們就可以用一些變數解釋這些以些微差距當選的現任者，然後把預測值變成倒數，加權給這些現任者，然後觀察他們的得票率有沒有比上次提高，等同於現任者的優勢。相關的討論請點Jens Hainmueller的文章。

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因為D的有無造成的平均變動程度。

\[ Y=\beta_{0} +\beta_{1}D+\beta_{2}X \]

4.1 實際操作：小班制與成績

DSS裡面的STAR.csv這筆資料是分組成小班與大班之後的成績資料，我們用來練習估計因果關係。首先讀取資料：

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         0

如果每一個班級的學生人數隨機分為兩種，一種是小班(small)，另一種是一般人數(regular)，我們想估計班級人數與閱讀(reading)成績的因果關係，也就是應用difference-in-means的估計式，我們先轉換一下班級人數這個變數的性質為類別，並且給定一個新變數叫做D。

star <- 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"

然後分別計算\(E[Y|D=\text{regular}]\)與\(E[Y|D=\text{small}]\)：

mean(star$reading[star$D=='small'])

## [1] 632.7

mean(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

我們可以用gtsummary::tbl_summary()這個功能呈現每個連續與類別變數對應特定變數的平均值或者比例：

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

結果是7.21，代表班級人數變少，學生平均閱讀分數會提高7.21分。

我們可以進一步計算數學成績的差異，以及畢業比例的差異。前提是這兩組學生是隨機分配。

• 用廻歸模型估計實驗效果：

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"))

Table 4.1: OLS and Causal Inference

	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

4.2 超過一個以上的實驗組

我們可以估計超過一個以上的實驗效果。例如Gerber et al. (2008)的投票實驗中，把18萬的家戶其中的10萬設定為控制組，另外8萬分成四組分別寄送信件。第一種信件是提到投票是公民責任，第二種是提到公民責任之外，告訴受訪者他們的投票紀錄是公開的，第三種除了提到公民責任還提示受訪者過去的投票紀錄，第四種除了提到公民責任還提示鄰居的投票紀錄。Gerber等人的全文在此。

預測受訪者是否投票\(Y_{i}\)的固定效果模型如下：

\[ Y_{i}=\beta_{0}+\beta_{1}D_{\text{1}i}+\beta_{2}D_{\text{2}i}+\beta_{3}D_{\text{3}i}+\beta_{4}D_{\text{4}i}+\sum_{k=1}^{K-1}\gamma_{k}C_{ki} \]

• \(\sum_{k=1}^{K-1}\gamma_{k}C_{ki}\)指的是1萬個地理區的固定效果。

5 次團體(subgroup): CATE

我們可以觀察實驗組中的某一個次團體跟控制組相對應的次團體之間的平均效果差異，這個差異稱為conditional average treatment effect (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}] \]

假設有許多自變數X，CATE等於：

\[\alpha_{x}=E[Y_{1}-Y_{0}|X = x]\]

• 當模型中有非常多自變數時，必須要用機器學習的模型，降低變數的維度或者面向(dimensions)，才能估計CATE。

CATE需要假設隨機分配的機制與觀察到的自變數之間互相獨立。寫成：

\[ Y_{1},Y_{0}\mathrel{\unicode{x2AEB}} T|X \]

5.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"))

Table 5.1: Interaction Model

	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的影響，圖5.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 5.1: Interacton Plot

5.2 估計高維度的CATE

hdcate使用機器學習降低自變數的維度，然後估計LATE:

library(hdcate)
# get simulation data
n_obs <- 300 # Num of observations
n_var <- 50 # Num of observed variables
n_rel_var <- 4 # Num of relevant variables
data <- HDCATE.get_sim_data(n_obs, n_var, n_rel_var)
# conditional expectation model is misspecified
x_formula <- paste(paste0('X', c(2:n_var)), collapse ='+')
# for example, and alternatively, the propensity score model is misspecified
# x_formula <- paste(paste0('X', c(1:(n_var-1))), collapse ='+')
# Example 1: full-sample estimator
# create a new HDCATE model
model <- HDCATE(data=data, y_name='Y', d_name='D', x_formula=x_formula)
# estimate HDCATE function, inference, and plot
HDCATE.set_condition_var(model, 'X2', min=-1, max=1, step=0.01)
HDCATE.fit(model)
HDCATE.inference(model)
HDCATE.plot(model)

6 Local Average Treatment Effect (LATE)

在實驗組中每一個個案應該都接收到刺激，控制組則不應該接收到刺激。有時候有些受試者不確實接受刺激，有些則跑去接受刺激，但是研究者無法禁止。

除了刺激\(D\)之外，還有決定給予\(D_{1}\)或者\(D_{0}\)的\(Z\)，也就是treatment assignment。

• 例如：政府無法強迫所有長輩都來接受流感疫苗。為了觀察長輩接種疫苗的實驗效果，政府決定送1,000元超商禮券給來打疫苗的長輩。但是為了確認這個政策可行，政府抽出1,000位民眾，其中又分成400位通知來疫苗可以獲得禮券，600位只通知要來打疫苗。被分到實驗組的民眾可能為了禮券而來打疫苗，有些沒被分到實驗組的民眾則是無論有沒有禮券都會來打疫苗，也有些民眾被分到控制組就不會去打疫苗。LATE等於實驗組且有打疫苗跟控制組而沒來打疫苗的比例。

• 另一個例子是Gerber et al.(2010)隨機打電話給兩群民眾，一群是通知要做好回收，一群是通知要記得去投票。因為有些民眾經常在家而且會接電話，有些民眾則常常不在家，就算在家也不會接電話，如果假設接電話與否是隨機，雖然實務上可能不是，通知去投票而且有接電話以及通知做好回收而且接電話是compliers，這兩群人的投票率差異等於LATE。

\[ \text{LATE}=\frac{E[Y_{1}-Y_{0}|D = 1]}{\text{Pr}(D=1)} \]

• 如果我們可以假設刺激的效果是因為某些外在的因素而有差異，也就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)。

  • 例如民眾按照住家附近一定的投票人口數被分配到附近投票所投票，但可能因為各種因素，造成有的投票所很快就投完票，有的投票所要排很久，\(Z\)可能是隨機，與\(D\)獨立，也可能是因為投票所人員難找，最後只找到不熟練的人員，造成投票所速度緩慢，因此投票所的分派不再是隨機，而是受到其他因素影響。被分到大排長龍的投票所的民眾，也許不想去排隊投票。反之被分到沒怎麼排隊的投票所的民眾，可能就完成投票，LATE就等於這兩種投票所之間的投票率差異。

6.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"))

Table 6.1: Instrumental Variable Model

	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”

7 因果關係需要注意的議題

7.1 選樣偏誤(selection bias)

有時候我們無法控制所有的變項。例如有人有參加職業訓練，有人沒參加，有參加職業訓練的人不止學到新的技能，而且參加的人原本就比較有企圖心，所以職業訓練對於未來收入的效果可能被高估。

出現選樣偏誤的原因是我們無法觀察\(E[Y_{0}|D=1]\)。也就是如果我們請民眾來參加職業訓練，但是他們沒有參加也沒有辦法聯絡，我們觀察不到這些人的未來工作收入。又或者像是在紐澤西的餐廳，如果有餐廳剛好在實施最低工資因為各種原因歇業，那麼我們無法統計到這些餐廳在實施最低工資之後的員工人數。

• 如果是調查實驗，隨機分配到實驗組與控制組的觀察對象最好可以代表母體。假設有一組有9成男性、1成女性，那麼實驗的結果無法推論到母體。
• 我們也要小心看不見的偏誤，例如小班制的老師跟一般班級的老師表現出來對學生的態度應該沒有太大差異，以免影響實驗結果。

7.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")

7.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] 1396

v.regular <- var(star$reading[star$D == 'regular']); v.regular

## [1] 1287

n.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"

• 兩個公式算出的標準差很接近。當個案數越小，差別越大。

回到小班制實驗，\(\texttt{gtsummary}\)可以列出使用到標準誤的t檢定的結果，結果是t<0.001，代表發生的機率很小，兩個樣本之間的差異不是偶然發生。

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

另一個計算聯合標準誤的方式，先計算聯合標準差，再計算標準誤，聯合標準差的公式：

\[ S_{p}=\sqrt{\frac{(n_{1}-1)s_{1}^2+(n_{2}-1)s_{2}^2}{n_{1}+n_{2}-2}} \]

聯合標準誤的公式變成：

\[ 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] 1396

v.regular <- var(star$reading[star$D == 'regular']); v.regular

## [1] 1287

n.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)，可以結論實驗效果並不是隨機發生。換句話說，小班制真的會提高成績。
• 套用以上的標準誤公式，可以確定平均數的差異是否大到不是隨機發生。

8 配對（matching）

• 在實驗中，對對照组和實驗组的個體進行匹配，使它们在一些重要的特徵上相似，以消除混雜變量(confounder)的影響，提高研究的内部有效性。
• 將觀察值按照一些關鍵特徵如年龄、性别、教育程度等進行配對，使得實驗组和對照组在這些特徵上盡可能相似。
• 而在非實驗的觀察研究(observational study)中，因為各種原因，我們無法進行隨機分配的實驗。而且在給予刺激之前，實驗組與控制組的特徵有很大的差異，或者有沒有被觀察到的差異。透過匹配，我們可以讓實驗組與控制組接近一致。
• 匹配有兩個方式，第一個是讓每一對的觀察值在各項特徵盡量接近。第二個作法是在已經配對的樣本當中，各項特徵盡量相同。也許有一對樣本的特徵差別很大，但是要是其他對樣本的差異很小，整體的差異還是小。

8.1 傾向分數(propensity score)

• 我們可以透過傾向分數(propensity score)進行配對。傾向分數的作法是用二元勝算對數迴歸模型(logit model)估計觀察到的資料。

0. 假設選舉制度越接近比例代表制，政黨數目越多。以往我們用迴歸模型描述選舉制度與政黨數目的關係，但是無法排除文化、其他制度等因素的影響，政黨數目也可能強化選舉制度。現在我們想估計兩者之間的因果關係，也就是如果採取新的制度而且比較接近比例代表制，會不會增加政黨數目？透過傾向分數找出因果關係的方式如下：
1. 收集觀察對象（例如100多個國家）的各項特徵以及行為，例如最近有無採用新的選舉制度。
2. 透過勝算對數迴歸模型，用這些特徵估計每一個國家是否改變選舉制度的機率。換句話說，我們觀察到每一個國家實際上有沒有改變選舉制度以及改變選舉制度的傾向(0-100%)。
3. 然後我們根據傾向分數匹配有無採用新選舉制度的國家，例如有一個國家有高達90%的機率採用新選舉制度而且真的換了選舉制度，另一個國家雖然有85%的機率採用新選舉制度但是沒有換選舉制度，這兩個國家可以配成一對。
4. 以此類推，我們可以得到50對左右的國家，然後我們統計有採用新選舉制度的50個國家的平均政黨數目，以及沒有採用新選舉制度的50個國家的平均政黨數目，兩者相減，就得到新選舉制度的平均效果(ATE)。
5. 另一個作法是把每一個國家的傾向分數倒過來，乘以每一個國家的政黨數目，然後分別求出有採用新選舉制度以及沒有採用新選舉制度國家的平均政黨數目，相減得到ATE。

8.1.1 實例

模擬200個學生的資料：

library(dplyr)

# Simulate data
set.seed(02138)

N <- 200  # Total students

df <- data.frame(
  ID = 1:N,
  prior_gpa = rnorm(N, mean = 3.0, sd = 0.5),  # Prior performance
  treatment = sample(c(0,1), N, replace = TRUE),  # Randomly assigned treatment (1 = new method, 0 = traditional)
  motivation = rnorm(N, mean = 50, sd = 10)  # Random motivation scores
)

# Generate final scores based on treatment effect
df$final_score <- with(df, 75 + 5*treatment + 2*prior_gpa + 0.5*motivation + rnorm(N, mean=0, sd=5))

假設在一定的成績條件以及動機以上給予刺激：

# Simulate observational study (self-selection)
df_obs <- df %>%
  mutate(treatment = ifelse(prior_gpa > 3.2 & motivation > 50, 1, sample(c(0,1), N, replace = TRUE, prob = c(0.5, 0.5))))

• 估計接收到刺激的機率成為配對分數：

# Estimate propensity scores (probability of receiving treatment)
ps_model <- glm(treatment ~ prior_gpa + motivation, data = df_obs, family = binomial)
df_obs$pscore <- predict(ps_model, type = "response")

• 用MatchIt套件，進行比對：

# Load libraries for causal inference
library(MatchIt)

# Perform propensity score matching
matched <- matchit(treatment ~ prior_gpa + motivation, data = df_obs, method = "nearest", distance = "logit")

• 用t檢定確認實驗刺激的效果是否在統計上顯著：

# Extract matched dataset
df_matched <- match.data(matched)

# Compare final scores in matched dataset
t.test(final_score ~ treatment, data = df_matched)

## 
##  Welch Two Sample t-test
## 
## data:  final_score by treatment
## t = -4.1, df = 164, p-value = 6e-05
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -6.918 -2.431
## sample estimates:
## mean in group 0 mean in group 1 
##           106.3           110.9

• 連同其他的變數進行迴歸分析如表 8.1 ：

Table 8.1: 配對後迴歸模型

	Dependent variable:
	final_score
	(1)	(2)
treatment	4.674***	0.172
	(1.136)	(1.013)
prior_gpa		2.288**
		(0.994)
motivation		0.539***
		(0.048)
Constant	106.300***	73.920***
	(0.803)	(4.078)
Observations	166	166
R2	0.094	0.495
Adjusted R2	0.088	0.486
Residual Std. Error	7.320 (df = 164)	5.498 (df = 162)
F Statistic	16.920*** (df = 1; 164)	52.910*** (df = 3; 162)
Note:	p<0.1; p<0.05; p<0.01

• 在沒有控制以前的gpa以及動機之前，分組的確影響成績，但是控制之後，分組的影響變得不顯著。

9 實驗設計的假設與類型

9.1 實驗假設

• 實驗設計需要隨機分配受試者到實驗組與控制組，並且確保給予刺激時，控制組不會受到影響，也就是SUTVA(stable unit treatment value assumption, SUTVA)的假設必須成立。
  • 把學生分成兩組，一組要求每天運動30分鐘，另一組不做任何要求，1個月後測驗體能。如果控制組的受試者意識到他們體能越來越差也去做運動，那麼運動的效果就不會明顯。
    • 建議徵求學生同意，一組每天運動30分鐘，另一組告知不要做任何運動，1個月後測驗體能，然後兩組互換，這樣對每個學生比較公平。

9.2 實驗設計類型

• 實驗設計可分以下幾種：
  1. 研究者可以隨機分配，例如實驗室的實驗、調查實驗、田野實驗等等，研究者可以操弄刺激以估計因果關係，但是要遵守研究倫理的規範。
  2. 自然實驗：例如8月出生的比較晚入學，9月出生比同年級的同學大幾個月，可以觀察學業的表現。又或者有人被抽到到金門服役，有人在本島服役，到金門服役有可能因為比較接近敵人，所以比較擔心台灣會遭到攻擊。自然實驗雖然可以確定隨機給予刺激，但是要注意研究問題是不是具有社會科學的重要性。
  3. 統計控制：用匹配或者迴歸模型的控制變項，估計因果關係。
  4. 注意避免選樣偏誤，也就是讓受試者自行選擇加入實驗組或者控制組，因為可能有不同的因素影響選擇。例如徵求志願者來打疫苗，然後比較有打疫苗跟沒打疫苗的人得到感冒的比例。志願接種疫苗者可能不是因為打疫苗所以比較健康，而是因為他們平常就比較重視自己的身體狀況。

10 作業

1. 請打開哈佛大學的民調資料(HXC23014 Harvard Poll Data.sav)，這筆資料有一個民調實驗，也就是把受訪者分成兩組，其中一組的受訪者先被提示中國在AI的專利申請數領先美國，然後問他們贊不贊成政府花費更多預算在AI上面(Q39)，另一個則是提示受訪者美國在AI專利申請數領先中國，然後問他們贊不贊成政府花費更多預算在AI上面(Q40)。資料中分組的變數是Q39_Q40Split。

(1) 請以圖形表示受試組與實驗組的次數分佈。
(2) 請問這兩組受訪者的贊成分數的平均數分別是多少？
(3) 請問這兩組受訪者的贊成分數的平均數差異為何？
(4) 請問聯合標準誤為何？

2. 請用以下模擬的資料，用age, education, motivation建立傾向分數模型，配對之後，計算income的差異，並以迴歸模型顯示。

# 1. Simulate Data (Replace with your actual data)
set.seed(02138) # For reproducibility

n <- 1000 # Sample size
age <- sample(20: 70, n, replace = TRUE)
education <- sample(1:4, n, replace = TRUE) # 1=low, 4=high
motivation <- rnorm(n, 0, 1) # Higher values indicate higher motivation
D <- rbinom(n, 1, prob = 0.3 + 0.2 * (education > 2) + 0.1 * motivation) # Treatment assignment depends on covariates

D[is.na(D)]<-1

income <- 20000 + 5000 * education + 3000 * motivation + 8000 * D + rnorm(n, 0, 10000)

df_obs <- data.frame(age, education, motivation, D, income)

11 更新講義時間

最後更新時間: 2025-03-14 22:00:42