イベントスタディ分析
UPDATE: 2024-12-13 16:46:19.719486
はじめに
今回は、イベントスタディ分析の内容と実行方法をまとめておく。Pythonで学ぶ効果検証入門のDIDの章で紹介されているイベントスタディ分析を再現する。fixestパッケージについては下記が詳しい。fixestパッケージの使い方は末尾にまとめておく。
- Fast Fixed-Effects Estimation: Short Introduction
- 【R】固定効果モデルの推定がめっちゃ速いパッケージ { fixest}
- 【R】差の差法でイベントスタディやるときのコード
クラスター標準誤差については下記が詳しい。1つ目はfixestパッケージの作者が、fixestパッケージの標準誤差の計算方法について、小サンプル補正などを絡めて、他のパッケージとの違いを説明してくれているもの。
イベントスタディ分析
データはPythonで学ぶ効果検証入門のDIDの章で利用されているものを利用する。ここでもまずは書籍に従って、DIDによる施策の分析をRでも再現しつつ、イベントスタディ分析の再現まで行う。
このデータセットは、アメリカの臓器提供登録率が変化するかどうかをカルフォルニア州を介入群として実験したデータセット。
State: 27州が記録されているQuarter: Q42010 ~ Q12012の6期間Quarter_Num:QuarterのQ42010を1としたものRate: 登録率IsTreatmentGroup: 介入群であれば1、対照群であれば0AfterTreatment:Quarter_Numが4以降であれば1が立つ。これは期間を識別するフラグなので、介入の有無は関係ない。IsTreatment:Quarter_Numが4以降&Californiaであれば1。これは介入期間かつ介入群を識別するフラグ。
library(tidyverse)
library(fixest)
organ_full <- read_csv('~/Desktop/ch4_organ_donations_full.csv')
organ_full %>%
filter(State == 'Alaska' | State == 'California') %>%
print(n = 12)## # A tibble: 12 × 7
## State Quarter Rate Quarter_Num IsTreatmentGroup AfterTreatment IsTreatment
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Alaska Q42010 0.75 1 0 0 0
## 2 Alaska Q12011 0.77 2 0 0 0
## 3 Alaska Q22011 0.77 3 0 0 0
## 4 Alaska Q32011 0.78 4 0 1 0
## 5 Alaska Q42011 0.78 5 0 1 0
## 6 Alaska Q12012 0.79 6 0 1 0
## 7 Califo… Q42010 0.267 1 1 0 0
## 8 Califo… Q12011 0.273 2 1 0 0
## 9 Califo… Q22011 0.274 3 1 0 0
## 10 Califo… Q32011 0.264 4 1 1 1
## 11 Califo… Q42011 0.261 5 1 1 1
## 12 Califo… Q12012 0.264 6 1 1 1集計するとこのような結果になる。介入したカルフォルニア州では登録率が低下しているように見える
organ_full %>%
group_by(IsTreatmentGroup, Quarter_Num) %>%
summarise(avg_rate = mean(Rate)) %>%
ggplot(., aes(x = factor(Quarter_Num), y = avg_rate, color = factor(IsTreatmentGroup), group = factor(IsTreatmentGroup))) +
geom_line(size = 1.2) +
geom_point(size = 3) +
geom_vline(xintercept = 3.5) +
labs(
title = "Difference-in-Differences (DID) Plot",
x = "After Treatment (0 = Before, 1 = After)",
y = "Average Rate",
color = "Group"
) +
ylim(c(0.2, 0.5)) +
scale_color_manual(values = c("0" = "royalblue", "1" = "tomato"),
labels = c("0" = "Control Group", "1" = "Treatment Group")) +
theme_bw()
書籍のモデルを再現する。知りたいのは\(\tau\)である。
\[
Rate_{it} = \sum_{i=1}^{27} State_{i} + \sum_{t=1}^{6}Quarter\_Num_{t} +
\tau IsTreatment_{it} + \epsilon_{it}
\]
書籍の結果とクラスター標準誤差がずれるのは、末尾に背景をまとめている。結果をみると、-0.022となっており、「臓器提供フォームの文言を変えたことで登録率が下がった」ことを意味する。
f_did <- feols(
Rate ~ IsTreatment | State + Quarter_Num,
data = organ_full,
cluster = 'State'
)
etable(f_did)## f_did
## Dependent Var.: Rate
##
## IsTreatment -0.0225** (0.0061)
## Fixed-Effects: ------------------
## State Yes
## Quarter_Num Yes
## _______________ __________________
## S.E.: Clustered by: State
## Observations 162
## R2 0.97932
## Within R2 0.00922
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ここまでの再現で、データセットやDIDの結果への理解が深まったところで、イベントスタディ分析に移る。イベントスタディ分析の主たる目的は「介入後の効果は時間によって上下に変化する」ことを分析するところにある。介入後は効果が小さくても、時間が経過すると効果が大きくなる場合もある。
イベントスタディのモデルは下記の通り。
\[ Y_{it} = \sum_{i=1}^{27} State_{i} + \sum_{t=1}^{6}Quarter\_Num_{t} + \sum_{t=1}^{m-2} \nu_{t} W_{i} 1_{t} + \sum_{t=m}^{T} \tau_{t} W_{i} 1_{t} + \epsilon_{it} \] イベントスタディ分析では、さきほどの分析とは異なり\(\tau\)は1つではなく、\(m\)期の施策効果であれば、\(\tau_m\)、\(m+1\)期施策効果であれば\(\tau_{m+1}\)と複数存在します。総じて、\(m\)期から\(T\)期までの\(T-m\)期の施策効果が係数\(\tau_t\)として表現される。そのため、この係数\(\tau_t\)を推定することで、イベントの効果を推定する。つまり、\(T=6,m=4\)であれば、\(t=1,2\)と\(t=4,5,6\)となり、\(t=3\)は基準点なので、モデルからは除外する。簡略化したイベントスタディ分析のモデルは下記の通り。
\[ Rate_{it} = State_{i} + Quarter\_Num_{t} \\ + \nu_{1} IsTreatmenrGroup_{i} * Quarter\_Num_{1} \\ + \nu_{2} IsTreatmenrGroup_{i} * Quarter\_Num_{2} \\ + \nu_{4} IsTreatmenrGroup_{i} * Quarter\_Num_{4} \\ + \nu_{5} IsTreatmenrGroup_{i} * Quarter\_Num_{5} \\ + \nu_{6} IsTreatmenrGroup_{i} * Quarter\_Num_{6} \\ + \epsilon_{it} \\ \]
イベントスタディ分析を行う前に、少しデータを前処理しておく。ここではパッケージによる便利な関数の使用は控え、分析方法への理解を深める。
organ_event <- organ_full %>%
mutate(
# case文でもよい
QuarterNum_1 = if_else(Quarter_Num == 1, 1, 0),
QuarterNum_2 = if_else(Quarter_Num == 2, 1, 0),
QuarterNum_4 = if_else(Quarter_Num == 4, 1, 0),
QuarterNum_5 = if_else(Quarter_Num == 5, 1, 0),
QuarterNum_6 = if_else(Quarter_Num == 6, 1, 0),
)
organ_event %>%
filter(State == 'Alaska' | State == 'California') %>%
as.data.frame()## State Quarter Rate Quarter_Num IsTreatmentGroup AfterTreatment
## 1 Alaska Q42010 0.7500 1 0 0
## 2 Alaska Q12011 0.7700 2 0 0
## 3 Alaska Q22011 0.7700 3 0 0
## 4 Alaska Q32011 0.7800 4 0 1
## 5 Alaska Q42011 0.7800 5 0 1
## 6 Alaska Q12012 0.7900 6 0 1
## 7 California Q42010 0.2666 1 1 0
## 8 California Q12011 0.2731 2 1 0
## 9 California Q22011 0.2743 3 1 0
## 10 California Q32011 0.2636 4 1 1
## 11 California Q42011 0.2607 5 1 1
## 12 California Q12012 0.2641 6 1 1
## IsTreatment QuarterNum_1 QuarterNum_2 QuarterNum_4 QuarterNum_5 QuarterNum_6
## 1 0 1 0 0 0 0
## 2 0 0 1 0 0 0
## 3 0 0 0 0 0 0
## 4 0 0 0 1 0 0
## 5 0 0 0 0 1 0
## 6 0 0 0 0 0 1
## 7 0 1 0 0 0 0
## 8 0 0 1 0 0 0
## 9 0 0 0 0 0 0
## 10 1 0 0 1 0 0
## 11 1 0 0 0 1 0
## 12 1 0 0 0 0 1実行結果がこちら。
# fl = 'Rate ~
# QuarterNum_1*IsTreatmentGroup +
# QuarterNum_2*IsTreatmentGroup +
# QuarterNum_4*IsTreatmentGroup +
# QuarterNum_5*IsTreatmentGroup +
# QuarterNum_6*IsTreatmentGroup | State + Quarter_Num'
f_pre <- paste0("IsTreatmentGroup", "*QuarterNum_", c(1:2, 4:6), collapse = " + ")
fl <- paste0('Rate ~ ', f_pre, ' | State + Quarter_Num')
f_event <- feols(
fml = as.formula(fl),
data = organ_event,
cluster = 'State'
)
etable(f_event)## f_event
## Dependent Var.: Rate
##
## IsTreatmentGroup x QuarterNum_1 -0.0029 (0.0051)
## IsTreatmentGroup x QuarterNum_2 0.0063** (0.0023)
## IsTreatmentGroup x QuarterNum_4 -0.0216*** (0.0050)
## IsTreatmentGroup x QuarterNum_5 -0.0203*** (0.0045)
## IsTreatmentGroup x QuarterNum_6 -0.0222* (0.0100)
## Fixed-Effects: -------------------
## State Yes
## Quarter_Num Yes
## _______________________________ ___________________
## S.E.: Clustered by: State
## Observations 162
## R2 0.97934
## Within R2 0.00979
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1あとはこれを可視化する。
f_event_df <- tibble(
label = c(names(f_event$coefficients), 't=3'),
time = c(1:2, 4:6, 3),
coef = c(f_event$coefficients, 0),
se = c(f_event$se, 0),
n = nrow(organ_event)) %>%
mutate(
ci_lwr = coef - qt(0.975, df = n)*se,
ci_upr = coef + qt(0.975, df = n)*se
)
f_event_df## # A tibble: 6 × 7
## label time coef se n ci_lwr ci_upr
## <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 IsTreatmentGroup:QuarterNum_1 1 -0.00294 0.00508 162 -0.0130 0.00710
## 2 IsTreatmentGroup:QuarterNum_2 2 0.00630 0.00227 162 0.00182 0.0108
## 3 IsTreatmentGroup:QuarterNum_4 4 -0.0216 0.00503 162 -0.0315 -0.0116
## 4 IsTreatmentGroup:QuarterNum_5 5 -0.0203 0.00447 162 -0.0291 -0.0115
## 5 IsTreatmentGroup:QuarterNum_6 6 -0.0222 0.0100 162 -0.0419 -0.00239
## 6 t=3 3 0 0 162 0 0プロットを見る限り、施策介入後は-2%近辺を推移しているため、施策実行後から時間が経過しても、この期間においては変化してないことがわかる。
ggplot(f_event_df, aes(x = time, y = coef)) +
geom_line(color = 'royalblue') +
geom_point(color = 'royalblue') +
geom_ribbon(aes(ymin = ci_lwr, ymax = ci_upr), alpha = 0.2, fill = 'royalblue') +
geom_vline(xintercept = 3.5, linetype = 'dotted', color = 'gray50') +
geom_hline(yintercept = 0, linetype = 'dotted', color = 'gray50') +
scale_x_continuous(breaks = 1:6, labels = paste("T", 1:6, sep = "")) +
coord_cartesian(ylim = c(-0.05, 0.05)) +
labs(
title = "Event Study Plot",
x = "Time",
y = "Coefficient"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
axis.title = element_text(size = 14),
axis.text = element_text(size = 12)
)
DIDを用いて介入効果を計測するためにはパラレルトレンド仮定が成立していることが前提になる。イベントスタディを実行することで、パラレルトレンド仮定を検証できる。先程の回帰モデルでは、介入前の時点でもパラメタを推定している。このパラメタは、施策を実施していないにも関わらず、介入群であることが、アウトカムへ与える影響を意味している。さきほどの可視化では、介入前の部分では0を中心に推移しているため、この結果から「介入前には施策効果はなかった」と結論を導けます。ただパラレルトレンド仮定は根本的に実証不可な点は注意したい。
介入時期が一律ではない場合
さきほどの例では、介入時期がすべての州で同一だったが、実際は介入時期がずれる可能性がある。例えば、前回のノートで扱ったMixtapeの「Castle-doctrine statutes and homicides」の章のデータセットは、異なるタイミングで行われた銃の規制に関するデータセットになっている。
その場合は、下記の画像の通り、分析前に介入時点を0として、前後のレコードを識別するフラグを作成する処理を加える必要がある。他の国では、lead7,8,9やlag5が存在しているが、ここでの説明で使用している国は該当してない。
Alabamaに介入したのは2006年なので、2006年を基準に前後のダミーフラグを作成する。2005年がlead1で、2007年がlag1になる。Texasに介入したのは2007年なので、2007年を基準に前後のダミーフラグを作成する。2006年がlead1で、2008年がlag1になる。Wyomingに介入していないのでNAしかない。説明のために｀NA`にしているが、実際には0で埋める。
castle <- read_csv('~/Desktop/castle.csv')
castle2 <- castle %>%
select(l_homicide, popwt, cdl, state, year, treatment_date) %>%
mutate(
time_til = year - treatment_date,
lead9 = case_when(time_til == -9 ~ 1, TRUE ~ 0),
lead8 = case_when(time_til == -8 ~ 1, TRUE ~ 0),
lead7 = case_when(time_til == -7 ~ 1, TRUE ~ 0),
lead6 = case_when(time_til == -6 ~ 1, TRUE ~ 0),
lead5 = case_when(time_til == -5 ~ 1, TRUE ~ 0),
lead4 = case_when(time_til == -4 ~ 1, TRUE ~ 0),
lead3 = case_when(time_til == -3 ~ 1, TRUE ~ 0),
lead2 = case_when(time_til == -2 ~ 1, TRUE ~ 0),
lead1 = case_when(time_til == -1 ~ 1, TRUE ~ 0),
lag0 = case_when(time_til == 0 ~ 1, TRUE ~ 0),
lag1 = case_when(time_til == 1 ~ 1, TRUE ~ 0),
lag2 = case_when(time_til == 2 ~ 1, TRUE ~ 0),
lag3 = case_when(time_til == 3 ~ 1, TRUE ~ 0),
lag4 = case_when(time_til == 4 ~ 1, TRUE ~ 0),
lag5 = case_when(time_til == 5 ~ 1, TRUE ~ 0)
)
event_study_formula <- as.formula(
paste("l_homicide ~ ",
paste(
paste(paste("lead", 1:9, sep = ""), collapse = " + "),
paste(paste("lag", 1:5, sep = ""), collapse = " + "), sep = " + "),
"| year + state"
),
)
event_study_reg <- lfe::felm(
event_study_formula,
weights = castle2$popwt,
data = castle2)
summary(event_study_reg)##
## Call:
## lfe::felm(formula = event_study_formula, data = castle2, weights = castle2$popwt)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -787.73 -166.65 6.99 150.48 884.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## lead1 -0.038693 0.034802 -1.112 0.2668
## lead2 0.002261 0.035967 0.063 0.9499
## lead3 0.006643 0.036762 0.181 0.8567
## lead4 -0.006975 0.037114 -0.188 0.8510
## lead5 0.018607 0.037330 0.498 0.6184
## lead6 0.053712 0.040069 1.340 0.1807
## lead7 -0.063573 0.050582 -1.257 0.2094
## lead8 -0.203023 0.082432 -2.463 0.0141 *
## lead9 -0.285209 0.296400 -0.962 0.3364
## lag1 0.079512 0.034807 2.284 0.0228 *
## lag2 0.085099 0.036101 2.357 0.0188 *
## lag3 0.079865 0.038332 2.084 0.0377 *
## lag4 0.041934 0.043856 0.956 0.3395
## lag5 0.089586 0.066579 1.346 0.1791
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.3 on 476 degrees of freedom
## Multiple R-squared(full model): 0.9367 Adjusted R-squared: 0.927
## Multiple R-squared(proj model): 0.06774 Adjusted R-squared: -0.07523
## F-statistic(full model):96.52 on 73 and 476 DF, p-value: < 2.2e-16
## F-statistic(proj model): 2.471 on 14 and 476 DF, p-value: 0.002202後は先程と同様に可視化すれば良い。
# order of the coefficients for the plot
plot_order <- c("lead9", "lead8", "lead7",
"lead6", "lead5", "lead4", "lead3",
"lead2", "lead1", "lag1",
"lag2", "lag3", "lag4", "lag5")
leadslags_plot <- tibble(
sd = c(event_study_reg$se[plot_order], 0),
mean = c(coef(event_study_reg)[plot_order], 0),
label = c(-9,-8,-7,-6, -5, -4, -3, -2, -1, 1,2,3,4,5, 0)
) %>%
mutate(
ymin = mean-1.96*sd,
ymax = mean+1.96*sd
)
ggplot(leadslags_plot, aes(x = label, y = mean)) +
geom_line(color = 'royalblue') +
geom_point(color = 'royalblue') +
geom_ribbon(aes(ymin = ymin, ymax = ymax), alpha = 0.2, fill = 'royalblue') +
geom_vline(xintercept = 0, linetype = 'dotted', color = 'gray50') +
geom_hline(yintercept = 0, linetype = 'dotted', color = 'gray50') +
scale_x_continuous(breaks = -9:5) +
coord_cartesian(ylim = c(-0.5, 0.5)) +
xlab("Years before and after castle doctrine expansion") +
ylab("log(Homicide Rate)") +
theme_bw()
fixestパッケージの使用方法
作者のイントロダクションを見ながらfixestパッケージの使用方法を簡単にまとめておく。
基本的にはRで回帰を実行するときのようにフォーミュラを記載すれば良い。固定効果を追加したい場合は、|で区切って変数を追加していく。
feols(y ~ x1 + x2 + x3 + x4 | a1 + a2 + a3 + a4, data)organ_full %>%
group_by(IsTreatmentGroup, AfterTreatment) %>%
summarise(avg_rate = mean(Rate))## # A tibble: 4 × 3
## # Groups: IsTreatmentGroup [2]
## IsTreatmentGroup AfterTreatment avg_rate
## <dbl> <dbl> <dbl>
## 1 0 0 0.445
## 2 0 1 0.459
## 3 1 0 0.271
## 4 1 1 0.2632x2差分の差分法を回帰モデルで実行する方法は下記の通り。
feols(Rate ~ IsTreatmentGroup*AfterTreatment, data = organ_full) ## OLS estimation, Dep. Var.: Rate
## Observations: 162Standard-errors: IID
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.444931 0.017047 26.100751 < 2.2e-16 ***
## IsTreatmentGroup -0.173597 0.088577 -1.959846 0.051773 .
## AfterTreatment 0.013926 0.024108 0.577645 0.564326
## IsTreatmentGroup:AfterTreatment -0.022459 0.125267 -0.179289 0.857941
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.148682 Adj. R2: 0.036196固定効果モデルを実行する方法は下記の通り。クラスター標準誤差ももちろん利用可能。cluster引数に指定すれば良い。
f1 <- feols(
Rate ~ IsTreatment | State + Quarter_Num,
data = organ_full,
cluster = 'State'
)
summary(f1)## OLS estimation, Dep. Var.: Rate
## Observations: 162Fixed-effects: State: 27, Quarter_Num: 6Standard-errors: Clustered (State)
## Estimate Std. Error t value Pr(>|t|)
## IsTreatment -0.022459 0.006131 -3.66304 0.0011185 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.021982 Adj. R2: 0.974196
## Within R2: 0.009221他の書き方として、summary()で渡しても良いとのこと。
f2 <- feols(
Rate ~ IsTreatment | State + Quarter_Num,
data = organ_full
)
summary(f2, cluster = 'State')## OLS estimation, Dep. Var.: Rate
## Observations: 162Fixed-effects: State: 27, Quarter_Num: 6Standard-errors: Clustered (State)
## Estimate Std. Error t value Pr(>|t|)
## IsTreatment -0.022459 0.006131 -3.66304 0.0011185 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.021982 Adj. R2: 0.974196
## Within R2: 0.009221固定効果の2つのカテゴリについてクラスター標準誤差を計算したい場合、summary()の中でse = "twoway"を利用する。1つの場合はcluster、2つの場合はtwoway、3つの場合はthreeway、4つの場合はfourwayを利用する。他にも、standard、heteroも利用できる。
f3 <- feols(
Rate ~ IsTreatment | State + Quarter_Num,
data = organ_full
)
summary(f3, se = 'twoway')## OLS estimation, Dep. Var.: Rate
## Observations: 162Fixed-effects: State: 27, Quarter_Num: 6Standard-errors: Clustered (State & Quarter_Num)
## Estimate Std. Error t value Pr(>|t|)
## IsTreatment -0.022459 0.005197 -4.32173 0.0075582 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.021982 Adj. R2: 0.974196
## Within R2: 0.009221もしくは、vcov引数を利用しても同じ。オプションはcluster、twoway、newey_west、driscoll_kraay、conleyが用意されている。
f4 <- feols(
Rate ~ IsTreatment | State + Quarter_Num,
data = organ_full,
vcov = 'twoway'
)
summary(f4)## OLS estimation, Dep. Var.: Rate
## Observations: 162Fixed-effects: State: 27, Quarter_Num: 6Standard-errors: Clustered (State & Quarter_Num)
## Estimate Std. Error t value Pr(>|t|)
## IsTreatment -0.022459 0.005197 -4.32173 0.0075582 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.021982 Adj. R2: 0.974196
## Within R2: 0.009221クラスター標準誤差の計算については、パッケージの作者が下記の通りnoteを用意してくれている。一部抜粋すると、
標準誤差は推定の重要な要素であると言うのは婉曲表現です。文字通り、論文の結果は標準誤差に依存します。したがって、標準誤差を計算する従来の「最良」の方法が存在しないのは残念です。 たとえば、さまざまなソフトウェアでまったく同じ推定を実行すると、異なる標準誤差が得られることは珍しくありません。最初に「バグがあるに違いない」と思ったとしても、その考えは脇に置いてください。バグなど存在しないのですから。多くの場合、それは開発者が小規模サンプル補正に関して行った選択に帰着しますが、驚くべきことに、実装に関しては多くの自由度があります。 複数の定義があると混乱が生じる可能性があるため、このドキュメントの目的は、このパッケージの標準誤差計算の面倒な詳細を明らかにすることです。
例えば、estimatrパッケージの下記の実行結果と同じもの(詳細までは調べていない)を得るには、小規模サンプル補正やクラスター標準誤差の計算方法を調整する必要がある。
estimatr::lm_robust(Rate ~ IsTreatment,
fixed_effects = State + Quarter_Num,
data = organ_full,
clusters = State,
se_type = 'CR0'
)## Estimate Std. Error t value Pr(>|t|) CI Lower
## IsTreatment -0.02245897 0.005903444 -3.804385 0.0007770557 -0.03459368
## CI Upper DF
## IsTreatment -0.01032427 26sscは小規模サンプル補正のオプション。
feols(
Rate ~ IsTreatment | State + Quarter_Num,
data = organ_full,
cluster = 'State',
ssc = ssc(adj = FALSE, cluster.adj = FALSE)
)## OLS estimation, Dep. Var.: Rate
## Observations: 162Fixed-effects: State: 27, Quarter_Num: 6Standard-errors: Clustered (State)
## Estimate Std. Error t value Pr(>|t|)
## IsTreatment -0.022459 0.005903 -3.80439 0.00077706 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.021982 Adj. R2: 0.974196
## Within R2: 0.009221fixef()を利用すれば、固定効果のパラメタを確認できる。いくつかのサンプルについてのみ、plot(fixef(f5))で可視化できる。
f5 <- feols(
Rate ~ IsTreatment | State,
data = organ_full,
cluster = 'State'
)
fixef(f5)## $State
## Alaska Arizona California
## 0.7733333 0.2404667 0.2713333
## Colorado Connecticut District of Columbia
## 0.6678500 0.3931500 0.3432833
## Florida Hawaii Louisiana
## 0.3963833 0.4209167 0.5566000
## Maryland Michigan Minnesota
## 0.4657167 0.2324667 0.5276167
## Missouri Montana Nebraska
## 0.4073500 0.6548833 0.4425167
## New Hampshire New Jersey New York
## 0.5441333 0.3125833 0.1297833
## North Carolina Ohio Pennsylvania
## 0.5249833 0.5656000 0.4542833
## South Carolina Tennessee Virginia
## 0.2690500 0.3352000 0.3409500
## Washington Wisconsin Wyoming
## 0.5870167 0.5721167 0.5910000
##
## attr(,"class")
## [1] "fixest.fixef" "list"
## attr(,"exponential")
## [1] FALSEetable()を利用すれば、複数のモデルを説明と共にレイアウトして表示してくれる。
etable(f1, f5)## f1 f5
## Dependent Var.: Rate Rate
##
## IsTreatment -0.0225** (0.0061) -0.0085*** (1.77e-18)
## Fixed-Effects: ------------------ ---------------------
## State Yes Yes
## Quarter_Num Yes No
## _______________ __________________ _____________________
## S.E.: Clustered by: State by: State
## Observations 162 162
## R2 0.97932 0.97702
## Within R2 0.00922 0.00125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1本ノートでも説明しているイベントスタディ分析も簡単にできる。
f_event <- feols(
Rate ~ i(Quarter_Num, IsTreatmentGroup, ref = 3) | State + Quarter_Num,
data = organ_full,
cluster = 'State'
)
f_event## OLS estimation, Dep. Var.: Rate
## Observations: 162Fixed-effects: State: 27, Quarter_Num: 6Standard-errors: Clustered (State)
## Estimate Std. Error t value Pr(>|t|)
## Quarter_Num::1:IsTreatmentGroup -0.002942 0.005084 -0.578719 0.56775872
## Quarter_Num::2:IsTreatmentGroup 0.006296 0.002266 2.778832 0.00999724 **
## Quarter_Num::4:IsTreatmentGroup -0.021565 0.005034 -4.284177 0.00022209 ***
## Quarter_Num::5:IsTreatmentGroup -0.020292 0.004473 -4.536282 0.00011432 ***
## Quarter_Num::6:IsTreatmentGroup -0.022165 0.010013 -2.213610 0.03583451 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.021976 Adj. R2: 0.973385
## Within R2: 0.009787iplot()で可視化できる。
iplot(f_event)