日々のつれづれ

不惑をむかえ戸惑いを隠せない男性の独り言

回帰モデルの係数の分散分析のまとめ・その２

R

共分散分析で回帰分析の分散分析が出てきたので、少し整理します。

回帰分析を分散分析する目的

回帰分析の分散分析は、データを意思決定に利用するための方法です。

これを使えば、回帰モデルに用いた異なる要因（説明変数、因子など）が目的変数に与える影響を統計的に評価できます。

一般的な回帰モデルは以下の式で表されます。

$$ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_n x_n + \varepsilon $$

  • $y$ : 目的変数
  • $x_1, x_2, \ldots, x_n$ : 説明変数、要因
  • $\beta_0, \beta_1, \beta_2, \ldots, \beta_n$ : 回帰係数（各説明変数の影響度）
  • $\varepsilon$ : 誤差項

この回帰モデルにおいて、異なる説明変数が目的変数に与える影響を評価するために分散分析を利用します。そして、全ての$x_n$が$y$に与える影響を見る「回帰モデル全体の有意性評価」と、個々の$x_n$が$y$に与える影響を個別にみる「特定の説明変数や切片に対する検定」の 2 つのアプローチがあります。

R で計算する

まずは、説明変数が 1 種類の単回帰モデルを考えます。単回帰モデルの分散分析は、説明変数が目的変数をどれだけ説明できるかを見ています。

$$y = \beta_0 + \beta_1 x + \varepsilon$$

なお、単回帰モデルの場合、説明変数が 1 つしかないので、「回帰モデル全体の有意性評価」と「特定の説明変数や切片に対する検定」は同じです。

ここからは、mtcarsデータセットを使い、目的変数を mpg（ガロンあたりのマイル数）、説明変数を wt（車両重量）にした回帰モデルで進めます。

# データセットの読み込み
data(mtcars)

# モデルの適合
model <- lm(mpg ~ wt, data = mtcars)

model

Call:
lm(formula = mpg ~ wt, data = mtcars)

Coefficients:
(Intercept)           wt
     37.285       -5.344

係数と切片の推定値が、$\hat{\beta}0 = 37.285$ 、$\hat{\beta}1 = -5.344$になりました。

続いて、回帰モデル全体の有意性を評価します。計算の流れは以下の通りです。

  1. 目的変数の総変動 $SS_y$ を回帰変動 $SS_R$ と残差変動 $SS_E$ に分解する
  2. それぞれの自由度で割って、それぞれの平均二乗 $MS_R$ と $MS_E$ を求める
  3. 回帰の平均二乗を残差の平均二乗で割って、F 統計量を計算する
  4. 自由度が（$df_R$, $df_E$）の F 分布から p 値を求めて有意水準で判断する

目的変数の総変動 $SS_y$

目的変数の総変動（$SS_y$）は、目的変数 $y$ の値がその平均値 $\bar{y}$ からどれだけずれているかを表す指標です。

目的変数とその平均値の差の二乗和を求めます。

$$SS_y = \sum_{i=1}^{n}{(y_i - \bar{y})2}$$

# 目的変数とその平均値
y <- mtcars$mpg # 目的変数
bar_y <- mean(y) # 平均値

# 総変動
SS_y <- sum((y - bar_y)^2)

回帰変動 $SS_R$

回帰変動（$SS_R$）は、回帰モデルが目的変数の平均値からどれだけずれているかを表す指標です。

回帰モデルによる目的変数の予測値 $\hat{y}$ とその平均値 $\bar{y}$ の差の二乗和を求めます。

$$SS_R = \sum{i=1}^{n}{(\hat{y}i - \bar{y})2}$$

# モデルの予測値
hat_y <- fitted(model)
# hat_y <- predict(model) # こちらでも可

# 回帰変動
SS_R <- sum((hat_y - bar_y)^2)

残差変動 $SS_E$

残差変動（$SS_E$）は、回帰モデルで説明できない目的変数の残りの変動を表す指標です。

目的変数 $y$ と回帰モデルによる予測値 $\hat{y}$ の差の二乗和を求めます。また、残差変動は総変動から回帰変動を引いた値になっています。

$$SS_E = SS_y - SS_R = \sum{i=1}^{n}{(y_i - \hat{y}i)2}$$

# 残差変動
SS_E <- sum((y - hat_y)^2)
# SS_E <- sum(residuals(model)^2) # こちらでも可

目的変数の総変動 $SS_y$

目的変数の総変動は、目的変数 $y$ の値がその平均値 $\bar{y}$ からどれだけずれているかを表す指標です。総変動は、目的変数とその平均値の差の二乗和で求められます。

$$SS_y = \sum_{i=1}^{n}{(y_i - \bar{y})2}$$

# 目的変数とその平均値
y <- mtcars$mpg # 目的変数
bar_y <- mean(y) # 平均値

# 総変動
SS_y <- sum((y - bar_y)^2)

回帰変動の自由度 $df_R$、残差変動の自由度 $df_E$

回帰変動と残差変動の自由度は、それぞれの平方和を求める際に使ったデータの個数や制約の数によって決まります。

  • 全体の自由度はサンプル数 - 1（$df_y = n - 1$）です。
  • 回帰変動の自由度は回帰係数の個数から 1 を引いた値です。単回帰モデルは回帰係数は $\beta_0$ と $\beta_1$ の 2 つなので、$df_R = 2 - 1 = 1$ です。
  • 残差変動の自由度は観測値の個数から回帰係数の個数を引いた値で、全体の自由度から回帰変動の自由度を引いた値です。今回、観測値は 32 個、回帰係数は 2 つなので、$df_E = df_y - df_R = 32 - 2 = 30$ です。
# 自由度
df_y <- nrow(mtcars) - 1 # 全体の自由度
df_R <- length(coefficients(model)) - 1 # 回帰変動の自由度
df_E <- df_y - df_R # 残差変動の自由度

回帰変動の平均平方 $MS_R$、残差変動の平均平方 $MS_E$

回帰変動と残差変動の平均平方は、それぞれの変動がどれだけ大きいかを表す指標です。

それぞれの平方和を対応する自由度で割って求めます。

$$MS_R = \frac{SS_R}{df_R}$$

$$MS_E = \frac{SS_E}{df_E}$$

# 平均平方
MS_y <- SS_y / df_y # 総変動の平均平方
MS_R <- SS_R / df_R # 回帰変動の平均平方
MS_E <- SS_E / df_E # 残差変動の平均平方

F 統計量と p 値

F 統計量は、回帰モデルが目的変数をどれだけ説明できるかを表す指標です。回帰モデルの当てはまりが良ければ、F 統計量は大きくなります。

回帰変動と残差変動の平均平方の比を求めます。

$$F = \frac{MS_R}{MS_E}$$

p 値は、F 統計量が得られる確率を表す指標です。回帰モデルが目的変数に有意な影響を与えていれば、p 値が小さくなります。

なお、このときの帰無仮説は次になります。

  • 帰無仮説 ( $H_0$ ): $\beta_1 = 0$（説明変数は目的変数に影響しない）
  • 対立仮説 ( $H_1$ ): $\beta_1 \neq 0$（説明変数は目的変数に影響する）
# F 統計量と p 値
F_R <- MS_R / MS_E # F統計量
p_R <- pf(F, df_R, df_E, lower.tail = FALSE) # p 値

分散分析表

回帰変動と誤差変動の平方和（SS）、自由度（DF）、平均平方（MS）、F 統計量（F_R）、および p 値（p_R）を分散分析表にまとて、単回帰モデルを評価します。

res <- data.frame(
  ID = c("Regression", "Residual", "Total"),
  DF = c(df_R, df_E, df_y),
  SS = c(SS_R, SS_E, SS_y),
  MS = c(MS_R, MS_E, MS_y),
  F = c(F_R, NA, NA),
  p = c(p_R, NA, NA)
)
cat("ANOVA table:")
res

ANOVA table:          ID DF        SS         MS        F            p
1 Regression  1  847.7252 847.725250 91.37533 1.293959e-10
2   Residual 30  278.3219   9.277398       NA           NA
3      Total 31 1126.0472  36.324103       NA           NA

F 統計量が 91.37533、p 値が 1.293959e-10 になりました。

p 値は十分に小さいので、この回帰モデルは単なる平均値よりも優れた予測力を持つとなります。

なお、これはanova(model)で同じ結果を取得できます。

anova(model)

Analysis of Variance Table

Response: mpg
          Df Sum Sq Mean Sq F value    Pr(>F)
wt         1 847.73  847.73  91.375 1.294e-10 ***
Residuals 30 278.32    9.28
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1