Data Analysis Using Statistical Packages: Card and Krueger (AER 1994)
Sho Kuroda / 黒田翔
2024年3月 (Last update: 2024年4月)
1 Overview
Card, David and Alan B. Krueger. 1994. Minimum Wages and Employment: A Case Study of the Fast-Food Industry in New Jersey and Pennsylvania, American Economic Review 84(4): 772-793.
Question. 最低賃金の引き上げは雇用を引き下げるか？
1992年4月に米国ニュージャージー州で最低賃金が引き上げられた（$4.25 → $5.05）一方で，隣接するペンシルベニア州ではとどめ置かれた． これを自然実験とみなし，DIDの枠組みを利用する． ファストフード店の雇用に関する調査データを用いて，最低賃金の引き上げによって雇用がほとんど変化しなかったことを明らかにした．
1.1 Reproduction/Replication by others
- Aaron Mamula (2021) Replicating Card and Krueger’s Minimum Wage Study in R
- Paul C. Bauer (2020) Section 10.5 in Applied Causal Analysis (with R)
- 矢内勇生 (2023) 差分の差分法（計量経済学応用）
- アウトカムの定義が元論文とは異なる点に注意
- 浅野正彦 (2022) 差分の差分法:
DID
- アウトカムの定義が元論文とは異なる点に注意．矢内先生の資料と重複する箇所が多い
1.2 Data
データは 著者 Card 本人が提供している data sets 一覧 または Mostly Harmless Econometrics の Data Archive からダウンロードできる．
Readme （または read.me）ファイルを読むと，「public.dat」がデータセットで，「codebook」がその名の通り対応する codebook （変数の名前やラベルなどの情報）である．
.dat は汎用のファイル拡張子で，テキストエディタで開くと
46 4 0 0 0 0 0 1 0 0 0 30.00 15.00 …
のような数字が入力されている．
Moodle には CardKrueger1994.dat というファイル名でアップロードしている．
- 1994年2月（最低賃金の引き上げ前；Wave 1）の変数名が
xyzの場合，対応する1994年11-12月（最低賃金の引き上げ後；Wave 2）の変数名はxyz2で定義されている - p. 775 によれば，Full-time-equivalent [FTE] employment = the number of full-time workers [incl. managers] + 0.5 × the number of part-time workers
- DID
分析の都合上，2月と11月の両方で賃金データが観測されていることを表す
logical 変数
wage_availableを作成しておく
.dat を読み込むには read.table 関数が使える．
library(tidyverse)
setwd("c://ws_stat")
ck_dat <- read.table("CardKrueger1994.dat", header = FALSE, dec = ".", na.strings = ".")
ck <- ck_dat %>%
rename(sheet = V1,
chain = V2,
co_owned = V3,
state = V4,
# location
south_NJ = V5,
central_NJ = V6,
north_NJ = V7,
PA1 = V8,
PA2 = V9,
shore = V10,
# first interview (February 1994)
ncalls = V11,
empft = V12,
emppt = V13,
nmgrs = V14,
wage_st = V15,
inctime = V16,
firstinc = V17,
bonus = V18,
pctaff = V19,
meals = V20,
open = V21,
hrsopen = V22,
psoda = V23,
pfry = V24,
pentree = V25,
nregs = V26,
nregs11 = V27,
# second interview (November 1994)
type2 = V28,
status2 = V29,
date2 = V30,
ncalls2 = V31,
empft2 = V32,
emppt2 = V33,
nmgrs2 = V34,
wage_st2 = V35,
inctime2 = V36,
firstinc2 = V37,
special2 = V38,
meals2 = V39,
open2r = V40,
hrsopen2 = V41,
psoda2 = V42,
pfry2 = V43,
pentree2 = V44,
nregs2 = V45,
nregs112 = V46) %>%
mutate(NJ = south_NJ + central_NJ + north_NJ,
state = ifelse(NJ == 1, "NJ", "PA"),
fte = empft + nmgrs + emppt / 2,
fte2 = empft2 + nmgrs2 + emppt2 / 2,
wage_available = !is.na(wage_st) & !is.na(wage_st2))2 Descriptive statistics
2.1 The number of stores interviewed (Table 1)
州ごとの店舗数．
table(ck$state)##
## NJ PA
## 331 79閉店すると FTE が 0 で記録されている (Table 3 の Notes より)．
ck %>%
group_by(state) %>%
summarise(closed = sum(fte2 == 0, na.rm = TRUE))## # A tibble: 2 × 2
## state closed
## <chr> <int>
## 1 NJ 5
## 2 PA 12.2 Means of key variables (Table 2)
Panel 1. Distribution of Store Types (%)
ck %>%
group_by(state) %>%
summarise(BurgerKing = mean(chain == 1),
KFC = mean(chain == 2),
RoyRogers = mean(chain == 3),
Wendys = mean(chain == 4),
CompnayOwned = mean(co_owned))## # A tibble: 2 × 6
## state BurgerKing KFC RoyRogers Wendys CompnayOwned
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NJ 0.411 0.205 0.248 0.136 0.341
## 2 PA 0.443 0.152 0.215 0.190 0.354一番右の列は「母平均が2つの州で等しい」を帰無仮説とする t 値．
ck %>%
mutate(BurgerKing = ifelse(chain == 1, 1, 0)) %>%
t.test(BurgerKing ~ state, data = .)##
## Welch Two Sample t-test
##
## data: BurgerKing by state
## t = -0.5152, df = 116.88, p-value = 0.6074
## alternative hypothesis: true difference in means between group NJ and group PA is not equal to 0
## 95 percent confidence interval:
## -0.1557951 0.0914714
## sample estimates:
## mean in group NJ mean in group PA
## 0.4108761 0.4430380t 値 (Welch’s t) の定義から計算することもできる．
\[ t = \frac{\bar{x}_{NJ} - \bar{x}_{PA}}{\sqrt{s_{NJ}^2 / n_{NJ} + s_{PA}^2 / n_{PA}}} \]
ck %>%
mutate(BurgerKing = ifelse(chain == 1, 1, 0)) %>%
group_by(state) %>%
summarise(
n = n(), # count (sample size)
x_bar = mean(BurgerKing),
x_var = var(BurgerKing)) %>%
print.data.frame(., digits = 5)## state n x_bar x_var
## 1 NJ 331 0.41088 0.24279
## 2 PA 79 0.44304 0.24992(.41088 - .44304) / sqrt(.24279/331 + .24992/79)## [1] -0.5151671または
ck_BK <- ck # copy
ck_BK$BurgerKing <- ifelse(ck$chain == 1, 1, 0)
x1 <- ck_BK$BurgerKing[ck_BK$state == "NJ"]
x2 <- ck_BK$BurgerKing[ck_BK$state == "PA"]
(mean(x1) - mean(x2)) / sqrt(var(x1) / length(x1) + var(x2) / length(x2))## [1] -0.5151975Panel 2. Means in Wave 1
カッコの中は標準偏差 (standard deviation; データのばらつき) ではなく標準誤差 (standard error; 推定量のばらつき)．
\[ SE = \frac{\sigma}{\sqrt{n}} \]
ck %>%
group_by(state) %>%
summarise(sample_size = sum(!is.na(fte)),
FTE = mean(fte, na.rm = TRUE),
FTE_se = sd(fte, na.rm = TRUE) / sqrt(sample_size),
StartingWage = mean(wage_st, na.rm = TRUE),
StartingWage_se = sd(wage_st, na.rm = TRUE) / sqrt(sample_size),
Wage425 = mean(wage_st == 4.25, na.rm = TRUE),
Wage425_se = sd(wage_st == 4.25, na.rm = TRUE) / sqrt(sample_size))## # A tibble: 2 × 8
## state sample_size FTE FTE_se StartingWage StartingWage_se Wage425 Wage425_se
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NJ 321 20.4 0.508 4.61 0.0193 0.322 0.0261
## 2 PA 77 23.3 1.35 4.63 0.0401 0.342 0.0544「Wage = $4.25 (percentage)」変数は再現できない．要確認．
t 値も先ほどと同様に以下のように計算できる．
ck %>%
t.test(fte ~ state, data = .)##
## Welch Two Sample t-test
##
## data: fte by state
## t = -2.0032, df = 98.562, p-value = 0.0479
## alternative hypothesis: true difference in means between group NJ and group PA is not equal to 0
## 95 percent confidence interval:
## -5.75630044 -0.02722102
## sample estimates:
## mean in group NJ mean in group PA
## 20.43941 23.33117余裕がある履修者向けの宿題：上で計算していない他の変数やWave2についても同様に計算する．
2.3 Distribution of starting wage rates (Figure 1)
最低賃金の引き上げ前と引き上げ後のそれぞれにおいて，NJとPAで starting wage （初任給に相当）の時間給の分布がどのように異なるかをプロットする．
NJ と PA を1つのグラフに同時に描画するには以下のように
fill と position="dodge"
を指定すればよい．
ggplot(data = ck, mapping = aes(x = wage_st, fill = state)) +
geom_histogram(position = "dodge") +
scale_fill_grey()
しかしながら，NJとPAで調査対象数が大きく異なるためうまく比較できない． そのため，NJとPAそれぞれでの割合を縦軸にしたい． 「ggplot2 histogram percentage by group」でGoogle検索すると以下のページがヒットする．
Stack Overflow: Let ggplot2 histogram show classwise percentages on y axis
ggplot(data = ck, mapping = aes(x = wage_st, fill = state)) +
geom_histogram(
aes(y = after_stat(c(count[group==1]/sum(count[group==1]), count[group==2]/sum(count[group==2])) * 100) ),
position = "dodge") +
ylab("% of stores") +
scale_fill_grey()
ggplot(data = ck, mapping = aes(x = wage_st2, fill = state)) +
geom_histogram(
aes(y = after_stat(c(count[group==1]/sum(count[group==1]), count[group==2]/sum(count[group==2])) * 100) ),
position = "dodge") +
ylab("% of stores") +
scale_fill_grey()
よく見ると < 5.05 がある． 僅かに法令不遵守の店舗があったという意味かもしれない．
ck %>% filter(state == "NJ") %>% select(wage_st2) %>% table## wage_st2
## 5 5.05 5.1 5.14 5.15 5.2 5.25 5.28 5.3 5.4 5.5 5.67 5.75
## 1 283 1 1 1 3 14 1 1 1 8 2 13 Employment effect of the minimum-wage increase (Section III)
3.1 III.A. DID (Table 3)
After-Before の difference は次のように計算できる． Difference の標準誤差が論文でレポートされている値と一致しない．
ck %>%
group_by(state) %>%
summarise(sample_size = sum(!is.na(fte)),
FTE_Feb = mean(fte, na.rm = TRUE),
FTE_Feb_se = sd(fte, na.rm = TRUE) / sqrt(sample_size),
FTE_Nov = mean(fte2, na.rm = TRUE),
FTE_Nov_se = sd(fte2, na.rm = TRUE) / sqrt(sample_size)) %>%
mutate(FTE_diff = FTE_Nov - FTE_Feb,
FTE_diff_se = sqrt(FTE_Feb_se^2 + FTE_Nov_se^2)) %>%
print.data.frame(., digits = 5)## state sample_size FTE_Feb FTE_Feb_se FTE_Nov FTE_Nov_se FTE_diff FTE_diff_se
## 1 NJ 321 20.439 0.50826 21.027 0.51869 0.58802 0.7262
## 2 PA 77 23.331 1.35115 21.166 0.94322 -2.16558 1.6478NJ-PA の difference は以下の通り．
20.439 - 23.331 # diff in February ## [1] -2.892sqrt(0.50826^2 + 1.35115^2) # standard error of diff in February ## [1] 1.44358421.027 - 21.166 # diff in November ## [1] -0.139sqrt(0.51869^2 + 0.94322^2) # standard error of diff in November ## [1] 1.076431DID 推定量．
(21.027 - 21.166) - (20.439 - 23.331) # DID estimator ## [1] 2.753sqrt(1.443584^2 + 1.076431^2) # standard error of DID ## [1] 1.800733要約
| Variable | PA | NJ | Difference: NJ-PA |
|---|---|---|---|
| FTE employment: Before (February) | 23.33 (1.35) | 20.44 (0.51) | -2.89 (1.44) |
| FTE employment: After (November) | 21.17 (0.94) | 21.03 (0.52) | -0.14 (1.08) |
| Change in FTE employment: After-Before | -2.17 (1.65) | 0.59 (0.73) | 2.75 (1.80) |
Note: Standard errors in parentheses.
標準誤差は After-Before の difference から計算しても一緒． ただし，論文の値とは一致しない（原因不明）．
sqrt(1.6478^2 + 0.7262^2)## [1] 1.8007253.2 III.B. Regression-adjusted models (Table 4)
ここまで使ってきたデータは wide 型．
Wide 形式：1行＝1店舗（すなわち，1行に複数時点のデータが含まれる）
| Store ID | Chain | NJ | FTE1 | FTE2 |
|---|---|---|---|---|
| 46 | 1 | 0 | 40.05 | 24.00 |
| 49 | 2 | 0 | 13.75 | 11.50 |
| … | … | … | … | … |
Long 形式：1行＝1観測（すなわち，1行には1時点のデータしか含まれない）
| Store ID | Chain | NJ | Wave | FTE |
|---|---|---|---|---|
| 46 | 1 | 0 | 1 | 40.05 |
| 46 | 1 | 0 | 2 | 24.00 |
| 49 | 2 | 0 | 1 | 13.75 |
| 49 | 2 | 0 | 2 | 11.50 |
| … | … | … | … | … |
3.2.1 Wide format
回帰分析でDID推定量を求めることもできる． ここで，一階の差分をとる（11月の値から2月の値を引く）と個体の固定効果が打ち消されるため，アウトカムの差分を treatment dummy （つまり NJ ダミー）に回帰すればよい．
回帰分析の枠組みで DID を推定することで，employment growth に影響を与えうる要因を条件づけることができる． たとえば，チェーンを説明変数として加えることで，チェーン間での growth の違いをコントロールしたもとで DID 推定ができる．
いずれも Table 4 でレポートされている推定値とは少し異なる．
Model (i)
\[ \Delta E_i = a + c NJ_i + \epsilon_i \]
この式は \(E_{i, Feb} = \mbox{baseline} + \gamma NJ_i\) と \(E_{i, Nov} = \mbox{baseline} + \gamma NJ_i + c NJ_i + a\) の差分をとることで得られる． \(\gamma NJ_i\) は処置の前後で共通する NJ 固有の効果であり（差分を取ることで相殺される），\(c NJ_i\) は NJ で処置後にのみ生じる部分を表している（差分をとった後に残る）． \(c\) が DID 推定量である． \(a\) は NJ と PA の両方で共通する処置後にのみ生じる部分で，time trend として解釈できる．
did_reg_i <- lm(formula = fte2 - fte ~ NJ, data = ck, subset = wage_available)
# summary(did_reg_i)Model (ii)
\[ \Delta E_i = a + b X_i + c NJ_i + \epsilon_i \]
did_reg_ii <- lm(formula = fte2 - fte ~ NJ + factor(chain) + co_owned, data = ck, subset = wage_available)Model (iii)
NJ の店舗でも starting wage が処置前から既に $5.05 を上回っている店舗では雇用に対する影響がほとんどないと考えられる一方，処置前には最低賃金 $4.25 に近かった店舗では雇用への影響が大きいと考えられる． そこで，このような within-NJ の変分を表す変数を下のように定義する．
\[ \Delta E_i = a' + c' GAP_i + \epsilon_i, \quad GAP = \begin{cases} 0 \quad \mbox{if PA} \\ 0 \quad \mbox{if NJ with } W_{i, Feb} \ge 5.05 \\ \frac{5.05 - W_{i, Feb}}{W_{i, Feb}} \quad \mbox{if NJ with } W_{i, Feb} < 5.05 \end{cases} \]
つまり，GAP 変数は
- 次の2つの変分を含む： (i) NJ と PA の違い，(ii) NJ の処置前の starting wage の違い
- DID の枠組みに位置付けると，Model (i), (ii) での NJ 変数に対応する．つまり，GAP 変数の回帰係数は DID 推定量のようなものに対応する
ck_with_gap <- ck %>%
mutate(gap = ifelse(state == "NJ" & wage_st < 5.05, (5.05-wage_st) / wage_st, 0))
ck_with_gap %>% filter(state == "NJ") %>% pull(gap) %>% hist(xlab = "GAP")
did_reg_iii <- lm(formula = fte2 - fte ~ gap, data = ck_with_gap, subset = wage_available)Model (iv) and (v)
\[ \Delta E_i = a' + b' X_i + c' GAP_i + \epsilon_i \]
did_reg_iv <- lm(formula = fte2 - fte ~ gap + factor(chain) + co_owned,
data = ck_with_gap, subset = wage_available)
did_reg_v <- lm(formula = fte2 - fte ~ gap + factor(chain) + co_owned + south_NJ + north_NJ + PA2 + shore,
data = ck_with_gap, subset = wage_available)まとめて出力．
library(modelsummary)
modelsummary(list(did_reg_i, did_reg_ii, did_reg_iii, did_reg_iv, did_reg_v),
gof_omit = "AIC|BIC|Log.Lik.|F|RMSE", stars = TRUE)| (1) | (2) | (3) | (4) | (5) | |
|---|---|---|---|---|---|
| (Intercept) | -1.879+ | -1.450 | -1.478* | -1.232 | -2.614* |
| (1.073) | (1.210) | (0.694) | (0.958) | (1.253) | |
| NJ | 2.277+ | 2.282+ | |||
| (1.190) | (1.197) | ||||
| factor(chain)2 | 0.235 | 0.394 | 0.638 | ||
| (1.297) | (1.287) | (1.298) | |||
| factor(chain)3 | -2.085 | -1.730 | -1.661 | ||
| (1.321) | (1.316) | (1.313) | |||
| factor(chain)4 | -0.757 | -0.171 | 0.160 | ||
| (1.491) | (1.505) | (1.520) | |||
| co_owned | 0.373 | 0.510 | 0.106 | ||
| (1.099) | (1.096) | (1.113) | |||
| gap | 17.052** | 16.363** | 16.534* | ||
| (6.090) | (6.237) | (6.994) | |||
| south_NJ | 2.200 | ||||
| (1.431) | |||||
| north_NJ | 1.876 | ||||
| (1.230) | |||||
| PA2 | 2.691 | ||||
| (1.800) | |||||
| shore | -2.355 | ||||
| (1.679) | |||||
| Num.Obs. | 351 | 351 | 351 | 351 | 351 |
| R2 | 0.010 | 0.020 | 0.022 | 0.029 | 0.046 |
| R2 Adj. | 0.008 | 0.006 | 0.019 | 0.015 | 0.021 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
注：論文中ではレポートされていないが，決定係数は非常に低い． この文脈で重要なことは DID 推定量がバイアスを持つかどうか（すなわち，興味のある因果効果が識別されているか否か）であって，アウトカムの変分が説明変数によってどれだけ説明されるかではない．
Probability value for controls の行は joint F test for exclusion of all control variables を表すとのこと． 論文でレポートされている値を全く再現できないが…
library(car)
ftest_ii <- linearHypothesis(did_reg_ii, c("factor(chain)2 = 0", "factor(chain)3 = 0",
"factor(chain)4 = 0", "co_owned = 0"))
ftest_iv <- linearHypothesis(did_reg_iv, c("factor(chain)2 = 0", "factor(chain)3 = 0",
"factor(chain)4 = 0", "co_owned = 0"))
ftest_v <- linearHypothesis(did_reg_v, c("factor(chain)2 = 0", "factor(chain)3 = 0",
"factor(chain)4 = 0", "co_owned = 0",
"south_NJ = 0", "north_NJ = 0", "PA2 = 0", "shore = 0"))
ftest_ii## Linear hypothesis test
##
## Hypothesis:
## factor(chain)2 = 0
## factor(chain)3 = 0
## factor(chain)4 = 0
## co_owned = 0
##
## Model 1: restricted model
## Model 2: fte2 - fte ~ NJ + factor(chain) + co_owned
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 349 26505
## 2 345 26236 4 268.08 0.8813 0.4752ftest_iv$`Pr(>F)`[2]; ftest_v$`Pr(>F)`[2]## [1] 0.6180603## [1] 0.3855093.2.2 Long format
Wide 形式のデータを long 形式に変換するには pivot_longer
を利用する． ちなみに以下のコードは ChatGPT が教えてくれた．
ポイント：
wage_st(starting wage) とfte(full-time equivalent employment) の2変数について，wave 1 と wave 2 それぞれが1つの行を構成するように変換する- Wave 1 の変数は
変数名1，wave 2 は変数名2で統一する wage_availableという変数名との混同を避けるため，wage_st*変数はwagest*に予め変更しておく
- Wave 1 の変数は
gap変数は before treatment の starting wage で定義するため，wide 形式の状態のうちに計算するstarts_with(文字列)で，当該文字列から始まる変数名の変数がすべて抽出されるwaveは文字列のため，Before/After ダミーを作成するためにas.numeric関数で数値型に変換する（wave= 1 ならば Before，wage= 2 ならば After）
ck_long <- ck %>%
mutate(gap = ifelse(state == "NJ" & wage_st < 5.05, (5.05-wage_st)/wage_st, 0)) %>%
rename(wagest1 = wage_st, wagest2 = wage_st2, fte1 = fte) %>%
dplyr::select(chain, co_owned, ends_with("NJ"), PA1, PA2, shore, state,
starts_with("wage"), fte1, fte2, gap) %>%
pivot_longer(cols = c(wagest1, wagest2, fte1, fte2),
names_to = c(".value", "wave"),
names_pattern = "([a-z]+)(\\d*)") %>%
mutate(after = as.numeric(wave) - 1) # = 1 if After (November)処置群ダミー \(NJ\) と処置後ダミー \(\mbox{After}\) の交差項の係数 \(\delta\) が DID 推定量．
\[ E_{it} = a + b X_{it} + c NJ_i + d \mbox{After}_t + \delta NJ_i \times \mbox{After}_t + \epsilon_i \]
did_reg_long_i <- lm(fte ~ NJ * after, data = ck_long, subset = wage_available)
# summary(did_reg_long_i)
did_reg_long_ii <- lm(fte ~ NJ * after + factor(chain) + co_owned, data = ck_long, subset = wage_available)
did_reg_long_iii <- lm(fte ~ gap * after, data = ck_long, subset = wage_available)
did_reg_long_iv <- lm(fte ~ gap * after + factor(chain) + co_owned, data = ck_long, subset = wage_available)
did_reg_long_v <- lm(fte ~ gap * after + factor(chain) + co_owned + south_NJ + north_NJ + PA2 + shore,
data = ck_long, subset = wage_available)まとめて出力．
- GAP 変数を利用する場合は「gap × after」の係数が DID 推定値
- 欠損値の関係で分析に使用される標本が wide 形式のデータを用いた場合と僅かに異なるため，推定値も若干変わっている
modelsummary(list(did_reg_long_i, did_reg_long_ii, did_reg_long_iii, did_reg_long_iv, did_reg_long_v),
gof_omit = "AIC|BIC|Log.Lik.|F|RMSE", stars = TRUE)| (1) | (2) | (3) | (4) | (5) | |
|---|---|---|---|---|---|
| (Intercept) | 23.705*** | 26.277*** | 22.892*** | 26.087*** | 26.872*** |
| (1.139) | (1.083) | (0.731) | (0.784) | (0.943) | |
| NJ | -3.039* | -2.373* | |||
| (1.260) | (1.129) | ||||
| after | -1.822 | -1.898 | -1.381 | -1.401 | -1.412 |
| (1.598) | (1.429) | (1.031) | (0.920) | (0.908) | |
| NJ × after | 2.361 | 2.466 | |||
| (1.771) | (1.583) | ||||
| factor(chain)2 | -10.607*** | -10.684*** | -10.257*** | ||
| (0.862) | (0.856) | (0.855) | |||
| factor(chain)3 | -1.473+ | -1.684+ | -1.655+ | ||
| (0.877) | (0.874) | (0.863) | |||
| factor(chain)4 | -1.364 | -1.753+ | -1.733+ | ||
| (0.975) | (0.985) | (0.985) | |||
| co_owned | -1.108 | -1.183 | -0.685 | ||
| (0.731) | (0.729) | (0.732) | |||
| gap | -19.681** | -18.909** | -15.895** | ||
| (6.431) | (5.804) | (6.094) | |||
| gap × after | 17.587+ | 17.889* | 17.782* | ||
| (9.054) | (8.076) | (7.972) | |||
| south_NJ | -3.837*** | ||||
| (0.941) | |||||
| north_NJ | -0.505 | ||||
| (0.805) | |||||
| PA2 | -1.452 | ||||
| (1.185) | |||||
| shore | -0.448 | ||||
| (1.097) | |||||
| Num.Obs. | 721 | 721 | 721 | 721 | 721 |
| R2 | 0.008 | 0.212 | 0.013 | 0.219 | 0.243 |
| R2 Adj. | 0.004 | 0.205 | 0.009 | 0.211 | 0.232 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |