研究の設計

使用するデータ：福祉国家に関するOECD諸国のデータ

• 研究の問い：社会保障費を規定する要因は何か？
• 検証仮説：高齢者の人口割合が増加することで、社会保障費は増加する。
• 利用するデータ：OECD諸国の比較福祉国家データ（Comparative Welfare State Dataset）データの入手先はこちら
• 変数の設定
  • 独立変数：高齢者の人口割合
  • 従属変数：社会保障支出の対GDP比

データフレームと変数の設定

まずはいつものように、データを読み込む。

df_cws<-read.csv("cws.csv")

次に、分析に必要と思われる変数を選択し、データフレームを構築しておこう。そして、head()コマンドによって、データフレームが意図したように構築できたかを確認する。

df_ex<- df_cws %>%　dplyr::select(year,id,socx_pub, pop,po65,postgini_2559,lisrpr_tot,
                    hunemr,leftseat,leftvot,rtseat,rtvot,vturn, pres)
head(df_ex)

##   year  id socx_pub     pop  po65 postgini_2559 lisrpr_tot hunemr leftseat
## 1 1960 AUL       NA 10275.0 874.9            NA         NA     NA     36.9
## 2 1961 AUL       NA 10508.2 894.6            NA         NA     NA     37.9
## 3 1962 AUL       NA 10700.5 913.6            NA         NA     NA     49.2
## 4 1963 AUL       NA 10906.9 933.0            NA         NA     NA     48.5
## 5 1964 AUL       NA 11121.6 948.1            NA         NA     NA     41.0
## 6 1965 AUL       NA 11340.9 966.3            NA         NA     NA     41.0
##   leftvot rtseat rtvot vturn pres
## 1    42.8   63.1  46.5 95.48    0
## 2    43.2   62.1  46.1 95.46    0
## 3    47.9   50.8  42.1 95.27    0
## 4    47.7   51.5  42.4 95.31    0
## 5    45.5   59.0  46.0 95.73    0
## 6    45.5   59.0  46.0 95.73    0

必要な変数の記述統計を確かめる。

library(psych)

df_ex %>% 
  with(describe(df_ex))

##               vars    n     mean       sd   median  trimmed     mad     min
## year             1 1210  1987.00    15.88  1987.00  1987.00   20.76 1960.00
## id*              2 1210    11.50     6.35    11.50    11.50    8.15    1.00
## socx_pub         3  651    20.97     5.14    20.71    20.97    5.68    9.90
## pop              4 1144 34795.58 54515.97 10133.05 21884.97 9536.45  314.90
## po65             5 1138  4534.67  6982.47  1434.06  2903.01 1423.69   33.90
## postgini_2559    6  123     0.28     0.05     0.28     0.28    0.05    0.18
## lisrpr_tot       7  137     9.90     3.77     9.12     9.71    4.32    3.73
## hunemr           8  670     6.74     3.63     6.37     6.49    3.20    0.00
## leftseat         9 1105    37.17    15.65    40.80    38.65   12.45    0.00
## leftvot         10 1105    37.06    13.42    40.40    38.74   10.08    0.00
## rtseat          11 1105    32.14    19.19    32.60    32.17   21.20    0.00
## rtvot           12 1105    29.82    15.96    32.60    30.69   18.24    0.00
## vturn           13 1105    78.58    12.13    80.42    79.71   12.50   42.25
## pres            14 1102     0.19     0.39     0.00     0.11    0.00    0.00
##                     max     range  skew kurtosis      se
## year            2014.00     54.00  0.00    -1.20    0.46
## id*               22.00     21.00  0.00    -1.21    0.18
## socx_pub          35.68     25.78  0.07    -0.58    0.20
## pop           311591.90 311277.00  2.91     9.10 1611.80
## po65           41394.14  41360.24  2.72     8.02  206.98
## postgini_2559      0.36      0.18 -0.15    -0.94    0.00
## lisrpr_tot        17.94     14.21  0.33    -0.98    0.32
## hunemr            21.64     21.64  0.90     1.59    0.14
## leftseat          68.60     68.60 -0.81     0.00    0.47
## leftvot           60.70     60.70 -1.13     0.86    0.40
## rtseat            85.10     85.10  0.00    -0.77    0.58
## rtvot             60.70     60.70 -0.38    -0.96    0.48
## vturn             95.77     53.52 -0.76     0.05    0.36
## pres               1.00      1.00  1.59     0.52    0.01

OLS推定の実践

ここまでの散布図も検討も踏まえて、高齢者人口割合で、社会保障支出を説明する単回帰モデルを実行する。

推定するモデルは、次のように表記される；

\(y_{welfare}=\beta_{0}+\beta_{1}x_{aged}+e_{i}.\)

これを行列形式で表記すると、以下のように書き直すことができる。

\(Y_{i}=X_{i}\beta+e.\)

この行列形式をもとにすると、OLSは、

\(\sum^n_{i=1}(Y_{i}-X_{i}\beta)^2\)

を最小化するものとして定義できる。これを解くと\(\beta_{1}\)と\(\beta_{0}\)は、

\(\hat{\beta}_{1}=Cor(Y,X)\frac{sd(Y)}{sd(X)},\)

\(\hat{\beta}_{0}=\hat{Y}-\hat{\beta}_{1}\hat{X},\)

として計算できる。この\(\hat{\beta}_{1}\)と\(\hat{\beta}_{0}\)を計算するためにRでは、次のコマンドが利用されるが・・・

result_ols <- lm(socx_pub ~ peraged, data = na.omit(df_ex)) 
summary(result_ols)

## 
## Call:
## lm(formula = socx_pub ~ peraged, data = na.omit(df_ex))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9707 -3.5175 -0.3949  3.0415  9.0872 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.624      2.592   0.627    0.532    
## peraged      141.266     17.675   7.992 3.97e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.945 on 91 degrees of freedom
## Multiple R-squared:  0.4124, Adjusted R-squared:  0.406 
## F-statistic: 63.88 on 1 and 91 DF,  p-value: 3.97e-12

上記の推定結果を、先ほどの散布図に反映しておく。

# 回帰モデルを推定：定数項と係数の値を抽出
reg_mod <- lm(socx_pub ~ peraged, data = na.omit(df_ex))
int_val <- coef(reg_mod)[1]
slope_val <- coef(reg_mod)[2]

# 注釈テキストを生成 (R の expression 形式に合わせる)
# round() により小数点以下3桁に丸めた値を表示するように指示
label_expr <- paste0("y[welfare] == ", round(int_val, 3), " + ", round(slope_val, 3), "*x[aged] + epsilon[t]")

# x_pos, y_pos は注釈の表示位置
x_pos <- 0.9 * max(df_ex$peraged, na.rm = TRUE)
y_pos <- 0.9 * max(df_ex$socx_pub, na.rm = TRUE)

# 散布図と回帰直線の作成
g5 <- df_ex %>%
  ggplot(aes(x = peraged, y = socx_pub)) +
    geom_point(size = 3, alpha = 0.2, colour = "black") +
    geom_point(size = 2, alpha = 0.2, colour = "red") +
    geom_smooth(method = "lm", se = TRUE) +  
    xlab("65%人口割合") +
    ylab("社会保障費の対GDP比") +
    ggtitle("高齢者人口と社会保障費の関係") +
    annotate(
      "text",
      x = x_pos,
      y = y_pos,
      label = label_expr, #先ほどの数式を書き入れると指示
      parse = TRUE,  
      size = 3,
      color = "black"
    ) +
    theme_bw()

plot(g5)

この結果と整合的なものが、上記の推定モデルの背景となる数式をもとに計算した場合にも得られるか確認しておこう。

results1 <- df_ex %>%
  # 欠損値を除いたデータ処理
  na.omit() %>%
  { 
    # 手計算による傾きと切片の計算
    beta1_manual <- cor(.$socx_pub, .$peraged, use = "complete.obs") * 
                     (sd(.$socx_pub, na.rm = TRUE) / sd(.$peraged, na.rm = TRUE))
    beta0_manual <- mean(.$socx_pub, na.rm = TRUE) - beta1_manual * mean(.$peraged, na.rm = TRUE)
    manual_vec <- c(beta0_manual, beta1_manual)
    
    # lm() を用いた線形回帰モデルの推定結果を再度計算
    lm_vec <- lm(socx_pub ~ peraged, data = .) %>% coef()
    
    # 手計算結果と lm() の結果を行ごとに結合
    rbind("manual" = manual_vec, "lm(df_ex$socx_pub ~ df_ex$peraged)" = lm_vec)
  } %>%
  { 
    #  列名の設定
    colnames(.) <- c("切片（定数項）", "傾き（係数）")
    .
  }

results1

##                                    切片（定数項） 傾き（係数）
## manual                                   1.624476     141.2655
## lm(df_ex$socx_pub ~ df_ex$peraged)       1.624476     141.2655

では、分析結果を解釈していく。それぞれの項目が何を意味するか、下図を参考に検討する。

• 赤四角：係数
• 青四角：標準誤差
• 緑四角：\(t\)値 (\(t\)検定の\(t\)統計量)
• 黄四角：\(p\)値 (\(t\)検定の有意確率値)

推定量を解釈する上でのポイントは以下の通りである。

• 各変数についての\(t\)検定の結果はどうか？

• 各係数の符号は、仮説にける符号条件を満たしているか？

• 各係数の値の大きさは、実質的に意味のあるものか？

次に解釈する項目は、モデルの妥当性についてである。

• グレー四角：決定係数（\(R^2\)）

• ピンク四角：\(F\)検定の結果

中心化(centering)

分析結果の意味について考えてみよう。まずは定数項に注目してみほしい。

定数項の意味：独立変数が0の時に予測されるYの値

→独立変数が0の時に意味があるのか？高齢者人口割合が0になる、という状態が存在するのか？

【中心化(centering)】：データの原点から平均値を差し引いた値をもとに各パラメターを算出

\(Y-\bar{Y}=\beta_{1}(X-\bar{X})+e.\)

results2 <- df_ex %>%
  na.omit() %>%
  mutate(
    centered_socx_pub = socx_pub - mean(socx_pub, na.rm = TRUE),# 従属変数の中心化
    centered_peraged  =peraged  - mean(peraged,  na.rm = TRUE)# 独立変数の中心化
  ) %>%
  { 
    orig_cor <- cor(.$socx_pub, .$peraged, use = "complete.obs")
    centered_cor <- cor(.$centered_socx_pub, .$centered_peraged, use = "complete.obs")
    lm_model <- lm(centered_socx_pub ~ centered_peraged, data = .)
    slope <- coef(lm_model)[2]
    intercept <- coef(lm_model)[1]
    tibble(
      統計量 = c("相関係数(社保, 高齢)", "相関係数(中心化・社保,中心化・高齢)", "中心化後の係数", "中心化後の定数項"),
      推定結果 = round(c(orig_cor, centered_cor, slope, intercept), 3)
    )# 元の相関と、中心化後の相関および回帰モデルの係数をまとめる
  }

results2

## # A tibble: 4 × 2
##   統計量                              推定結果
##   <chr>                                  <dbl>
## 1 相関係数(社保, 高齢)                   0.642
## 2 相関係数(中心化・社保,中心化・高齢)    0.642
## 3 中心化後の係数                       141.   
## 4 中心化後の定数項                       0

g6 <- df_ex %>%
  # 中心化：socx_pub と peraged の各値から平均を引く、中心化済みの値を利用
  mutate(
    centered_socx_pub = socx_pub - mean(socx_pub, na.rm = TRUE),
    centered_peraged  = peraged  - mean(peraged,  na.rm = TRUE)
  ) %>%
  ggplot(aes(x = centered_peraged, y = centered_socx_pub)) +
    geom_point(size = 3, alpha = 0.2, colour = "black") +
    geom_point(size = 2, alpha = 0.2, colour = "red") +
    geom_smooth(method = "lm") +
    labs(
      x = "中心化された65%人口割合",
      y = "中心化された社会保障費の対GDP比",
      title = "中心化後の高齢者人口と社会保障費の関係"
    ) +
    theme_bw()

g6

上記の結果において、係数（傾き）は変わらず、定数項（切片）の値のみが変化していることに注目しよう。

正規化(normalize)

【正規化(normalize)】:特徴量の値の範囲を一定に抑える処理

results3 <- df_ex %>%
  # 正規化：平均を引き、標準偏差で割る
  mutate(
    norm_socx_pub = (socx_pub - mean(socx_pub, na.rm = TRUE)) / sd(socx_pub, na.rm = TRUE),# 従属変数の正規化
    norm_peraged  = (peraged  - mean(peraged,  na.rm = TRUE)) / sd(peraged, na.rm = TRUE)# 独立変数の正規化
  ) %>%
  { 
    orig_cor <- cor(.$socx_pub, .$peraged, use = "complete.obs")
    norm_cor <- cor(.$norm_socx_pub, .$norm_peraged, use = "complete.obs")
    lm_model <- lm(norm_socx_pub ~ norm_peraged, data = .)
    slope <- coef(lm_model)[2]
    intercept <- coef(lm_model)[1]
    tibble(
      統計量 = c("相関係数(社保, 高齢)", "相関係数(正規化・社保,正規化・高齢)", "正規化後の係数", "正規化後の定数項"),
      推定結果 = round(c(orig_cor, norm_cor, slope, intercept), 3)
    )# 元の相関と、正規化後の相関および回帰モデルの係数をまとめる
  }

results3

## # A tibble: 4 × 2
##   統計量                              推定結果
##   <chr>                                  <dbl>
## 1 相関係数(社保, 高齢)                   0.585
## 2 相関係数(正規化・社保,正規化・高齢)    0.585
## 3 正規化後の係数                         0.721
## 4 正規化後の定数項                      -0.296

正規化後の係数・定数項が変化していることがわかるはず。

残差(residuals)の検討

次に算出された残差(residuals)について検討する。以下の残差と予測値の計算を示す。

library(gridExtra) #複数の図を1枚に図にアレンジするためのパッケージ

#------- データ前処理 
# df_exから欠損値を除去し、単回帰分析を実施して予測値と残差を計算
df_ex_clean <- df_ex %>%
  na.omit() %>%
  mutate(
    pred  = predict(lm(socx_pub ~ peraged, data = .)),
    resid = resid(lm(socx_pub ~ peraged, data = .))
  )

#------- 第1パネル：実測値＆回帰直線に対する残差プロット 
p1 <- ggplot(df_ex_clean, aes(x = peraged, y = socx_pub)) +
  geom_point(size = 3, alpha = 0.2, shape = 21, fill = "lightblue", color = "black") +
  geom_point(size = 2, alpha = 0.2, color = "red") +
  # 回帰直線を追加（信頼区間は非表示）
  geom_smooth(method = "lm", se = FALSE, color = "black", size = 1.2) +
  # 各観測点から予測値 (pred) への垂直赤線を追加
  geom_segment(aes(x = peraged, xend = peraged, y = socx_pub, yend = pred),
               color = "red", size = 1.2) +
  labs(x = "X: 高齢者人口割合", y = "Y: 社会保障費の対GDP比",
       title = "回帰直線に対する残差") +
  theme_bw()

#------- 第2パネル：残差プロット 
p2 <- ggplot(df_ex_clean, aes(x = peraged, y = resid)) +
  geom_point(size = 3, alpha = 0.2, shape = 21, fill = "lightblue", color = "black") +
  geom_point(size = 2, alpha = 0.2, color = "red") +
  # y=0 の水平線（参考線）を点線で描画
  geom_hline(yintercept = 0, linetype = "dashed", size = 1.2) +
  # 各観測点の残差から0へ向かう垂直赤線を描画
  geom_segment(aes(x = peraged, xend = peraged, y = resid, yend = 0),
               color = "red", size = 1.2) +
  labs(x = "X: 高齢者人口割合", y = "残差", title = "残差 vs X") +
  theme_bw()

# --- プロットの配置 ---
# grid.arrange() を用いて2つのプロットを横並び (1行2列) に配置
grid.arrange(p1, p2, ncol = 2)

このように見ると、残差がそれぞれのXに対して不均一に分散しているように見える・・・これは何を意味しているのか、どのような問題があるのかについては、のちほど詳しく検討しよう。

Goodness of fitを考える

「モデルの当てはまり(Goodness of fit)」についても検討しよう。

• モデルの当てはまりを測る指標：\(R^2\)(アール・スクエア)

\(R^2=\frac{\sum^n_i=1(\hat{Y}-\bar{Y})^2}{\sum^n_{i=1}(Y_{i}-\bar{Y})^2}\)

\(R^2=1-\frac{\sum^n_{i=1}\hat{\mu}_i^2}{\sum^n_{i=1}(Y_{i}-\bar{Y})^2}\)

• \(\sum^n_{i=1}(Y_{i}-\bar{Y})^2\)：\(Y_i\)の全変動

• \(\sum^n_{i=1}(\hat{Y}_{i}-\bar{Y})^2\):\(\hat{Y}_i\)で説明される部分。回帰直線によって説明される部分。

• \(\sum^n_{i=1}\hat{\mu}_i^2\):\(\hat{Y}_i\)で説明されない部分。回帰直線によって説明されない部分

• 更に確認しておこう。\(y\)と\(x\)の相関係数\(r\)は・・・

df_ex %>%
  na.omit() %>%
    with(cor(peraged,socx_pub))

## [1] 0.6422189

相関係数\(r\)の２乗値は・・・

df_ex %>%
  na.omit() %>%
    with((cor(peraged,socx_pub))^2)

## [1] 0.4124451

相関係数\(r\)の2乗値は\(R^2\)に一致する。

• \(R^2\)の評価：一概にどの値であればいいとは言い難い。いま単回帰分析で、「0.412」という値であれば、十分なモデルの当てはまりであると判断できる

①線形性：推定量\(\beta_1\)は\(Y_1,...,Y_n\)に対する線形関数である

\(\hat{\beta}=\sum^n_{i=1}a_iY_i.\)

※仮定が侵害される場合：従属変数の形状が離散型確率変数であるために。独立変数と従属変数の間にU字・逆U字などの関係がある場合

②誤差項の期待値：誤差項の期待値はゼロである。全ての\(X\)に対して、誤差項の平均は0になる。

\(E[(e_i|X_{i})]=0\) （→\(E[Y_i|X_{i}]=\beta_{0}+\beta_{1}X_{i}.\)）

③誤差間の相関：誤差どうしの相関はゼロである

\(Cov(e_i,e_j)=0,\) \(\forall i\neq j.\)

※仮定が侵害される場合：時系列データや空間データで、系列相関（serial correlation）がある場合。今期(\(t\))の誤差が、前期(\(t-1\))の誤差に影響を受けている場合。

⇒推定量の有効性の侵害

④誤差の分散：誤差の分散は均一(homoskedasticity)である

\(Var(e_i|X_{i})=\sigma^2,\) \(\forall i.\)

※仮定が侵害される場合：従属変数の形状が離散型確率変数であるために、誤差の分散が不均一（heteroskedasticity）である場合。

⇒推定量の有効性の侵害

⑤誤差項と独立変数の相関：誤差項と独立変数の相関はゼロである(=誤差項は独立変数に対して外生的・exoneousである)

\(Cov(e_i,X_{i})=0\)

※仮定が侵害される場合：モデルの特定化(specification)で変数の除去があるために、誤差項と独立変数が相関を持つ場合（内生性・endogeneity）。

⇒推定量の不偏性の侵害

⑥誤差項の正規性：誤差は正規分布に従う（分布解析（\(t\)定・\(F\)検定）を正しく実行するための追加仮定）

\(e_i \sim i.i.d.,\quad N(0,\sigma^2)\)

※仮定が侵害される場合：\(t\)検定や\(F\)検定における誤謬が生じる可能性

OLS推定量がBLUEであるために

①から⑤の仮定が満たされるとき、OLS推定量は…

Best Linear Unbiased Estimator (BLUE): 最良線形不偏推定量

しかし、①から⑤が満たされる、OLS推定量がBLUEとなる場合はかなり限られているのではないかと推察される。どのようにすれば、OLS推定量がBLUEであるかを確かめられるだろうか。そのために、どのような検定が必要だろうか。

「①線形性」をどのように検定するか？

線形の関係を仮定して行った推定だったが、ひょっとしたら両者の関係は線形のモデルに当てはめるよりも、非線形(non-linear)のモデルに当てはめる方がよいのかもしれない。両者の関係が、U字や逆U字となる場合もあるだろう。それを検討するために、Ramsey-RESET検定を行う。2乗項、3乗項を加えたときに、それらの変数が統計的に有意である可能性はどうか検定する。

#------ Ramsey-RESET検定のためのパッケージ
library(lmtest)

#------ 単回帰モデル result_ols を対象に Ramsey RESET test を実施。2乗項、3乗項を追加したときに追加した項が有意かを確認。追加した項が有意ならば、線形でない可能性
reset_result <- resettest(result_ols, power = 2:3)
print(reset_result)

## 
##  RESET test
## 
## data:  result_ols
## RESET = 4.0635, df1 = 2, df2 = 89, p-value = 0.02048

検定結果から、帰無仮説は棄却される。すなわち、2乗項や3乗項を加えた場合に、それらの項が統計的に有意であるの可能性を否定できないという結果。実際にそうなのか確かめてみよう。

以下で2乗項I(peraged^2)を含んだモデルを推定する。^がべき乗の処理を意味する。

result_sq <- lm(socx_pub ~ peraged + I(peraged^2), data = na.omit(df_ex))
summary(result_sq)

## 
## Call:
## lm(formula = socx_pub ~ peraged + I(peraged^2), data = na.omit(df_ex))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2288 -3.0375 -0.1404  2.6462  9.4724 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -34.04      12.69  -2.682 0.008698 ** 
## peraged        643.22     175.95   3.656 0.000431 ***
## I(peraged^2) -1721.59     600.63  -2.866 0.005170 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.797 on 90 degrees of freedom
## Multiple R-squared:  0.4616, Adjusted R-squared:  0.4496 
## F-statistic: 38.58 on 2 and 90 DF,  p-value: 7.942e-13

2乗項が統計的に有意に0ではないとは、どういうことを意味しているのだろうか。図示によって確かめてみよう。

df_clean <- df_ex %>% 
  na.omit()

#------ 散布図と二次回帰曲線（U字型の関係）のプロット
g8 <- df_clean %>%
  ggplot(aes(x = peraged, y = socx_pub)) +
    geom_point(color = "blue", alpha = 0.5) +
    geom_smooth(method = "lm", formula = y ~ x + I(x^2), color = "red", se = TRUE) + #2次式に当てはめることを指示
    labs(
      x = "65%人口割合",
      y = "社会保障費の対GDP比",
      title = "両者のlinkは線形？非線形？"
    ) +
    theme_bw()

plot(g8)

この図示から明らかなように、高齢者人口割合が20%に近づく（5人に1人！）になってくると、社会保障費の対GDP比割合は緩やかになり、やや漸減に転じる可能性。1ケタ台から10％台にかけては増えていくのに対して、徐々にカーブしてくことを考えると、線形よりも非線形のモデルの方が望ましい可能性がある。線形性のチェックがいかに重要かわかってもらえるだろう。

「②誤差項の期待値はゼロ」をどのように検定するか？

#------- モデルの残差を抽出
resid_ols <- resid(result_ols)

#------- 残差の平均が0であるかを1標本t検定を用いて検定
t_test_res <- t.test(resid_ols, mu = 0)
print(t_test_res)

## 
##  One Sample t-test
## 
## data:  resid_ols
## t = -1.4896e-15, df = 92, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.8080578  0.8080578
## sample estimates:
##     mean of x 
## -6.060766e-16

「帰無仮説：true meanが0である」=誤差項の期待値（平均）は0であるを棄却できない。よって、この推定結果からは、誤差項の期待値は0であるとの仮定が満たされることが示唆される

「③誤差項間の相関はゼロ」をどのように検定するか？

誤差項間の相関、系列相関の有無を確かめるために、時系列データに対するDurbin-Watson testとより一般的なBreusch-Godfrey testを行って、系列相関の有無を確かめる。

• DW検定の帰無仮説：自己相関は0である

• BG検定の帰無仮説：（指定したラグまでの）自己相関は0である

#------ パッケージの読み込み
library(lmtest)     

#------　Durbin-Watson検定
dwtest_result <- dwtest(result_ols)
print(dwtest_result)

## 
##  Durbin-Watson test
## 
## data:  result_ols
## DW = 0.55457, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

#------、Breusch-Godfrey 検定
bgtest_result <- bgtest(result_ols)
print(bgtest_result)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  result_ols
## LM test = 47.556, df = 1, p-value = 5.345e-12

DW検定とBG検定で異なる結果が得られている。DW検定の結果にもとづくと、帰無仮説は棄却されず、系列相関を否定できる。それに対して、BG検定では、帰無仮説を棄却できず、系列相関を否定できない。推定量は有効であるとはいえず、誤差の再考が必要と考えられる。

「④誤差項の均一分散」をどのように検定するか？

誤差項の不均一分散(heteroskedasiticity)が生じていないかを確認するために、Breusch-Pagan Cook-Weisberg検定を行う。

• BPCW検定の帰無仮説: 誤差項の分散は均一である。

bptest_result <- bptest(result_ols)
print(bptest_result)

## 
##  studentized Breusch-Pagan test
## 
## data:  result_ols
## BP = 0.38949, df = 1, p-value = 0.5326

推定の結果、誤差項の不均一分散は生じていないことが明らかである。もし誤差項の不均一分散が認められるならば、誤差項の不均一分散に対処し、誤差の過小評価を避けるためにホワイト修正済み頑健誤差（White-adjusted robust error）を計算し報告する。

library(sandwich)

robust_coefs <- coeftest(result_ols, vcov = vcovHC(result_ols, type = "HC1"))
print(robust_coefs)

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   1.6245     2.6276  0.6182     0.538    
## peraged     141.2655    18.1232  7.7947 1.017e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

「⑤誤差項と独立変数間の無相関」をどのように検定するか？

誤差項のと独立変数間の無相関を検定することは、モデルに除外変数(omitted varibles)が存在しないかを検定することである。除外変数無視バイアスにより、推定量は不偏性を持たないことから、除外変数バイアスのために、再びRamsey-RESET検定が必要に。

• Ramsey-RESET検定の帰無仮説: モデルに除外変数はない。

#------ Ramsey-RESET検定のためのパッケージ
library(lmtest)

reset_result_simple <- resettest(result_ols)
print(reset_result_simple)

## 
##  RESET test
## 
## data:  result_ols
## RESET = 4.0635, df1 = 2, df2 = 89, p-value = 0.02048

Ramsey-RESET検定の結果からは、モデルに除外変数がある可能性を棄却できない。現段階で、これは当然だろう。1つの要因で、社会保障費の変動を説明できないのは自明である。従って、複数の変数を含めた重回帰分析（multiple regression analysis）を行っていくことになる。

「⑥誤差項の正規性」をどのように検定するか？

まずはqqプロットによって、視覚的に確認を試みる。

library(car)
qqPlot(result_ols,             
       main = "正規性についてのQQプロット",
       id = FALSE)     

直線状に各プロットがより多く一致しているほど、正規性が高いことを意味する。ここでの判断から、かなりの程度誤差の正規性が保たれていることが見て取れるが、やや外れる観察もあるようで判別が難しい。

続いて、正規性検定を行う。まずシャピロ・ウィルク検定(Shapiro-Wilk test)から。

shapiro.test(resid_ols)

## 
##  Shapiro-Wilk normality test
## 
## data:  resid_ols
## W = 0.96412, p-value = 0.01172

帰無仮説「正規分布に従う」が5％水準で棄却されるので、誤差の正規性は否定される。

続いて、Jarque Bera testの結果からは、帰無仮説が棄却されないことが示されている。これらの結果を総合すると、極めて境界的なケースであることが示唆される。

library(tseries)
jarque.bera.test(resid_ols)

## 
##  Jarque Bera Test
## 
## data:  resid_ols
## X-squared = 4.0549, df = 2, p-value = 0.1317

よって総合的に、この推定結果をジャッジしてみよう。回帰診断（regression diagnostics）を、gvlma()関数を使って行う。

library(gvlma)
reg_diag <- gvlma(result_ols)
summary(reg_diag)

## 
## Call:
## lm(formula = socx_pub ~ peraged, data = na.omit(df_ex))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9707 -3.5175 -0.3949  3.0415  9.0872 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.624      2.592   0.627    0.532    
## peraged      141.266     17.675   7.992 3.97e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.945 on 91 degrees of freedom
## Multiple R-squared:  0.4124, Adjusted R-squared:  0.406 
## F-statistic: 63.88 on 1 and 91 DF,  p-value: 3.97e-12
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = result_ols) 
## 
##                      Value  p-value                   Decision
## Global Stat        12.6870 0.012911 Assumptions NOT satisfied!
## Skewness            1.7258 0.188950    Assumptions acceptable.
## Kurtosis            2.3291 0.126974    Assumptions acceptable.
## Link Function       7.7794 0.005285 Assumptions NOT satisfied!
## Heteroscedasticity  0.8527 0.355781    Assumptions acceptable.

これらの分析結果からすると、特に\(\texttt{Gloval-test}\)の結果から、本推定結果はガウス・マルコフ仮定を全体として満たさない可能性が高いと判断される。とりわけ問題となるのは、「線形性」のでようである。既述の通り、問題になるのは逆Ｕ字性である。

選択バイアスの問題

調査時点の母集団の中でデータの収集が可能、回答が可能なサンプルは常に限られる（例：各国の中にはデータが取れない国もある／高齢者がアンケートに答えてくれない／不真面目な人はきちんと答えてくれないなどたくさん！）

⇒データの収集には観察されたサンプリングに対して、母集団を反映できない選択の偏り（selection bias）が不可避。また、観察できない要素≒データから漏れる要素も多く存在する

パネル・データの特性を考える①：固定効果モデル

パネル・データを使うことの利点を固定効果モデル（Fixed Effect model: FE model）の説明をもとに考えていこう。

通常の最小二乗法(OLS)のモデルは下記の通り；

\[y_{it}=\beta_{0}+\alpha z_{i}+\beta x_{it}+\upsilon_{it}.\]

\(y_{it}\)はサンプル\(i\)の\(t\)期のの生活費、\(z_{i}\)は職種などのいったん決まったら時間的に変化しない変数、\(x_{it}\)は\(t\)期の給料のように時期ごとに変化する変数、\(\upsilon_{it}\)は観察できない誤差であるとする。この誤差の中には、「真の誤差」に加えて、個人\(i\)が個人\(i\)であることによる属性・要因（個性、性質、性向、特性）などが含まれる。その特性を\(\mu_{i}\)、真の誤差を\(\epsilon_{it}\)とすると、上式は以下のように書き直せる；

\[y_{it}=\beta_{0}+\alpha z_{i}+\beta x_{it}+\mu_{i}+\epsilon_{it}.\]

\(\upsilon_{it}\)は\(\mu_{i}+\epsilon_{it}\)として書き直されていることになるが、ここで\(\mu_{i}\)は各国の歴史的文脈のようなものでデータ化が困難であったり、アンケートでも質問が困難であったりする観察不可能な個体固有成分の異質性（heterogeneity）である。ここで説明変数と\(\mu_{i}\)の間に相関があるならば、\(\beta\)の推定値はバイアスをもつ。

ここで、個人間（between）で要因が変動すると考えると、変化する変数の個人の平均をもとに、以下の式を設定できる；

\[\bar{y_{i}}=\alpha z_{i}+\beta\bar{x}_{i}+\mu_{i}+\epsilon_{i}\]

バイアスとなっている\(\mu_{i}\)を除去するために、上2本の式の差をとると；

\[(y_{it}-\bar{y}_{i})=\beta(x_{it}-\bar{x}_{i})+(\epsilon_{it}-\bar{\epsilon}_{i})\]

が導かれ、観測できな要因である\(\mu_{i}\)を除いた推定量を導き、個人内での独立変数の変動によって、従属変数の変動を説明することが可能になる。パネル推定の固定効果推定量が一致推定量ならば、

• 推定される\(\beta\)は「同じ個体内で\(x\)が変わった時の平均的な影響」となる

• 時点で不変(time-invariant)な\(z_i\)も消えてしまうため、係数\(\alpha\)は推定できない。

• 説明変数と\(\mu_i\)の相関が許容されることで内生性バイアスを回避可能

といった利点がある。一方で、時変しない変数（執政制度、選挙制度、民主主義の程度など）の効果は測定できないため注意を要する。

パネル・データの特性を考える②：ランダム効果モデル（Random Effect model: RE model）

固定効果モデルにおいては「\(\mu_i\)」をダミー変数で消したが、ランダム効果モデルは\(\mu_i\)を確率的にランダムな係数とみなす；

\[ y_{it} = \beta_0 + \alpha\,z_i + \beta\,x_{it} + \upsilon_i + \epsilon_{it},\quad \upsilon_i\sim \mathrm{i.i.d.}(0,\sigma_\upsilon^2),\ \epsilon_{it}\sim \mathrm{i.i.d.}(0,\sigma_\epsilon^2).\]

上式においては、個体固有成分\(\upsilon_i\)はすべての説明変数と無相関であり、\(\epsilon_{it}\)は個体・時点ともに独立であることが仮定されることになる。

ランダム効果推定量の特徴として、以下の点が挙げられる。

• 時点で不変な変数\(z_i\)の係数\(\alpha\)も推定可能

• 無相関の仮定が満たされるときにFEより分散が小さく有効性が高い

• 無相関の仮定が侵害されているならば（\(\upsilon_i\)と説明変数が相関）しているならば不偏推定量にならない

\(\mu_i\)または\(\upsilon_i\)に関係して、Wu-Hausman 検定によりFEとREの推定値の差に関する統計的有意性の判別が必要に。

パネル・データ分析の設定

ここまでの分析でも、ある程度リサーチ・デザインは意識してきたが、以下では改めてリサーチ・デザインを明示する。

研究の設計：各国の社会保障費は、何によって規定されているのか。

メインの仮説の設定：

• メインの仮説1：左派政党の議席率が高いほど、社会保障費が減る

対立仮説の検討：

• 対立仮説1：高齢者人口費が高まるほど、社会保障費が増える

• 対立仮説2：比例代表制の場合ほど、社会保障費が増える

• 対立仮説3：投票率が高いほど、社会保障費が増える

利用するデータ：Comparative Welfare State dataset

変数の設定：

• 従属変数：社会保障費の対GDP比

• カギとなる独立変数(key independent variables)：左派議席比率

• 制御変数(control variables)：高齢者人口比、 選挙制度ダミー、失業率、インフレ率、投票率

パネル・データの推定

パネル・データのプロット。社会保障費の対GDP比に対する仮説に対応した各要因について、国別にパネルを折りたたみながら図示してみよう。

まずはメインの独立変数である左派議席率について検討する。左派議席率が増えることによって、社会保障支出が増えている国はポーランド、スウェーデンぐらいで、左派議席率の増加に伴って、社会保障支出が減る国、あるいは影響が認められないと思われる国が多いようである。

#------ データ読み込みと整形
df_cws <- read.csv("cws.csv", fileEncoding = "SJIS") %>%
  transmute(
    id,
    year,
    socx_pub,
    popold = po65 / pop,
    leftseat,
    vturn,
    singmemd,
    hunemr,
    cpi
  )

#------ 国別のパネル図
df_cws %>%
  ggplot(aes(x = leftseat, y = socx_pub)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE, color = "steelblue", fill = "lightblue", alpha = 0.3) +
  facet_wrap(~ id, scales = "free") +
  theme_bw() +
  labs(
    title = "国別：左派議席比率と社会保障支出",
    x     = "左派議席比率",
    y     = "社会保障費対GDP比"
  )

高齢者人口割合については以下の通り。多くの国で、高齢者人口割合の増加が社会保障支出を押し上げていることが分かる。I-2内の各種散布図を国別に仕分けたものだが、I-2内の各散布図が各国の結果を集積したものであるとわかるように、各国で明確な正の結びつきが見てとれる。

#------　国別のプロット：高齢者人口割合

df_cws %>%
  ggplot(aes(x = popold, y = socx_pub)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE, color = "steelblue", fill = "lightblue", alpha = 0.3) +
  facet_wrap(~ id, scales = "free") +
  theme_bw() +
  labs(
    title = "国別：65歳以上人口比と社会保障支出",
    x     = "65歳以上人口比",
    y     = "社会保障費対GDP比"
  )

投票率については、左派議席率の結果とよく似て、明確な傾向は確かめられない。

#------ 国別のプロット：投票率

df_cws %>%
  ggplot(aes(x = vturn, y = socx_pub)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE, color = "steelblue", fill = "lightblue", alpha = 0.3) +
  facet_wrap(~ id, scales = "free") +
  theme_bw() +
  labs(
    title = "国別：左派議席比率と社会保障支出",
    x     = "投票率",
    y     = "投票率"
  )

続いて、既述の固定効果モデルとランダム効果モデルの推定を順に行い、どちらが一致推定量(consistent estimator)となるかを評価する。手順は、

①固定効果モデル推定の実施

②ランダム効果モデル推定の実施

③どちらが一致推定量かについてのハウスマン検定(Wu-Hausman test)の実施

\(H_0\): 固定効果モデルとの差が偶然である（ランダム効果モデルが、一致推定量として望ましい

\(H_1\): 固定効果モデルとの差が偶然ではない（固定効果モデルが、一致推定量として望ましい）

となる。上記の推定から検定までを、効率化されたコードをもとに実行する。まずはなるべく、以下のコードを見ずに、パネル・データの推定モデルはどのようにすれば実施できるか、どのようにすればハウスマン検定ができるか、それらのコードを効率化させるにはどうすればよいか自身で考えて、多くのエラーを出しながら頑張ってみてほしい。

library(plm) #パネル・データ分析用のパッケージ
library(lmtest) #Wu-Hausman testのためのパッケージ

hausman_result <- 
  read.csv("cws.csv", fileEncoding = "SJIS") %>%
  # データ整形
  transmute(
    id,
    year,
    socx_pub,
    popold = po65 / pop,
    leftseat,
    vturn,
    prsys = case_when(
      singmemd == 0        ~ 1,
      singmemd %in% c(1,2) ~ 0,
      TRUE                 ~ NA_real_
    ),
    hunemr,
    cpi
  ) %>%
  { 
    # 固定効果モデルの推定
    model_FE <- plm(socx_pub ~ popold + leftseat + vturn + prsys + hunemr + cpi,
                    data  = .,
                    index = c("id", "year"),
                    model = "within")
    # ランダム効果モデル
    model_RE <- plm(socx_pub ~ popold + leftseat + vturn + prsys + hunemr + cpi,
                    data  = .,
                    index = c("id", "year"),
                    model = "random")
    # モデル要約を出力
    print(summary(model_FE))
    print(summary(model_RE))
    # Wu-Hausman検定
    phtest(model_FE, model_RE)
  }

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = socx_pub ~ popold + leftseat + vturn + prsys + 
##     hunemr + cpi, data = ., model = "within", index = c("id", 
##     "year"))
## 
## Unbalanced Panel: n = 21, T = 11-31, N = 537
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -4.171415 -0.937344 -0.032683  0.835710  6.086397 
## 
## Coefficients:
##            Estimate Std. Error t-value  Pr(>|t|)    
## popold   52.5319548  6.1969363  8.4771 2.499e-16 ***
## leftseat -0.0107755  0.0099337 -1.0847    0.2785    
## vturn     0.0300136  0.0221282  1.3564    0.1756    
## prsys    -0.4616750  0.4052846 -1.1391    0.2552    
## hunemr    0.4261277  0.0325927 13.0743 < 2.2e-16 ***
## cpi       0.0529909  0.0054992  9.6362 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    2646.2
## Residual Sum of Squares: 1232.1
## R-Squared:      0.53437
## Adj. R-Squared: 0.51064
## F-statistic: 97.5497 on 6 and 510 DF, p-value: < 2.22e-16
## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = socx_pub ~ popold + leftseat + vturn + prsys + 
##     hunemr + cpi, data = ., model = "random", index = c("id", 
##     "year"))
## 
## Unbalanced Panel: n = 21, T = 11-31, N = 537
## 
## Effects:
##                  var std.dev share
## idiosyncratic  2.416   1.554 0.136
## individual    15.292   3.911 0.864
## theta:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8810  0.9237  0.9237  0.9221  0.9264  0.9288 
## 
## Residuals:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -4.4805 -1.0368 -0.1162 -0.0049  0.9080  6.4994 
## 
## Coefficients:
##               Estimate Std. Error z-value Pr(>|z|)    
## (Intercept)  3.5801631  2.3291024  1.5371  0.12426    
## popold      55.0774946  6.0982487  9.0317  < 2e-16 ***
## leftseat    -0.0075766  0.0097797 -0.7747  0.43850    
## vturn        0.0396441  0.0213939  1.8531  0.06387 .  
## prsys       -0.3019353  0.3987817 -0.7571  0.44896    
## hunemr       0.4221173  0.0324805 12.9960  < 2e-16 ***
## cpi          0.0526473  0.0054927  9.5849  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    2760.2
## Residual Sum of Squares: 1287.4
## R-Squared:      0.53361
## Adj. R-Squared: 0.52833
## Chisq: 588.888 on 6 DF, p-value: < 2.22e-16

print(hausman_result)

## 
##  Hausman Test
## 
## data:  socx_pub ~ popold + leftseat + vturn + prsys + hunemr + cpi
## chisq = 9.1342, df = 6, p-value = 0.1662
## alternative hypothesis: one model is inconsistent

本推定モデルの場合、帰無仮説は棄却されないことがわかる。ランダム効果推定量を採用することになる。

DIDとは何か？

パネル・データを用いて行う、有力な自然実験（natural experiment）の一つに「差分の差法(Difference in Difference: DID)」がある。

下記図のように、何らかの政策介入、構造変化に実際に効果があったのか、そのもとで関心のある変数に効果差が生じていたのかを知りたいケースは、公共政策の評価において多々あるだろう。

政策介入が行われるグループとそうでないグループが分かれ、時間経過のもとに政策介入によって変化が生じたのであれば、それを自然に生じた実験の結果とみなせる。DIDでは、

• 政策介入が起こった群と起こっていない群の時間的経過

• 構造変化が起こった群と起こっていない群の時間的経過

これらについて、処置群と対照群との変化の差をとって、さらにその差をとることにより、共通トレンド以外の要因を取り除くことにより、政策介入の因果効果の特定を目指す。

DIDの実践

先ほどから取り上げている、AS (2004: 505)には、構造変化(structural break)について、次のような記述がある。

さらに、AS (2004: FN20, 505)では、構造変化としてこれらの時期が想定されているという。

上記のもとに、以下の研究設計を検討し、DIDとして実践してみたい（政策介入とは異なるので、自然実験のデザインとしては疑問もあるが…）。

研究の設計：1980年代に民主主義体制の各国で生じた不況の構造変化は、各国の社会保障費を押し上げたのか。また、そこで左派政党が果たした役割は、どのようなものだったと考えられるか。

メインの仮説の設定：

• メインの仮説1：1980年代に民主主義体制の各国で生じた不況の構造変化によって、各国の社会保障は拡充された。

利用するデータ：Comparative Welfare State dataset

変数の設定：

• 従属変数：社会保障費の対GDP比

• 処置群：構造変化が発生した国

• 統制群：構造変化が発生していない国

• メインの独立変数：左派政党の議席比率

満たされるべき仮定：処置が起きなかった場合の平均アウトカムの継時的変化は、処置群と対照群で一致する（平行トレンド仮定（common trend assumption））。

→[処置効果]=[処置後に観測される処置群の実際の変化]―[処置前と同じなら観測されたであろう変化]

二元変数の固定効果モデルによって、DIDは以下のように表される。

\[ Y_{it}=\alpha_i+\gamma_t+\tau_tD_{it}+\beta_1(\mathrm{leftseat}\times D_{it})+\epsilon_{it}. \]

各項のnotationは以下の通り。

• \(Y_{it}\): 国\(_i\)の年\(t\)における社会保障費対GDP比

• \(\alpha_i\): 国ごとの時間不変（time invariant）の固定効果

• \(\gamma_{t}\): 全体の年のショック

• \(D_{it}\): 構造変化後ダミー（各不況開始後=1、それより前=0）

• \(\beta_1\): 左派議席率の構造変化のもとでの係数の変化（差分の差分項）

#------ データの成形
df_cws <- df_cws %>%
  # ID変数とyear変数を因子化
  mutate(
    id_f   = factor(id),
    year_f = factor(year),
    # AS(2004)の記述に沿って構造変化ダミーを設定
    structural_break = case_when(
      id == "NZL" & year > 1978                             ~ 1,
      id == "NZL"                                          ~ 0,
      id %in% c("AUT","BEL","DNK","FIN","SWE","GBR") & year > 1981 ~ 1,
      id %in% c("AUT","BEL","DNK","FIN","SWE","GBR")               ~ 0,
      id %in% c("CAN","DEU","ITA","NLD","NOR","CHE","USA") & year > 1982 ~ 1,
      id %in% c("CAN","DEU","ITA","NLD","NOR","CHE","USA")               ~ 0,
      id %in% c("AUL","FRA","IRL","JPN") & year > 1983     ~ 1,
      id %in% c("AUL","FRA","IRL","JPN")                    ~ 0,
      TRUE                                                ~ NA_real_
    ),
    # DIDとそれとメインの変数との交互作用項の設定
    D      = structural_break,
    D_left = structural_break * leftseat
  ) %>%
  filter(!is.na(D)) 


#------ モデルの当てはめ
model_did <- lm(
  socx_pub ~ 
    D      +         # τ：構造変化による切片の変化
    leftseat         # 構造変化前のスロープ
    + D_left         # β：構造変化によるスロープ差分の差分
    + popold + vturn + singmemd + hunemr + cpi  # 共変量
    + id_f + year_f,                          # 国・年固定効果
  data = df_cws
)

#------　モデルの推定結果の解釈
coef(summary(model_did))

##                  Estimate Std. Error     t value     Pr(>|t|)
## (Intercept) -15.134395844 3.42856894 -4.41420199 1.500230e-05
## D             2.379308526 1.01711110  2.33928087 2.009650e-02
## leftseat     -0.032569130 0.02238309 -1.45507750 1.468825e-01
## D_left       -0.001254188 0.02015729 -0.06222009 9.504365e-01
## popold       83.918208377 8.54993926  9.81506486 1.735706e-19
## vturn         0.148101326 0.02936073  5.04419721 8.676883e-07
## singmemd      0.271207264 0.22242624  1.21931326 2.238564e-01
## hunemr        0.307168327 0.04818821  6.37434546 8.601061e-10
## cpi           0.101908097 0.01934047  5.26916217 2.926627e-07
## id_fBEL       6.776985786 0.60819985 11.14269553 9.697727e-24
## id_fCAN       4.437048133 0.92196380  4.81260557 2.558583e-06
## id_fFIN      11.930170931 0.94271114 12.65517120 8.653048e-29
## id_fFRA      12.655394316 0.92065980 13.74600514 1.615659e-32
## id_fITA       4.167613159 0.58572724  7.11527968 1.137403e-11
## id_fJPN      -0.946149242 1.11457942 -0.84888454 3.967449e-01
## id_fNOR       7.495076111 0.72016016 10.40751286 2.328732e-21
## id_fNZL       5.978635938 0.45709627 13.07959905 3.110915e-30
## id_fSWE      11.934212806 0.76669915 15.56570499 8.097942e-39
## id_fUSA       3.650778329 1.36799636  2.66870470 8.104736e-03
## year_f1981   -0.515455521 0.80427810 -0.64089215 5.221706e-01
## year_f1982   -0.504894369 0.82368117 -0.61297306 5.404426e-01
## year_f1983   -1.841326055 0.95765442 -1.92274584 5.563016e-02
## year_f1984   -3.140681533 1.05515107 -2.97652311 3.196802e-03
## year_f1985   -2.721446636 1.08624509 -2.50537072 1.285896e-02
## year_f1986   -3.319707673 1.12148323 -2.96010459 3.366101e-03
## year_f1987   -3.763216620 1.15390103 -3.26129930 1.260842e-03
## year_f1988   -4.006269671 1.17644230 -3.40541111 7.677331e-04
## year_f1989   -4.443019188 1.22762819 -3.61918962 3.565878e-04
## year_f1990   -4.453359409 1.28169476 -3.47458659 6.014404e-04
## year_f1991   -4.031205852 1.33679646 -3.01557190 2.825045e-03
## year_f1992   -3.622763516 1.37560921 -2.63357028 8.967627e-03
## year_f1993   -4.026557177 1.41609777 -2.84341751 4.826300e-03
## year_f1994   -4.487723108 1.44558714 -3.10442933 2.122343e-03
## year_f1995   -4.963861112 1.47725348 -3.36019591 8.986637e-04
## year_f1996   -5.139766614 1.50604469 -3.41275836 7.482207e-04
## year_f1997   -5.547583234 1.53104078 -3.62340657 3.511046e-04
## year_f1998   -5.339600867 1.54991584 -3.44509085 6.677332e-04
## year_f1999   -5.473870324 1.56901858 -3.48872243 5.719019e-04
## year_f2000   -6.039912123 1.60249110 -3.76907686 2.037874e-04
## year_f2001   -6.111385871 1.64372087 -3.71801926 2.470706e-04
## year_f2002   -5.773398349 1.67295461 -3.45101913 6.538837e-04
## year_f2003   -5.631633735 1.70719283 -3.29876839 1.110068e-03
## year_f2004   -6.040658138 1.73920266 -3.47323419 6.043400e-04
## year_f2005   -6.538758838 1.77460300 -3.68463191 2.799178e-04
## year_f2006   -7.071467980 1.81854440 -3.88853194 1.288203e-04
## year_f2007   -7.413688769 1.85675906 -3.99281142 8.552935e-05
## year_f2008   -7.182198939 1.92632992 -3.72843657 2.375900e-04
## year_f2009   -5.481725605 1.94316626 -2.82102758 5.165120e-03
## year_f2010   -6.328417406 2.13464711 -2.96461995 3.318740e-03

現時点で、\(D\)項が5％水準で統計的に有意であることが明らか。\(D\)項の係数がDIDにおけるAverage Treatment Effect on the Treated (ATT)で、構造変化（処置）を受けた後に、社会保障費対GDP比が平均でどれだけシフトしたかを表す。いまそれが5％水準で統計的に有意となっている。

さらに、構造変化前後の結果の平均の比較を図によって示し対比させてみよう。シンプルな平均の差分の描画だが、推定結果との関係でどのようなことを読み取れるだろうか。

library(grid) 

#------- 構造変化前後をもとにグループ化、そのもとでの平均と信頼区間を計算
df_summary <- df_cws %>%
  group_by(D) %>%
  summarise(
    mean_socx = mean(socx_pub, na.rm = TRUE),
    n         = sum(!is.na(socx_pub)),
    sd        = sd(socx_pub, na.rm = TRUE),
    se        = sd / sqrt(n),
    ci_mult   = qt(0.975, df = n - 1),
    lower     = mean_socx - ci_mult * se,
    upper     = mean_socx + ci_mult * se
  ) %>%
  # 
  mutate(
    label = paste0(
      round(mean_socx, 2), "\n[",
      round(lower, 2), ", ",
      round(upper, 2), "]"
    )
  )

#-------  図示
ggplot(df_summary, aes(x = factor(D, labels = c("前", "後")), y = mean_socx)) +
  geom_col(fill = "gray70", width = 0.6) +
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = 0.2, size = 0.8) +
  # ラベルをバーの右側に配置。hjust=0 で左揃え、nudge_x で少し右へ
  geom_label(aes(label = label),
             hjust       = 0,
             nudge_x     = 0.25,
             nudge_y     = 0, 
             fill         = "white",
             label.r      = unit(0.15, "lines"),
             size         = 3) +
  labs(
    x     = "構造変化",
    y     = "社会保障費対GDP比 (平均)",
    title = "構造変化前後の平均比較"
  ) +
  scale_x_discrete(expand = expansion(add = c(0.5, 1.5))) +
  theme_bw()

次に推定結果をもとに、左派議席率が社会保障費対GDP比に与える効果の違いを検討しよう。先の結果からも、\(\beta_1(\mathrm{leftseat}\times D_{it})\)は統計的に有意ではなく、構造変化のもとで左派議席率の影響が変化していたとは解釈できない結果になっている。これを、図示すると以下のようになる。

#----- 共変量の平均と固定効果のベースレベルを準備
cov_means <- df_cws %>% summarise(
  popold   = mean(popold,   na.rm=TRUE),
  vturn    = mean(vturn,    na.rm=TRUE),
  singmemd = mean(singmemd, na.rm=TRUE),
  hunemr   = mean(hunemr,   na.rm=TRUE),
  cpi      = mean(cpi,      na.rm=TRUE)
)

id_levels   <- levels(df_cws$id_f)
year_levels <- levels(df_cws$year_f)

base_id     <- id_levels[1]
base_year   <- year_levels[1]

#-----  model_did から係数を取り出す
co <- coef(model_did)
beta0      <- co["(Intercept)"]
beta1      <- co["leftseat"]
tau        <- co["D"]
delta_beta <- co["D_left"]

#----- x の範囲設定
x_min <- min(df_cws$leftseat, na.rm=TRUE)
x_max <- max(df_cws$leftseat, na.rm=TRUE)
x_seq <- seq(x_min, x_max, length.out = 50)

#-----  予測用データフレーム作成し、他の変数についても mutate() で追加
df_line <- tibble(
  leftseat = rep(x_seq, 2),
  period   = rep(c("前","後"), each = length(x_seq))
) %>%
  mutate(
    D       = if_else(period == "前", 0L, 1L),
    D_left  = D * leftseat,
    # 共変量は平均値
    popold   = cov_means$popold,
    vturn    = cov_means$vturn,
    singmemd = cov_means$singmemd,
    hunemr   = cov_means$hunemr,
    cpi      = cov_means$cpi,
    # 固定効果因子はベースレベルかつ元の因子レベルを指定
    id_f     = factor(base_id,   levels = id_levels),
    year_f   = factor(base_year, levels = year_levels)
  )

#----- predict() で標準誤差を取って、直線に関する95% CI 計算
pred <- predict(model_did, newdata = df_line, se.fit = TRUE)
df_line <- df_line %>%
  mutate(
    fit   = pred$fit,
    se    = pred$se.fit,
    lower = fit - 1.96 * se,
    upper = fit + 1.96 * se
  )

#-----  図におけるyとxの最小値、最大値を設定
y_pre_min  <- df_line %>% filter(period=="前",  leftseat==x_min) %>% pull(fit)
y_post_min <- df_line %>% filter(period=="後",  leftseat==x_min) %>% pull(fit)
y_pre_max  <- df_line %>% filter(period=="前",  leftseat==x_max) %>% pull(fit)
y_post_max <- df_line %>% filter(period=="後",  leftseat==x_max) %>% pull(fit)

#-----　図示
ggplot(df_line, aes(x = leftseat, color = period, fill = period)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, color = NA) +
  geom_line(aes(y = fit), size = 1) +
  
  # 切片の変化に関するτを矢印も含めて描画、説明を付記
  geom_segment(x = x_min, xend = x_min,
               y = y_pre_min, yend = y_post_min,
               arrow = arrow(length = unit(0.2, "cm")),
               color = "blue") +
  annotate("text",
           x     = x_min,
           y     = (y_pre_min + y_post_min)/2,
           label = "切片シフト (τ)",
           hjust = -0.1,
           color = "blue") +
  
  # 傾きの変化に関するδを矢印も含めて描画、説明を付記 
  geom_segment(x = x_max, xend = x_max,
               y = y_pre_max, yend = y_post_max,
               arrow = arrow(length = unit(0.2, "cm")),
               color = "darkgreen") +
  annotate("text",
           x     = x_max,
           y     = (y_pre_max + y_post_max)/2,
           label = "直線の傾きの変化 (δ)",
           hjust =  1.1,
           color = "darkgreen") +
  
  scale_color_manual(values = c("前"="steelblue","後"="tomato")) +
  scale_fill_manual(values  = c("前"="steelblue","後"="tomato")) +
  labs(
    x      = "左派議席率",
    y      = "社会保障費対GDP比（予測値）",
    color  = "時期",
    fill   = "時期",
    title  = "構造変化前後の回帰直線\n（左派議席率の効果への注目)"
  ) +
  theme_bw()

ではこれらの結果は、平行トレンド仮定を満たしているのだろうか。もし構造変化が起こらなければ、処置群と対照群の個体は時間変化のもとで平行した変化を示すという平行トレンド仮定(common trend assumption)を満たさなくてはならない。仮定が満たされているかを、チェックしておく。

この検定において、

• \(H_0\): 処置前の各国における年に対する回帰係数の平均は0

• \(H_1\): 平均の傾きは0ではない（国ごとに異なるトレンドがある）

df_cws %>%
  #各国ごとに「最初の構造変化年」を計算
  group_by(id) %>%
  mutate(
    first_break_year = min(year[structural_break == 1], na.rm = TRUE)
  ) %>%
  ungroup() %>%
  na.omit() %>%
  #その年より前だけを残す
  filter(year < first_break_year) %>%
  #国ごとにsocx_pub ~ year の傾きを回帰：トレンドの検出を行うため！
  group_by(id) %>%
  summarise(
    trend = coef(lm(socx_pub ~ year, data = cur_data()))[2]
  ) %>%
  ungroup() %>%
  #傾きの平均がゼロかどうかの t-test
  pull(trend) %>%
  t.test(mu = 0)

## 
##  One Sample t-test
## 
## data:  .
## t = 2.87, df = 3, p-value = 0.06405
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.07170661  1.38913831
## sample estimates:
## mean of x 
## 0.6587159

推定結果から、構造変化前における各国の年次トレンドは0.659[-0.0717,1.3891]となっており、信頼区間はゼロをまたいでいることから、真の傾きがゼロである可能性を排除できない\(t=2.87(p-value=0.064)\)。

帰無仮説「処置前のトレンドの傾きの平均はゼロ」＝平行トレンドが成り立っているを棄却できないことから、平行トレンド仮定を満たしていると解釈することになる。