データ解析： さまざまな種類の重回帰分析
I.質的変数（離散型確率変数）を従属変数とする回帰分析
質的変数とは何か？：ダミー変数（二値変数、二項選択型変数）、名義変数、順序変数を含み、値と値との距離が任意に定まるもの。元のデータは、ほぼ文字情報。
質的変数の変動を説明するために：質的変数を従属変数にする場合、第４回授業「ガウス・マルコフ定理」に従うと、OLSを使った場合にBLUEを得られない。
⇒質的変数／離散型確率変数のための推定法を知る。OLSが線形モデルであったのに対して、質的変数の分析は非線形モデルとなる。
- 質的変数と推定モデル
- 線形モデルと非線形モデルの対照：線形モデルでは量的変動を説明・予測するのに対して、非線形モデルでは発生確率を説明・予測することに主眼が置かれる
(1) 二値変数を従属変数とする分析
ロジスティック回帰分析：従属変数の形態がダミー変数（二値変数、二項選択型変数）の場合
\(Y\)が0または1からなる、ロジスティック回帰モデルは以下の通り；
\[log(\frac{p(Y=1)}{1-p(Y=1)})=\beta_{o}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{k}x_{k}, i=1,...,n.\] \(Y=\{1,0\}\)のもと、\(Y=1\)になる確率\(p\)の説明が目的。
分析の設定
研究の設計：有権者の投票参加は何によって規定されているのか？
仮説の設定：
仮説1：政治的有効性感覚が高い有権者ほど、投票に行く
仮説2：政治的関心が高い有権者ほど、投票に行く
仮説3：政治的信頼が高い有権者ほど、投票に行く
仮説4：教育程度が高い有権者ほど、投票に行く
仮説5：年齢が高い有権者ほど、投票に行く
利用するデータ：Japan Election Study 2019
変数の設定
従属変数：投票参加(Q1)
カギとなる独立変数(key independent variables)： 政治的有効性感覚(aQ49_1)、政治関心(aQ29)、政治的信頼(aQ40_1)、教育(aQ59)、年齢(AGE)
制御変数(control variables)：党派性(aQ7)、性別(SEX)、職業(aQ60)、居住年数(aQ58)、収入(aQ68)
Rで実践するロジスティック回帰分析(1)：投票参加確率の図示とオッズ比の解釈
①データの読み込み
②変数とデータフレームの設定
#投票参加
votep <- Q1
votep[Q1 == 1] <- 1
votep[Q1 == 2] <- 1
votep[Q1 == 3] <- 0
#政治的有効性感覚
efficacy <- aQ49_1
#政党支持
psu_rul2019<-jes2019$aQ7
psu_rul2019[jes2019$aQ7==1|jes2019$aQ7==3]<-1
psu_rul2019[jes2019$aQ7==2|jes2019$aQ7==4|jes2019$aQ7==5|jes2019$aQ7==6|jes2019$aQ7==7|jes2019$aQ7==8|jes2019$aQ7==11|jes2019$aQ7==12]<-0
psu_rul2019[jes2019$aQ7==9]<-0
#性別
gender2019 <- jes2019$SEX
gender2019[jes2019$SEX == 2] <- 1
gender2019[jes2019$SEX == 1] <- 0
#年齢
age2019 = jes2019$AGE
#教育
educ2019 <- jes2019$aQ59
educ2019[jes2019$aQ59 == 4] <- 1
educ2019[jes2019$aQ59 == 3 | jes2019$Q61 == 2 |jes2019$Q61 == 1] <- 0
educ2019[jes2019$aQ59 == 5] <- NA
#職業
employment2019 <- jes2019$aQ60
employment2019[jes2019$aQ60 ==1 | jes2019$aQ60 == 2 | jes2019$aQ60 ==3 ] <- 1
employment2019[jes2019$aQ60 == 6 ] <- 0
employment2019[jes2019$aQ60 == 4 | jes2019$aQ60 == 5 | jes2019$aQ60 ==7 | jes2019$aQ60 == 8 ] <- NA
#居住年数
residy2019 <- jes2019$aQ58
residy2019[jes2019$aQ58 == 6] <- NA
#収入
income2019 <- jes2019$aQ68
income2019[jes2019$aQ68 == 88 | jes2019$aQ68 ==99] <- NA
#政治関心
interest2019 <- jes2019$aQ29
interest2019[jes2019$aQ29 ==1] <- 4
interest2019[jes2019$aQ29 ==2] <- 3
interest2019[jes2019$aQ29 ==3] <- 2
interest2019[jes2019$aQ29 ==4] <- 1
interest2019[jes2019$aQ29 ==5 | jes2019$aQ29 ==6] <- NA
#政権信頼
trust2019 <- jes2019$aQ40_1
trust2019[jes2019$aQ40_1 ==1] <- 4
trust2019[jes2019$aQ40_1 ==2] <- 3
trust2019[jes2019$aQ40_1 ==3] <- 2
trust2019[jes2019$aQ40_1 ==4] <- 1
trust2019[jes2019$aQ40_1 ==5 | jes2019$aQ40_1 ==6] <- NA
#データフレームの設定
df_part <- na.omit(data.frame(votep, efficacy, psu_rul2019, gender2019, age2019, educ2019,
employment2019, residy2019, income2019, interest2019, trust2019))③ロジスティック回帰分析の実行
#glm関数を使ったロジスティック回帰分析
model_full <- glm(data = df_part,
votep ~ psu_rul2019 + efficacy+ interest2019 + trust2019 + age2019 +
gender2019 + educ2019 + employment2019 + residy2019 + income2019 ,
family = binomial("logit"))
summary(model_full)##
## Call:
## glm(formula = votep ~ psu_rul2019 + efficacy + interest2019 +
## trust2019 + age2019 + gender2019 + educ2019 + employment2019 +
## residy2019 + income2019, family = binomial("logit"), data = df_part)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4923 -0.6332 0.4348 0.7042 2.0874
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.233032 0.717822 -5.897 3.70e-09 ***
## psu_rul2019 -0.117697 0.226605 -0.519 0.6035
## efficacy 0.648238 0.087010 7.450 9.33e-14 ***
## interest2019 0.961681 0.141312 6.805 1.01e-11 ***
## trust2019 0.306615 0.135392 2.265 0.0235 *
## age2019 -0.009388 0.009252 -1.015 0.3103
## gender2019 -0.002987 0.201568 -0.015 0.9882
## educ2019 -0.068594 0.121258 -0.566 0.5716
## employment2019 -0.088952 0.296303 -0.300 0.7640
## residy2019 0.106547 0.095824 1.112 0.2662
## income2019 0.027849 0.035652 0.781 0.4347
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 922.21 on 791 degrees of freedom
## Residual deviance: 726.21 on 781 degrees of freedom
## AIC: 748.21
##
## Number of Fisher Scoring iterations: 5分析結果から何がわかるだろうか？政治的有効性感覚、政治関心、政治的信頼、年齢の各要素の効果はどのように解釈できるだろうか？
上記の設定した仮説はどのように評価できるだろうか？
③ロジスティック回帰分析の精査。リンク関数であるロジスティック関数に当てはめるとはどういうことだろうか？
#政治的有効性感覚を例にロジスティック分布への当てはめを図示するために政治的有効性感覚のみからなるモデルを推定
model_eff <- glm(data = df_part,
votep ~ efficacy,
family = binomial("logit"))
summary(model_eff)##
## Call:
## glm(formula = votep ~ efficacy, family = binomial("logit"), data = df_part)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1299 -1.0795 0.4675 0.6756 1.6354
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.83047 0.27924 -6.555 5.56e-11 ***
## efficacy 0.79787 0.07873 10.135 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 922.21 on 791 degrees of freedom
## Residual deviance: 801.82 on 790 degrees of freedom
## AIC: 805.82
##
## Number of Fisher Scoring iterations: 4efficacyの係数値は「0.79787」。これは先ほどのロジスティック関数の式に当てはめると…
\[log(\frac{p(Y=1)}{1-p(Y=1)})=-1.830+0.7979x\]
\[\Leftrightarrow p=\frac{exp(-1.830+0.7979x)}{1+exp(-1.830+0.7979x)}\]
上式を利用してモデル化された曲線を作図すると、以下のようになる。
calc.prob <- function(xbeta) {
p <- exp(xbeta) / (1 + exp(xbeta))
return (p)
}
plot(efficacy, votep, xlim = range(efficacy), xlab = '政治的有効性感覚', ylab = '投票参加確率')
x.new <- seq(1, 5, 1)
lines(x.new, calc.prob(cbind(1, x.new) %*% coef(model_eff)))
y.predicted <- predict(model_eff, newdata = data.frame(efficacy = x.new), type = 'link', se.fit = TRUE)
alpha <- 0.05
ci.lower <- y.predicted$fit - (qnorm(1 - alpha / 2) * y.predicted$se.fit)
ci.upper <- y.predicted$fit + (qnorm(1 - alpha / 2) * y.predicted$se.fit)
lines(x.new, calc.prob(ci.lower), col = 'darkgray')
lines(x.new, calc.prob(ci.upper), col = 'darkgray')
さらに全ての変数に関しても投票確率を示す曲線を描くと…
library(broom)
# Predict the probability (p) of diabete positivity
probabilities <- predict(model_full, type = "response")
#量的変数のみを抽出
df_part.m <- df_part %>%
dplyr::select_if(is.numeric)
predictors <- colnames(df_part.m)
# ロジット分析の結果を統合し、図示のために成形
df_part.m <- df_part.m %>%
mutate(logit = log(probabilities/(1-probabilities))) %>%
gather(key = "predictors", value = "predictor.value", -logit)
#図示
ggplot(df_part.m, aes(logit, predictor.value))+
geom_point(size = 0.5, alpha = 0.5) +
geom_smooth(method = "loess") +
theme_bw() +
facet_wrap(~predictors, scales = "free_y")
④分析結果を実質的に解釈するために、オッズ比(odds ratio)を算出。
## OR 2.5 % 97.5 %
## (Intercept) 0.01450834 0.003459372 0.05790615
## psu_rul2019 0.88896525 0.569784487 1.38690132
## efficacy 1.91216901 1.616273991 2.27441389
## interest2019 2.61608965 1.992214317 3.46931762
## trust2019 1.35881806 1.043210888 1.77505978
## age2019 0.99065639 0.972749127 1.00873535
## gender2019 0.99701760 0.673068410 1.48474216
## educ2019 0.93370537 0.736926541 1.18614308
## employment2019 0.91488979 0.507993534 1.62709529
## residy2019 1.11243061 0.921359832 1.34221355
## income2019 1.02824024 0.959131411 1.10320432政治的有効性感覚のオッズ比が「1.912」という時の意味：政治的有効性感覚が1ポイント増えるにしたがって、投票に行く確率が1.912倍上がるという解釈
オッズ比の結果も踏まえて、上記に設定した仮説はどのように評価できるだろうか？
⑤モデルの決定係数の確認のためにマックファデンの疑似決定係数（Pseudo\(R^2\)）を計算。
## McFadden
## 0.2125329Rで実践するロジスティック回帰分析(2)：交差項(interaction terms)を含んだモデルの分析
- 分析を進めていく中で、仮説1と仮説2の組み合わせが重要ではないかと考えたとする
⇒仮説1×仮説2：選挙への関心が高いもとで、政治的有効性感覚が高い場合に、投票に行く確率が高くなる
関心のある2つ以上の変数\(X\)、\(Z\)を乗した項を加えることで、\(Z\)のもとでの\(X\)の効果の推移（=限界効果(marginal effect)）を測ることができる。
交差項を含んだモデル：
\[log(\frac{p(Y=1)}{1-p(Y=1)})=\beta_{o}+\beta_{1}x_{efficacy}+\beta_{2}x_{interest}+\beta_{3}x_{efficacy\times interest}+...+\beta_{k}x_{k}, i=1,...,n.\]
①交差項を含んだモデルの推定
- 交差項がダミー変数の時：男女（女性=1／男性=0）という性別のもとで、政治的有効性感覚が与える影響の違いを測る
まず下記は交差項を含まないモデル
model_dum <- glm(data = df_part,
votep ~ efficacy + gender2019,
family = binomial("logit"))
summary(model_dum)##
## Call:
## glm(formula = votep ~ efficacy + gender2019, family = binomial("logit"),
## data = df_part)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1839 -0.9962 0.5210 0.7455 1.7254
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.66359 0.29030 -5.731 1e-08 ***
## efficacy 0.79034 0.07892 10.015 <2e-16 ***
## gender2019 -0.35941 0.17786 -2.021 0.0433 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 922.21 on 791 degrees of freedom
## Residual deviance: 797.76 on 789 degrees of freedom
## AIC: 803.76
##
## Number of Fisher Scoring iterations: 4この分析結果は、性別のもとでどのように解釈できるか？
→下記表のように、定数項が変わって係数（＝接線の傾き）は変わらないと解釈できる
これを図で表すと…
plot(votep ~ efficacy, data = df_part,
col = as.factor(gender2019),
pch = 20,
xlab = "政治的有効性感覚",
ylab = "投票参加",
main = "性別ダミー変数に関する視覚化")
model_coef <- coefficients(model_dum)
model_coef## (Intercept) efficacy gender2019
## -1.6635874 0.7903439 -0.3594131# 性別が男性の時の回帰直線（性別=0）
abline(model_coef[1], model_coef[2], col = "black")
# 性別が女性の時の回帰直線(性別=1)
abline(model_coef[1] + model_coef[3], model_coef[2], col = "red")
# add a legend to the plot
legend("bottomright",
c("gender2019 = 0", "gender2019 = 1"),
col=c("black","red"),
lwd = 1, # line width = 1 for adding a line to the legend
pch = 20)
図からは、性別ダミー変数を考慮しても傾きは変わらず、定数項（切片）の位置だけが変わることがわかる。
では次に、交差項を含むモデルを推定
model_int_dum <- glm(data = df_part,
votep ~ efficacy + gender2019 + efficacy:gender2019,
family = binomial("logit"))
summary(model_int_dum)##
## Call:
## glm(formula = votep ~ efficacy + gender2019 + efficacy:gender2019,
## family = binomial("logit"), data = df_part)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2180 -1.0372 0.5522 0.7616 1.6447
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.83441 0.36459 -5.031 4.87e-07 ***
## efficacy 0.84096 0.10285 8.177 2.92e-16 ***
## gender2019 0.06683 0.56795 0.118 0.906
## efficacy:gender2019 -0.12674 0.16034 -0.790 0.429
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 922.21 on 791 degrees of freedom
## Residual deviance: 797.14 on 788 degrees of freedom
## AIC: 805.14
##
## Number of Fisher Scoring iterations: 4交差項の結果を図で表すために、Zeligパッケージを使う。
library(Zelig)
library(HH)
#パッケージZeligを使って、ロジテスティック回帰分析を実行
zelig_int <- zelig(data = df_part,
votep ~ efficacy + gender2019 + efficacy:gender2019,
model = "logit")## How to cite this model in Zelig:
## R Core Team. 2007.
## logit: Logistic Regression for Dichotomous Dependent Variables
## in Christine Choirat, Christopher Gandrud, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
## "Zelig: Everyone's Statistical Software," https://zeligproject.org/#性別が女性の時と男性の時でのシミュレーションの設定
x.out1 <- setx(zelig_int, gender2019 = 0, efficacy = 1:5)
x.out2 <- setx(zelig_int, gender2019 = 1, efficacy = 1:5)
# シミュレーションの実行
s.outcome <- sim(zelig_int, x = x.out1, x1 = x.out2)
#結果の図示
Zelig::ci.plot(s.outcome,
ci = 95,
xlab = "政治的有効性感覚",
ylab = "投票参加の期待値",
main = "性別ごとでの政治的有効性感覚の効果"
) 
女性と男性とによって、政治的有効性感覚が投票参加にもたらす効果は異ならないという結果
- 交差項がダミー変数ではない時：政治関心のもとで、政治的有効性感覚が与える影響の違いを測る
library(interplot)
#交差項を含むモデルの推定
model_int_dum <- glm(data = df_part,
votep ~ efficacy + interest2019 + efficacy:interest2019,
family = binomial("logit"))
summary(model_int_dum)##
## Call:
## glm(formula = votep ~ efficacy + interest2019 + efficacy:interest2019,
## family = binomial("logit"), data = df_part)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7862 -0.7615 0.3995 0.7845 1.8159
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.00440 0.99567 -2.013 0.0441 *
## efficacy 0.07491 0.29138 0.257 0.7971
## interest2019 0.27473 0.38916 0.706 0.4802
## efficacy:interest2019 0.21958 0.11163 1.967 0.0492 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 922.21 on 791 degrees of freedom
## Residual deviance: 731.13 on 788 degrees of freedom
## AIC: 739.13
##
## Number of Fisher Scoring iterations: 5#政治関心のもとでの政治的有効性感覚が与える影響の図示
interplot(m = model_int_dum, var1 = "efficacy", var2 = "interest2019")+
xlab("政治的関心の推移")+
ylab("投票参加")
分析結果の解釈：政治的関心が増加するもとで、政治的有効性感覚の効果も大きくなるが、その統計的有意性が認められるのは、政治的関心が中程度の時（2または3）。
プロビットモデルを使った分析
もう1つの二値変数を従属変数とする推定はプロビット・モデル。プロビットモデルでは、推定式の誤差項\(\mu_{i}\)が正規分布に従うと仮定し、確率を推定する。そのための確率密度関数は、以下の通り；
\[f(x)=\frac{1}{\sqrt{2\pi\sigma}}exp\{-\frac{(x-\mu)^2}{2\sigma^2} \}.\]
推定式は上式を積分したもので求められる；
\[P(Y_{i}=1)=\frac{1}{\sqrt{2\pi\sigma}}\int^{\beta_{0}+\beta_{1}X_{i}}_{-\infty}exp\{-\frac{(x-\mu)^2}{2\sigma^2}\}dx\]
Rでのプロビット分析の実行:分析結果は、上記のロジット分析の場合と大きく変わらないことを確かめてみよう。
#glm関数を使ったロジスティック回帰分析
model_probit <- glm(data = df_part,
votep ~ psu_rul2019 + efficacy+ interest2019 + trust2019 + age2019 +
gender2019 + educ2019 + employment2019 + residy2019 + income2019 ,
family = binomial("probit"))
summary(model_probit)##
## Call:
## glm(formula = votep ~ psu_rul2019 + efficacy + interest2019 +
## trust2019 + age2019 + gender2019 + educ2019 + employment2019 +
## residy2019 + income2019, family = binomial("probit"), data = df_part)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5452 -0.6577 0.4277 0.7156 2.0678
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.473611 0.409598 -6.039 1.55e-09 ***
## psu_rul2019 -0.083429 0.130431 -0.640 0.5224
## efficacy 0.377275 0.049715 7.589 3.23e-14 ***
## interest2019 0.555050 0.080697 6.878 6.06e-12 ***
## trust2019 0.181979 0.077234 2.356 0.0185 *
## age2019 -0.004586 0.005340 -0.859 0.3904
## gender2019 0.006353 0.116556 0.055 0.9565
## educ2019 -0.039214 0.070360 -0.557 0.5773
## employment2019 -0.053470 0.169517 -0.315 0.7524
## residy2019 0.057554 0.055629 1.035 0.3008
## income2019 0.015032 0.020610 0.729 0.4658
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 922.21 on 791 degrees of freedom
## Residual deviance: 725.87 on 781 degrees of freedom
## AIC: 747.87
##
## Number of Fisher Scoring iterations: 5(2) 順序変数を従属変数とする分析
分析の設定
研究の設計：有権者の政治的有効性感覚は何によって規定されているのか？
仮説の設定：
仮説1：政治的関心が高い有権者ほど、政治的有効性感覚が高い
仮説2：政治的信頼が高い有権者ほど、政治的有効性感覚が高い
仮説3：教育程度が高い有権者ほど、政治的有効性感覚が高い
仮説4：年齢が高い有権者ほど、政治的有効性感覚が高い
利用するデータ：Japan Election Study 2019
変数の設定
従属変数：政治的有効性感覚(aQ49_1)
カギとなる独立変数(key independent variables)： 政治関心(aQ29)、政治的信頼(aQ40_1)、教育(aQ59)、年齢(AGE)
制御変数(control variables)：党派性(aQ7)、性別(SEX)、職業(aQ60)、居住年数(aQ58)、収入(aQ68)
Rで実践する順序ロジスティック回帰分析：従属変数の形態が順序変数の場合
①順序ロジット分析の実行
#plor関数を使ったロジスティック回帰分析
library(MASS)
#順序ロジット推定
model_ologit <- polr(data = df_part,
as.factor(efficacy) ~ psu_rul2019 + interest2019 + trust2019 + age2019 +
gender2019 + educ2019 + employment2019 + residy2019 + income2019 )
summary(model_ologit)## Call:
## polr(formula = as.factor(efficacy) ~ psu_rul2019 + interest2019 +
## trust2019 + age2019 + gender2019 + educ2019 + employment2019 +
## residy2019 + income2019, data = df_part)
##
## Coefficients:
## Value Std. Error t value
## psu_rul2019 0.072022 0.159028 0.45289
## interest2019 0.931623 0.097887 9.51736
## trust2019 -0.097631 0.093882 -1.03993
## age2019 0.014772 0.006666 2.21623
## gender2019 -0.006869 0.144299 -0.04760
## educ2019 0.048917 0.088503 0.55272
## employment2019 -0.006281 0.210548 -0.02983
## residy2019 0.148999 0.070860 2.10273
## income2019 -0.015239 0.025454 -0.59867
##
## Intercepts:
## Value Std. Error t value
## 1|2 -0.3290 0.5159 -0.6377
## 2|3 1.6340 0.4887 3.3434
## 3|4 3.0686 0.4978 6.1637
## 4|5 4.2466 0.5088 8.3463
##
## Residual Deviance: 2094.377
## AIC: 2120.377#従属変数のefficacyをas.factor()にしていることに注意しよう！classを因子に変換している。numeric(数値)ではerrorが出る。
#オッズ比の算出
exp(coef(model_ologit))## psu_rul2019 interest2019 trust2019 age2019 gender2019
## 1.0746785 2.5386256 0.9069835 1.0148821 0.9931548
## educ2019 employment2019 residy2019 income2019
## 1.0501337 0.9937388 1.1606723 0.9848769② 分析結果の解釈：分析結果はどのように解釈できるだろうか？interceptが段階ごとに異なる値になっている。1から2へ、2から3へと推移するときそれぞれに、定数項（切片）が異なることに注意。
③ 交差項を含んだモデルについても検討してみよう。所得の推移のもとでの政治的有効性感覚の効果を測る。連続変数である所得のもとでの交差項の解釈でもある。下記の分析結果はどのように解釈できるだろうか。
library(interplot)
#交差項を含むモデルの推定
model_ologit_int <- polr(data = df_part,
as.factor(efficacy) ~ psu_rul2019 + interest2019 + trust2019 + age2019 +
gender2019 + educ2019 + employment2019 + residy2019 + income2019+income2019:interest2019 )
#政治関心のもとでの政治的有効性感覚が与える影響の図示
interplot(m = model_ologit_int, var1 = "interest2019", var2 = "income2019")+
xlab("所得の推移の推移")+
ylab("投票参加")
II. 時系列データの分析
時系列データの分析とは何か？
- 時系列データ(time-series data):ある特定の集団、地域、国に関する長期間のデータが集積されたもの。時間とともに変動する現象が記録されたものであることから、代表的な社会科学が扱うデータ例として、経済成長率、日経平均株価、支持率などのデータがある
- 時系列データ分析ですること
可視化と特性の理解：時系列データを図で示し、記述統計により推移・動態の特徴を明らかにする
モデリング：時系列データの特徴、扱うべき問題を踏まえた上でモデルを作る
予測：作成したモデル（＝過去と現実の説明）をもとに、将来の推移・動態を予測する
⇒時系列データの分析において重要なことは何か？：モデリングと分析をするために、時系列データの特性を踏まえた「前処理」をすること。では、問題となる時系列データの特性とは何なのか？
時系列データを扱うに際して注意すべき特性とは何か？
- 時系列データは時間的に連続している：現在の値は過去の値に規定される。
⇒自己相関（autoregression）、系列相関(serial correlation )の問題：\(y_{t}\)と\(y_{t-1}\)が相関しているとき、原系列（original series）は自己相関があり、系列相関がある判断する。系列相関がある変数をもとにした最小二乗法は有効性をもたない
- 時系列データの中には頻繁に平均に戻るものもあれば、増えたら増えたまま、減ったら減ったままのものもある
⇒系列の定常性(stationary)、非定常性(non-stationary)の問題：\(y_{t}\)と\(y_{t-j}\)の同時分布が\(t\)に依存しないならば（期待値が\(t\)に依存しない／分散が\(t\)に依存しない／共分散が\(t\)に依存しない）、系列は弱定常で、平均が0で1次以上の自己共分散が全て0の弱定常な時系列をホワイトノイズと呼ぶ。これに対して、上昇／下降傾向といったトレンドを持つ場合や構造変化がある場合、系列は非定常になる
- データが非定常な場合に何が問題になるのか？
⇒系列が非定常で、\(y_{t+1}-y_{t}\)の差分系列が定常になるとき、データは単位根(unit root)をもつ（1次の和分過程(\(I(1)\))）。また、ある系列が、
\[y_{t}=c+y_{t-1}+\mu_{t}\] で表現されるとき、その非定常時系列モデルは、ランダム・ウオーク（random walk process）過程に従うという。なお、\(c\)はドリフト項である。 単位根の有無を確認し、単位根を持つ系列を含むモデルならば、単位根を考慮した推定モデルを試す必要がある。
- 非定常なデータどうし(\(x_{t}\)、\(y_{t}\))の線形和(\(x_{t}+\beta y_{t}=z_{t}\))が定常過程となるならばどうか？
⇒非定常なデータの線形和は多くの場合、非定常。しかし定常になるならば、両者のかく乱項の和が定常になっていることを意味する。この時、\(x_{t}\)と\(y_{t}\)は共和分(cointegration)の関係にあるとする。
- 時系列のデータは互いに相互に入り組んだ関係にある
⇒変数間の内生性（endogeneity）を考慮しつつ、多変量時系列データを、あたかも構造方程式モデリングのような発想に立って分析できる推定方法があればよい。ベクトル自己回帰モデル(Vector Autoregressive model)の利用。
Rを使った時系列データの分析
時系列データの特性の確認
① データセットの読み込みと設定
#データの読み込み
df_app <- read.csv("approve_econ.csv", fileEncoding = "SJIS")
attach(df_app)
#変数の設定
#主観的社会志向の経済評価
df_app$sociotropic_gd<-(q5c1+q5c2) #肯定的方向
df_app$sociotropic_bd<-(q5c3+q5c4) #否定的方向
#主観的個人志向の経済評価
df_app$liv_gd<-(q3c1+q3c2) #肯定的方向
df_app$liv_bd<-(q3c4+q3c5) #否定的方向② 時系列データの描画
## 内閣支持率と与党支持率
library(tidyverse)
library(lubridate)
#グラフ用のデータの設定
df_aprul <- data.frame(yearmon, approve, ruling)
colnames(df_aprul)<-c("yearmon","内閣支持率","与党支持率")
tb_long1 <- df_aprul %>% tidyr::gather(支持率,value,-yearmon)
# 2つの系列を併記
p1 <- ggplot2::ggplot(tb_long1, aes(x = as.Date(yearmon,"%Y-%m-%d"), y = value)) +
geom_line(aes(linetype= 支持率), size = 1) +
theme_minimal()+
xlab("年")+
ylab("支持率")+
ggtitle("与党支持率と内閣支持率")
p1
#上記のプロットに経済イベントをテキストとして挿入
windows(20, 10)
ggplot2::ggplot(tb_long1, aes(x = as.Date(yearmon,"%Y-%m-%d"), y = value)) +
geom_line(aes(linetype= 支持率), size = 1) +
theme_minimal()+
xlab("年")+
ylab("支持率")+
annotate("text", x = as.Date("1976-02-01", "%Y-%m-%d"),
y = c(0.35), label = c("第1次石油ショック(1973)"))+ #テキストの描画
geom_vline(xintercept = as.Date("1973-10-01"))+ #指定した年月に垂線を言描画
annotate("text", x = as.Date("1986-01-01", "%Y-%m-%d"),
y = c(0.35), label = c("消費税導入(1989)"))+
geom_vline(xintercept = as.Date("1989-04-01"))+
annotate("text", x = as.Date("1992-12-01", "%Y-%m-%d"),
y = c(0.35), label = c("バブル経済崩壊(1990)"))+
geom_vline(xintercept = as.Date("1990-01-01"))+
annotate("text", x = as.Date("2000-09-01", "%Y-%m-%d"),
y = c(0.35), label = c("アジア通貨危機(1997)"))+
geom_vline(xintercept = as.Date("1997-07-01"))+
annotate("text", x = as.Date("2004-11-01", "%Y-%m-%d"),
y = c(0.35), label = c("グローバル経済危機(2008)"))+
geom_vline(xintercept = as.Date("2008-09-01"))+
annotate("text", x = as.Date("2016-03-01", "%Y-%m-%d"),
y = c(0.35), label = c("東日本大震災(2011)"))+
geom_vline(xintercept = as.Date("2011-03-01"))+
annotate("text", x = as.Date("2015-03-01", "%Y-%m-%d"),
y = c(0.35), label = c("アベノミクス導入(2012)"))+
geom_vline(xintercept = as.Date("2012-12-01"))+
ggtitle("経済的出来事と内閣支持率・与党支持率")
## 内閣支持率と社会志向の経済評価
#グラフ用のデータの設定
df_aprul <- data.frame(yearmon, approve, df_app$sociotropic_gd)
colnames(df_aprul)<-c("yearmon","内閣支持率","社会志向")
tb_long1 <- df_aprul %>% tidyr::gather(系列,value,-yearmon)
# 2つの系列を併記
p2 <- ggplot2::ggplot(tb_long1, aes(x = as.Date(yearmon,"%Y-%m-%d"), y = value)) +
geom_line(aes(linetype= 系列), size = 1) +
theme_minimal()+
xlab("年")+
ylab("割合")+
ggtitle("社会志向の経済評価と内閣支持率")
## 内閣支持率と個人志向の経済評価
#グラフ用のデータの設定
df_aprul <- data.frame(yearmon, approve, df_app$liv_gd)
colnames(df_aprul)<-c("yearmon","内閣支持率","個人志向")
tb_long1 <- df_aprul %>% tidyr::gather(系列,value,-yearmon)
# 2つの系列を併記
p3 <- ggplot2::ggplot(tb_long1, aes(x = as.Date(yearmon,"%Y-%m-%d"), y = value)) +
geom_line(aes(linetype= 系列), size = 1) +
theme_minimal()+
xlab("年")+
ylab("割合")+
ggtitle("個人志向の経済評価と内閣支持率")
##支持率と経済評価に関する図の統合
library(gridExtra)
windows(30,40)
grid.arrange(p1, p2, p3, ncol = 1)
経済評価は何らかのトレンドを持っているように見えるが、内閣支持率は頻繁な後進に特徴づけられているように見える。平均回帰的な性質を持っていることがうかがえるが…
③ 各系列の性質について考える：自己相関(autocorrelation)
コレログラムにより自己相関を確認すると、10次以上の自己相関が確認できる。
#自己回帰の確認：Autoregression(AR)モデルへの当てはめ
result_ar <- ar(approve, aic = T)
#自己回帰のモデルの次数の特定
result_ar$order## [1] 13## [1] 0.850 -0.043 0.094 -0.097 0.017 0.001 0.039 0.030 -0.115 0.007
## [11] -0.028 0.182 -0.111この結果から、ARモデルに内閣支持率を当てはめると次式で表すことができる；
\[ \begin{aligned} y_{t}=0.850y_{t-1}-0.043y_{t-2}+0.094y_{t-3}-0.097y_{t-4}+0.017y_{t-5}+0.001y_{t-6}+0.039y_{t-7}\\ +0.030y_{t-8}-0.115y_{t-9}+0.007y_{t-10}-0.028y_{t-11}+0.182y_{t-12}-0.111y_{t-13}+\epsilon_{t}. \end{aligned} \]
④単位根の確認のために、Augmented Dicky-Fueler testを実行。
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.694 -3.219 -0.430 2.442 56.909
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.108225 1.176846 5.190 2.79e-07 ***
## z.lag.1 -0.166702 0.031519 -5.289 1.67e-07 ***
## z.diff.lag1 0.023775 0.044401 0.535 0.59251
## z.diff.lag2 -0.023986 0.044343 -0.541 0.58874
## z.diff.lag3 0.065981 0.043377 1.521 0.12871
## z.diff.lag4 -0.024010 0.042686 -0.562 0.57398
## z.diff.lag5 -0.009595 0.041808 -0.230 0.81854
## z.diff.lag6 -0.004766 0.041625 -0.115 0.90887
## z.diff.lag7 0.036325 0.041279 0.880 0.37918
## z.diff.lag8 0.059768 0.040767 1.466 0.14309
## z.diff.lag9 -0.057824 0.040241 -1.437 0.15120
## z.diff.lag10 -0.048237 0.039607 -1.218 0.22370
## z.diff.lag11 -0.078776 0.039455 -1.997 0.04627 *
## z.diff.lag12 0.109466 0.038854 2.817 0.00498 **
## z.diff.lag13 -0.029250 0.038677 -0.756 0.44976
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.201 on 671 degrees of freedom
## Multiple R-squared: 0.1169, Adjusted R-squared: 0.09843
## F-statistic: 6.342 on 14 and 671 DF, p-value: 4.829e-12
##
##
## Value of test-statistic is: -5.289 13.9879
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1 6.43 4.59 3.78-5.289がADF検定の\(t\)統計量。臨界値であるcritical values for test statisticsのtau2の値と\(t\)統計量の比較のもとに、単位根を確認。いま1％の値が－3.43で\(t\)統計量よりも十分に小さい。よって、単位根の存在は棄却できる。
⑤共和分の確認のために、内閣支持率と社会志向の経済評価の金木を評価する。内閣支持率を従属変数とする回帰モデルの残差が定常か否かを調べればよい。
#共和分の確認
result_coint <- lm(data = df_app, approve ~ sociotropic_gd)
resid_coint <- resid(result_coint)
#残差のプロットから定常性／非定常性を視覚的に確認
plot(resid_coint, type="l") #図示からは、残差は定常であると考えられるが…
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.044 -3.146 -0.399 2.214 57.882
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.009465 0.236982 0.040 0.968
## z.lag.1 -0.166918 0.025942 -6.434 2.33e-10 ***
## z.diff.lag1 0.007519 0.040501 0.186 0.853
## z.diff.lag2 -0.029748 0.039812 -0.747 0.455
## z.diff.lag3 0.037212 0.039538 0.941 0.347
## z.diff.lag4 -0.025138 0.038832 -0.647 0.518
## z.diff.lag5 -0.009184 0.038375 -0.239 0.811
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.229 on 684 degrees of freedom
## Multiple R-squared: 0.08692, Adjusted R-squared: 0.07891
## F-statistic: 10.85 on 6 and 684 DF, p-value: 1.48e-11
##
##
## Value of test-statistic is: -6.4342 20.7049
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1 6.43 4.59 3.78図示からも残差の定常性は肯定的だが、残差に対するADF検定の結果から、1％水準で残差の非定常性は棄却できる=単位根は認められないという結果。すなわち、内閣支持率と社会志向の経済評価の間には、共和分が成立していると考えられる。では、共和分が成立している変数間のモデルとして、何を採用するのがよいか？
⇒いくつか方法が考えられるが、代表的な推定モデルとして誤差修正モデル(Error Correction Model: ECM)を採用
時系列データの推定方法の例①：グレンジャー因果性の確認とベクトル自己回帰(VAR)モデル
グレンジャー因果性(Granger causality)とは：現在と過去の\(x_{t}\)の値だけに基づく\(x\)の予測値と、現在と過去の\(x_{t}\)と\(y_{t}\) に基づいた将来の\(\bar{x}\)の予測を比較した場合に、後者の平均平方誤差(mean squared error: MSE)の方が小さくなる場合、\(y_t\)から\(x_t\)へのグレンジャー因果性があると判断
グレンジャー因果性はどのようにすれば確かめられるのか？：VARモデルの推定のもと、その\(F\)検定の結果を参照すればよい。
①VARモデルの推定
#VARモデルの推定
library(vars)
#推定に含めるすべての変数をモデルに投入
df_var <- na.omit(data.frame(approve, ruling, df_app$sociotropic_gd, df_app$liv_gd,
cpi, unemp, nav_high))
#var modelの次数の設定
VARselect(df_var)## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 2 1 3
##
## $criteria
## 1 2 3 4 5
## AIC(n) 13.48232 13.27584 13.15963 13.19292 13.24357
## HQ(n) 13.62525 13.54385 13.55270 13.71106 13.88679
## SC(n) 13.85176 13.96855 14.17561 14.53217 14.90609
## FPE(n) 716636.49105 582968.67244 519063.15490 536740.44305 564810.20665
## 6 7 8 9 10
## AIC(n) 13.30416 13.34746 13.40183 13.44779 13.51434
## HQ(n) 14.07244 14.24082 14.42025 14.59129 14.78290
## SC(n) 15.28994 15.65651 16.03414 16.40338 16.79319
## FPE(n) 600370.81401 627355.79537 662982.27125 694944.50028 743797.62097## $approve
## Estimate Std. Error t value Pr(>|t|)
## approve.l1 0.8095409995 0.0438876878 18.445742763 8.782464e-62
## ruling.l1 0.0956525945 0.1038000207 0.921508434 3.571157e-01
## df_app.sociotropic_gd.l1 0.1691479560 0.1174900803 1.439678615 1.504240e-01
## df_app.liv_gd.l1 0.0816239496 0.2714499411 0.300696140 7.637392e-01
## cpi.l1 0.3922015813 0.6571331391 0.596837319 5.508170e-01
## unemp.l1 -0.3845487522 2.4136531654 -0.159322291 8.734628e-01
## nav_high.l1 0.0000640893 0.0003480576 0.184134186 8.539637e-01
## approve.l2 -0.0309833889 0.0554228907 -0.559035961 5.763234e-01
## ruling.l2 -0.0050115324 0.1202830476 -0.041664495 9.667785e-01
## df_app.sociotropic_gd.l2 0.0004706630 0.1399221712 0.003363749 9.973171e-01
## df_app.liv_gd.l2 -0.2293627656 0.2873182804 -0.798288105 4.249853e-01
## cpi.l2 -0.7788589947 0.9874343300 -0.788770424 4.305244e-01
## unemp.l2 3.1764224377 3.0327894533 1.047360025 2.953101e-01
## nav_high.l2 0.0001031753 0.0005283767 0.195268531 8.452418e-01
## approve.l3 0.0231565463 0.0435094315 0.532219005 5.947502e-01
## ruling.l3 -0.0365861188 0.1028866470 -0.355596376 7.222545e-01
## df_app.sociotropic_gd.l3 -0.0983670713 0.1176789911 -0.835893225 4.035121e-01
## df_app.liv_gd.l3 0.2552928862 0.2416845179 1.056306330 2.912079e-01
## cpi.l3 0.4109332940 0.6663011048 0.616738125 5.376164e-01
## unemp.l3 -2.4629825448 2.3942050186 -1.028726665 3.039783e-01
## nav_high.l3 -0.0001681264 0.0003501777 -0.480117226 6.313003e-01
## const 1.6854661573 2.8155320485 0.598631494 5.496204e-01
##
## $ruling
## Estimate Std. Error t value Pr(>|t|)
## approve.l1 3.153979e-02 0.0183412474 1.7196099 8.596382e-02
## ruling.l1 6.151033e-01 0.0433794067 14.1796152 4.072274e-40
## df_app.sociotropic_gd.l1 1.106000e-01 0.0491006644 2.2525154 2.461159e-02
## df_app.liv_gd.l1 5.403209e-02 0.1134425342 0.4762948 6.340192e-01
## cpi.l1 -2.583929e-01 0.2746246631 -0.9408947 3.470969e-01
## unemp.l1 8.837530e-01 1.0086977020 0.8761326 3.812710e-01
## nav_high.l1 -1.713844e-04 0.0001454579 -1.1782404 2.391179e-01
## approve.l2 -2.213472e-02 0.0231619618 -0.9556498 3.395929e-01
## ruling.l2 1.544348e-01 0.0502678825 3.0722352 2.210480e-03
## df_app.sociotropic_gd.l2 -2.665601e-02 0.0584753330 -0.4558505 6.486448e-01
## df_app.liv_gd.l2 -7.917950e-02 0.1200741239 -0.6594218 5.098508e-01
## cpi.l2 -3.934119e-01 0.4126619159 -0.9533517 3.407547e-01
## unemp.l2 -9.298156e-01 1.2674429765 -0.7336153 4.634391e-01
## nav_high.l2 2.663633e-04 0.0002208156 1.2062699 2.281379e-01
## approve.l3 -6.842779e-03 0.0181831691 -0.3763249 7.067942e-01
## ruling.l3 1.063750e-01 0.0429976957 2.4739691 1.360807e-02
## df_app.sociotropic_gd.l3 -4.723860e-02 0.0491796127 -0.9605323 3.371329e-01
## df_app.liv_gd.l3 -3.900559e-02 0.1010031687 -0.3861819 6.994842e-01
## cpi.l3 6.517718e-01 0.2784560776 2.3406628 1.954108e-02
## unemp.l3 -3.736520e-01 1.0005700632 -0.3734391 7.089394e-01
## nav_high.l3 -9.055145e-05 0.0001463439 -0.6187578 5.362856e-01
## const 4.776809e+00 1.1766482227 4.0596746 5.492623e-05
##
## $df_app.sociotropic_gd
## Estimate Std. Error t value Pr(>|t|)
## approve.l1 2.445962e-03 0.0155717052 0.15707734 8.752311e-01
## ruling.l1 -1.096740e-02 0.0368290836 -0.29779180 7.659541e-01
## df_app.sociotropic_gd.l1 7.530477e-01 0.0416864270 18.06457874 9.326144e-60
## df_app.liv_gd.l1 7.185493e-02 0.0963126260 0.74605926 4.558925e-01
## cpi.l1 -1.810569e-01 0.2331561319 -0.77654775 4.376989e-01
## unemp.l1 -5.419674e-02 0.8563835884 -0.06328559 9.495579e-01
## nav_high.l1 5.067598e-04 0.0001234936 4.10352919 4.568878e-05
## approve.l2 -1.772703e-02 0.0196644881 -0.90147442 3.676591e-01
## ruling.l2 5.066340e-02 0.0426773943 1.18712489 2.355978e-01
## df_app.sociotropic_gd.l2 1.702868e-01 0.0496455136 3.43005437 6.403692e-04
## df_app.liv_gd.l2 -2.933513e-01 0.1019428406 -2.87760565 4.134234e-03
## cpi.l2 1.892882e-01 0.3503496554 0.54028381 5.891803e-01
## unemp.l2 -5.100300e-01 1.0760581313 -0.47398000 6.356681e-01
## nav_high.l2 -4.222213e-04 0.0001874723 -2.25217935 2.463293e-02
## approve.l3 2.506624e-02 0.0154374968 1.62372431 1.049036e-01
## ruling.l3 -6.105520e-02 0.0365050113 -1.67251545 9.488819e-02
## df_app.sociotropic_gd.l3 -7.364843e-02 0.0417534541 -1.76388838 7.820511e-02
## df_app.liv_gd.l3 2.036154e-01 0.0857516140 2.37447899 1.785361e-02
## cpi.l3 -4.545446e-02 0.2364089999 -0.19227045 8.475885e-01
## unemp.l3 7.600978e-01 0.8494832292 0.89477673 3.712268e-01
## nav_high.l3 -2.570735e-05 0.0001242459 -0.20690703 8.361451e-01
## const 3.202093e+00 0.9989734538 3.20538305 1.412692e-03
##
## $df_app.liv_gd
## Estimate Std. Error t value Pr(>|t|)
## approve.l1 6.031081e-03 6.474060e-03 0.93157635 3.518901e-01
## ruling.l1 -1.201574e-02 1.531198e-02 -0.78472781 4.328897e-01
## df_app.sociotropic_gd.l1 9.686309e-02 1.733146e-02 5.58885786 3.324032e-08
## df_app.liv_gd.l1 4.953774e-01 4.004274e-02 12.37121740 8.135207e-32
## cpi.l1 -4.772404e-01 9.693652e-02 -4.92322625 1.072376e-06
## unemp.l1 3.780889e-01 3.560483e-01 1.06190339 2.886610e-01
## nav_high.l1 7.269541e-05 5.134346e-05 1.41586485 1.572784e-01
## approve.l2 2.823068e-03 8.175667e-03 0.34530129 7.299759e-01
## ruling.l2 3.881440e-03 1.774347e-02 0.21875321 8.269087e-01
## df_app.sociotropic_gd.l2 -5.191971e-02 2.064052e-02 -2.51542660 1.212097e-02
## df_app.liv_gd.l2 4.172583e-02 4.238355e-02 0.98448168 3.252331e-01
## cpi.l2 4.861476e-01 1.456607e-01 3.33753512 8.918510e-04
## unemp.l2 -3.640685e-01 4.473797e-01 -0.81377962 4.160595e-01
## nav_high.l2 -1.639114e-04 7.794310e-05 -2.10296220 3.584035e-02
## approve.l3 8.881837e-05 6.418262e-03 0.01383838 9.889630e-01
## ruling.l3 1.021475e-02 1.517725e-02 0.67303054 5.011592e-01
## df_app.sociotropic_gd.l3 -8.883555e-03 1.735933e-02 -0.51174525 6.089974e-01
## df_app.liv_gd.l3 2.169887e-01 3.565192e-02 6.08631174 1.941489e-09
## cpi.l3 -2.987657e-02 9.828892e-02 -0.30396686 7.612472e-01
## unemp.l3 -3.373635e-02 3.531794e-01 -0.09552185 9.239288e-01
## nav_high.l3 1.079137e-04 5.165622e-05 2.08907576 3.707652e-02
## const 2.154914e+00 4.153312e-01 5.18842405 2.811610e-07
##
## $cpi
## Estimate Std. Error t value Pr(>|t|)
## approve.l1 -3.620165e-03 2.517604e-03 -1.43794073 1.509164e-01
## ruling.l1 -3.691909e-03 5.954456e-03 -0.62002452 5.354518e-01
## df_app.sociotropic_gd.l1 4.281918e-03 6.739783e-03 0.63531995 5.254360e-01
## df_app.liv_gd.l1 -2.184063e-02 1.557164e-02 -1.40258972 1.612010e-01
## cpi.l1 1.156366e+00 3.769624e-02 30.67589237 6.833548e-130
## unemp.l1 -7.998084e-02 1.384585e-01 -0.57765214 5.636925e-01
## nav_high.l1 -1.796805e-05 1.996622e-05 -0.89992242 3.684839e-01
## approve.l2 -2.646923e-03 3.179317e-03 -0.83254460 4.053974e-01
## ruling.l2 3.385378e-03 6.900000e-03 0.49063446 6.238452e-01
## df_app.sociotropic_gd.l2 1.836357e-03 8.026593e-03 0.22878409 8.191063e-01
## df_app.liv_gd.l2 -2.342172e-02 1.648193e-02 -1.42105506 1.557646e-01
## cpi.l2 -3.979681e-01 5.664387e-02 -7.02579242 5.232809e-12
## unemp.l2 -8.101073e-02 1.739750e-01 -0.46564571 6.416200e-01
## nav_high.l2 1.714794e-05 3.031017e-05 0.56574862 5.717535e-01
## approve.l3 5.582911e-03 2.495905e-03 2.23682800 2.562480e-02
## ruling.l3 -6.455183e-03 5.902061e-03 -1.09371673 2.744711e-01
## df_app.sociotropic_gd.l3 7.338356e-03 6.750619e-03 1.08706404 2.773983e-01
## df_app.liv_gd.l3 6.978790e-04 1.386416e-02 0.05033691 9.598689e-01
## cpi.l3 2.412480e-01 3.822216e-02 6.31173010 5.010598e-10
## unemp.l3 3.125009e-02 1.373428e-01 0.22753346 8.200782e-01
## nav_high.l3 -3.691608e-06 2.008784e-05 -0.18377328 8.542468e-01
## const 9.303984e-01 1.615121e-01 5.76054783 1.276617e-08
##
## $unemp
## Estimate Std. Error t value Pr(>|t|)
## approve.l1 5.301962e-04 6.965949e-04 0.7611256 4.468491e-01
## ruling.l1 1.530313e-03 1.647536e-03 0.9288494 3.533007e-01
## df_app.sociotropic_gd.l1 -1.687870e-03 1.864828e-03 -0.9051076 3.657328e-01
## df_app.liv_gd.l1 -2.093440e-03 4.308512e-03 -0.4858846 6.272074e-01
## cpi.l1 2.318928e-03 1.043016e-02 0.2223291 8.241253e-01
## unemp.l1 7.575111e-01 3.831002e-02 19.7731826 6.137546e-69
## nav_high.l1 -8.029671e-06 5.524445e-06 -1.4534800 1.465573e-01
## approve.l2 -1.746541e-04 8.796841e-04 -0.1985418 8.426812e-01
## ruling.l2 -3.115068e-03 1.909159e-03 -1.6316447 1.032229e-01
## df_app.sociotropic_gd.l2 -1.083249e-03 2.220875e-03 -0.4877577 6.258806e-01
## df_app.liv_gd.l2 -1.879393e-03 4.560378e-03 -0.4121135 6.803877e-01
## cpi.l2 -1.485592e-02 1.567277e-02 -0.9478806 3.435310e-01
## unemp.l2 6.281845e-02 4.813709e-02 1.3049907 1.923427e-01
## nav_high.l2 4.547635e-07 8.386509e-06 0.0542256 9.567715e-01
## approve.l3 -1.099289e-03 6.905911e-04 -1.5918092 1.118979e-01
## ruling.l3 2.794307e-03 1.633039e-03 1.7111082 8.752262e-02
## df_app.sociotropic_gd.l3 -1.770235e-03 1.867826e-03 -0.9477514 3.435967e-01
## df_app.liv_gd.l3 3.543316e-03 3.836069e-03 0.9236841 3.559823e-01
## cpi.l3 1.295415e-02 1.057568e-02 1.2249009 2.210417e-01
## unemp.l3 1.664042e-01 3.800134e-02 4.3789040 1.383098e-05
## nav_high.l3 6.598520e-06 5.558097e-06 1.1871905 2.355720e-01
## const 5.402786e-02 4.468873e-02 1.2089817 2.270951e-01
##
## $nav_high
## Estimate Std. Error t value Pr(>|t|)
## approve.l1 -7.58186073 4.95078986 -1.531444666 1.261301e-01
## ruling.l1 -4.52098809 11.70925415 -0.386103849 6.995420e-01
## df_app.sociotropic_gd.l1 9.72121827 13.25357356 0.733479029 4.635221e-01
## df_app.liv_gd.l1 1.84620122 30.62115333 0.060291694 9.519412e-01
## cpi.l1 -19.56041598 74.12849135 -0.263871767 7.919596e-01
## unemp.l1 0.32241638 272.27430354 0.001184160 9.990555e-01
## nav_high.l1 1.18582195 0.03926295 30.202060141 2.955783e-127
## approve.l2 15.71597370 6.25202874 2.513739836 1.217854e-02
## ruling.l2 11.21425599 13.56863673 0.826483619 4.088232e-01
## df_app.sociotropic_gd.l2 -39.00688667 15.78404564 -2.471285726 1.370967e-02
## df_app.liv_gd.l2 7.80756836 32.41119553 0.240891094 8.097130e-01
## cpi.l2 -18.05033231 111.38841253 -0.162048564 8.713163e-01
## unemp.l2 -170.19285230 342.11652612 -0.497470421 6.190201e-01
## nav_high.l2 -0.12453987 0.05960401 -2.089454704 3.704231e-02
## approve.l3 -8.59439832 4.90812032 -1.751056973 8.039234e-02
## ruling.l3 -4.56379930 11.60622022 -0.393220119 6.942816e-01
## df_app.sociotropic_gd.l3 29.20914536 13.27488381 2.200331526 2.812293e-02
## df_app.liv_gd.l3 -16.66025367 27.26343815 -0.611084104 5.413505e-01
## cpi.l3 41.70594317 75.16269191 0.554875592 5.791644e-01
## unemp.l3 121.11684406 270.08043795 0.448447303 6.539750e-01
## nav_high.l3 -0.07348354 0.03950211 -1.860243331 6.328777e-02
## const 2.73970265 317.60861030 0.008626034 9.931201e-01②グランジャー因果性の確認：以下いずれの場合にも、グランジャーの意味での因果性が確認されるという結果
##与党支持率の場合
df_var_rul <- na.omit(data.frame(approve, ruling))
#var modelの次数の設定
VARselect(df_var_rul)## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 2 2 3
##
## $criteria
## 1 2 3 4 5 6
## AIC(n) 5.460657 5.403611 5.394391 5.399361 5.409761 5.410420
## HQ(n) 5.475917 5.429043 5.429997 5.445140 5.465713 5.476546
## SC(n) 5.500107 5.469360 5.486440 5.517709 5.554409 5.581368
## FPE(n) 235.252047 222.207386 220.168435 221.265582 223.579407 223.727592
## 7 8 9 10
## AIC(n) 5.420664 5.430153 5.432287 5.439652
## HQ(n) 5.496962 5.516625 5.528932 5.546470
## SC(n) 5.617911 5.653700 5.682134 5.715799
## FPE(n) 226.032176 228.188786 228.677983 230.370752#VARモデルの推定
result_var_rul <- VAR(df_var_rul, p = 3, type = "const")
#グランジャー因果性の検定
causality(result_var_rul, cause = "ruling")## $Granger
##
## Granger causality H0: ruling do not Granger-cause approve
##
## data: VAR object result_var_rul
## F-Test = 0.3506, df1 = 3, df2 = 1380, p-value = 0.7887
##
##
## $Instant
##
## H0: No instantaneous causality between: ruling and approve
##
## data: VAR object result_var_rul
## Chi-squared = 121.77, df = 1, p-value < 2.2e-16##社会志向の経済評価の場合
df_var_soc <- na.omit(data.frame(approve, df_app$sociotropic_gd))
#var modelの次数の設定
VARselect(df_var_soc)## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 1 1 1 1
##
## $criteria
## 1 2 3 4 5 6
## AIC(n) 5.311681 5.312279 5.317035 5.320288 5.329999 5.341176
## HQ(n) 5.326995 5.337803 5.352769 5.366232 5.386153 5.407539
## SC(n) 5.351265 5.378252 5.409397 5.439039 5.475139 5.512705
## FPE(n) 202.690654 202.811964 203.778973 204.443296 206.438898 208.759907
## 7 8 9 10
## AIC(n) 5.349969 5.357882 5.363247 5.372771
## HQ(n) 5.426542 5.444665 5.460240 5.479973
## SC(n) 5.547887 5.582190 5.613944 5.649857
## FPE(n) 210.604595 212.279225 213.422798 215.467200#VARモデルの推定
result_var_soc <- VAR(df_var_soc, p = 1, type = "const")
#グランジャー因果性の検定
causality(result_var_soc, cause = "df_app.sociotropic_gd")## $Granger
##
## Granger causality H0: df_app.sociotropic_gd do not Granger-cause
## approve
##
## data: VAR object result_var_soc
## F-Test = 1.1485, df1 = 1, df2 = 1386, p-value = 0.2841
##
##
## $Instant
##
## H0: No instantaneous causality between: df_app.sociotropic_gd and
## approve
##
## data: VAR object result_var_soc
## Chi-squared = 7.3533, df = 1, p-value = 0.006694③インパルス・レスポンス関数(Impulse Response Function: IRF)の図示
#レスポンスが内閣支持率、インパルスが与党支持率の場合のIRF
irf_rul <- irf(result_var, response = "approve", impulse = "ruling")
plot(irf_rul)
#レスポンスが内閣支持率、インパルスが社会志向の経済評価の場合のIRF
irf_socio <- irf(result_var, response = "approve", impulse = "df_app.sociotropic_gd")
plot(irf_socio)
#レスポンスが内閣支持率、インパルスが個人志向の経済評価の場合のIRF
irf_ego <- irf(result_var, response = "approve", impulse = "df_app.liv_gd")
plot(irf_ego)
時系列データの推定方法の例②：誤差修正モデル
- 誤差修正モデルは以下のように表される；
\[\Delta y_{t}=(\alpha_{1}-1)(y_{t-1}-\frac{\alpha_{0}}{1-\alpha_{1}}-\frac{\beta_{0}+\beta_{1}}{1-\alpha_{1}}x_{t-1})+\beta_{0}\Delta x_{t}+\epsilon_{t}\]
ここで、\(y_{t-1}-\frac{\alpha_{0}}{1-\alpha_{1}}-\frac{\beta_{0}+\beta_{1}}{1-\alpha_1}x_{t-1}\)が誤差修正項(error correction term)。
- 誤差修正項をめぐる評価：
\(\alpha_{1}-1<0\)：各変数は共和分の関係にある
\(\alpha_{1}-1=0\)：誤差修正はない
\(\alpha_{1}-1>0\)：不均衡（disequilibrium）
①誤差修正モデルの推定
#誤差修正モデルの設定
library(ecm)
#全ての予測変数を長期均衡(equilibrium)としても、一過性の効果としても設定(transient effect)
df_ecm <- na.omit(data.frame(df_app$approve, df_app$sociotropic_gd, df_app$liv_gd, df_app$ruling, df_app$unemp, df_app$cpi, df_app$nav_high))
xeq <- xtr <- df_ecm[c('df_app.sociotropic_gd', 'df_app.liv_gd', "df_app.ruling", 'df_app.unemp', "df_app.cpi", "df_app.nav_high")]
#誤差修正モデルの推定
model_ecm <- ecm(df_ecm$df_app.approve, xeq, xtr, includeIntercept=TRUE)
summary(model_ecm)##
## Call:
## lm(formula = dy ~ ., data = x, weights = weights)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.012 -2.720 -0.354 1.937 52.325
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.829e+00 2.337e+00 -2.066 0.0392 *
## deltadf_app.sociotropic_gd 1.351e-01 1.010e-01 1.337 0.1815
## deltadf_app.liv_gd 3.285e-01 2.117e-01 1.552 0.1211
## deltadf_app.ruling 1.008e+00 7.978e-02 12.641 < 2e-16 ***
## deltadf_app.unemp -3.598e+00 2.071e+00 -1.738 0.0827 .
## deltadf_app.cpi -2.787e-01 5.739e-01 -0.486 0.6274
## deltadf_app.nav_high -4.526e-05 2.999e-04 -0.151 0.8801
## df_app.sociotropic_gdLag1 7.284e-02 6.130e-02 1.188 0.2352
## df_app.liv_gdLag1 1.280e-01 1.647e-01 0.777 0.4374
## df_app.rulingLag1 2.425e-01 5.219e-02 4.648 4.03e-06 ***
## df_app.unempLag1 6.016e-01 3.544e-01 1.697 0.0901 .
## df_app.cpiLag1 4.156e-02 2.660e-02 1.563 0.1186
## df_app.nav_highLag1 -4.128e-05 5.155e-05 -0.801 0.4236
## yLag1 -2.155e-01 2.364e-02 -9.118 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.537 on 682 degrees of freedom
## Multiple R-squared: 0.288, Adjusted R-squared: 0.2744
## F-statistic: 21.22 on 13 and 682 DF, p-value: < 2.2e-16②分析結果の解釈：内閣支持率に対する長期均衡となっている変数は何か？短期的影響をもつ変数は何か？誤差修正項の解釈は？
III. パネル・データの分析
パネル・データの分析とは何か？
- 個体（国／地域／グループ／個人など）に関する複数期間（年／月／日／時間）からなるデータ。
第1回授業のデータの再掲
パネル・データはなぜ望ましいのか？
- 選択バイアスの問題：調査時点の母集団の中でデータの収集が可能、回答が可能なサンプルは常に限られる（例：各国の中にはデータが取れない国もある／高齢者がアンケートに答えてくれない／不真面目な人はきちんと答えてくれないなどたくさん！）
⇒データの収集には観察されたサンプリングに対して、母集団を反映できない選択の偏り（selection bias）が不可避
- パネル・データの特性を考える：パネル・データを使うことの利点を固定効果モデル（fixed effect model: FE model）の説明をもとに考えていこう。
通常の最小二乗法の式は下記の通り；
\[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_{i}\)は\(\mu_{i}+\epsilon_{it}\)として書き直されていることになるが、ここで\(\mu_{i}\)は各国の歴史的文脈のようなものでデータ化が困難であったり、アンケートでも質問が困難であったりする観察不可能な異質性（heterogeneity）である。ここで説明変数と\(mu_{i}\)の間に相関があるならば、\(\alpha\)の推定値はバイアスをもつ。
しかしここで、個人間（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}\)を除いた推定量を導き、個人内での独立変数の変動によって、従属変数の変動を説明することが可能になる。
Rで実践するパネル・データ分析
パネル・データ分析の設定
研究の設計：各国の社会保障費は、何によって規定されているのか。
仮説の設定：
仮説1：高齢者人口費が高まるほど、社会保障費が増える
仮説2：小選挙区制の国ほど、社会保障費が減る
仮説3：左派政党の議席比率が高いほど、社会保障費が増える
仮説4：投票率が高いほど、社会保障費が増える
利用するデータ：Comparative Welfare State dataset
変数の設定
従属変数：社会保障費の対GDP比
カギとなる独立変数(key independent variables)： 高齢者人口比、高齢者人口比、左派議席比率、比例代表制ダミー、投票率
制御変数(control variables)：失業率、インフレ率
パネル・データの推定
①パネル・データのプロット。社会保障費の対GDP比に対する仮説に対応した要素のプロット
#データの読み込み
cws <- read.csv("cws.csv", fileEncoding = "SJIS")
attach(cws)
#データフレームの設定
df_cws <- data.frame(id, year, socx_pub, popold = po65/cws$pop, leftseat, vturn, singmemd, hunemr, cpi)
#国別のプロット：65歳以上人口比
library(lattice)
xyplot(socx_pub ~ popold | id,
group = id, data = df_cws,
type = c("p", "smooth"),
scales = "free")
#国別のプロット：左派議席比率
library(lattice)
xyplot(socx_pub ~ leftseat | id,
group = id, data = df_cws,
type = c("p", "smooth"),
scales = "free")
#国別のプロット：投票率
library(lattice)
xyplot(socx_pub ~ vturn | id,
group = id, data = df_cws,
type = c("p", "smooth"),
scales = "free")
② 固定効果モデルの推定
#データの読み込み
cws <- read.csv("cws.csv", fileEncoding = "SJIS")
attach(cws)
#データフレームの設定
df_cws <- data.frame(id, year, socx_pub, popold = po65/cws$pop, leftseat, vturn, singmemd, hunemr, cpi)
df_cws$prsys <- singmemd
df_cws$prsys[singmemd == 0] <- 1
df_cws$prsys[singmemd == 1 | singmemd == 2] <- 0
library(plm)
#パネル固定効果モデルの推定
model_FE <- plm(data= df_cws, socx_pub ~ popold + leftseat + vturn +prsys + hunemr + cpi, index=c("id", "year"), model = "within")
summary(model_FE)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = socx_pub ~ popold + leftseat + vturn + prsys +
## hunemr + cpi, data = df_cws, 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③ ランダム効果モデルの推定
#ランダム効果モデルの推定
model_RE <- plm(data= df_cws, socx_pub ~ popold + leftseat + vturn +prsys + hunemr + cpi, index=c("id", "year"), model = "random")
summary(model_RE)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = socx_pub ~ popold + leftseat + vturn + prsys +
## hunemr + cpi, data = df_cws, 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④固定効果モデルとランダム効果モデルのどちらを採用するか？：ハウスマン検定による判断
\(H_0\): ランダム効果モデルが望ましい
\(H_1\): 固定効果でモルが望ましい
##
## 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本推定モデルの場合、帰無仮説は棄却されず。ランダム効果推定量が採用される。
⑤国ごとに有意な差があるかをBreusch-Pagan Lagrange Multiplier testにより確認
##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: socx_pub ~ popold + leftseat + vturn + prsys + hunemr + cpi
## chisq = 3682.4, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects帰無仮説が棄却され、国ごとに有意な効果差があることが明らかに。ランラム効果モデルが妥当な推定量として追認