以下は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("ch7ex.csv")
data1
## agecl age exp sal bon size
## 1 -19 19.1 0.9 198.5 158.7 large
## 2 20-24 23.1 2.0 245.5 422.3 large
## 3 25-29 27.4 4.1 292.0 789.4 large
## 4 30-34 32.5 7.3 334.4 992.0 large
## 5 35-39 37.5 10.3 374.4 1237.6 large
## 6 40-44 42.6 13.5 401.6 1364.4 large
## 7 45-49 47.5 17.1 419.2 1497.0 large
## 8 50-54 52.4 20.5 445.2 1673.5 large
## 9 55-59 57.4 23.4 440.2 1628.2 large
## 10 -19 19.0 0.9 189.4 121.7 medium
## 11 20-24 23.0 2.1 225.4 395.0 medium
## 12 25-29 27.5 4.3 263.9 680.6 medium
## 13 30-34 32.5 7.0 294.0 788.3 medium
## 14 35-39 37.5 9.5 325.3 895.3 medium
## 15 40-44 42.5 12.0 348.9 1025.6 medium
## 16 45-49 47.5 15.2 365.6 1128.7 medium
## 17 50-54 52.3 17.2 382.9 1173.9 medium
## 18 55-59 57.4 19.2 384.1 1157.5 medium
## 19 -19 19.1 0.9 186.1 94.2 small
## 20 20-24 22.9 2.2 214.0 285.9 small
## 21 25-29 27.5 4.0 242.2 451.6 small
## 22 30-34 32.5 6.0 268.6 543.7 small
## 23 35-39 37.5 8.4 293.4 616.2 small
## 24 40-44 42.6 10.8 312.4 678.0 small
## 25 45-49 47.5 12.6 327.3 713.2 small
## 26 50-54 52.4 14.1 331.4 703.9 small
## 27 55-59 57.4 15.9 330.2 714.2 small
使用するパッケージを読み込みますが、パッケージはインストールしておく必要があります。 以下のようにコマンドでインストールするか、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
年収(ainc)と企業規模ダミー(large, med, small)を作成し、表示します。
data1$ainc <- (data1$sal*12+data1$bon)/10
data1$large <- as.integer(data1$size=="large")
data1$med <- as.integer(data1$size=="medium")
data1$small <- as.integer(data1$size=="small")
data1
## agecl age exp sal bon size ainc large med small
## 1 -19 19.1 0.9 198.5 158.7 large 254.07 1 0 0
## 2 20-24 23.1 2.0 245.5 422.3 large 336.83 1 0 0
## 3 25-29 27.4 4.1 292.0 789.4 large 429.34 1 0 0
## 4 30-34 32.5 7.3 334.4 992.0 large 500.48 1 0 0
## 5 35-39 37.5 10.3 374.4 1237.6 large 573.04 1 0 0
## 6 40-44 42.6 13.5 401.6 1364.4 large 618.36 1 0 0
## 7 45-49 47.5 17.1 419.2 1497.0 large 652.74 1 0 0
## 8 50-54 52.4 20.5 445.2 1673.5 large 701.59 1 0 0
## 9 55-59 57.4 23.4 440.2 1628.2 large 691.06 1 0 0
## 10 -19 19.0 0.9 189.4 121.7 medium 239.45 0 1 0
## 11 20-24 23.0 2.1 225.4 395.0 medium 309.98 0 1 0
## 12 25-29 27.5 4.3 263.9 680.6 medium 384.74 0 1 0
## 13 30-34 32.5 7.0 294.0 788.3 medium 431.63 0 1 0
## 14 35-39 37.5 9.5 325.3 895.3 medium 479.89 0 1 0
## 15 40-44 42.5 12.0 348.9 1025.6 medium 521.24 0 1 0
## 16 45-49 47.5 15.2 365.6 1128.7 medium 551.59 0 1 0
## 17 50-54 52.3 17.2 382.9 1173.9 medium 576.87 0 1 0
## 18 55-59 57.4 19.2 384.1 1157.5 medium 576.67 0 1 0
## 19 -19 19.1 0.9 186.1 94.2 small 232.74 0 0 1
## 20 20-24 22.9 2.2 214.0 285.9 small 285.39 0 0 1
## 21 25-29 27.5 4.0 242.2 451.6 small 335.80 0 0 1
## 22 30-34 32.5 6.0 268.6 543.7 small 376.69 0 0 1
## 23 35-39 37.5 8.4 293.4 616.2 small 413.70 0 0 1
## 24 40-44 42.6 10.8 312.4 678.0 small 442.68 0 0 1
## 25 45-49 47.5 12.6 327.3 713.2 small 464.08 0 0 1
## 26 50-54 52.4 14.1 331.4 703.9 small 468.07 0 0 1
## 27 55-59 57.4 15.9 330.2 714.2 small 467.66 0 0 1
企業規模がわかるように、勤続年数と年収の散布図を作成します。
ggplot(data1,aes(x=exp,y=ainc, shape=size))+geom_point()+xlab("勤続年数")+ylab("年収(万円)")+theme_classic()
reg1 <-lm(ainc~exp+large+med, data=data1)
summary(reg1)
##
## Call:
## lm(formula = ainc ~ exp + large + med, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -96.830 -22.567 6.202 27.112 56.926
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 241.15 16.41 14.692 3.53e-13 ***
## exp 17.58 1.18 14.894 2.65e-13 ***
## large 93.93 18.87 4.978 4.92e-05 ***
## med 40.62 18.67 2.175 0.0401 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 39.46 on 23 degrees of freedom
## Multiple R-squared: 0.924, Adjusted R-squared: 0.9141
## F-statistic: 93.19 on 3 and 23 DF, p-value: 5.142e-13
勤続年数の係数の推定値は17.6となっており、勤続年数が1年増加すると、年収が17.6万円増加することがわかります。また、推定値は1%水準で有意であり(0.1%水準でも有意です)、勤続年数は年収に影響していると言えます。
大企業ダミーの係数の推定値は93.9で、1%水準で有意となっています。したがって大企業では小企業に比べて年収が93.9万円高く、統計的に差があると言えます。
中企業ダミーの係数の推定値は40.6で、5%水準で有意となっています。したがって中企業では小企業に比べて年収が40.6万円高く、統計的に差があると言えます。
決定係数は0.924で、年収の92.4%は勤続年数と企業規模によって説明されることがわかります。
7-3の結果を使って予測値を計算し、その散布図を作成します。
data1$painc1 <- predict(reg1)
ggplot(data1,aes(x=exp,y=painc1, shape=size))+geom_point()+geom_line()+xlab("勤続年数")+ylab("年収(万円)")+theme_classic()
reg2 <-lm(ainc~exp*large+exp*med, data=data1)
summary(reg2)
##
## Call:
## lm(formula = ainc ~ exp * large + exp * med, data = data1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -86.55 -16.76 10.92 26.80 57.65
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 259.603 25.800 10.062 1.74e-09 ***
## exp 15.359 2.647 5.803 9.25e-06 ***
## large 64.288 34.874 1.843 0.0794 .
## med 22.904 35.766 0.640 0.5289
## exp:large 3.233 3.171 1.020 0.3195
## exp:med 2.141 3.405 0.629 0.5363
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.3 on 21 degrees of freedom
## Multiple R-squared: 0.9276, Adjusted R-squared: 0.9103
## F-statistic: 53.79 on 5 and 21 DF, p-value: 2.919e-11
係数ダミーの推定値を見ると、大企業は3.2、中企業は2.1と小企業に比べて、勤続年数の限界効果が大きいことが示されていますが、ともに有意ではありません。したがって、小企業と比べた大企業と中企業の勤続年数の限界効果に、統計的な差があるとは言えないという結果になっています。
定数項ダミーについては、大企業ダミーは10%水準で有意で、小企業に比べて64.3万円高いと言えますが、中企業ダミーは有意ではなく、小企業に比べて統計的に高いとは言えないという結果となっています。 なお、この回帰モデルの定数項ダミーの推定値は、勤続年数が0年の場合の年収の差であり、小企業と中企業では、勤続開始時点の年収に統計的な差があるとは言えないことが示されています。
決定係数は0.928で、年収の92.8%はこれらの変数で説明されることがわかります。係数ダミーを導入しましたが、説明力はそれほど上昇しませんでした。
7-5の結果を使って予測値を計算し、その散布図を作成します。
data1$painc2 <- predict(reg2)
ggplot(data1,aes(x=exp,y=painc2, shape=size))+geom_point()+geom_line()+xlab("勤続年数")+ylab("年収(万円)")+theme_classic()
これまでの結果をまとめて示します。
stargazer(reg1,reg2,type="text",star.cutoffs=c(0.1,0.05,0.01),keep.stat = c("n","rsq","adj.rsq"))
##
## =========================================
## Dependent variable:
## ----------------------------
## ainc
## (1) (2)
## -----------------------------------------
## exp 17.576*** 15.359***
## (1.180) (2.647)
##
## large 93.929*** 64.288*
## (18.868) (34.874)
##
## med 40.617** 22.904
## (18.671) (35.766)
##
## exp:large 3.233
## (3.171)
##
## exp:med 2.141
## (3.405)
##
## Constant 241.153*** 259.603***
## (16.414) (25.800)
##
## -----------------------------------------
## Observations 27 27
## R2 0.924 0.928
## Adjusted R2 0.914 0.910
## =========================================
## Note: *p<0.1; **p<0.05; ***p<0.01