第８章 練習問題
以下はWindows10、R version 4.2.1で実行しています。
データを読み込むにあたって、作業場所の指定とそこにデータが置いてあることが前提になります。 R本体を操作している場合は、「ファイル」→「ディレクトリの変更」でデータの置いてある場所を指定するのが簡単です。 RStudioを操作している場合は、例えばデスクトップの「R」という名前のフォルダにデータがあるとすると
setwd("C:/Users/ユーザー名/Desktop/R")を最初に実行するのがよいでしょう。ユーザー名は各自異なるので注意です。
使用するデータを読み込み、表示します。なお、R version 4.1以前で読み込む際には
data1 <- read.csv("ファイル名.csv", fileEncoding = "UTF-8-BOM")のように、エンコードのオプションを指定する必要があります。version 4.2以降では必要ありませんので、以下ではオプションを指定していません。
ワイド形式
data2 <- read.csv("ch8ex_wide.csv")
data2## pref flp00 tfr00 flp10 tfr10
## 1 北海道 57.4 1.23 62.7 1.26
## 2 青森 61.2 1.47 65.7 1.38
## 3 岩手 63.4 1.56 66.1 1.46
## 4 宮城 59.2 1.39 63.5 1.30
## 5 秋田 64.6 1.45 67.5 1.31
## 6 山形 68.2 1.62 69.3 1.48
## 7 福島 62.4 1.65 65.0 1.52
## 8 茨城 56.8 1.47 62.5 1.44
## 9 栃木 59.5 1.48 63.6 1.44
## 10 群馬 59.2 1.51 63.8 1.46
## 11 埼玉 54.9 1.30 61.3 1.32
## 12 千葉 55.3 1.30 61.4 1.34
## 13 東京 57.6 1.07 64.3 1.12
## 14 神奈川 54.2 1.28 60.3 1.31
## 15 新潟 65.4 1.51 68.0 1.43
## 16 富山 66.7 1.45 69.3 1.42
## 17 石川 65.5 1.45 68.9 1.44
## 18 福井 66.4 1.60 69.3 1.61
## 19 山梨 58.9 1.51 62.5 1.46
## 20 長野 62.8 1.59 64.3 1.53
## 21 岐阜 60.7 1.47 64.2 1.48
## 22 静岡 61.1 1.47 64.4 1.54
## 23 愛知 58.2 1.44 62.6 1.52
## 24 三重 59.4 1.48 63.9 1.51
## 25 滋賀 57.1 1.53 62.0 1.54
## 26 京都 55.3 1.28 62.3 1.28
## 27 大阪 53.6 1.31 60.9 1.33
## 28 兵庫 53.4 1.38 59.9 1.41
## 29 奈良 50.3 1.30 58.0 1.29
## 30 和歌山 54.7 1.45 61.5 1.47
## 31 鳥取 65.6 1.62 68.9 1.54
## 32 島根 65.5 1.65 69.4 1.68
## 33 岡山 58.7 1.51 63.5 1.50
## 34 広島 58.2 1.41 62.8 1.55
## 35 山口 58.6 1.47 63.0 1.56
## 36 徳島 58.0 1.45 64.3 1.42
## 37 香川 60.7 1.53 65.0 1.57
## 38 愛媛 57.8 1.45 62.9 1.50
## 39 高知 63.6 1.45 67.3 1.42
## 40 福岡 57.6 1.36 63.2 1.44
## 41 佐賀 62.7 1.67 66.7 1.61
## 42 長崎 59.6 1.57 64.2 1.61
## 43 熊本 62.4 1.56 66.3 1.62
## 44 大分 59.6 1.51 63.9 1.56
## 45 宮崎 62.9 1.62 66.7 1.68
## 46 鹿児島 58.7 1.58 64.1 1.62
## 47 沖縄 56.5 1.82 63.4 1.87ロング形式
data1 <- read.csv("ch8ex_long.csv")
data1## pref flp tfr year
## 1 北海道 57.4 1.23 2000
## 2 青森 61.2 1.47 2000
## 3 岩手 63.4 1.56 2000
## 4 宮城 59.2 1.39 2000
## 5 秋田 64.6 1.45 2000
## 6 山形 68.2 1.62 2000
## 7 福島 62.4 1.65 2000
## 8 茨城 56.8 1.47 2000
## 9 栃木 59.5 1.48 2000
## 10 群馬 59.2 1.51 2000
## 11 埼玉 54.9 1.30 2000
## 12 千葉 55.3 1.30 2000
## 13 東京 57.6 1.07 2000
## 14 神奈川 54.2 1.28 2000
## 15 新潟 65.4 1.51 2000
## 16 富山 66.7 1.45 2000
## 17 石川 65.5 1.45 2000
## 18 福井 66.4 1.60 2000
## 19 山梨 58.9 1.51 2000
## 20 長野 62.8 1.59 2000
## 21 岐阜 60.7 1.47 2000
## 22 静岡 61.1 1.47 2000
## 23 愛知 58.2 1.44 2000
## 24 三重 59.4 1.48 2000
## 25 滋賀 57.1 1.53 2000
## 26 京都 55.3 1.28 2000
## 27 大阪 53.6 1.31 2000
## 28 兵庫 53.4 1.38 2000
## 29 奈良 50.3 1.30 2000
## 30 和歌山 54.7 1.45 2000
## 31 鳥取 65.6 1.62 2000
## 32 島根 65.5 1.65 2000
## 33 岡山 58.7 1.51 2000
## 34 広島 58.2 1.41 2000
## 35 山口 58.6 1.47 2000
## 36 徳島 58.0 1.45 2000
## 37 香川 60.7 1.53 2000
## 38 愛媛 57.8 1.45 2000
## 39 高知 63.6 1.45 2000
## 40 福岡 57.6 1.36 2000
## 41 佐賀 62.7 1.67 2000
## 42 長崎 59.6 1.57 2000
## 43 熊本 62.4 1.56 2000
## 44 大分 59.6 1.51 2000
## 45 宮崎 62.9 1.62 2000
## 46 鹿児島 58.7 1.58 2000
## 47 沖縄 56.5 1.82 2000
## 48 北海道 62.7 1.26 2010
## 49 青森 65.7 1.38 2010
## 50 岩手 66.1 1.46 2010
## 51 宮城 63.5 1.30 2010
## 52 秋田 67.5 1.31 2010
## 53 山形 69.3 1.48 2010
## 54 福島 65.0 1.52 2010
## 55 茨城 62.5 1.44 2010
## 56 栃木 63.6 1.44 2010
## 57 群馬 63.8 1.46 2010
## 58 埼玉 61.3 1.32 2010
## 59 千葉 61.4 1.34 2010
## 60 東京 64.3 1.12 2010
## 61 神奈川 60.3 1.31 2010
## 62 新潟 68.0 1.43 2010
## 63 富山 69.3 1.42 2010
## 64 石川 68.9 1.44 2010
## 65 福井 69.3 1.61 2010
## 66 山梨 62.5 1.46 2010
## 67 長野 64.3 1.53 2010
## 68 岐阜 64.2 1.48 2010
## 69 静岡 64.4 1.54 2010
## 70 愛知 62.6 1.52 2010
## 71 三重 63.9 1.51 2010
## 72 滋賀 62.0 1.54 2010
## 73 京都 62.3 1.28 2010
## 74 大阪 60.9 1.33 2010
## 75 兵庫 59.9 1.41 2010
## 76 奈良 58.0 1.29 2010
## 77 和歌山 61.5 1.47 2010
## 78 鳥取 68.9 1.54 2010
## 79 島根 69.4 1.68 2010
## 80 岡山 63.5 1.50 2010
## 81 広島 62.8 1.55 2010
## 82 山口 63.0 1.56 2010
## 83 徳島 64.3 1.42 2010
## 84 香川 65.0 1.57 2010
## 85 愛媛 62.9 1.50 2010
## 86 高知 67.3 1.42 2010
## 87 福岡 63.2 1.44 2010
## 88 佐賀 66.7 1.61 2010
## 89 長崎 64.2 1.61 2010
## 90 熊本 66.3 1.62 2010
## 91 大分 63.9 1.56 2010
## 92 宮崎 66.7 1.68 2010
## 93 鹿児島 64.1 1.62 2010
## 94 沖縄 63.4 1.87 2010使用するパッケージを読み込みますが、パッケージはインストールしておく必要があります。 以下のようにコマンドでインストールするか、RやRStudioからクリックでインストールすることもできます。
install.packages(“stargazer”, dependencies = TRUE)パッケージの機能が使えるように読み込みます。
library(ggplot2)
library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer８－１
データ年の違いがわかるように散布図を作成します。
data1$year <- as.factor(data1$year)
ggplot(data1,aes(x=flp,y=tfr, shape=year))+geom_point()+xlab("女性労働力率(%)")+ylab("合計特殊出生率")+labs(shape="データ年")+
theme_classic()+scale_shape_manual(values = c(16,2))
2000年、2010年ともに正の相関があるように見えます。 2時点まとめた場合でも正の相関になりそうです。
８－２
年ごと、そして２年まとめて推定します。
reg00 <- lm(tfr~flp, data=data1, subset = year == 2000)
summary(reg00)##
## Call:
## lm(formula = tfr ~ flp, data = data1, subset = year == 2000)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36307 -0.06874 0.01230 0.05237 0.40698
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.383186 0.247764 1.547 0.129
## flp 0.018227 0.004135 4.408 6.4e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1125 on 45 degrees of freedom
## Multiple R-squared: 0.3016, Adjusted R-squared: 0.2861
## F-statistic: 19.43 on 1 and 45 DF, p-value: 6.403e-05reg10 <- lm(tfr~flp, data=data1, subset = year == 2010)
summary(reg10)##
## Call:
## lm(formula = tfr ~ flp, data = data1, subset = year == 2010)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35044 -0.09224 0.00877 0.08243 0.41372
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.458973 0.442843 1.036 0.3055
## flp 0.015730 0.006875 2.288 0.0269 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.127 on 45 degrees of freedom
## Multiple R-squared: 0.1042, Adjusted R-squared: 0.0843
## F-statistic: 5.235 on 1 and 45 DF, p-value: 0.02689data1$y2010 <- as.integer(data1$year==2010)
reg_all1 <- lm(tfr~flp+y2010, data=data1)
summary(reg_all1)##
## Call:
## lm(formula = tfr ~ flp + y2010, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36480 -0.07922 0.01173 0.06330 0.41535
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.43029 0.21774 1.976 0.05116 .
## flp 0.01744 0.00363 4.804 6.09e-06 ***
## y2010 -0.08129 0.02968 -2.739 0.00741 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1194 on 91 degrees of freedom
## Multiple R-squared: 0.2023, Adjusted R-squared: 0.1848
## F-statistic: 11.54 on 2 and 91 DF, p-value: 3.409e-05reg_all2 <- lm(tfr~flp*y2010, data=data1)
summary(reg_all2)##
## Call:
## lm(formula = tfr ~ flp * y2010, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36307 -0.08323 0.01038 0.06184 0.41372
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.383186 0.264230 1.450 0.150
## flp 0.018227 0.004410 4.133 8e-05 ***
## y2010 0.075786 0.494771 0.153 0.879
## flp:y2010 -0.002497 0.007850 -0.318 0.751
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.12 on 90 degrees of freedom
## Multiple R-squared: 0.2032, Adjusted R-squared: 0.1767
## F-statistic: 7.652 on 3 and 90 DF, p-value: 0.00013stargazer(reg00, reg10, reg_all1, reg_all2, type="text", digits = 4, omit.stat = c("ser","f"))##
## ====================================================
## Dependent variable:
## ---------------------------------------
## tfr
## (1) (2) (3) (4)
## ----------------------------------------------------
## flp 0.0182*** 0.0157** 0.0174*** 0.0182***
## (0.0041) (0.0069) (0.0036) (0.0044)
##
## y2010 -0.0813*** 0.0758
## (0.0297) (0.4948)
##
## flp:y2010 -0.0025
## (0.0079)
##
## Constant 0.3832 0.4590 0.4303* 0.3832
## (0.2478) (0.4428) (0.2177) (0.2642)
##
## ----------------------------------------------------
## Observations 47 47 94 94
## R2 0.3016 0.1042 0.2023 0.2032
## Adjusted R2 0.2861 0.0843 0.1848 0.1767
## ====================================================
## Note: *p<0.1; **p<0.05; ***p<0.01いずれのモデルにおいても女性労働力率の係数は正で、有意になっています。 つまり、女性労働力率の上昇は合計特殊出生率を上昇させる効果があることわかります。 限界効果としては、女性労働力率が1％ポイント上昇すると0.016～0.018くらいの範囲で合計特殊出生率が増加しています。 （3）の2010年ダミーが-0.081で1％水準で有意になっています。 これは、女性労働力率に変化がなかった場合でも2000年から2010年にかけて合計特殊出生率が0.081減少していることを示しています。 また（4）の推定で、女性労働力率と2010年ダミーとの交差項の係数は有意ではありませんので、この10年間で女性労働力率の影響に変化はないようです。
８－３
差分データを作成し、散布図を作成します。
data2$dflp <- data2$flp10 - data2$flp00
data2$dtfr <- data2$tfr10 - data2$tfr00
ggplot(data2,aes(x=dflp,y=dtfr))+geom_point()+xlab("女性労働力率(%)の差分")+ylab("合計特殊出生率の差分")+
theme_classic()
差分データの散布図を見ると、やや右上がりで正の相関があるように見えます。
８－４
差分データで回帰分析します。
dreg1 <- lm(dtfr~1+dflp, data=data2)
dreg2 <- lm(dtfr~0+dflp, data=data2)
summary(dreg1)##
## Call:
## lm(formula = dtfr ~ 1 + dflp, data = data2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.105689 -0.046646 0.001749 0.041094 0.140993
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.091147 0.026261 -3.471 0.001156 **
## dflp 0.019599 0.005459 3.590 0.000813 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05694 on 45 degrees of freedom
## Multiple R-squared: 0.2227, Adjusted R-squared: 0.2054
## F-statistic: 12.89 on 1 and 45 DF, p-value: 0.0008126summary(dreg2)##
## Call:
## lm(formula = dtfr ~ 0 + dflp, data = data2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.144711 -0.051254 0.005289 0.032771 0.132527
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## dflp 0.001624 0.001922 0.845 0.402
##
## Residual standard error: 0.06341 on 46 degrees of freedom
## Multiple R-squared: 0.01528, Adjusted R-squared: -0.006122
## F-statistic: 0.714 on 1 and 46 DF, p-value: 0.4025stargazer(dreg2, dreg1, type="text", digits = 4, omit.stat = c("ser","f"))##
## =========================================
## Dependent variable:
## ----------------------------
## dtfr
## (1) (2)
## -----------------------------------------
## dflp 0.0016 0.0196***
## (0.0019) (0.0055)
##
## Constant -0.0911***
## (0.0263)
##
## -----------------------------------------
## Observations 47 47
## R2 0.0153 0.2227
## Adjusted R2 -0.0061 0.2054
## =========================================
## Note: *p<0.1; **p<0.05; ***p<0.01トレンド項があるモデル（2）では、女性労働力率の係数が正で有意になっています。 2000年から2010年にかけて、女性労働力率の増加は合計特殊出生率の増加をもたらしたことがわかります。 限界効果で見ると、女性労働力率が1％ポイント増加すると、合計特殊出生率が0.02増加していることになります。 また、トレンド項の係数が負で有意ですので、女性労働力率が変化しなくても、合計特殊出生率が趨勢的に減少していることもわかります。
付録
8-4の結果から、回帰直線を示します。
ggplot(data2,aes(x=dflp,y=dtfr))+geom_point()+xlab("女性労働力率(%)の差分")+ylab("合計特殊出生率の差分")+
theme_classic()+
geom_abline(aes(intercept = dreg1$coefficients[[1]],slope = dreg1$coefficients[[2]],linetype="line1"))+
geom_abline(aes(intercept = 0,slope = dreg2$coefficients[[1]],linetype="line2"))+
scale_linetype(name="回帰",labels=c(line1="トレンド項あり",line2="トレンド項なし"))