重回帰分析

公開

2026年7月1日

更新日

2026年7月2日

複数の説明変数を用いる重回帰分析について学びます。

パッケージの読み込み

library(tidyverse)
library(car)      # VIF の計算

データの準備

単回帰分析で用いた学習時間とテスト得点の例に，睡眠時間と模試の事前得点を説明変数として追加する（n = 80）。

set.seed(123)
n <- 80

df <- tibble(
  study_hours  = runif(n, min = 0, max = 8),
  sleep_hours  = runif(n, min = 4, max = 9),
  prior_score  = rnorm(n, mean = 55, sd = 10),
  noise        = rnorm(n, mean = 0, sd = 6)
) |>
  mutate(
    test_score = 20 + 5 * study_hours + 2 * sleep_hours +
      0.4 * prior_score + noise,
    test_score = pmin(test_score, 100)
  ) |>
  select(-noise)

summary(df)

  study_hours        sleep_hours     prior_score      test_score    
 Min.   :0.004998   Min.   :4.052   Min.   :34.47   Min.   : 44.02  
 1st Qu.:2.088011   1st Qu.:5.212   1st Qu.:48.56   1st Qu.: 63.63  
 Median :3.765116   Median :6.279   Median :54.05   Median : 75.84  
 Mean   :4.011478   Mean   :6.465   Mean   :54.51   Mean   : 74.64  
 3rd Qu.:6.043770   3rd Qu.:7.617   3rd Qu.:60.74   3rd Qu.: 84.99  
 Max.   :7.954158   Max.   :8.925   Max.   :76.87   Max.   :100.00

単回帰から重回帰へ

モデル

説明変数が 1 つの単回帰モデルを拡張し，$k$ 個の説明変数 $x_1, \dots, x_k$ で目的変数 $y$ を説明する。

\[ y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \cdots + \beta_k x_{ki} + \varepsilon_i, \quad \varepsilon_i \sim N(0, \sigma^2) \]

$\beta_j$ は「他の説明変数を一定に保ったときの，$x_j$ が 1 単位増えることによる $y$ の変化量」を表す偏回帰係数である。

最小二乗法（OLS）

単回帰と同様に，残差二乗和を最小にする係数を求める。

\[ \hat{\boldsymbol{\beta}} = \arg\min_{\boldsymbol{\beta}} \sum_{i=1}^{n} \varepsilon_i^2 = \arg\min_{\boldsymbol{\beta}} \sum_{i=1}^{n} \left( y_i - \beta_0 - \sum_{j=1}^{k} \beta_j x_{ji} \right)^2 \]

行列表記では，計画行列 $X$（1 列目は切片用の 1）を用いて次の閉じた解が得られる。

\[ \hat{\boldsymbol{\beta}} = (X^\top X)^{-1} X^\top \mathbf{y} \]

説明変数間に強い相関がある（多重共線性が高い）場合，$X^\top X$ がほぼ特異となり係数の推定が不安定になる点に注意する。

`lm()` による推定

model <- lm(test_score ~ study_hours + sleep_hours + prior_score, data = df)
summary(model)


Call:
lm(formula = test_score ~ study_hours + sleep_hours + prior_score, 
    data = df)

Residuals:
   Min     1Q Median     3Q    Max 
-8.140 -4.047 -0.688  3.178 13.537 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  18.8086     4.3744   4.300 5.02e-05 ***
study_hours   4.8953     0.2628  18.626  < 2e-16 ***
sleep_hours   1.7347     0.4246   4.085 0.000108 ***
prior_score   0.4582     0.0628   7.296 2.42e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.277 on 76 degrees of freedom
Multiple R-squared:  0.8543,    Adjusted R-squared:  0.8486 
F-statistic: 148.6 on 3 and 76 DF,  p-value: < 2.2e-16

出力の読み方

表 1: 回帰係数の推定結果

項目	Estimate	Std. Error	t value	Pr(>\|t\|)
切片 β̂₀	18.8086	4.3744	4.2997	1e-04
学習時間 β̂₁	4.8953	0.2628	18.6260	0e+00
睡眠時間 β̂₂	1.7347	0.4246	4.0852	1e-04
事前得点 β̂₃	0.4582	0.0628	7.2963	0e+00

単回帰の出力との違いは，各説明変数の係数が「他の変数を一定に保ったときの」効果として推定される点である。

自由度調整済み決定係数

説明変数を増やすと $R^2$ は単調に増加してしまうため，変数の数で罰則をかけた自由度調整済み決定係数で当てはまりを比較する。

\[ R^2_{\text{adj}} = 1 - (1 - R^2) \frac{n - 1}{n - k - 1} \]

r2     <- summary(model)$r.squared
r2_adj <- summary(model)$adj.r.squared
cat(sprintf("R²       = %.3f\n", r2))

R²       = 0.854

cat(sprintf("調整済み R² = %.3f\n", r2_adj))

調整済み R² = 0.849

多重共線性の確認（VIF）

説明変数間の相関が強いと係数の標準誤差が膨らみ，推定が不安定になる。VIF（分散拡大係数）で確認する。

\[ \text{VIF}_j = \frac{1}{1 - R_j^2} \]

ここで $R_j^2$ は，$x_j$ を他のすべての説明変数で回帰したときの決定係数である。

表 2: 各説明変数の VIF

vif(model)

study_hours sleep_hours prior_score 
   1.026433    1.027571    1.003329

VIF の目安	解釈
< 5	問題なし
5〜10	やや注意
> 10	多重共線性が強く疑われる

モデル比較とモデル選択

説明変数の組み合わせを変えたモデルを AIC（赤池情報量規準） で比較する。値が小さいほど良いモデルとされる。

model_simple <- lm(test_score ~ study_hours, data = df)
model_full   <- lm(test_score ~ study_hours + sleep_hours + prior_score, data = df)

tibble(
  モデル = c("単回帰（study_hours のみ）", "重回帰（全変数）"),
  AIC    = c(AIC(model_simple), AIC(model_full)),
  `調整済みR²` = c(summary(model_simple)$adj.r.squared,
                summary(model_full)$adj.r.squared)
) |>
  knitr::kable(digits = 3)

モデル	AIC	調整済みR²
単回帰（study_hours のみ）	548.958	0.711
重回帰（全変数）	499.069	0.849

step() を使うと，AIC を基準に変数増減法で自動的にモデルを選択できる。

model_step <- step(model_full, direction = "both")
formula(model_step)

予測

new_data <- tibble(
  study_hours = c(2, 4, 6),
  sleep_hours = c(6, 7, 7.5),
  prior_score = c(50, 55, 60)
)
predict(model, newdata = new_data, interval = "confidence")

       fit      lwr      upr
1 61.91856 60.22936 63.60777
2 75.73503 74.47514 76.99493
3 88.68413 86.83826 90.53000

predict(model, newdata = new_data, interval = "confidence") |>
  round(2) |>
  cbind(new_data) |>
  select(study_hours, sleep_hours, prior_score, fit, lwr, upr) |>
  knitr::kable(col.names = c("学習時間", "睡眠時間", "事前得点",
                       "予測得点", "下限 (95%)", "上限 (95%)"))

学習時間	睡眠時間	事前得点	予測得点	下限 (95%)	上限 (95%)
2	6.0	50	61.92	60.23	63.61
4	7.0	55	75.74	74.48	76.99
6	7.5	60	88.68	86.84	90.53

まとめ

表 3: 重回帰分析で使用した主な関数

関数	役割
lm(y ~ x1 + x2 + …, data)	重回帰モデルの OLS 推定
summary(model)	偏回帰係数・R²・F 統計量の出力
summary(model)$adj.r.squared	自由度調整済み決定係数の取得
vif(model)	多重共線性（VIF）の確認
AIC(model)	モデルの AIC の計算
step(model, direction = ‘both’)	AIC によるステップワイズ変数選択
predict(model, newdata, interval = ‘confidence’)	予測値と信頼区間の計算

単回帰と同様に，重回帰でも残差診断（plot(model)）で線形性・等分散性・正規性・独立性を確認することが望ましい。