以下は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以降では必要ありませんので、以下ではオプションを指定していません。
data1 <- read.csv("ch6ex.csv")
data1## pref theft rhinc unemp police lpopd
## 1 北海道 4.42 524.3 4.6 1.98 5.48
## 2 青森 2.76 435.5 5.3 1.77 6.00
## 3 岩手 2.67 502.2 4.0 1.68 5.84
## 4 宮城 5.45 402.2 4.9 1.64 6.61
## 5 秋田 2.14 459.3 4.3 1.93 5.77
## 6 山形 3.05 547.9 3.6 1.77 5.97
## 7 福島 4.72 627.1 4.4 1.85 6.12
## 8 茨城 7.63 602.3 4.5 1.64 6.60
## 9 栃木 5.47 570.1 4.3 1.71 6.50
## 10 群馬 5.83 503.8 4.3 1.71 6.76
## 11 埼玉 7.69 596.1 4.3 1.57 7.94
## 12 千葉 7.47 558.4 4.1 1.88 7.47
## 13 東京 8.01 549.1 3.9 3.22 9.16
## 14 神奈川 5.13 501.4 3.9 1.70 8.73
## 15 新潟 4.51 515.4 3.7 1.80 6.23
## 16 富山 3.91 636.6 3.1 1.82 6.36
## 17 石川 5.10 589.2 3.4 1.71 6.72
## 18 福井 3.49 545.9 3.3 2.20 6.59
## 19 山梨 6.20 557.1 4.4 1.98 6.77
## 20 長野 3.86 570.0 3.4 1.65 6.48
## 21 岐阜 6.47 553.8 3.4 1.72 6.82
## 22 静岡 4.49 556.2 4.0 1.68 7.20
## 23 愛知 7.31 558.0 3.4 1.79 7.83
## 24 三重 6.30 495.9 3.4 1.68 6.78
## 25 滋賀 5.81 561.5 3.5 1.62 6.99
## 26 京都 6.79 490.4 4.4 2.51 7.71
## 27 大阪 11.78 488.7 5.3 2.43 8.80
## 28 兵庫 7.62 411.0 4.6 2.11 7.60
## 29 奈良 5.10 570.2 4.9 1.81 7.37
## 30 和歌山 5.34 529.5 4.5 2.27 6.76
## 31 鳥取 4.59 499.2 3.9 2.14 6.46
## 32 島根 3.51 551.4 2.9 2.17 6.28
## 33 岡山 5.50 517.0 4.1 1.83 6.76
## 34 広島 4.42 533.6 3.7 1.83 7.12
## 35 山口 3.53 577.6 4.0 2.22 6.71
## 36 徳島 4.12 522.6 5.0 2.04 6.62
## 37 香川 5.14 614.6 4.0 1.87 6.88
## 38 愛媛 6.17 494.8 4.4 1.76 6.72
## 39 高知 5.98 514.7 4.9 2.25 6.44
## 40 福岡 7.93 509.8 5.3 2.11 7.52
## 41 佐賀 4.73 566.0 4.1 2.03 6.44
## 42 長崎 2.40 446.5 4.4 2.24 6.71
## 43 熊本 4.19 494.7 4.5 1.72 6.46
## 44 大分 3.07 543.0 4.5 1.76 6.47
## 45 宮崎 4.56 460.1 4.6 1.82 6.39
## 46 鹿児島 3.57 559.7 4.7 1.83 6.21
## 47 沖縄 4.84 429.7 6.3 1.82 7.11使用するパッケージを読み込みますが、パッケージはインストールしておく必要があります。 以下のようにコマンドでインストールするか、RやRStudioからクリックでインストールすることもできます。
install.packages(“stargazer”, dependencies = TRUE)パッケージの機能が使えるように読み込みます。
library(psych)
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推定に先立って記述統計を示します。
describe(data1, fast=TRUE)## Warning in FUN(newX[, i], ...): min の引数に有限な値がありません: Inf を返します## Warning in FUN(newX[, i], ...): max の引数に有限な値がありません: -Inf を返しま
## す## vars n mean sd min max range se
## pref 1 47 NaN NA Inf -Inf -Inf NA
## theft 2 47 5.21 1.83 2.14 11.78 9.64 0.27
## rhinc 3 47 528.60 53.97 402.20 636.60 234.40 7.87
## unemp 4 47 4.22 0.66 2.90 6.30 3.40 0.10
## police 5 47 1.92 0.30 1.57 3.22 1.65 0.04
## lpopd 6 47 6.84 0.76 5.48 9.16 3.68 0.11都道府県名（pref）は文字情報なので、おかしな表示になっています。
reg1 <- lm(theft~rhinc+unemp+police, data=data1)
summary(reg1)##
## Call:
## lm(formula = theft ~ rhinc + unemp + police, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9081 -1.1305 -0.3197 1.0332 5.2612
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.210062 4.619218 -0.911 0.3672
## rhinc 0.006428 0.005594 1.149 0.2568
## unemp 0.712821 0.455155 1.566 0.1247
## police 1.567631 0.883825 1.774 0.0832 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.776 on 43 degrees of freedom
## Multiple R-squared: 0.1214, Adjusted R-squared: 0.06014
## F-statistic: 1.981 on 3 and 43 DF, p-value: 0.1311実質世帯収入の係数の推定値は0.006であり、実質世帯収入が千円増加すると窃盗犯認知件数が人口千人当たり0.006件増えることを示しています。ただし、統計的には有意ではなく、実質世帯収入が窃盗犯認知件数に影響するかどうかは明らかではありません。
失業率の係数の推定値は0.713で、失業率が１％ポイント増加すると窃盗犯認知件数が人口千人当たり0.713件増えることを示しています。ただし、統計的には有意ではなく、失業率が窃盗犯認知件数に影響するかどうかは明らかではありません。
警察官数の係数の推定値は1.568で、警察官数が人口千人当たり１人増えると窃盗犯認知件数が人口千人当たり1.568件増えることを示しています。統計的には10%水準で有意で、警察官数は窃盗犯認知件数に影響していると言えます。
決定係数は0.121で、窃盗犯認知件数の12％はこれらの変数で説明されることがわかりました。
data1$unemp2 <- data1$unemp^2
reg2 <- lm(theft~rhinc+unemp+unemp2+police, data=data1)
summary(reg2)##
## Call:
## lm(formula = theft ~ rhinc + unemp + unemp2 + police, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9332 -1.1480 -0.3432 1.0131 5.2832
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.073432 8.308676 -0.611 0.5447
## rhinc 0.006369 0.005679 1.122 0.2684
## unemp 1.138214 3.416120 0.333 0.7406
## unemp2 -0.048858 0.388773 -0.126 0.9006
## police 1.562808 0.894940 1.746 0.0881 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.796 on 42 degrees of freedom
## Multiple R-squared: 0.1218, Adjusted R-squared: 0.03812
## F-statistic: 1.456 on 4 and 42 DF, p-value: 0.2328失業率の1次項、2次項ともに有意ではなく、失業率と窃盗犯認知件数に2次関数的な関係があるとは言えないという結果となりました。
失業率は有意ではありませんでしたが、推定結果を失業率で偏微分すると、11.6%のあたりで窃盗犯認知件数が最も高くなるという、上に凸の関係性があることが示されており、失業率が上がると窃盗犯認知件数は増えるものの増えにくくなるという逓減的な関係が示唆されています。
また、失業率の2次項をモデルに入れた場合、１次項のみのモデルに比べて自由度修正済み決定係数が低下していて、２次項の追加はモデルを改善してはいないこともわかります。
reg3 <- lm(theft~rhinc+unemp+police+lpopd, data=data1)
summary(reg3)##
## Call:
## lm(formula = theft ~ rhinc + unemp + police + lpopd, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1276 -0.8993 0.1609 0.7981 2.8305
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.336461 3.480684 -3.257 0.00223 **
## rhinc 0.004935 0.004006 1.232 0.22484
## unemp 0.568340 0.326210 1.742 0.08878 .
## police -0.224118 0.689653 -0.325 0.74682
## lpopd 1.750737 0.269808 6.489 7.86e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.27 on 42 degrees of freedom
## Multiple R-squared: 0.5613, Adjusted R-squared: 0.5195
## F-statistic: 13.43 on 4 and 42 DF, p-value: 3.918e-07stargazer(reg1, reg3, type="text", keep.stat=c("n","rsq","adj.rsq"))##
## =========================================
## Dependent variable:
## ----------------------------
## theft
## (1) (2)
## -----------------------------------------
## rhinc 0.006 0.005
## (0.006) (0.004)
##
## unemp 0.713 0.568*
## (0.455) (0.326)
##
## police 1.568* -0.224
## (0.884) (0.690)
##
## lpopd 1.751***
## (0.270)
##
## Constant -4.210 -11.336***
## (4.619) (3.481)
##
## -----------------------------------------
## Observations 47 47
## R2 0.121 0.561
## Adjusted R2 0.060 0.519
## =========================================
## Note: *p<0.1; **p<0.05; ***p<0.01対数化した人口密度を説明変数に加えたところ、失業率の係数の推定値が正で10％水準で有意になった一方、警察官数の係数の推定値は負で有意ではなくなりました。 人口密度の係数の推定値は正で1％水準で有意で、人口密度が高い都道府県ほど、窃盗犯認知件数が高くなることが示されています。 都市部で犯罪が多いことをコントロールしたところ、失業率と警察官数が窃盗犯認知件数に与える影響が、より理解しやすいものとなりました。
attach(data1)
data2 <- data.frame(theft, rhinc, unemp, police, lpopd)
cor(data2)## theft rhinc unemp police lpopd
## theft 1.00000000 0.03384438 0.18905377 0.25953353 0.72550227
## rhinc 0.03384438 1.00000000 -0.49396219 -0.11059542 -0.01940897
## unemp 0.18905377 -0.49396219 1.00000000 0.09727683 0.08110936
## police 0.25953353 -0.11059542 0.09727683 1.00000000 0.40220564
## lpopd 0.72550227 -0.01940897 0.08110936 0.40220564 1.00000000警察官数と人口密度には一定程度の正の相関関係があり似通った変数であることがわかります。また、それぞれの窃盗犯認知件数との相関は人口密度のほうが高くなっています。結果として、関係性の強い人口密度の方が有意に、関係性の弱い方の警察官数は有意ではなくなったと考えられます。
失業率と人口密度には弱い正の相関関係があるため、失業率の係数の推定値は小さくなりましたが、人口密度によって窃盗犯発生に対する都市的な影響がコントロールされた結果、失業率と窃盗犯の関係性がクリアになったと言えます。