第4章 線形単回帰モデルの推定と検定
実証例 4.1 （線形単回帰）
(pp.128-129)
lm()関数を用いれば簡単に線形回帰モデルのパラメータを推定することができます．回帰モデルの推定に用いる関数は多くの場合，はじめの引数formulaで回帰式を指定します．たとえば，\(y = \beta_0 + \beta_1 x + \varepsilon\)というモデルなら，formula = y ~ xのように式を表現します．その次にdataという引数で，どのデータフレームを用いているかを指定します．
mydata <- read.csv("ch04_wage.csv")
result <- lm(formula = wage ~ productivity, data = mydata)
result
##
## Call:
## lm(formula = wage ~ productivity, data = mydata)
##
## Coefficients:
## (Intercept) productivity
## 276.1296 0.5468さらに，summary()関数を用いることで簡単に\(t\)検定を行うことができます．
summary(result)
##
## Call:
## lm(formula = wage ~ productivity, data = mydata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.618 -17.612 4.186 21.946 37.052
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 276.12961 87.61057 3.152 0.00525 **
## productivity 0.54682 0.02442 22.395 4.04e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.77 on 19 degrees of freedom
## Multiple R-squared: 0.9635, Adjusted R-squared: 0.9616
## F-statistic: 501.5 on 1 and 19 DF, p-value: 4.037e-15教科書では不均一分散に頑健な標準誤差を用いています．estimatrパッケージを用いれば簡単に不均一分散に頑健な標準誤差を用いた検定が行えます．formula =という箇所は毎回書くのが面倒ならば，省略することができます．
library(estimatr)
result_robust <- lm_robust(wage ~ productivity, data = mydata, se_type = "stata")
summary(result_robust)
##
## Call:
## lm_robust(formula = wage ~ productivity, data = mydata, se_type = "stata")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 276.1296 71.25559 3.875 1.019e-03 126.990 425.2693 19
## productivity 0.5468 0.02046 26.722 1.553e-16 0.504 0.5896 19
##
## Multiple R-squared: 0.9635 , Adjusted R-squared: 0.9616
## F-statistic: 714.1 on 1 and 19 DF, p-value: < 2.2e-16\(\beta_1 = 1\)という帰無仮説の下での\(t\)統計量は次のように求められます．
beta1 <- result_robust$coefficients["productivity"]
se <- result_robust$std.error["productivity"]
tstat <- (beta1 - 1) / se
tstat
## productivity
## -22.14628区間推定は次の通り行えます．1.96をqnorm(0.975)としても構いません．qnorm(p)関数は標準正規分布の\(100p\%\)分位点となる値を求める関数です．
upper <- beta1 + 1.96 * se
lower <- beta1 - 1.96 * se
lower
## productivity
## 0.506712
upper
## productivity
## 0.5869272練習問題 4-2
この練習問題のデータセットのln(2013人口）という変数は最後の丸括弧が全角になっています．また，いくつかの変数名は数字から始まっていますが，Rでは変数の名前を数字から始めることは推奨されません．さらに，変数名に日本語を含めることもコーディングの効率性の観点から推奨されません．したがって，はじめにデータセットの変数名を変更します．後は実証例4.1と同じような手続きです．
library(readxl)
mydata2 <- read_excel("data for chap 4 exercise 2.xlsx")
colnames(mydata2) <- c("prefecture", "pop2013", "gdp2013", "log_pop2013", "log_gdp2013")