このページに対応するRmdファイル：GitHub

1 Data

このページで使用するデータは，カーセンサーnetという中古車情報サイトに掲載されている人気車種ランキング Top 80 の車両に関する車種レベルの属性．

ランキングはカーセンサーに掲載されている自動車への問い合わせデータに基づいて車種別に作成されている．
- このページでは，車種とは「プリウス」（トヨタのハイブリッド車）とか「N-BOX」（ホンダのトールワゴン型軽自動車）のようなカテゴリーを指す語として用いる．
一般的に一つの車種には様々なバリエーションがあるが，そのうち一つを選択しその車種を代表する属性とした．
- たとえば，トヨタのプリウスを例にとると，まず世代と呼ばれる区分けがあり（初代 1997～，2代目 2003～，3代目 2009～，4代目 2015～），さらに「Z」・「G」・「X」・「U」など幾つかのグレードに分かれ，そこからさらにハイブリッドかプラグインハイブリッドか，駆動方式が前輪駆動（FF）か四輪駆動（4WD）か，などの基準で枝分かれする．
- 幾つかあるバリエーションのうちデータ入力者（黒田）がいずれを選択したか（の規則性）は分析結果に影響を及ぼすが，この演習ではそれは気にしないことにしよう．
2025年6月のランキングに基づいてデータを作成した．データセットにランクを表す列はないが，ランキングのとおり昇順で並んでいる（列名を表す行を除いて1行目が1位で，80行目が80位）．
変数
- model: 車種
- maker: 製造メーカー
- year: 世代の開始年（一つ前の世代からモデルチェンジした年）［西暦］
- price: 新車時価格（税込）［万円］
- seat: 乗車定員
- length: 全長［m］
- width: 全幅［m］
- height: 全高［m］
- weight: 車両重量［kg］
- WLTC: WLTCモード燃費［km/L］
- ps: 最高出力［ps］（注：ドイツ語の Pferdestärke の略で馬力の意．1 PS ≒ 0.74 kW）
- disp: 排気量［cc］
- hybrid: ハイブリッド車ダミー
- electric: 電気自動車ダミー

car <- read.csv("https://raw.githubusercontent.com/kurodaecon/dasp/refs/heads/main/data/carsensor_ranking_of_most_popular_car_models_asof202506.csv")
str(car)

## 'data.frame':    80 obs. of  14 variables:
##  $ model   : chr  "N-BOX" "Serena" "Prius" "Tanto" ...
##  $ maker   : chr  "Honda" "Nissan" "Toyota" "Daihatsu" ...
##  $ year    : int  2023 2022 2022 2019 2023 2024 2022 2022 2021 2019 ...
##  $ price   : num  174 362 299 148 555 ...
##  $ seat    : int  4 7 5 4 7 6 7 5 7 4 ...
##  $ length  : num  3.4 4.69 4.6 3.4 5 4.31 4.8 4.26 4.7 3.4 ...
##  $ width   : num  1.48 1.7 1.78 1.48 1.85 1.7 1.75 1.7 1.73 1.48 ...
##  $ height  : num  1.79 1.89 1.42 1.76 1.94 1.76 1.84 1.7 1.9 1.68 ...
##  $ weight  : int  910 1900 1360 880 2060 1460 1710 1280 1640 830 ...
##  $ WLTC    : num  21.6 17 32.6 22.7 10.6 25.6 13.2 18.4 23 25 ...
##  $ ps      : int  58 98 98 52 182 106 150 120 98 49 ...
##  $ disp    : int  658 1433 1797 658 2493 1496 1496 1490 1797 657 ...
##  $ hybrid  : int  0 1 1 0 0 1 0 0 1 0 ...
##  $ electric: int  0 0 0 0 0 0 0 0 0 0 ...

このページ全体を通した分析の目的は「車両の価格と属性の関係を定量的に記述すること」としよう．因果推論は目的ではない．

1.1 Summary stats

長くなるのでこのページには結果を出力しない設定にしているが，実際にはこれを見ながらデータの分布と特徴を理解するのもデータ分析の重要な工程．

summary(car[, c(-1, -2)])
table(car$maker)
car[car$maker == "Nissan", ]
car[car$maker %in% c("BMW", "Jeep", "Mercedes-Benz", "Porsche"), ]
hist(car$price)
hist(car$price[car$price <= 1000])
car[car$price > 1000, ]
table(car$seat)
table(car$hybrid)
table(car$electric)
car[car$electric == 1, ]
car[, c(-1, -2, -14)] |> cor(use = "complete.obs") |> round(2)  # "|>" = pipe in base R 
# round(cor(car[, c(-1, -2, -14)], use = "complete.obs"), 2)  # same as above 
plot(car[, c("price", "weight", "height", "WLTC", "ps", "year")])

1.2 Create variables

car$accel <- car$ps / car$weight  # 加速性能の大まかな指標 (power-to-weight ratio) 
car$kei <- 1*(car$disp <= 660)  # 軽自動車ダミー

1.3 Trim dataset

次の車両は除外する：乗車定員が3以下（スポーツカー・軽トラック等）または6以上（ミニバン・ワゴン車），海外メーカー，1,000万円を超える自動車，電気自動車．

car <- car[car$seat %in% 4:5, ]
car <- car[!(car$maker %in% c("BMW", "Jeep", "Mercedes-Benz", "Porsche")), ]
car <- car[car$price <= 1000, ]
car <- car[car$electric == 0, ]
# which(colnames(car) == "electric")  # 14
car <- car[, -14]  # remove "electric" variable (column) 
nrow(car)  # sample size

## [1] 50

2 Simple regression analysis / 単回帰分析

1つの説明変数をもつ次のような回帰モデル（単回帰モデル）を考える．

\[ \mbox{Price} = \alpha + \beta \ \mbox{Weight} + u \]

plot(price ~ weight, car, cex = 1.5, xlab = "Curb Weight [kg]", ylab = "Price [10,000 JPY]")

上の plot 関数は plot(formula = y ~ x, data = dataset) の形式で変数を指定しており，plot(x = dataset$x, y = dataset$y) と同じ意味．

2.1 Estimation with `lm` function

lm (linear model) 関数で推定する．

lm(formula = price ~ weight, data = car)

## 
## Call:
## lm(formula = price ~ weight, data = car)
## 
## Coefficients:
## (Intercept)       weight  
##   -169.6602       0.3651

次のように formula = や data = を省略してもよい（出力は省略）．

lm(price ~ weight, car)

以上の推定結果は以下のように回帰式が推定されたことを意味している．

\[ \mbox{Price} = -170 + 0.37 \ \mbox{Weight} + u \]

weight の回帰係数は「weight が1単位高いと price は 0.37 単位高い傾向にある」＝「車両重量が 1kg 高い車の価格は3,700円高い」と解釈される．

ある特定の車種について「設計を変更して Weight を1単位引き上げると Price を 0.37 単位引き上げることができる（市場適正価格が3,700円上昇する）」あるいは「政策的介入によって重量が1kg増加したら均衡価格が3,700円上昇する」などのように解釈することはできない．この回帰係数はあくまで両変数の相関的な関係を記述しているだけであり（注：相関係数とは異なる），因果関係を示すものではない．

次のように考えてもよい．

weight = 1,000kg の平均価格は $-170 + 0.37 \cdot 1000 = 200$ 万円
weight = 1,100kg の平均価格は $-170 + 0.37 \cdot 1100 = 237$ 万円
- 100kg で37万円差が付くので，1kg 当たり3,700円違う

説明変数が1つの場合，回帰式は下図のように散布図に書き足すことができる．大雑把に言えば，回帰式の推定（＝回帰係数パラメタを計算すること）とはデータに当てはまる直線の切片 $\alpha$ と傾き $\beta$ を決める手続きである．

plot(price ~ weight, car, xlim = c(-100, 2100), ylim = c(-200, 800))
abline(v = 0, col = "grey")  # Y-axis
abline(lm(price ~ weight, car), col = "red")
points(0, -169.66, pch = 15, col = "red")
text(250, -170, expression(alpha == -170), cex = 1.5)
text(600, -70, expression(beta == 0.37), cex = 1.5)

lm 関数で計算した回帰オブジェクトを summary 関数に渡すと決定係数なども表示される．

summary(lm(price ~ weight, car))

## 
## Call:
## lm(formula = price ~ weight, data = car)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -212.197  -30.221   -5.594   27.261  170.713 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -169.66019   29.72521  -5.708 6.99e-07 ***
## weight         0.36510    0.02422  15.075  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 64.15 on 48 degrees of freedom
## Multiple R-squared:  0.8256, Adjusted R-squared:  0.822 
## F-statistic: 227.3 on 1 and 48 DF,  p-value: < 2.2e-16

Multiple R-squared が決定係数，Adjusted R-squared が自由度調整済み決定係数を表す．

2.2 Estimation by hand

回帰モデル

\[ Y_i = \alpha + \beta X_i + u_i \]

の係数パラメタ（$\alpha, \beta$）を最小二乗法（OLS）によって求めよう．

最小二乗法とは，誤差 $u = Y - (\alpha + \beta X_i)$ （残差 $\hat{u}_i = Y_i - \hat{Y}_i = Y_i - (\hat{\alpha}+\hat{\beta}X_i)$ で読み替えてもよい）の二乗和を最小化するように係数パラメタを求める方法である．

誤差の二乗和 $Q$ は以下で与えられる．

\[ Q = \sum_{i} u_i^2 = \sum_{i} (Y_i - \alpha - \beta X_i)^2 \]

$Q$ は $\hat \alpha, \hat \beta$ のいずれについても下に凸の関数であるので，$Q$ をそれぞれで微分して 0 とおく．

\[ \frac{\partial Q}{\partial \alpha} = -2 \sum_i (Y_i - \alpha - \beta X_i) = 0 \]

\[ \frac{\partial Q}{\partial \beta} = -2 \sum_i X_i (Y_i - \alpha - \beta X_i) = 0 \]

上の連立方程式を解くことによって，最小二乗推定量は次のように与えられる．

\[ \hat{\alpha} = \bar{Y} - \hat{\beta} \bar{X}, \quad \hat{\beta} = \frac{S_{XY}}{S_{XX}} \] \[ S_{XX} = \sum_i (X_i - \bar{X})^2, \quad S_{XY} = \sum_i (X_i - \bar{X}) (Y_i - \bar{Y}) \]

$\bar{Y}, \bar{X}, S_{XX}, S_{XY}$ という4つの統計量は次のように計算できる．

x <- car$weight
y <- car$price
y_bar <- mean(y)  # Y bar
x_bar <- mean(x)  # X bar
s_xx <- sum((x - x_bar)^2)  # S_XX
s_xy <- sum((x - x_bar) * (y - y_bar))  # S_XY
c(y_bar, x_bar, s_xx, s_xy)

## [1]     257.07    1168.80 7015928.00 2561523.20

計算した統計量を使って回帰係数を計算する．

beta_hat <- s_xy / s_xx
alpha_hat <- y_bar - beta_hat * x_bar
c(alpha_hat, beta_hat)

## [1] -169.6601939    0.3651011

2.3 Estimation by numerical calculation / 数値計算による推定

OLSは残差二乗和の最小化によって回帰係数パラメタを推定する方法（モデルをデータに当てはめる原理）であるが，これを解析的に解くのではなく数値計算によって解いてみる．

残差二乗和の関数（目的関数） $Q$ の最小化は R に built-in されている汎用の最適化関数 optim によって行う．

2.3.1 Example of `optim`

例：次の二次関数は $x=2$ で最小値 3 をとる．

\[ f(x) = x^2 -4x + 7 = (x-2)^2 + 3 \]

最小値を取るような $x$ の値を optim 関数を使って求めてみる．

optim(par = パラメタの初期値ベクトル, fn = 目的関数) の形で指定
- $par ：最適化されたパラメタ（＝目的関数の引数；この例では $x$）
- $value ：最小化された目的関数の値
- $convergence ：最適化計算が収束していれば 0

obj_function <- function (x) x^2 - 4*x + 7
optim(par = 0, fn = obj_function, method = "BFGS")  # ニュートン法の亜種を利用

## $par
## [1] 2
## 
## $value
## [1] 3
## 
## $counts
## function gradient 
##        6        3 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

なお，ニュートン法は次の更新式を用いて実装できる．

\[ x_{\mbox{new}} = x - \frac{f'(x)}{f''(x)} = x - \frac{2x-4}{2} = 2 \]

今扱っている目的関数は単純な二次式なので初期値によらず1回の更新で収束するが，通常は何度も更新を繰り返す．

grad_f <- function (x) 2 * x - 4  # 目的関数の1階微分（勾配）
newton_method <- function (initial) {
  x <- initial
  for (i in 1:10) {  # 10 = max iteration 
    step <- grad_f(x) / 2  # 分母の2 = 目的関数の2階微分（Hessian）
    x_new <- x - step
    cat("i=", i, ", f'=", grad_f(x), ", f'/f''=", step, ", x=", x, 
        ", x_new=", x_new, ", Δx=", x_new - x, "\n", sep = "")
    if (abs(x_new - x) < 0.01) return(x_new)  # 0.01 = tolerance 
    x <- x_new
  }
  cat("Not converged")
}
newton_method(initial = 10)

## i=1, f'=16, f'/f''=8, x=10, x_new=2, Δx=-8
## i=2, f'=0, f'/f''=0, x=2, x_new=2, Δx=0

## [1] 2

2.3.2 In the case of regression

回帰分析の場合，目的関数は残差二乗和 $Q$ であり，その引数は推定したい回帰係数パラメタ $(\alpha, \beta)$ となる．

\[ (\hat{\alpha}, \hat{\beta}) = \arg \min_{\alpha, \beta} Q(\alpha, \beta) \] \[ Q (\alpha, \beta) := \sum_i \hat{u}_i^2 = \sum_i (Y_i - \alpha - \beta X_i)^2 \]

sum_of_squared_residuals <- function (coef_vector) {
  # coef_vector[1]: alpha (intercept)
  # coef_vector[2]: beta (slope of X)
  sum((y - coef_vector[1] - coef_vector[2] * x)^2)
}
sum_of_squared_residuals(coef_vector = c(0, 1))  # alpha=1, beta=2

## [1] 44588206

sum_of_squared_residuals(coef_vector = c(0, 0.5))  # alpha=1, beta=3

## [1] 5682450

optim(par = c(0, 0), fn = sum_of_squared_residuals)

## $par
## [1] -169.6753516    0.3651259
## 
## $value
## [1] 197529.7
## 
## $counts
## function gradient 
##      107       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

$\hat{\alpha} = -170$ に固定して $\hat{\beta}$ の値によって残差二乗和がどのように変化するかをプロットすると，$\hat{\beta} = 0.37$ で最小化されていることが分かる．

2.4 Coefficient of determination / 決定係数

\[ R^2 = \frac{\mbox{Variation in } Y \mbox{ explained by } X}{\mbox{Total variation in } Y} = \frac{S_{\hat{Y}\hat{Y}}}{S_{YY}} = 1 - \frac{\sum \hat{u}_i^2}{S_{YY}} \]

\[ S_{YY} = \sum_i (Y_i - \bar{Y})^2, \quad S_{\hat{Y} \hat{Y}} = \sum_i (\hat{Y}_i - \bar{\hat{Y}})^2 \]

s_yy <- sum((y - y_bar)^2)
y_pred <- alpha_hat + beta_hat * x  # y hat
s_yhyh <- sum((y_pred - y_bar)^2)
c(s_yhyh, s_yy)

## [1]  935215 1132745

s_yhyh / s_yy  # R^2

## [1] 0.8256185

res <- y - y_pred  # residuals
1 - sum(res^2) / s_yy  # R^2 ... should be the same as above

## [1] 0.8256185

決定係数が 0.83 であることは「説明変数 (weight) によって被説明変数 (price) の変動のうち 83% を説明することができる」という意味．

余裕のある履修者向けの宿題：自由度調整済み決定係数を計算．

\[ \bar{R}^2 = 1 - \frac{n - 1}{n - K - 1} \frac{\sum \hat{u}_i^2}{S_{YY}} \]

$K$: 説明変数の数（定数項は含めない）

ヒント：$n$ は nrow(car) で与えられる（欠損値があればその分の調整も必要）．

2.5 Standard error and test / 標準誤差と検定

帰無仮説「$H_0: \beta = 0$」の検定を行うために，$\hat{\beta}$ の標準誤差（s.e.）を求める．

\[ V(\hat{\beta}) = \frac{\sigma^2}{S_{XX}}, \quad s.e.(\hat{\beta}) = \sqrt{V(\hat{\beta})} \]

誤差分散 $\sigma^2$ は残差 $\hat{u}$ を利用して次のように推定する．

\[ \hat{\sigma}^2 = s^2 = \frac{1}{n - 2} \sum \hat{u}_i^2 \]

n <- nrow(car)  # sample size 
sum(res^2) / (n - 2)  # estimated sigma^2

## [1] 4115.202

beta_se <- sqrt(sum(res^2)/(n - 2) / s_xx)  # standard error
beta_se

## [1] 0.02421882

\[ t = \frac{\hat{\beta}}{\sqrt{V(\hat{\beta})}} \]

beta_hat / beta_se

## [1] 15.0751

\[ p = 2 \times [ 1 - \Pr(T \le |t|) ] \]

curve(dt(x, df = n - 2), from = -20, to = 20, main = "Density of t-distribution")
abline(v = beta_hat / beta_se * c(-1, 1), col = "red")
text(c(-17, 17), 0.02, labels = c("Pr(T < -15.1)", "Pr(T > 15.1)"))

2 * (1 - pt(q = beta_hat / beta_se, df = n - 2))

## [1] 0

注：あまりにも小さいので 0 と表示されているが，厳密には 0 ではない．このような場合，論文・レポート等では「$p < 0.001$」（表中では「$< 0.001$」）のように書くことで，厳密には0ではないものの 0.1% 水準で有意であることが報告できる．

2.5.1 Critical value approach / 臨界値を用いて検定を行う場合

例として，有意水準 $\alpha = 0.05$ と設定した場合を考えよう．

t分布の97.5%分位点（両側5%に対応する臨界値）は qt 関数で計算できる．

qt(p = 0.975, df = n - 2)  # df: degree of freedom (自由度)

## [1] 2.010635

qt(p = c(0.975, 0.995, 0.9995), df = n - 2)  # 5%, 1%, and 0.1%

## [1] 2.010635 2.682204 3.505068

curve(dt(x, df = n - 2), from = -4, to = 4, main = "Density of t-distribution")
abline(v = c(-2.01, 2.01), col = "red")
text(c(-3, 3), 0.04, labels = "2.5%"); text(0, 0.15, labels = "95%")

よって，$t = 15.1$ は有意水準5%に対応する臨界値 (2.01) よりも大きいので，5%水準で有意と言える．

同様に考えて $t = 15.1$ であれば1%でも0.1%でも有意なので，論文・レポートでは「変数 $X_j$ の回帰係数は0.1%水準で有意（に0と異なる）」と書く場合もある．これは，単に帰無仮説（$\beta_j = 0$）を棄却するかどうかの二値的な情報だけでなく，その結論の確実性がどの程度高いか（← 非常にラフな表現）という情報を報告する意図があるものと思われる．ただし，このことは「有意になるように（あるいは有意にならないように），$p$ 値を計算した後に有意水準を設定してよい」ことを意味しない．
Rでは非常に小さい $p$ 値が 2e-16 のように表示される場合があるが，これをそのまま「$p = 2\times10^{-16}$」のようにレポートするのはNG．2e-16 という厳密な数値が計算されているわけではなくて，数値計算上の精度の限界によりこれより小さい値である（$p < 2\times10^{-16}$）ことしか分からないという意味である．そして，この場合の $2$ あるいは $10^{-16}$ という数値の厳密な値は実務上ほとんど意味をなさない．たとえば，これが $2.5 \times 10^{-15}$ であったとしても仮説検定の結論は何も変わらない．よって，このような場合は「$p < 0.001$」と書けば十分である．

自由度が大きければ標準正規分布の場合とほとんど同じ臨界値になる．

c(qt(p = 0.975, df = 10000), qnorm(p = 0.975))

## [1] 1.960201 1.959964

2.5.2 MC Simulation

$\beta = 0$ の設定でデータを生成し，その回帰係数の $t$ 値の分布と第一種の過誤を確認してみよう．

t_list <- NULL
for (i in 1:1000) {
  set.seed(100*i)  # ちょっとごまかす
  x_i <- 1:1000
  y_i <- 1 + 0 * x_i + rnorm(1000)  # x_i の係数を 0 に設定
  t_i <- summary(lm(y_i ~ x_i))$coef["x_i", "t value"]
  t_list <- c(t_list, t_i)
}
# qt(p = 0.975, df = 1000-1)  # 1.962341
hist(t_list, breaks = 100)
abline(v = qt(p = c(0.025, 0.975), df = 1000-1), col = "red")

table(t_list <= -1.962 | t_list >= 1.962)  # TRUE = 第一種の過誤

## 
## FALSE  TRUE 
##   953    47

2.5.3 Robust standard error / 頑健な標準誤差

詳しい議論はこの授業では行わないが，標準誤差は誤差項の不均一分散を前提として計算した方がよい．

fixest パッケージの feols 関数が便利．下のようにvcov オプションを指定する．

# install.packages("fixest")
library(fixest)
fixest::feols(price ~ weight, car, vcov = "hetero")  # HC1

## OLS estimation, Dep. Var.: price
## Observations: 50 
## Standard-errors: Heteroskedasticity-robust 
##                Estimate Std. Error  t value   Pr(>|t|)    
## (Intercept) -169.660194  43.287229 -3.91941 2.8103e-04 ***
## weight         0.365101   0.042941  8.50237 3.8751e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 62.9   Adj. R2: 0.821986

$\beta = 0$ かつ不均一な誤差分散をもつモデルに基づいてデータを生成し，第一種の過誤を通常のSE vs. 頑健なSEで比較してみよう．

plot(x = 1:100, y = 1 + rnorm(100, sd = 1:100), xlab = "x", ylab = "y")

t_list <- t_r_list <- NULL
for (i in 1:1000) {
  set.seed(i)
  data_i <- data.frame(
    x_i = 1:1000,
    y_i = 1 + 0 * (1:1000) + rnorm(1000, sd = 1:1000)
  )
  t_i <- fixest::feols(y_i ~ x_i, data_i)$coeftable["x_i", "t value"]
  t_r_i <- fixest::feols(y_i ~ x_i, data_i, vcov = "hetero")$coeftable["x_i", "t value"]
  t_list <- c(t_list, t_i)
  t_r_list <- c(t_r_list, t_r_i)
}
qt(p = 0.975, df = 1000-1)

## [1] 1.962341

table(t_list <= -1.962 | t_list >= 1.962)  # 第一種の過誤が5%を超える

## 
## FALSE  TRUE 
##   928    72

table(t_r_list <= -1.962 | t_r_list >= 1.962)  # 5%がTRUEになるはず -> 概ねOK

## 
## FALSE  TRUE 
##   951    49

2.5.4 Cluster-robust standard error / クラスター頑健標準誤差

vcov = ~ cluster の形式でオプションを指定する．

fixest::feols(price ~ weight, car, vcov = ~ maker)

## OLS estimation, Dep. Var.: price
## Observations: 50 
## Standard-errors: Clustered (maker) 
##                Estimate Std. Error  t value   Pr(>|t|)    
## (Intercept) -169.660194  31.124525 -5.45101 1.5859e-03 ** 
## weight         0.365101   0.027755 13.15432 1.1913e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 62.9   Adj. R2: 0.821986

上記を regression table の形式で比較するには huxtable パッケージを使うと便利．

推定値の下段の丸カッコが標準誤差を示す．
modelsummary パッケージが使えればそちらでもよい．黒田の環境ではうまくいかなかったのでこのページでは huxtable で代用．

# install.packages("huxtable")
library(huxtable)
lm0 <- fixest::feols(price ~ weight, car)
lm1 <- fixest::feols(price ~ weight, car, vcov = "hetero")
lm2 <- fixest::feols(price ~ weight, car, vcov = ~ maker)
huxtable::huxreg("Classical SE" = lm0, "Robust SE" = lm1, "Clustered SE" = lm2)

	Classical SE	Robust SE	Clustered SE
(Intercept)	-169.660 ***	-169.660 ***	-169.660 **
	(29.725)	(43.287)	(31.125)
weight	0.365 ***	0.365 ***	0.365 ***
	(0.024)	(0.043)	(0.028)
N	50	50	50
R2	0.826	0.826	0.826
logLik	-277.987	-277.987	-277.987
AIC	559.975	559.975	559.975
* p < 0.001; p < 0.01; * p < 0.05.

3 Multiple regression analysis / 重回帰分析

2つの説明変数をもつ次のような回帰モデル（重回帰モデル）を考える．

\[ \mbox{Price} = \alpha + \beta_1 \ \mbox{Weight} + \beta_2 \ \mbox{PS} + u \]

summary(lm(price ~ weight + ps, car))

## 
## Call:
## lm(formula = price ~ weight + ps, data = car)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -141.669  -26.642    2.186   17.320  211.752 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -103.71044   31.12576  -3.332 0.001708 ** 
## weight         0.23679    0.03844   6.160 1.67e-07 ***
## ps             0.77649    0.19347   4.014 0.000218 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 56.39 on 46 degrees of freedom
##   ( 1 個の観測値が欠損のため削除されました )
## Multiple R-squared:  0.8686, Adjusted R-squared:  0.8629 
## F-statistic: 152.1 on 2 and 46 DF,  p-value: < 2.2e-16

weight の係数は「ps の水準を固定したもとで weight が1単位高い値を取ると price は 0.24 単位高い値を取る傾向にある」≒「馬力が同じ自動車同士で比較すると，車両が 1kg 重い自動車の価格は 2,400 円高い」と解釈される．

なお，「単回帰分析では因果関係が推定できないが，重回帰分析では因果関係が推定できる」という主張は一般には誤り．

重回帰分析によって因果効果を推定するためには，処置変数の条件付き独立（かなり大雑把に言うと，交絡変数がすべてコントロール変数として正確な関数形でモデルに含められていること）が成立していなければならず，この仮定を社会経済データを用いた実証分析において満たすことは容易ではない．

3.1 Multicollinearity / 多重共線性

3.1.1 Example of perfect multico. / 完全な多重共線性の例

lm(price ~ weight + I(2 * weight), car)$coef

##   (Intercept)        weight I(2 * weight) 
##  -169.6601939     0.3651011            NA

cor(car$weight, 2 * car$weight)

## [1] 1

R では，該当する変数（のうち一つ）が自動的に drop される．

3.1.2 Example of imperfect multico. / 不完全な多重共線性の例

互いに強い相関関係にある weight, length, width という3つの変数に着目しよう．いずれの変数も price と正の相関を持ち，単回帰分析では正の回帰係数を得る．

しかしながら，3つの変数を一つの回帰式に説明変数として同時に含めると，いずれの変数も標準誤差が大きくなる．ラフに言えば，これは推定が安定しなくなることを意味する．実際，width の回帰係数は有意に負の値をとっているが，これは解釈しにくい結果である．

car[, c("price", "weight", "length", "width")] |> cor() |> round(2)

##        price weight length width
## price   1.00   0.91   0.86  0.80
## weight  0.91   1.00   0.94  0.91
## length  0.86   0.94   1.00  0.97
## width   0.80   0.91   0.97  1.00

# round(cor(car[, c("price", "weight", "length", "width")]), 2)  # same as above 
lm0 <- lm(price ~ weight, car)
lm1 <- lm(price ~ length, car)
lm2 <- lm(price ~ width, car)
lm3 <- lm(price ~ weight + length + width, car)
huxtable::huxreg(lm0, lm1, lm2, lm3)

	(1)	(2)	(3)	(4)
(Intercept)	-169.660 ***	-635.980 ***	-974.501 ***	75.468
	(29.725)	(77.217)	(135.563)	(164.080)
weight	0.365 ***			0.330 ***
	(0.024)			(0.073)
length		225.199 ***		160.922
		(19.270)		(86.590)
width			746.497 ***	-510.580 *
			(81.783)	(240.597)
N	50	50	50	50
R2	0.826	0.740	0.634	0.841
logLik	-277.987	-287.979	-296.490	-275.644
AIC	561.975	581.958	598.980	561.288
* p < 0.001; p < 0.01; * p < 0.05.

3.1.3 VIF for diagnostics / VIF による診断

変数 $j$ の分散拡大係数（variance inflation factor, VIF）は，変数 $j$ をそれ以外の変数に回帰したときの決定係数 $R_j^2$ を使って以下のように計算される．

\[ VIF_j = \frac{1}{1 - R_j^2} \]

経験則として「$VIF > 10$ ならば多重共線性が生じている」と判断されることが多い．

たとえば $j =$ length としてこの変数のVIFを計算してみよう．

lm_length <- lm(length ~ weight + width, car)
r2_length <- summary(lm_length)$r.squared  # cf. str(summary(lm_length))
1 / (1 - r2_length)  # VIF

## [1] 31.69281

car パッケージの vif 関数を使うと便利．なお，このパッケージ名がデータセットを格納しているオブジェクト名と同じなので少しだけややこしいが，パッケージの方の名称は An R Companion to Applied Regression という書籍名から来ている．

# install.packages("car")
library(car)
car::vif(lm(price ~ weight + length + width, car))

##    weight    length     width 
##  9.449551 31.692806 19.094157

多重共線性への対処：

深刻な多重共線性が生じている場合は該当する変数のうち（少なくとも）一つをモデルから外すことがある．とりわけ予測モデルとして回帰モデルを推定する際は変数の除去は珍しくない．予測が目的の場合は，相関が高い変数を合成しそれを説明変数としてモデルに加える場合もあるだろう．
一方で key variable $X_1$ がアウトカムに与える因果効果を測ろうとして回帰モデルを推定する場合に，$X_1$ と相関がある（かつアウトカムに影響を与える）変数 $X_2$ を除外すると除外変数バイアスが生じる可能性があるため，$X_2$ は除外するべきではないだろう．
- ただし，$X_2$ が $X_1$ の影響を受ける変数でありかつアウトカムに影響を与える変数である場合に，$X_1$ がアウトカムに与える効果（$X_1$ から直接与える効果と $X_2$ を経由して間接的に与える効果の両方を合計したもの）を推定したい場合は，むしろ $X_2$ は除外するべきである．これは多重共線性とは別の理由による（いわゆる bad control の問題）．

3.2 Comparing the size of the coefficients / 回帰係数の大きさを比較する

例：weight と ps の回帰係数の（絶対値の）大きさを比較する．

lm(price ~ weight + ps, car)$coef

##  (Intercept)       weight           ps 
## -103.7104391    0.2367862    0.7764891

上の結果から「$|0.24| < |0.78|$ なので ps の方が weight の分散をより説明する」と結論付けることはできない．

なぜなら，係数の大きさは変数のスケール（ばらつき）に依存するから．

たとえば，線形モデル $y = \beta_0 + \beta_1 x + u$ において金額ベースの変数 $x$ の単位を「円」から「千円」に変えると，推定される値そのもの ($\hat{\beta}_1$) は1千倍になるが，推定された回帰係数の「解釈」は全く変わらない（この変数の重要性が1千倍になるわけではない）．

3.2.1 Comparing the change per one SD / 1標準偏差当たりの変化量の比較

c(0.24 * sd(car$weight), 0.78 * sd(car$ps, na.rm = TRUE))

## [1] 90.81462 58.69636

「$|90.8| > |58.7|$ なので weight の方が price の分散に対する説明力がより高い」と言える．

3.2.2 Standardized partial regression coefficient / 標準偏回帰係数

より一般的な方法は，標準化（平均 0，分散 1 になるように変換）した変数を使って回帰分析をするもの．

car$price_s <- (car$price - mean(car$price)) / sd(car$price)
car$weight_s <- (car$weight - mean(car$weight)) / sd(car$weight)
car$ps_s <- (car$ps - mean(car$ps, na.rm = TRUE)) / sd(car$ps, na.rm = TRUE)
lm(price_s ~ weight_s + ps_s, car)$coef

## (Intercept)    weight_s        ps_s 
## 0.007155726 0.589295242 0.384312017

「$|0.59| > |0.38|$ なので weight の方が price の分散をより説明する」と言える．

なお，t 値と p 値は標準化の前後で不変．

スケーリング（分散を 1 にする）とセンタリング（平均を 0 にする）を同時に行ってくれる scale という関数を使えばより簡単．

lm(price ~ weight + ps, as.data.frame(scale(car[, c(-1, -2)])))

（出力は省略）

4 Variable transformation / 変数変換

4.1 Dummy variables / ダミー変数

hybrid 変数は，ハイブリッド車のときに 1 をとるダミー変数である．ハイブリッド車以外であるガソリン車やディーゼル車で 0 をとる．

\[ \mbox{Price} = \alpha + \beta \ \mbox{Hybrid} + u \]

table(car$hybrid)

## 
##  0  1 
## 45  5

lm(price ~ hybrid, car)

## 
## Call:
## lm(formula = price ~ hybrid, data = car)
## 
## Coefficients:
## (Intercept)       hybrid  
##      249.38        76.86

この回帰係数 77 は「ハイブリッド車はハイブリッド車以外と比較して 77 万円高い」ことを意味している（ただし有意ではない）．

mean(car$price[car$hybrid == 0])

## [1] 249.3844

mean(car$price[car$hybrid == 1])

## [1] 326.24

mean(car$price[car$hybrid == 1]) - mean(car$price[car$hybrid == 0])

## [1] 76.85556

下の図のように理解できるだろう．

「□」は hybrid の値ごとの平均価格を表しており，これは回帰式に基づいて計算される価格の予測値（当てはめ値）に等しい．
- hybrid = 0 のとき，価格の予測値は $249 + 77 \cdot 0 = 249$ 万円（すなわち，切片の推定値そのもの）
- hybrid = 1 のとき，価格の予測値は $249 + 77 \cdot 1 = 326$ 万円
傾きは $\Delta y / \Delta x = (326 - 249) / 1 = 77$ である．

条件分岐によってダミー変数を作ることもできる．

ifelse(論理演算, TRUE の場合の値, FALSE の場合の値) のように指定する．

car$recent <- ifelse(car$year >= 2022, 1, 0)  # = 1 if recent model 
table(car$recent)

## 
##  0  1 
## 39 11

summary(lm(price ~ weight + recent, car))

## 
## Call:
## lm(formula = price ~ weight + recent, data = car)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -198.196  -26.286   -7.163   28.501  150.919 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -167.99503   29.37240  -5.719 7.14e-07 ***
## weight         0.35746    0.02446  14.617  < 2e-16 ***
## recent        33.03056   22.11425   1.494    0.142    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 63.34 on 47 degrees of freedom
## Multiple R-squared:  0.8335, Adjusted R-squared:  0.8264 
## F-statistic: 117.7 on 2 and 47 DF,  p-value: < 2.2e-16

lm(price ~ weight + recent, car)$coef

##  (Intercept)       weight       recent 
## -167.9950335    0.3574592   33.0305631

4.2 Quardatic and interaction terms / 二次・交差項

\[ \mbox{Price} = \alpha + \beta_1 \ \mbox{PS} + \beta_2 \ \mbox{PS}^2 + u \]

car$ps2 <- car$ps^2
lm(price ~ ps + ps2, car)$coef

## (Intercept)          ps         ps2 
##  3.90668080  2.82808263 -0.00307522

上の回帰式は次のような非線形な関係（二乗項の係数が負なので上に凸の二次関数）を表す．

データセットを直接操作せずに，回帰モデルを指定する際に I 関数を使って二乗項を表現することもできる．

lm(price ~ ps + I(ps^2), car)$coef

## (Intercept)          ps     I(ps^2) 
##  3.90668080  2.82808263 -0.00307522

\[ \mbox{Price} = \alpha + \beta_1 \ \mbox{Weight} + \beta_2 \ \mbox{Hybrid} + \beta_3 \ \mbox{Weight} \times \mbox{Hybrid} + u \]

car$wh <- car$weight * car$hybrid
lm(price ~ weight + hybrid + wh, car)$coef

##  (Intercept)       weight       hybrid           wh 
## -148.6298253    0.3441707 -168.3087621    0.1583126

このように，交差項を用いることでカテゴリー（この場合はハイブリッド車とガソリン車（・ディーゼル車））による回帰係数の異質性を表現できる．

\[ \mbox{Price} = -149 + 0.34 \mbox{Weight} - 168 \mbox{Hybrid} + 0.16 \mbox{Weight} \times \mbox{Hybrid} + u \] \[ = \left\{ \begin{array}{l} -317 + 0.50 \mbox{Weight} + u \quad \mbox{if} \quad \mbox{Hybrid} \\ -149 + 0.34 \mbox{Weight} + u \quad \mbox{if} \quad \mbox{Not Hybrid} \end{array} \right. \]

「var1 : var2」は「var1 × var2」を表し，「var1 * var2」は「var1 + var2 + (var1 × var2)」を表す．

lm(price ~ weight + hybrid + weight:hybrid, car)$coef

##   (Intercept)        weight        hybrid weight:hybrid 
##  -148.6298253     0.3441707  -168.3087621     0.1583126

lm(price ~ weight*hybrid, car)$coef

##   (Intercept)        weight        hybrid weight:hybrid 
##  -148.6298253     0.3441707  -168.3087621     0.1583126

4.3 Log-linear model

\[ \log(\mbox{Price}) = \alpha + \beta \ \mbox{Weight} + u \]

説明変数や被説明変数に対数を取る場合は I なしで log を使ってよい．

lm(log(price) ~ weight, car)$coef

## (Intercept)      weight 
## 3.941721733 0.001258482

アウトカムに自然対数を取る場合，「$x$ が1単位大きな値を取る時に $y$ が $(100 \times \beta)$ % 高い値を取る」と解釈される．

\[ \beta = \frac{\partial \log(y)}{\partial x} = \frac{\partial \log(y)}{\partial y} \frac{\partial y}{\partial x} = \frac{1}{y} \frac{\partial y}{\partial x} = \frac{\partial y / y}{\partial x} \]

よって，「1 kg重い自動車の価格は 0.13 %高い」と解釈できる．

なお，上は $\beta$ が0に近い値（たとえば 0.01 や 0.001 くらいのオーダー）を取る場合にはよい近似となる．0から離れた値（たとえば 0.5 またはそれ以上）を取る場合は $\exp(\beta) - 1$ で厳密に計算するとよい．

\[ \Delta \% := \frac{y_{x+1} - y_{x}}{y_{x}} = \frac{\exp(\alpha + \beta (x+1)) - y_{x}}{y_{x}} = \frac{y_{x} \cdot \exp(\beta) - y_{x}}{y_{x}} = \exp(\beta) - 1 \]

4.4 Log-log model

\[ \log(\mbox{Price}) = \alpha + \beta \ \log(\mbox{Weight}) + u \]

lm(log(price) ~ log(weight), car)$coef

## (Intercept) log(weight) 
##   -5.342678    1.533081

両辺に対数を取る場合は「$x$ が1%大きな値を取る時に $y$ が $\beta$ % 高い値を取る」と解釈される．

$\beta$ は経済学では弾力性として知られる．

\[ \beta = \frac{\partial \log(y)}{\partial \log(x)} = \frac{\partial \log(y)}{\partial y} \frac{\partial x}{\partial \log(x)} \frac{\partial y}{\partial x} = \frac{x}{y} \frac{\partial y}{\partial x} = \frac{\partial y / y}{\partial x / x} \]

「1 %重い自動車の価格は 1.5 %高い」と解釈できる．

4.5 Categorical variables

複数の値を持つカテゴリー変数（質的変数）はダミー変数化する（0-1 の量的変数に変換する）ことで回帰式に含めることができる．

いずれか一つのカテゴリーをベースラインとする（すべてのカテゴリーをダミー変数として含めると多重共線性により推定できなくなるため，一つは除外しなければならない）．

maker 変数はカテゴリー変数なので，これをダミー変数にしてみよう．ダイハツをベースライン（reference group）とする．

table(car$maker)

## 
## Daihatsu    Honda    Mazda   Nissan   Subaru   Suzuki   Toyota 
##        7        6        1        7        3       10       16

car$Honda <- 1 * (car$maker == "Honda")
car$Mazda <- 1 * (car$maker == "Mazda")
car$Nissan <- 1 * (car$maker == "Nissan")
car$Subaru <- 1 * (car$maker == "Subaru")
car$Suzuki <- 1 * (car$maker == "Suzuki")
car$Toyota <- 1 * (car$maker == "Toyota")
head(car[, c("maker", "Honda", "Mazda", "Nissan", "Subaru", "Suzuki", "Toyota")])

maker	Honda	Suzuki	Toyota
Honda	1	0	0
Toyota	0	0	1
Daihatsu	0	0	0
Toyota	0	0	1
Suzuki	0	1	0
Suzuki	0	1	0

lm(price ~ weight + Honda + Mazda + Nissan + Subaru + Suzuki + Toyota, car)$coef

##  (Intercept)       weight        Honda        Mazda       Nissan       Subaru 
## -192.7228606    0.4035732  -16.8371837 -147.6797964  -28.2007247  -87.1606913 
##       Suzuki       Toyota 
##    3.0663751  -26.1403619

たとえば，Honda の回帰係数 -17 は，ホンダの車がベースラインであるダイハツの車と比較して 17 万円低いことを表している（ただし統計的に有意ではない）．

lm 関数にカテゴリー変数をそのまま含めると自動的に水準別ダミー変数が作成され推定される．

lm(price ~ weight + maker, car)$coef

##  (Intercept)       weight   makerHonda   makerMazda  makerNissan  makerSubaru 
## -192.7228606    0.4035732  -16.8371837 -147.6797964  -28.2007247  -87.1606913 
##  makerSuzuki  makerToyota 
##    3.0663751  -26.1403619

4.5.1 Change the reference group

ベースラインのカテゴリーを変更するには relevel 関数を使う．

car$maker <- relevel(as.factor(car$maker), ref = "Honda")
lm(price ~ weight + maker, car)$coef

##   (Intercept)        weight makerDaihatsu    makerMazda   makerNissan 
##  -209.5600443     0.4035732    16.8371837  -130.8426127   -11.3635410 
##   makerSubaru   makerSuzuki   makerToyota 
##   -70.3235076    19.9035588    -9.3031782

今度は，ダイハツがベースラインのホンダと比較して 17 万円高いと解釈できる．

5 Prediction / 予測

当てはめ値は定義通り次のように計算できる．

\[ \hat{Y}_i = \hat{\alpha} + \hat{\beta} X_i \]

lm(price ~ weight + ps, car)$coef

##  (Intercept)       weight           ps 
## -103.7104391    0.2367862    0.7764891

car$pred0 <- -103.71 + 0.2368 * car$weight + 0.7765 * car$ps
head(car[, c("model", "price", "pred0")], 3)

model	price	pred0
N-BOX	174	157
Prius	299	294
Tanto	148	145

lm(y ~ x, data)$fitted.value はアウトカム変数 and/or 説明変数の欠損により除外された個体を含まないベクトルのため，元のデータセットと行数が合わないことがある．

たとえば {1,000 kg, 100 ps} という仮想的な車種の価格を予測するには，上の計算式にこの属性値を代入すればよい．

-103.71 + 0.2368 * 1000 + 0.7765 * 100

## [1] 210.74

5.1 `predict` function

推定に用いたサンプルの当てはめ値も，out-of-sample の予測も，predict 関数で計算できる．

予測値を元のデータセットにマージしたい場合，変数に欠損値があると予測値ベクトルのサイズが元のデータセットの行数より小さくなってしまうため，次のように na.action = na.exclude オプションを追加する．

lm1 <- lm(price ~ weight + ps, car, na.action = na.exclude)
car$pred1 <- predict(lm1, newdata = car)
head(car[, c("model", "price", "pred0", "pred1")], 3)

model	price	pred0	pred1
N-BOX	174	157	157
Prius	299	294	294
Tanto	148	145	145

predict(lm1, newdata = data.frame(weight = 1000, ps = 100))

##        1 
## 210.7247

5.2 Inverse transformation of log outcomes

Log-linear model の場合，予測値を元の変数のスケールに戻すには次のようにする．

\[ \hat{y} = \exp \left( \hat{\log{y}} + \frac{\sigma^2}{2} \right) \]

Log-linear model の当てはめ値は $\hat{\log{y}}$ であって $\log{\hat{y}}$ ではないことに注意．
$\sigma^2$ は誤差分散．真の値は未知なので残差の分散（正確には，誤差分散の不偏推定量になるように調整したもの）で代用する．$K$ は定数項を除いた説明変数の数．

\[ \hat{\sigma}^2 = \frac{1}{n - K - 1} \sum_{i = 1}^n \hat{u}_i^2 \quad \mbox{where} \quad \hat{u}_i = Y_i - \hat{Y}_i \]

## correct way 
lm0 <- lm(log(price) ~ weight + ps, car, na.action = na.exclude)
sigma2 <- sum(lm0$residuals^2) / (length(lm0$residuals) - 3)
# sigma2 <- summary(lm2)$sigma^2  # same as above 
car$pred_correct <- exp(predict(lm0, newdata = car) + sigma2/2)
## incorrect way 
car$pred_incorrect <- exp(predict(lm0, newdata = car))
head(car[, c("model", "price", "pred_correct", "pred_incorrect")], 3)

model	price	pred_correct	pred_incorrect
N-BOX	174	162	160
Prius	299	266	263
Tanto	148	156	154

上のような補正を行わないと予測値が過小評価される．このことを次のように平均値で確認する．

car2 <- car[!is.na(car$ps), ]  # 欠損のある車種は推定に使われていないので除外
c(mean(car2$price), mean(car2$pred_correct), mean(car2$pred_incorrect))

## [1] 259.8918 260.2486 257.3103

6 Model selection / モデル選択

被説明変数の予測を目的として回帰分析を行う場合，複数のモデル候補があった場合に，予測力が最も優れたモデルを一つ選択したい場合があるだろう．しかし，説明変数を追加すれば決定係数は必ず増加するが，それは標本に特有のデータの変動（あるいはノイズ）に過適合しているだけかもしれない．この過適合によって説明変数過剰のモデルが選択されることを防ぐために，モデルの説明力とモデルの簡潔さのトレードオフの間でバランスを取るための何らかの指標が必要となる．

6.1 AIC

実務的に広く利用されるのはAIC（赤池情報量規準）と呼ばれる規準．これはモデルを最尤推定することを前提として「最尤推定量の対数尤度」と「パラメタの数（＝定数項を含む説明変数の数＝K+1）」のトレードオフとして記述される．

\[ AIC = -2 \log L + 2(K+1) \]

対数尤度は大きいほどモデルのデータへの当てはまりが良く，principle of parsimony の観点からパラメタ数は少ない方がよい．すなわち，複数のモデル間で AIC が最も小さいモデルが望ましい．

注：

因果推論を目的として回帰分析を行う場合には一般にモデル選択は行わない場合が多い．分析の目的に応じてモデル選択を行うかどうかを決めよう．
- 典型的な AIC の使いどころは時系列分析におけるラグの次数決定．
AICは「真のモデル」を選択する基準ではない．あくまで予測力の観点で（比較対象の中で相対的に最も）望ましいモデルを選択する基準である．
「AICは入れ子関係のモデル間の比較にしか使えない」という主張に稀に出くわすが，それは誤り．入れ子関係のモデル間でないと適用できない一部の検定と混同されているものと思われる．詳しくは Burnham and Anderson (2002) Model Selection and Multimodel Inference (Ch. 2) などを参照．
モデル選択後に回帰係数の検定は行わない．少なくとも，オーソドックスな t 検定をナイーブに適用するのはNG．

データセットが欠損値を含むため na.omit 関数を使って欠損値を含む行を予めすべて削除する．

lm1 <- lm(price ~ weight, na.omit(car))
lm2 <- lm(price ~ weight + ps, na.omit(car))
lm3 <- lm(price ~ weight + ps + length + width, na.omit(car))
AIC(lm1, lm2, lm3)

df	AIC
3	507
4	497
6	498

2つ目のモデルの値が最も小さいが，2つ目と3つ目の差は十分に小さい（2を下回っている）ため，2つ目と3つ目のモデルはほとんど同程度と判断できる．

6.2 Stepwise model selection with AIC

考え得る説明変数をすべて含めたモデル (full model) からスタートして stepwise で最適なモデルを選択することもできる．

step 関数を使うと AIC を規準としてモデル選択が行われる．

lm_full <- lm(price ~ weight + length + width + height + WLTC + 
                ps + disp + hybrid + accel + kei + maker, na.omit(car))
lm_step <- step(lm_full, trace = 0)  # trace = 1 にすると計算過程が出力される
# AIC(lm_step)
summary(lm_step)$coef

##                  Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept)  -587.5679793 4.571500e+02 -1.285285 0.2071329299
## weight          0.2067181 1.238871e-01  1.668601 0.1041196936
## length       -168.6160919 8.209700e+01 -2.053864 0.0475182032
## width         996.5137653 3.534776e+02  2.819171 0.0078717951
## height       -219.0329144 6.892895e+01 -3.177662 0.0030983378
## ps              3.1288058 8.913864e-01  3.510044 0.0012534436
## disp           -0.1232391 4.430405e-02 -2.781667 0.0086526134
## hybrid         92.0504882 2.474967e+01  3.719262 0.0006973619
## accel       -3186.4444056 1.496767e+03 -2.128884 0.0403765717
## kei           147.3452243 5.329504e+01  2.764708 0.0090288881

しかし符号が直感と反する変数が幾つかある． VIFをチェックすると案の定多重共線性が生じている．

car::vif(lm_step)

##     weight     length      width     height         ps       disp     hybrid 
##  55.871573  54.803133  81.432658   2.273168 110.198829  22.175713   1.508772 
##      accel        kei 
##  65.599128  17.280466

6.3 Cross validation / 交差検証

回帰モデルの予測性能（汎化性能）を評価するための手法として交差検証はよく用いられる．交差検証ではデータを「推定用」と「評価用」に分割し，推定用のサブサンプルで推定されたモデルの性能を評価用のサブサンプル（out-of-sample）で評価する．

機械学習の文脈では訓練データ・検証データ・テストデータと区別するが，ここでは検証データで最終評価を行う（テストデータのように扱う）場合を考える．

この演習で使っている car データはサンプルサイズが小さいため，評価用のサブサンプルはサイズ1としよう．

評価指標には MSE (mean squared error) を用いる． MSE の標本対応は次の通り． $\hat{Y}_i$ は $i$ 以外のデータ（すなわち $1,2,\ldots,i-2,i-1,i+1,i+2,\ldots,n$）によって推定されたモデルのパラメタ（$\hat{\alpha}_{-i}, \hat{\beta}_{-i}$）に基づいて予測された $i$ の値である．

\[ \mbox{MSE} = \frac{1}{n} \sum_{i = 1}^n (Y_i - \hat{Y}_i)^2 \quad \mbox{where} \quad \hat{Y}_i = \hat{\alpha}_{-i} + \hat{\beta}_{-i} \ X_i \]

つまり，以下のように Step 1, Step 2, …, Step n を実行し，二乗誤差の平均値を計算する．

Step 1. $i=2,\ldots,n$ を推定用，$i=1$ を評価用として，二乗誤差を計算
Step 2. $i=1,3,\ldots,n$ を推定用，$i=2$ を評価用として，二乗誤差を計算
Step 3. $i=1,2,4,\ldots,n$ を推定用，$i=3$ を評価用として，二乗誤差を計算
Step 4. $i=1,2,3,5,\ldots,n$ を推定用，$i=4$ を評価用として，二乗誤差を計算
…
Step n. $i=1,\ldots,n-1$ を推定用，$i=n$ を評価用として，二乗誤差を計算

このような交差検証を Leave-One-Out CV と呼ぶ．

説明変数が少ないモデル（alternative model）の方が予測誤差が小さい．これは full model で生じている深刻な多重共線性が幾つかの変数削除によって解消されたため（そして，それによる予測力の向上が説明変数を減らすことによる説明力の低下を相殺してなお上回るため）であると考えられる．

car2 <- na.omit(car)
se <- NULL
for (i in 1:nrow(car2)) {
  # Full model 
  lm_full_i <- lm(price ~ weight + length + width + height + WLTC + 
                    ps + disp + hybrid + accel + kei + year, car2[-i,])
  pred_full_i <- as.numeric(predict(lm_full_i, newdata = car2[i, ]))  # hat{y}_i
  se_full_i <- (car2$price[i] - pred_full_i)^2  # (squared error)_i
  # Alternative model 
  lm_alt_i <- lm(price ~ weight + height + ps + hybrid + kei + year, car2[-i,])
  pred_alt_i <- as.numeric(predict(lm_alt_i, newdata = car2[i, ]))
  se_alt_i <- (car2$price[i] - pred_alt_i)^2
  # stack
  se <- rbind(se, c(id = i, price = car2$price[i], 
                    pred_full = pred_full_i, se_full = se_full_i, 
                    pred_alt = pred_alt_i, se_alt = se_alt_i))
}
summary(se[, c(4, 6)])

##     se_full             se_alt        
##  Min.   :    4.12   Min.   :    1.46  
##  1st Qu.:  261.67   1st Qu.:  314.94  
##  Median : 1315.18   Median : 1281.30  
##  Mean   : 3356.59   Mean   : 2916.42  
##  3rd Qu.: 2707.38   3rd Qu.: 2758.49  
##  Max.   :41321.30   Max.   :34592.28

# plot(pred_full ~ price, se, cex = 1.5); abline(a = 0, b = 1, col = "red")
# plot(pred_alt ~ price, se, cex = 1.5); abline(a = 0, b = 1, col = "red")

単位を万円に戻すためにMSEの平方根をとると，Full model では 58 万円，Alternative model では 54 万円となる．

c(sqrt(3356.59), sqrt(2916.42))  # RMSE (Root MSE)

## [1] 57.93609 54.00389

Alternative model は以下のように推定される．

lm_alt <- lm(price ~ weight + height + ps + hybrid + kei + year, car2)
summary(lm_alt)$coef

##                  Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept) -1.711715e+04 4.178485e+03 -4.096495 2.118593e-04
## weight       2.949684e-01 3.916656e-02  7.531128 4.755765e-09
## height      -1.284230e+02 5.828417e+01 -2.203395 3.370089e-02
## ps           8.867450e-01 1.760835e-01  5.035935 1.186617e-05
## hybrid       3.924293e+01 2.183725e+01  1.797064 8.027690e-02
## kei          8.514798e+01 2.203004e+01  3.865086 4.204883e-04
## year         8.469391e+00 2.069523e+00  4.092435 2.144453e-04

car::vif(lm_alt)

##   weight   height       ps   hybrid      kei     year 
## 5.589411 1.626768 4.304064 1.175648 2.955349 1.189786

これでもなお不安な符号が幾つか残っているが（全高 height が負，軽自動車ダミー kei が正），何らかの除外変数によるものなのか，あるいはこの符号で実際に正しいのかは不明．

7 Take home messages

回帰分析：summary(lm(formula = Y ~ X1 + X2, data = dataset))
説明変数の二乗は I(X1^2)
2つの説明変数の掛け算（交差項）は X1:X2

Data Analysis Using Statistical Packages: Regression Analysis

Sho Kuroda / 黒田翔

最終更新：2025年6月