研究の設計

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

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

①データの読み込み

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

③必要なデータのデータフレーム化

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

df_ex<-data.frame(year,id,socx_pub, pop,po65,postgini_2559,lisrpr_tot,
                    hunemr,leftseat,leftvot,rtseat,rtvot,vturn)
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
## 1    42.8   63.1  46.5 95.48
## 2    43.2   62.1  46.1 95.46
## 3    47.9   50.8  42.1 95.27
## 4    47.7   51.5  42.4 95.31
## 5    45.5   59.0  46.0 95.73
## 6    45.5   59.0  46.0 95.73

④記述統計を確かめる。

library(psych)
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
##                     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

⑤変数の設定

社会保障支出の分布をみると・・・

hist(socx_pub)

hist(po65)

ヒストグラムからは、高齢者人口の分布が偏っていることが明らか。後に見る残差の分布の観点から、分析に利用する変数の分布に著しい偏りがあることは望ましくない。よって変数を割合値に変換することで調整

• 独立変数である高齢者の対全人口比

peraged=po65/pop

調整後の値を改めてヒストグラムで確認すると、分布の偏りが解消され、分析し安い形状に変換されていることがわかる。

hist(peraged)

散布図から見る最小二乗法（Ordinary Least Square: OLS）

①次に、高齢者の人口割合と社会保障支出の関係性を散布図によって確認しよう。

plot(socx_pub~peraged,ylab="社会保障支出割合(対GDP比)",xlab="高齢者人口割合")

②散布図を見ると、両者の間には何らかの線形（linear）の関係があることが示唆される。ここに回帰直線を引いて確かめる。

plot(socx_pub~peraged,ylab="社会保障支出割合(対GDP比)",xlab="高齢者人口割合")
x<-peraged #x（独立変数）として高齢者人口割合を指定
y<-socx_pub #y（従属変数）として高齢者人口割合を指定
z<-lsfit(x,y) #回帰直線を算出
abline(z) #図に回帰直線を重ね描き

この回帰直線は、\(y_{welfare}=\beta_{0}+\beta_{1}x_{aged}+\epsilon_{t}.\)で表される。 ③回帰直線をもとに、全変動、回帰変動、残差変動とは何であるかを確認する。 \(i\)番目の観察\(y_i\)の残差は、 \[|y_{i}-\beta_{0}-\beta_{1}x_{i}|\] こうした個々の残差の平方和、

\[\sum^{n}_{i=1}(y_{i}-\beta_{0}-\beta_{1}x_{i})^2\] を最小化するのが最小二乗法。

ここでもう少し違った形式の散布図も出力してみよう。

# plot the relationship between the two
g = ggplot(data.frame(x = peraged, y = socx_pub), aes(x = peraged, y = socx_pub))
g = g + geom_point(size = 3, alpha = .2, colour = "black")
g = g + geom_point(size = 2, alpha = .2, colour = "red")
g = g +  geom_smooth(method=lm)
g

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では、次のコマンドが利用されるが・・・

df_lm<-na.omit(data.frame(socx_pub,peraged))

result_sols<-lm(socx_pub~peraged,data=df_lm)
summary(result_sols)

## 
## Call:
## lm(formula = socx_pub ~ peraged, data = df_lm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.4541  -3.2005   0.1049   3.0495  13.0503 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.659      1.010   2.632  0.00869 ** 
## peraged      127.634      6.947  18.372  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.173 on 649 degrees of freedom
## Multiple R-squared:  0.3421, Adjusted R-squared:  0.3411 
## F-statistic: 337.5 on 1 and 649 DF,  p-value: < 2.2e-16

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

# slope
beta1 <- cor(x=df_lm$socx_pub,y=df_lm$peraged)*(SD(df_lm$socx_pub)/SD(df_lm$peraged))
# intercept
beta0 <- mean(df_lm$socx_pub)-beta1*mean(df_lm$peraged)
# results are the same as using the lm command
results <- rbind("manual" = c(beta0, beta1), "lm(df_lm$socx_pub ~ df_lm$peraged)" = coef(lm(df_lm$socx_pub ~ df_lm$peraged)))
# set column names
colnames(results) <- c("切片（定数項）", "傾き（係数）")
# print results
results

##                                    切片（定数項） 傾き（係数）
## manual                                   2.658863     127.6341
## lm(df_lm$socx_pub ~ df_lm$peraged)       2.658863     127.6341

②分析結果の解釈（1）：推定量の解釈と検討

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

推定量を解釈する上でのポイント

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

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

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

③分析結果の解釈（2）：モデルの妥当性に関する解釈と検討

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

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

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

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

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

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

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

# 従属変数の中心化
mean.socx_pub <- df_lm$socx_pub-mean(df_lm$socx_pub)
# 独立変数の中心化
mean.peraged <- df_lm$peraged - mean(df_lm$peraged)
# compare correlations
results2 <- rbind("cor(社保, 高齢)" = cor(df_lm$socx_pub, df_lm$peraged), "cor(中心化・社保,中心化・高齢 )" = cor(mean.socx_pub, mean.peraged),
"中心化後の係数" = coef(lm(mean.socx_pub ~ mean.peraged))[2],"中心化後の定数項" = coef(lm(mean.socx_pub ~ mean.peraged))[1])
# print results
results2

##                                 mean.peraged
## cor(社保, 高齢)                 5.849266e-01
## cor(中心化・社保,中心化・高齢 ) 5.849266e-01
## 中心化後の係数                  1.276341e+02
## 中心化後の定数項                1.565575e-15

# plot the relationship between the two
g2 = ggplot(data.frame(x = mean.peraged, y = mean.socx_pub), aes(x = mean.peraged, y = mean.socx_pub))
g2 = g2 + geom_point(size = 3, alpha = .2, colour = "black")
g2 = g2 + geom_point(size = 2, alpha = .2, colour = "red")
g2 = g2 +  geom_smooth(method=lm)
g2

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

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

# 従属変数の正規化
norm.socx_pub <- (df_lm$socx_pub-mean(df_lm$socx_pub))/SD(df_lm$socx_pub)
# 独立変数の正規化
norm.peraged <- (df_lm$peraged - mean(df_lm$peraged))/SD(df_lm$peraged)
# 正規化した結果とそうでない結果を比較する。
results3 <- rbind("cor(社保, 高齢)" = cor(df_lm$socx_pub, df_lm$peraged), "cor(正規化・社保,正規化・高齢 )" = cor(norm.socx_pub, norm.peraged),
"正規化後の係数" = coef(lm(norm.socx_pub ~ norm.peraged))[2],
"正規化後の定数項" = coef(lm(norm.socx_pub ~ norm.peraged))[1])
# print results
results3

##                                 norm.peraged
## cor(社保, 高齢)                 5.849266e-01
## cor(正規化・社保,正規化・高齢 ) 5.849266e-01
## 正規化後の係数                  5.849266e-01
## 正規化後の定数項                6.666458e-17

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

# plot the relationship between the two
g3 = ggplot(data.frame(x = norm.peraged, y = norm.socx_pub), aes(x = norm.peraged, y = norm.socx_pub))
g3 = g3 + geom_point(size = 3, alpha = .2, colour = "black")
g3 = g3 + geom_point(size = 2, alpha = .2, colour = "red")
g3 = g3 +  geom_smooth(method=lm)
g3

⑤次に算出された残差(residuals)について検討していこう。以下の残差と予測値の計算をしっかり理解しておこう。

x<-df_lm$peraged
y<-df_lm$socx_pub

# "multiplot.R"関数を読み込む
source("multiplot.R")
#yの長さをとる
n <- length(y)
#単回帰分析を実行
fit <- lm(y ~ x)
# 残差を計算
e <- resid(fit)
# 予測値を計算
yhat <- predict(fit)
# 縦1行、横2行の図のパネルを作成
par(mfrow=c(1, 2))
# 残差を回帰直線に対してプロット
plot(x, y, xlab = "X:高齢者人口割合", ylab = "残差", bg = "lightblue",
col = "black", cex = 2, pch = 21,frame = FALSE,main = "回帰直線に対する残差")
# 回帰直線をプロット
abline(fit, lwd = 2)
# 赤い直線を各データポイントの残差から回帰直線に対して描画
for (i in 1 : n){lines(c(x[i], x[i]), c(y[i], yhat[i]), col = "red" , lwd = 2)}
# xに対する残差を描画
plot(x, e, xlab = "X:高齢者人口割合", ylab = "残差", bg = "lightblue",
col = "black", cex = 2, pch = 21,frame = FALSE,main = "残差 vs X")
# y=0線を描画
abline(h = 0, lwd = 2)
# 赤い直線をX軸に対して描画
for (i in 1 : n){lines(c(x[i], x[i]), c(e[i], 0), col = "red" , lwd = 2)}

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

⑥次に、「推論(inference)」のプロセスをより詳しく検討していこう。

2つの計算するべきパラメター\(\beta_0\)と\(\beta1\)。これらの計算の手順は、次のようなプロセスになっている。

# データの長さをとる
n <- length(y)
# beta1を計算
beta1 <- cor(y, x) * sd(y) / sd(x)
# beta0を計算
beta0 <- mean(y) - beta1 * mean(x)
# 誤差項を計算
e <- y - beta0 - beta1 * x
# 分散の不偏推定量を計算
sigma <- sqrt(sum(e^2) / (n-2))
# (X_i - X Bar)
ssx <- sum((x - mean(x))^2)
# 標準誤差を計算
seBeta0 <- (1 / n + mean(x) ^ 2 / ssx) ^ .5 * sigma
seBeta1 <- sigma / sqrt(ssx)
# 仮説検定 H0: beta0 = 0 and beta1 = 0
tBeta0 <- beta0 / seBeta0; tBeta1 <- beta1 / seBeta1
# p値を計算 Ha: beta0 != 0 and beta1 != 0 (両側検定)
pBeta0 <- 2 * pt(abs(tBeta0), df = n - 2, lower.tail = FALSE)
pBeta1 <- 2 * pt(abs(tBeta1), df = n - 2, lower.tail = FALSE)
# 分析結果を表で出力
coefTable <- rbind(c(beta0, seBeta0, tBeta0, pBeta0), c(beta1, seBeta1, tBeta1, pBeta1))
colnames(coefTable) <- c("Estimate", "Std. Error", "t value", "P(>|t|)")
rownames(coefTable) <- c("(Intercept)", "x")
# 表を表示
coefTable

##               Estimate Std. Error   t value      P(>|t|)
## (Intercept)   2.658863   1.010178  2.632073 8.688354e-03
## x           127.634050   6.947201 18.372010 5.149670e-61

上記の分析結果が、OLSの単回帰の分析結果と一致するか確認してみよう。

result<-lm(y~x)
summary(result)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.4541  -3.2005   0.1049   3.0495  13.0503 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.659      1.010   2.632  0.00869 ** 
## x            127.634      6.947  18.372  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.173 on 649 degrees of freedom
## Multiple R-squared:  0.3421, Adjusted R-squared:  0.3411 
## F-statistic: 337.5 on 1 and 649 DF,  p-value: < 2.2e-16

推論に関して確認するべきポイント

• \(t\)値と\(p\)値に注目：帰無仮説は棄却できるか

• 係数に注目：係数の値は、\(x\)の変化量に対して実質的に意味のある\(y\)の変化量を示唆しているか。

• 定数項に注目：\(x\)がゼロの時の値に、実質的に意味があるのかどうか、あらかじめ確認し必要ならば中心化、正規化の処置をとる

⑦最後に、「モデルの当てはまり(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\)は・・・

cor(y,x)

## [1] 0.5849266

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

(cor(y,x))^2

## [1] 0.3421391

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

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

①線形性：推定量\(\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）。

⇒推定量は不偏性を持たない