第5章 線形重回帰モデルの推定と検定

はじめに，不足している変数をつくります．

mydata <- read.csv("youdou.csv")

mydata$lny80 <- log(mydata$y80)
mydata$lny99 <- log(mydata$y99)
mydata$lny90 <- log(mydata$y90)
mydata$growthrate8099 <- (mydata$lny99-mydata$lny80) / 19
mydata$growthrate8090 <- (mydata$lny90-mydata$lny80) / 10

mydata$growthrate8099 <-  mydata$growthrate8099 * 100

実証例 5.1 （線形重回帰）

(pp.151-152)

線形重回帰モデルの推定もlm()関数で行えます．lm(y ~ x1 + x2 + x3, data)のようにして，複数の説明変数を含めることができます．

model1 <- lm(growthrate8099 ~ trust80 + education80 + lny80, data = mydata)

model1
## 
## Call:
## lm(formula = growthrate8099 ~ trust80 + education80 + lny80, 
##     data = mydata)
## 
## Coefficients:
## (Intercept)      trust80  education80        lny80  
##     6.04885      0.02058      2.61208     -2.38309

model2 <- lm(growthrate8099 ~ norm80 + education80 + lny80, data = mydata)

model2
## 
## Call:
## lm(formula = growthrate8099 ~ norm80 + education80 + lny80, data = mydata)
## 
## Coefficients:
## (Intercept)       norm80  education80        lny80  
##      5.2909       0.3383       4.3872      -1.9911

実証例 5.2 （FWL定理１）

(p.154)

model1
## 
## Call:
## lm(formula = growthrate8099 ~ trust80 + education80 + lny80, 
##     data = mydata)
## 
## Coefficients:
## (Intercept)      trust80  education80        lny80  
##     6.04885      0.02058      2.61208     -2.38309

reg1 <- lm(trust80 ~ education80 + lny80, data = mydata)
reg1 # 回帰
## 
## Call:
## lm(formula = trust80 ~ education80 + lny80, data = mydata)
## 
## Coefficients:
## (Intercept)  education80        lny80  
##       2.674      -11.289       -1.025
res1 <- reg1$residuals # 残差

reg2 <- lm(growthrate8099 ~ education80 + lny80, data = mydata)
reg2
## 
## Call:
## lm(formula = growthrate8099 ~ education80 + lny80, data = mydata)
## 
## Coefficients:
## (Intercept)  education80        lny80  
##       6.104        2.380       -2.404
res2 <- reg2$residuals

model3 <- lm(res2 ~ 0 + res1) # 定数項を0として推定
model3
## 
## Call:
## lm(formula = res2 ~ 0 + res1)
## 
## Coefficients:
##    res1  
## 0.02058

実証例 5.3 （FWL定理２）

(pp.155-156)

reg1 <- lm(trust80 ~ education80 + lny80, data = mydata)
reg1 # 回帰
## 
## Call:
## lm(formula = trust80 ~ education80 + lny80, data = mydata)
## 
## Coefficients:
## (Intercept)  education80        lny80  
##       2.674      -11.289       -1.025
res1 <- reg1$residuals # 残差

model4 <- lm(mydata$growthrate8099 ~ 0 + res1) # 定数項を0として推定
model4
## 
## Call:
## lm(formula = mydata$growthrate8099 ~ 0 + res1)
## 
## Coefficients:
##    res1  
## 0.02058

実証例 5.4 （標準誤差）

(p.163)

library(estimatr)

result_robust1 <- lm_robust(growthrate8099 ~ trust80 + education80 + lny80, 
                    data = mydata, se_type = "stata")
summary(result_robust1)
## 
## Call:
## lm_robust(formula = growthrate8099 ~ trust80 + education80 + 
##     lny80, data = mydata, se_type = "stata")
## 
## Standard error type:  HC1 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
## (Intercept)  6.04885    0.42643 14.1849 8.041e-18    5.189   6.9088 43
## trust80      0.02058    0.07564  0.2721 7.868e-01   -0.132   0.1731 43
## education80  2.61208    2.70857  0.9644 3.403e-01   -2.850   8.0744 43
## lny80       -2.38309    0.49147 -4.8489 1.658e-05   -3.374  -1.3920 43
## 
## Multiple R-squared:  0.5619 ,    Adjusted R-squared:  0.5313 
## F-statistic: 20.21 on 3 and 43 DF,  p-value: 2.531e-08


result_robust2 <- lm_robust(growthrate8099 ~ norm80 + education80 + lny80, 
                            data = mydata, se_type = "stata")
summary(result_robust2)
## 
## Call:
## lm_robust(formula = growthrate8099 ~ norm80 + education80 + lny80, 
##     data = mydata, se_type = "stata")
## 
## Standard error type:  HC1 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
## (Intercept)   5.2909     0.6682   7.918 6.204e-10  3.94324   6.6385 43
## norm80        0.3383     0.1370   2.469 1.758e-02  0.06202   0.6145 43
## education80   4.3872     1.9611   2.237 3.051e-02  0.43233   8.3421 43
## lny80        -1.9911     0.5746  -3.465 1.213e-03 -3.14987  -0.8324 43
## 
## Multiple R-squared:  0.6391 ,    Adjusted R-squared:  0.614 
## F-statistic: 41.04 on 3 and 43 DF,  p-value: 1.11e-12

実証例 5.5 （多項式モデル）

(pp.170-171)

2乗項は，mydata$y80sq <- mydata$y80sq^2などとして新しい変数として定義して回帰式に含めるか，または，lm()関数の中の回帰式にI(y80^2)として含めることで，その係数を推定するように指示することができます．

library(estimatr)

result_robust3 <- lm_robust(growthrate8099 ~ y80 + I(y80^2), 
                    data = mydata, se_type = "stata") # 2乗項
summary(result_robust3)
## 
## Call:
## lm_robust(formula = growthrate8099 ~ y80 + I(y80^2), data = mydata, 
##     se_type = "stata")
## 
## Standard error type:  HC1 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
## (Intercept)  6.51866    1.38538   4.705 2.535e-05  3.72662   9.3107 44
## y80         -1.22615    0.70791  -1.732 9.027e-02 -2.65285   0.2005 44
## I(y80^2)     0.08935    0.08861   1.008 3.188e-01 -0.08923   0.2679 44
## 
## Multiple R-squared:  0.5503 ,    Adjusted R-squared:  0.5299 
## F-statistic: 27.39 on 2 and 44 DF,  p-value: 1.879e-08

実証例 5.6 （交差項１）

(pp.173-174)

D * lny80としてlm()関数の回帰式に含めることで，D，lny80，そしてそれらの交差項の係数を推定するように指示できます．

mydata$D <- 1 * (mydata$did > 0.4)

model6 <- lm(growthrate8099 ~ D * lny80, data = mydata)
model6
## 
## Call:
## lm(formula = growthrate8099 ~ D * lny80, data = mydata)
## 
## Coefficients:
## (Intercept)            D        lny80      D:lny80  
##     5.74905     -0.17551     -1.91120      0.06441


lm(growthrate8099 ~ lny80, data = subset(mydata, D == 0)) # 標本を分けて推定
## 
## Call:
## lm(formula = growthrate8099 ~ lny80, data = subset(mydata, D == 
##     0))
## 
## Coefficients:
## (Intercept)        lny80  
##       5.749       -1.911
lm(growthrate8099 ~ lny80, data = subset(mydata, D == 1)) # 標本を分けて推定
## 
## Call:
## lm(formula = growthrate8099 ~ lny80, data = subset(mydata, D == 
##     1))
## 
## Coefficients:
## (Intercept)        lny80  
##       5.574       -1.847

実証例 5.7 （交差項２）

(p.175)

mydata$D1 <- mydata$D
mydata$D2 <- 1 * (mydata$lny80 > 1.4)

model7 <- lm(growthrate8099 ~ D1 * D2, data = mydata)
model7
## 
## Call:
## lm(formula = growthrate8099 ~ D1 * D2, data = mydata)
## 
## Coefficients:
## (Intercept)           D1           D2        D1:D2  
##     3.45474     -0.23329     -0.58003      0.04725

実証例 5.8 （結合仮説）

(p.178)

出力結果の一番下に\(F\)統計量があります．

library(estimatr)

model5 <- lm_robust(growthrate8099 ~ y80 + I(y80^2), se_type = "stata", data = mydata)
summary(model5)
## 
## Call:
## lm_robust(formula = growthrate8099 ~ y80 + I(y80^2), data = mydata, 
##     se_type = "stata")
## 
## Standard error type:  HC1 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper DF
## (Intercept)  6.51866    1.38538   4.705 2.535e-05  3.72662   9.3107 44
## y80         -1.22615    0.70791  -1.732 9.027e-02 -2.65285   0.2005 44
## I(y80^2)     0.08935    0.08861   1.008 3.188e-01 -0.08923   0.2679 44
## 
## Multiple R-squared:  0.5503 ,    Adjusted R-squared:  0.5299 
## F-statistic: 27.39 on 2 and 44 DF,  p-value: 1.879e-08