分析結果の解釈については、書籍をご参照ください。以下は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("ch6.csv")
data1## pref her hinc univ stu popd rcpi
## 1 北海道 43.3 525.9 0.69 12.32 240.5 100.3
## 2 青森 43.6 433.3 0.76 12.09 405.1 99.5
## 3 岩手 44.2 500.2 0.39 11.44 344.5 99.6
## 4 宮城 49.5 395.8 0.60 13.35 739.8 98.4
## 5 秋田 44.6 454.7 0.68 11.71 319.3 99.0
## 6 山形 44.8 555.6 0.53 12.10 389.6 101.4
## 7 福島 45.7 631.5 0.42 12.57 453.9 100.7
## 8 茨城 50.5 592.7 0.31 13.84 733.9 98.4
## 9 栃木 52.0 567.2 0.46 14.51 661.9 99.5
## 10 群馬 52.6 489.7 0.66 13.81 865.6 97.2
## 11 埼玉 56.9 601.5 0.41 15.90 2811.4 100.9
## 12 千葉 56.0 557.3 0.43 15.82 1750.7 99.8
## 13 東京 66.5 560.6 1.01 16.47 9528.5 102.1
## 14 神奈川 61.5 513.9 0.33 16.41 6205.8 102.5
## 15 新潟 46.2 512.8 0.78 13.45 508.1 99.5
## 16 富山 52.0 629.6 0.47 12.42 578.6 98.9
## 17 石川 54.7 596.3 1.04 13.12 829.1 101.2
## 18 福井 56.0 547.0 0.64 13.46 730.3 100.2
## 19 山梨 56.3 552.1 0.84 13.15 874.8 99.1
## 20 長野 48.9 558.6 0.43 13.13 650.7 98.0
## 21 岐阜 55.1 541.1 0.59 13.87 918.9 97.7
## 22 静岡 53.0 547.9 0.38 14.48 1345.8 98.5
## 23 愛知 58.7 550.7 0.67 16.11 2504.6 98.7
## 24 三重 50.5 490.9 0.39 13.52 881.8 99.0
## 25 滋賀 55.0 564.3 0.57 13.86 1080.8 100.5
## 26 京都 66.4 495.3 1.30 13.66 2223.8 101.0
## 27 大阪 60.5 490.7 0.62 15.71 6643.3 100.4
## 28 兵庫 60.6 415.5 0.69 14.05 1988.8 101.1
## 29 奈良 58.9 557.1 0.81 14.12 1594.7 97.7
## 30 和歌山 49.5 533.2 0.31 12.57 864.1 100.7
## 31 鳥取 43.5 495.7 0.52 11.42 636.6 99.3
## 32 島根 47.0 555.8 0.29 10.51 534.6 100.8
## 33 岡山 50.5 511.8 0.88 13.60 866.0 99.0
## 34 広島 59.9 532.5 0.70 13.72 1230.7 99.8
## 35 山口 42.7 576.4 0.71 11.90 823.0 99.8
## 36 徳島 51.7 519.5 0.53 11.68 748.1 99.4
## 37 香川 50.6 609.7 0.41 12.24 970.9 99.2
## 38 愛媛 52.2 490.8 0.36 12.60 827.9 99.2
## 39 高知 47.4 515.2 0.41 9.48 626.1 100.1
## 40 福岡 54.4 502.7 0.67 15.37 1847.4 98.6
## 41 佐賀 43.0 555.8 0.24 11.94 623.6 98.2
## 42 長崎 44.6 449.2 0.73 12.27 821.6 100.6
## 43 熊本 46.1 490.7 0.50 12.71 638.8 99.2
## 44 大分 46.4 535.4 0.43 12.12 648.4 98.6
## 45 宮崎 45.1 449.1 0.63 12.25 596.8 97.6
## 46 鹿児島 42.7 545.1 0.36 10.91 497.5 97.4
## 47 沖縄 39.2 427.6 0.56 12.91 1226.2 99.5都道府県名（pref）、進学率（her）、世帯収入（hinc）、大学数（univ）、生徒数（stu）、人口密度（popd）、消費者物価地域差指数（rcpi）となっています。
本章で使用するパッケージを読み込みますが、パッケージはインストールしておく必要があります。 以下のようにコマンドでインストールするか、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実質世帯収入（rhinc）、人口密度の対数（lpopd）、1万円単位にした実質世帯収入（rhinc10）を作成します。
data1$rhinc <- round(data1$hinc*100/data1$rcpi,1)
data1$lpopd <- round(log(data1$popd),2)
data1$rhinc10 <- data1$rhinc/10記述統計の計算に使用するために、必要な変数だけでデータフレーム（data2）をつくります。
attach(data1)
data2 <- data.frame(her, rhinc, univ, stu, lpopd)表６－１ 記述統計
describe(data2, fast=TRUE)## vars n mean sd min max range se
## her 1 47 51.09 6.59 39.20 66.50 27.30 0.96
## rhinc 2 47 528.60 53.97 402.20 636.60 234.40 7.87
## univ 3 47 0.58 0.22 0.24 1.30 1.06 0.03
## stu 4 47 13.21 1.56 9.48 16.47 6.99 0.23
## lpopd 5 47 6.84 0.76 5.48 9.16 3.68 0.11記述統計を出すために、以下のようにsummary関数を使った場合は標準偏差が出ないので、sd関数も使う必要があります。
summary(data2)## her rhinc univ stu
## Min. :39.20 Min. :402.2 Min. :0.2400 Min. : 9.48
## 1st Qu.:45.40 1st Qu.:497.6 1st Qu.:0.4100 1st Qu.:12.18
## Median :50.50 Median :533.6 Median :0.5600 Median :13.13
## Mean :51.09 Mean :528.6 Mean :0.5774 Mean :13.21
## 3rd Qu.:55.55 3rd Qu.:560.6 3rd Qu.:0.6900 3rd Qu.:13.87
## Max. :66.50 Max. :636.6 Max. :1.3000 Max. :16.47
## lpopd
## Min. :5.480
## 1st Qu.:6.440
## Median :6.710
## Mean :6.835
## 3rd Qu.:7.115
## Max. :9.160sd(data2$her)## [1] 6.593639基本の回帰モデル（6.1）式
reg1 <- lm(her~rhinc+univ+stu, data=data1)
summary(reg1)##
## Call:
## lm(formula = her ~ rhinc + univ + stu, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.2216 -2.3877 -0.6028 2.5130 8.4340
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.09741 8.00144 0.012 0.99034
## rhinc 0.01704 0.01188 1.434 0.15877
## univ 8.67383 3.03392 2.859 0.00653 **
## stu 2.79971 0.41204 6.795 2.55e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.202 on 43 degrees of freedom
## Multiple R-squared: 0.6204, Adjusted R-squared: 0.5939
## F-statistic: 23.42 on 3 and 43 DF, p-value: 3.84e-09実質世帯収入を1万円単位にした場合
reg2 <- lm(her~rhinc10+univ+stu, data=data1)
summary(reg2)##
## Call:
## lm(formula = her ~ rhinc10 + univ + stu, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.2216 -2.3877 -0.6028 2.5130 8.4340
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.09741 8.00144 0.012 0.99034
## rhinc10 0.17041 0.11883 1.434 0.15877
## univ 8.67383 3.03392 2.859 0.00653 **
## stu 2.79971 0.41204 6.795 2.55e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.202 on 43 degrees of freedom
## Multiple R-squared: 0.6204, Adjusted R-squared: 0.5939
## F-statistic: 23.42 on 3 and 43 DF, p-value: 3.84e-09実質世帯収入の2次項を作成し、表6-2をつくるため、使用する変数だけでデータフレーム（data3）をつくります。
data1$rhinc2 <- data1$rhinc^2
attach(data1)## 以下のオブジェクトは data1 (pos = 3) からマスクされています:
##
## her, hinc, lpopd, popd, pref, rcpi, rhinc, rhinc10, stu, univdata3 <- data.frame(her, rhinc, rhinc2, univ, stu)表６－２ 実質世帯収入を2乗した変数の追加
head(data3)## her rhinc rhinc2 univ stu
## 1 43.3 524.3 274890.5 0.69 12.32
## 2 43.6 435.5 189660.2 0.76 12.09
## 3 44.2 502.2 252204.8 0.39 11.44
## 4 49.5 402.2 161764.8 0.60 13.35
## 5 44.6 459.3 210956.5 0.68 11.71
## 6 44.8 547.9 300194.4 0.53 12.10実質世帯収入の2次項（rhinc2）を使った回帰モデル（6.2）式
reg3 <- lm(her~rhinc+rhinc2+univ+stu, data=data1)
summary(reg3)##
## Call:
## lm(formula = her ~ rhinc + rhinc2 + univ + stu, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.6739 -2.5038 -0.7806 2.5165 9.1531
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.920e+01 4.443e+01 -0.657 0.51458
## rhinc 1.319e-01 1.717e-01 0.768 0.44665
## rhinc2 -1.107e-04 1.651e-04 -0.671 0.50615
## univ 8.593e+00 3.056e+00 2.812 0.00746 **
## stu 2.791e+00 4.149e-01 6.728 3.57e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.229 on 42 degrees of freedom
## Multiple R-squared: 0.6244, Adjusted R-squared: 0.5886
## F-statistic: 17.46 on 4 and 42 DF, p-value: 1.655e-08対数化した実質世帯収入(lrhinc)を使った回帰モデル（6.3）式
data1$lrhinc <- log(data1$rhinc)
reg3a <- lm(her~lrhinc+univ+stu, data=data1)
summary(reg3a)##
## Call:
## lm(formula = her ~ lrhinc + univ + stu, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.095 -2.410 -0.637 2.480 8.391
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -47.1372 38.1492 -1.236 0.22331
## lrhinc 8.9818 6.0786 1.478 0.14680
## univ 8.6888 3.0270 2.870 0.00634 **
## stu 2.7970 0.4115 6.797 2.53e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.196 on 43 degrees of freedom
## Multiple R-squared: 0.6214, Adjusted R-squared: 0.595
## F-statistic: 23.53 on 3 and 43 DF, p-value: 3.618e-09対数化した人口密度を加えた回帰モデル
reg4 <- lm(her~rhinc+univ+stu+lpopd, data=data1)
summary(reg4)##
## Call:
## lm(formula = her ~ rhinc + univ + stu + lpopd, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6631 -2.0735 -0.0494 1.9474 6.6751
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.61564 7.20729 -1.473 0.148232
## rhinc 0.02196 0.01010 2.174 0.035428 *
## univ 7.99297 2.56730 3.113 0.003325 **
## stu 0.89445 0.56537 1.582 0.121141
## lpopd 4.92562 1.15196 4.276 0.000107 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.549 on 42 degrees of freedom
## Multiple R-squared: 0.7355, Adjusted R-squared: 0.7103
## F-statistic: 29.2 on 4 and 42 DF, p-value: 1.221e-11表6－3 推定に使用した変数の相関行列（N=47）
round(cor(data2), digit=3)## her rhinc univ stu lpopd
## her 1.000 0.126 0.404 0.737 0.794
## rhinc 0.126 1.000 -0.223 0.076 -0.019
## univ 0.404 -0.223 1.000 0.228 0.236
## stu 0.737 0.076 0.228 1.000 0.797
## lpopd 0.794 -0.019 0.236 0.797 1.000ここでは、round関数で小数点以下第3位までの表示にしています。
付録
ここまでの結果をまとめています。
stargazer(reg1, reg2, reg3, reg4, type="text", keep.stat=c("n", "rsq", "adj.rsq"))##
## ================================================
## Dependent variable:
## -----------------------------------
## her
## (1) (2) (3) (4)
## ------------------------------------------------
## rhinc 0.017 0.132 0.022**
## (0.012) (0.172) (0.010)
##
## rhinc10 0.170
## (0.119)
##
## rhinc2 -0.0001
## (0.0002)
##
## univ 8.674*** 8.674*** 8.593*** 7.993***
## (3.034) (3.034) (3.056) (2.567)
##
## stu 2.800*** 2.800*** 2.791*** 0.894
## (0.412) (0.412) (0.415) (0.565)
##
## lpopd 4.926***
## (1.152)
##
## Constant 0.097 0.097 -29.205 -10.616
## (8.001) (8.001) (44.432) (7.207)
##
## ------------------------------------------------
## Observations 47 47 47 47
## R2 0.620 0.620 0.624 0.736
## Adjusted R2 0.594 0.594 0.589 0.710
## ================================================
## Note: *p<0.1; **p<0.05; ***p<0.01