日々のつれづれ

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

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

R 勉強会

R Advent Calender 2023として投稿することにしました

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

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

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

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

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

$$ 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$ : 誤差項

この回帰モデルにおいて、異なる説明変数が目的変数に与える影響を評価するために分散分析を利用します。そして、回帰モデルや説明変数により目的変数の総変動（ばらつき）が減少（説明）する量を見ることで、その影響が分かります。

これには、以下の 2 つのアプローチがあります。

  • 回帰モデル全体の有意性評価：回帰モデルが目的変数をどれだけ説明できるかを見るために、回帰モデルと単なる平均値（切片だけのモデル）を比較する。
  • 特定の説明変数や切片に対する検定：各説明変数や切片が目的変数にどれだけ影響を与えるかを見るために、その係数が 0 であるかどうかを検定する。

重回帰モデルの場合

重回帰モデルは説明変数が複数あるため、「回帰モデル全体の有意性評価」と「特定の説明変数や切片に対する検定」の 2 つのアプローチは異なります。

回帰モデル全体の有意性評価の計算手順

ここでは、回帰モデル全体の有意性評価を整理します。この評価では、回帰モデルが目的変数を説明できる程度を見ています。

次の計算手順は、サンプル数が $m$、説明変数の総数が $n$の例になります。

  1. 目的変数の総変動を回帰変動と残差変動に分解する
    • $SS_y$ : 目的変数の総変動、目的変数のばらつきの程度
      • $y_i$ : $i$ 番目の実測値
      • $\bar{y}$ : 目的変数の平均値 $$ SS_y = \sum_{i = 1}^{n}(y_i - \bar{y})2 $$
    • $SS_R$ : 回帰変動、回帰モデルが目的変数を説明できる程度
      • $\hat{y}i$ : $i$ 番目の予測値 $$ SS_R = \sum{i = 1}^{n}(\hat{y}_i - \bar{y})2 $$
    • $SS_E$ : 残差変動、回帰モデルが目的変数を説明できない程度 $$ SS_E = \sum{i = 1}^{n}(y_i - \hat{y}i)2 $$
  2. それぞれの変動を自由度で割って平均二乗を求める
    • $df_y$ : 全体の自由度、観測値の数 - 1 $$ df_y = m - 1 $$
    • $df_R$ : 回帰の自由度、説明変数の個数 $$ df_R = n $$
    • $df_E$ : 残差の自由度、全体の自由度 - 回帰の自由度 $$ df_E = df_y - df_R = m - n - 1 $$
    • $MS_y$ : 全体の平均二乗、総変動を全体の自由度で割った値 $$ MS_y = \frac{SS_y}{df_y} $$
    • $MS_R$ : 回帰の平均二乗、回帰変動を回帰の自由度で割った値 $$ MS_R = \frac{SS_R}{df_R} $$
    • $MS_E$ : 残差の平均二乗、残差変動を残差の自由度で割った値 $$ MS_E = \frac{SS_E}{df_E} $$
  3. 回帰の平均二乗を残差の平均二乗で割って、F 統計量$F_R$を求める $$ F_R = \frac{MS_R}{MS_E} $$
  4. 回帰と残差の自由度の F 分布から p 値$p_R$を求めて有意水準で判断する $$ p_R = P(F > F_R) = 1 - P(F \leq F_R) $$ なお、このときの帰無仮説は
    • 帰無仮説 ( $H_0$ ): $\beta_1 = 0$（説明変数は目的変数に影響しない）
    • 対立仮説 ( $H_1$ ): $\beta_1 \neq 0$（説明変数は目的変数に影響する）

R で計算する

回帰分析で分散分析すると、異なる要因（説明変数、因子など）が目的変数に与える影響を統計的に評価できます。そして、これはデータに基づいた合理的な意思決定に有効です。

説明変数が 2 つの重回帰モデルを例に確認します。式は次のようになります。

できるだけスクラッチ実装をしてみます。

$$ y{mpg} = \beta_0 + \beta_1 x{wt} + \beta{2} x{hp} + \varepsilon $$

これは、mtcarsデータセットを使った回帰モデルで、目的変数が mpg（ガロンあたりのマイル数）、説明変数が wt（車両重量）と hp（馬力）です。そして、wt と hp に交互作用はないと考えています。

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

# 回帰モデルの作成
model <- lm(mpg ~ wt + hp, data = mtcars)
model

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

Coefficients:
(Intercept)           wt           hp
   37.22727     -3.87783     -0.03177

係数と切片の推定値が、それぞれ $\hat{\beta}0 = 37.22727$ 、$\hat{\beta}1 = -3.87783$ 、$\hat{\beta}_2 = -0.03177$ になりました。

続いて、回帰モデル全体の有意性を評価します。計算の流れは単回帰モデルと同じです。

  1. 目的変数の総変動 $SS_y$ を回帰変動 $SS_R$ と残差変動 $SS_E$ に分解する 回帰変動は、回帰モデルの目的変数の平均値からのズレの程度です。良い回帰モデルほど値は大きくなります。
# 目的変数とその推定値
y <- mtcars$mpg # 目的変数
bar_y <- mean(y) # 平均値
hat_y <- fitted(model) # 推定値
# hat_y <- predict(model) # 推定値

# 総変動と回帰変動と残差変動
SS_y <- sum((y - bar_y)^2) # 総変動
SS_R <- sum((hat_y - bar_y)^2) # 回帰変動
SS_E <- sum((y - hat_y)^2) # 残差変動
  1. それぞれの自由度で割って、それぞれの平均二乗 $MS_R$ と $MS_E$ を求める
# 自由度
df_y <- nrow(mtcars) - 1 # 全体の自由度
df_R <- length(coef(model)) - 1 # 回帰の自由度（説明変数の数）
df_E <- df_y - df_R # 残差の自由度

# 平均二乗、分散
MS_y <- SS_y / df_y # 全体の平均二乗
MS_R <- SS_R / df_R # 回帰の平均二乗
MS_E <- SS_E / df_E # 残差の平均二乗
  1. 回帰の平均二乗を残差の平均二乗で割って、F 統計量を計算する F 統計量は、回帰モデルが目的変数をどれだけ説明できるかを表す指標です。回帰モデルの当てはまりが良ければ、F 統計量は大きくなります。
# F統計量の計算
F_R <- MS_R / MS_E
  1. 自由度が（$df_R$, $df_E$）の F 分布から p 値を求めて有意水準で判断する p 値は、F 統計量が得られる確率を表す指標です。回帰モデルが目的変数に有意な影響を与えていれば、p 値が小さくなります。
# p値の計算
p_R <- pf(F_R, df_R, df_E, lower.tail = FALSE)
# p_R <- 1 - pf(F_R, df_R, df_E)
  1. 分散分析表を作成する
# 分散分析表の作成
res <- data.frame(
   Df = c(df_R, df_E, df_y),
   `Sum Sq` = c(SS_R, SS_E, SS_y),
   `Mean Sq` = c(MS_R, MS_E, MS_y),
   `F value` = c(F_R, NA, NA),
   `Pr(>F)` = c(p_R, NA, NA),
   row.names = c("Regression", "Residuals", "Total"),
   check.names = FALSE
)
res

           Df    Sum Sq    Mean Sq  F value       Pr(>F)
Regression  2  930.9994 465.499716 69.21121 9.109054e-12
Residuals  29  195.0478   6.725785       NA           NA
Total      31 1126.0472  36.324103       NA           NA

なお、これはanova(model)では取得できません。

summary(model)の返り値を使って計算します。

f_stat = summary(model)$fstatistic
p_value = pf(f_stat[1], f_stat[2], f_stat[3], lower.tail=FALSE)
cat("F statistic: ", f_stat[1], "\n")
cat("p value: ", p_value, "\n")

F statistic:  69.21121
p value:  9.109054e-12