44の例題で学ぶ計量経済学（オーム社） - R入門

目次

最小二乗法による回帰直線の傾き、切片、決定係数、自由度調整済み決定係数を求める（pp.96-108）

データは以下のとおり（p.97）。xxiは説明変数、yyiは目的変数（被説明変数）。これをメモ帳に貼り付けて「table3_2.csv」というテキストファイルで保存をして、カレントディレクトリに置いておく。

i, xxi, yyi
1,   2,   3
2,   4,   5
3,   6,   6
4,   8,  10

以下、計算。

> dtf <- read.csv("table3_2.csv", header = TRUE)
> xxi <- dtf$xxi
> yyi <- dtf$yyi
> n <- length(xxi)
> ssxy <- sum((xxi - mean(xxi)) * (yyi - mean(yyi)))
> ssxx <- sum((xxi - mean(xxi)) ^ 2)
> b <- ssxy / ssxx
> xxm <- mean(xxi)
> yym <- mean(yyi)
> a <- yym - b * xxm
> yyei <- a + b * xxi
> ui <- yyi - yyei
> ssu2 <- sum(ui ^ 2)
> ssyy <- sum((yyi - yym) ^ 2)
> ssyeye <- sum((yyei - yym) ^ 2)
> rr2 <- ssyeye / ssyy
> sig2 <- ssu2 / (n - 2)
> sy2 <- ssyy / (n - 1)
> arr2 <- 1 - sig2 / sy2
> cat(sprintf("傾き β=%.2f\n", b))
傾き β=1.10
> cat(sprintf("切片 α=%.2f\n", a))
切片 α=0.50
> cat(sprintf("決定係数 R^2=%.3f\n", rr2))
決定係数 R^2=0.931
> cat(sprintf("自由度調整済み決定係数 adj.R^2=%.3f\n", arr2))
自由度調整済み決定係数 adj.R^2=0.896

推定した回帰モデルのパラメーターの分散と標準誤差（p.171）

データはp.97に掲載されている（表示省略）。xxiは説明変数、yyiは目的変数（被説明変数）。これをメモ帳に貼り付けて「table3_2.csv」というテキストファイルで保存をして、カレントディレクトリに置いておく。

> dtf <- read.csv("table3_2.csv", header = TRUE)
> n <- nrow(dtf)
> k <- 2
> xxm <- mean(dtf$xxi)
> yym <- mean(dtf$yyi)
> ssxy <- sum((dtf$xxi - xxm) * (dtf$yyi - yym))
> ssxx <- sum((dtf$xxi - xxm) ^ 2)
> bh <- ssxy / ssxx
> ah <- yym - bh * xxm
> yyhi <- ah + bh * dtf$xxi
> ui <- dtf$yyi - yyhi
> sh2 <- sum(ui ^ 2) / (n - k)
> sbh2 <- sh2 / ssxx
> sbh <- sqrt(sbh2)
> sah2 <- sh2 * sum(dtf$xxi ^ 2) / (n * ssxx)
> sah <- sqrt(sah2)
> # β^の分散
> print(sbh2)
[1] 0.045
> # β^の標準誤差
> print(sbh)
[1] 0.212132
> # α^の分散
> print(sah2)
[1] 1.35
> # α^の標準誤差
> print(sah)
[1] 1.161895

単回帰モデルにおけるt分布による傾きの仮説検定（両側検定、片側検定）を行う（pp.178-180）

データは以下のとおり（p.97）。xxiは説明変数、yyiは目的変数（被説明変数）。これをメモ帳に貼り付けて「table3_2.csv」というテキストファイルで保存をして、カレントディレクトリに置いておく。

i, xxi, yyi
1,   2,   3
2,   4,   5
3,   6,   6
4,   8,  10

以下、計算。

> dtf <- read.csv("table3_2.csv", header = TRUE)
> xxi <- dtf$xxi
> yyi <- dtf$yyi
> n <- nrow(dtf)
> k <- 2
> degf <- n - k
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1.0, n), xxi), ncol = 2)
> mxb <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> b <- as.vector(mxb)
> yyhi <- mxxx %*% mxb
> ssu2 <- sum((yyi - yyhi) ^ 2)
> sh2 <- ssu2 / (n - k)
> ssxx <- sum((xxi - mean(xxi)) ^ 2)
> sb2 <- sh2 / ssxx
> sb <- sqrt(sb2)
> tb <- b[2] / sb
> cat(sprintf("残差分散 σ^2 = %.1f\n", sh2))
残差分散 σ^2 = 0.9
> cat(sprintf("βの推定値の分散 sb^2 = %.3f\n", sb2))
βの推定値の分散 sb^2 = 0.045
> cat(sprintf("βの推定値の標準誤差 sβ = %.5f\n", sb))
βの推定値の標準誤差 sβ = 0.21213
> cat(sprintf("βの推定値のt値 tβ = %.3f\n", tb))
βの推定値のt値 tβ = 5.185
> cat(sprintf("両側検定 t0.05(2) = %.3f\n", abs(qt(0.05 / 2, n - k))))
両側検定 t0.05(2) = 4.303
> cat(sprintf("片側検定 t0.05(2) = %.3f\n", qt(0.05, n - k, lower.tail = FALSE)))
片側検定 t0.05(2) = 2.920

単回帰モデルにおけるt分布による切片の仮説検定（両側検定）を行う（p.181）

データは以下のとおり（p.97）。xxiは説明変数、yyiは目的変数（被説明変数）。これをメモ帳に貼り付けて「table3_2.csv」というテキストファイルで保存をして、カレントディレクトリに置いておく。

i, xxi, yyi
1,   2,   3
2,   4,   5
3,   6,   6
4,   8,  10

以下、計算。

> dtf <- read.csv("table3_2.csv", header = TRUE)
> xxi <- dtf$xxi
> yyi <- dtf$yyi
> n <- nrow(dtf)
> k <- 2
> degf <- n - k
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1.0, n), xxi), ncol = 2)
> mxb <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> b <- as.vector(mxb)
> yyhi <- mxxx %*% mxb
> ssu2 <- sum((yyi - yyhi) ^ 2)
> sh2 <- ssu2 / (n - k)
> ssxx <- sum((xxi - mean(xxi)) ^ 2)
> sa2 <- sh2 * sum(xxi ^ 2) / (n * ssxx)
> sa <- sqrt(sa2)
> ta <- b[1] / sa
> cat(sprintf("残差分散 σ^2 = %.1f\n", sh2))
残差分散 σ^2 = 0.9
> cat(sprintf("αの推定値の分散 sa^2 = %.2f\n", sa2))
αの推定値の分散 sa^2 = 1.35
> cat(sprintf("αの推定値の標準誤差 sα = %.5f\n", sa))
αの推定値の標準誤差 sα = 1.16190
> cat(sprintf("αの推定値のt値 tα = %.3f\n", ta))
αの推定値のt値 tα = 0.430
> cat(sprintf("両側検定 t0.05(2) = %.3f\n", abs(qt(0.05 / 2, n - k))))
両側検定 t0.05(2) = 4.303

単回帰モデルにおける回帰係数のp値を求める（p.181）

データは以下のとおり（p.97）。xxiは説明変数、yyiは目的変数（被説明変数）。これをメモ帳に貼り付けて「table3_2.csv」というテキストファイルで保存をして、カレントディレクトリに置いておく。

i, xxi, yyi
1,   2,   3
2,   4,   5
3,   6,   6
4,   8,  10

以下、計算。

> dtf <- read.csv("table3_2.csv", header = TRUE)
> xxi <- dtf$xxi
> yyi <- dtf$yyi
> n <- nrow(dtf)
> k <- 2
> degf <- n - k
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1.0, n), xxi), ncol = 2)
> mxb <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> b <- as.vector(mxb)
> yyhi <- mxxx %*% mxb
> sssse <- sum((yyi - yyhi) ^ 2)
> sh2 <- sssse / degf
> mxcc <- sh2 * solve(t(mxxx) %*% mxxx)
> seb <- sqrt(diag(mxcc))
> tb <- b / seb
> pb <- 2 * pt(abs(tb), degf, lower.tail = FALSE)
> rr2 <- sum((yyhi - mean(yyi)) ^ 2) / sum((yyi - mean(yyi)) ^ 2)
> arr2 <- 1 - (sssse / (n - k)) / (sum((yyi - mean(yyi)) ^ 2) / (n - 1))
> cat(sprintf("回帰式 Y^i =  %.1f  +  %.1f Xi\n", b[1], b[2]))
回帰式 Y^i =  0.5  +  1.1 Xi
> cat(sprintf("            (%.3f) (%.3f)\n", tb[1], tb[2]))
            (0.430) (5.185)
> cat(sprintf("αの推定値の t値 tα = %.3f\n", tb[1]))
αの推定値の t値 tα = 0.430
> cat(sprintf("αの推定値の p値 pα = %.3f\n", pb[1]))
αの推定値の p値 pα = 0.709
> cat(sprintf("βの推定値の t値 tβ = %.3f\n", tb[2]))
βの推定値の t値 tβ = 5.185
> cat(sprintf("βの推定値の p値 pβ = %.3f\n", pb[2]))
βの推定値の p値 pβ = 0.035

重回帰モデルにおけるt検定（pp.185-188）

以下の表4.5（p.186）の値をCSV形式で入力したtable4_5.csvを、カレントディレクトリに置いておく。

i, xx2i, xx3i, yyi
1,    2,   26,   3
2,    4,    5,   5
3,    6,   21,   6
4,    8,    8,  10
5,   10,   20,  11

計算する。

> dtf <- read.csv("table4_5.csv", header = TRUE)
> xx2i <- dtf$xx2i
> xx3i <- dtf$xx3i
> yyi <- dtf$yyi
> kk <- 3
> n <- nrow(dtf)
> degf <- n - kk
> qv <- 0.05
> tv <- qt(1 - 0.05 / 2, 2)
> r <- lm(yyi ~ 1 + xx2i + xx3i)
> bhk <- as.vector(r$coefficient)
> sh2 <- sum(r$residuals ^ 2) / degf
> sk2 <- sk <- tk <- double(kk)
> xx2m <- mean(xx2i)
> xx3m <- mean(xx3i)
> ss22 <- sum((xx2i - xx2m) * (xx2i - xx2m))
> ss23 <- sum((xx2i - xx2m) * (xx3i - xx3m))
> ss33 <- sum((xx3i - xx3m) * (xx3i - xx3m))
> d1 <- xx2m ^ 2 * ss33 - 2 * xx2m * xx3m * ss23 + xx3m ^ 2 * ss22
> d2 <- ss22 * ss33 - ss23 ^ 2
> sk2[1] <- sh2 * (1 / n + d1 / d2)
> sk2[2] <- sh2 * ss33 / d2
> sk2[3] <- sh2 * ss22 / d2
> sk <- sqrt(sk2)
> tk <- bhk / sk
> # 残差分散 σ^^2
> print(sh2)
[1] 0.7197232
> # 推定値 β^^k の分散
> print(sk2)
[1] 1.583391003 0.018451538 0.002263992
> # 推定値 β^^k の標準誤差
> print(sk)
[1] 1.25832865 0.13583644 0.04758143
> # 推定値 β^^k のt値
> print(tk)
[1]  1.1219364  7.6037916 -0.7999399

重回帰モデルにおける仮説検定（pp.205-206）

> yyi <- c(3, 6, 8, 12, 11)
> xx2i <- c(2, 4, 6, 8, 10)
> xx3i <- c(1, 1, 2, 2, 4)
> n <- length(yyi)
> kk <- 3
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- as.matrix(data.frame(rep(1, n), xx2i, xx3i))
> mxxxtxxi <- solve(t(mxxx) %*% mxxx)
> mxbh <- mxxxtxxi %*% t(mxxx) %*% mxy
> bhk <- as.vector(mxbh)
> mxyh <- mxxx %*% mxbh
> ui <- as.vector(mxy - mxyh)
> sh2 <- sum(ui ^ 2) / (n - kk)
> sk2 <- as.vector(diag(sh2 * mxxxtxxi))
> sk <- sqrt(sk2)
> tk <- bhk / sk
> print(n - kk)  # 自由度
[1] 2
> print(sh2)  # 残差分散
[1] 0.1818182
> print(sk)  # 推定した回帰モデルのパラメーターの標準誤差
[1] 0.4490578 0.1574592 0.4065578
> print(tk)  # t値
[1]  2.631773 11.835681 -5.366563

異常値を含むダミー変数を用いた回帰直線の推定（pp.226-227）

> dtf <- read.csv("table5_1.csv", header = TRUE)
> n <- nrow(dtf)
> xxi <- dtf$xxi
> ddi <- dtf$ddi
> yyi <- dtf$yyi
> k <- 3
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1.0, n), xxi, ddi), ncol = 3)
> mxbh <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> mxyh <- mxxx %*% mxbh
> mxu <- mxy - mxyh
> ssy2 <- sum(t(mxy) %*% mxy) - sum(yyi) ^ 2 / n
> ssyh2 <- sum(t(mxyh) %*% mxyh) - sum(yyi) ^ 2 / n
> sse2 <- as.vector(t(mxu) %*% mxu)
> s2 <- sse2 / (n - k)
> s <- sqrt(s2)
> mxvv <- s2 * solve(t(mxxx) %*% mxxx)
> mxs <- sqrt(diag(mxvv))
> rr2 <- ssyh2 / ssy2
> tval <- as.vector(mxbh) / as.vector(mxs)
> # 推定された回帰モデルのパラメーター α^,β^,γ^
> print(as.vector(mxbh))
[1]  0.5  1.1 15.8
> # 決定係数 R^2
> print(rr2)
[1] 0.9936886
> # 各パラメーターの t値
> print(tval)
[1]  0.4303315  5.1854497 13.8309443
> # 図の出力
> png("fig5_1.png", width = 512, height = 512)
> plot(xxi, yyi, asp = 1., type = "n", xlim = c(0, 12), ylim = c(0, 30), xlab = "X", ylab = "Y")
> x <- seq(0, 12, length = 100)
> lines(x, mxbh[1] + mxbh[2] * x, lty = "solid")
> lines(x, mxbh[1] + mxbh[2] * x + mxbh[3], lty = "dashed")
> points(xxi, yyi, pch = 20)
> dev.off()
null device 
          1

グループに対するダミー変数（性別）（pp.231-232）

以下の表5.4（p.232）の値をCSV形式で入力したtable5_4.csvを、カレントディレクトリに置いておく。

i,  sex,   xxi,  yyi
 1, 女性, 255.5, 15.8
 2, 女性, 259.7, 16.8
 3, 女性, 385.3, 16.4
 4, 女性, 286.5, 15.7
 5, 女性, 302.1, 16.3
 6, 女性, 317.8, 17.4
 7, 女性, 321.5, 18.3
 8, 女性, 393.0, 17.6
 9, 女性, 393.7, 17.4
10, 女性, 289.2, 15.3
11, 男性, 352.2,  9.7
12, 男性, 367.1,  5.2
13, 男性, 276.9,  7.1
14, 男性, 348.4,  9.6
15, 男性, 277.8,  6.9
16, 男性, 386.5,  8.2
17, 男性, 337.7,  7.3
18, 男性, 219.2,  6.6
19, 男性, 330.2,  8.7
20, 男性, 395.8, 10.1

計算する。

> dtf <- read.csv("table5_4.csv")
> r <- lm(yyi ~ 1 + xxi, dtf, trimws(sex) == "女性")
> print(r$coefficients)  # 女性だけ
 (Intercept)          xxi 
13.722088020  0.009293487 
> r <- lm(yyi ~ 1 + xxi, dtf, trimws(sex) == "男性")
> print(r$coefficients)  # 男性だけ
(Intercept)         xxi 
  3.5590101   0.0133088 
> ddi <- ifelse(trimws(dtf$sex) == "女性", 1, 0)
> r <- lm(dtf$yyi ~ 1 + dtf$xxi + ddi)
> print(summary(r))  # 女性と男性
Call:
lm(formula = dtf$yyi ~ 1 + dtf$xxi + ddi)
Residuals:
    Min      1Q  Median      3Q     Max 
-3.1720 -0.4944 -0.1475  0.7592  1.5878 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.19004    1.74965   2.395   0.0284 *  
dtf$xxi      0.01139    0.00519   2.195   0.0423 *  
ddi          8.85968    0.53539  16.548 6.45e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.193 on 17 degrees of freedom
Multiple R-squared:  0.9417,    Adjusted R-squared:  0.9348 
F-statistic: 137.2 on 2 and 17 DF,  p-value: 3.234e-11

---