単回帰分析では、目的変数と説明変数の間に直線的な関係があるかを調べる。Rのlm関数を使えば簡単に分析できるが、その裏側にある計算ロジックを理解することも重要だ。この記事では、lm関数に頼らずに回帰直線の主要な要素、つまり回帰係数の推定、その検定、そして95%信頼区間の計算をRスクリプトで実装する。さらに、これらの手動計算の結果をlm関数の出力と比較し、その精度を確認していく。
単回帰直線の基本モデル
単回帰分析のモデルは、以下に示す直線の方程式で表される。
$$Y_i = \beta_0 + \beta_1 X_i + \epsilon_i$$
ここで、
- $Y_i$: 目的変数（従属変数）
- $X_i$: 説明変数（独立変数）
- $\beta_0$: 切片。$X=0$のときの$Y$の期待値。
- $\beta_1$: 回帰係数。$X$が1単位増加したときの$Y$の変化量を示し、回帰直線の傾きを表す。
- $\epsilon_i$: 誤差項。観測値と回帰直線とのずれで、平均0、分散$\sigma^2$の正規分布に従うと仮定される。
このモデルから、私たちはデータのパターンに最もフィットする回帰直線を見つけ、その傾き（$\beta_1$）が統計的に意味のあるものか（つまり0ではないか）を判断し、その信頼できる範囲（信頼区間）を求める。
回帰直線の傾きと切片の推定（最小二乗法）
回帰直線をデータに最もよくフィットさせるには、最小二乗法を用いる。これは、実際の$Y$の値と回帰直線から予測される$\hat{Y}$の値との差（残差）の二乗和が最小になるように、$\beta_0$と$\beta_1$を推定する方法だ。
- 回帰係数($\beta_1$)の推定値 $\hat{\beta}_1$
$$\hat{\beta}_1 = \frac{\sum_{i=1}^{n}(X_i – \bar{X})(Y_i – \bar{Y})}{\sum_{i=1}^{n}(X_i – \bar{X})^2}$$
- 切片($\beta_0$)の推定値 $\hat{\beta}_0$
$$\hat{\beta}_0 = \bar{Y} – \hat{\beta}_1 \bar{X}$$
ここで、$\bar{X}$と$\bar{Y}$はそれぞれ$X$と$Y$の平均値である。
回帰係数に対する標準誤差の計算
回帰係数の推定値の精度を示すのが標準誤差だ。この値は、推定された回帰係数のばらつきの度合いを表し、検定や信頼区間の計算に不可欠となる。
- 残差平方和 (RSS)
$$RSS = \sum_{i=1}^{n}(Y_i – \hat{Y}_i)^2 = \sum_{i=1}^{n}(Y_i – (\hat{\beta}_0 + \hat{\beta}_1 X_i))^2$$
- 残差分散の推定値 ($\hat{\sigma}^2$)
$$\hat{\sigma}^2 = \frac{RSS}{n-2}$$
- 回帰係数($\hat{\beta}_1$)の標準誤差 $SE(\hat{\beta}_1)$
$$SE(\hat{\beta}_1) = \sqrt{\frac{\hat{\sigma}^2}{\sum_{i=1}^{n}(X_i – \bar{X})^2}}$$
回帰直線の傾きに関する検定（t検定）
回帰直線の傾き、つまり回帰係数$\beta_1$が統計的に0と異なるかどうかを判断するためには、t検定を用いる。帰無仮説は「$\beta_1 = 0$」（$X$と$Y$の間に線形な関係がない）、対立仮説は「$\beta_1 \ne 0$」（$X$と$Y$の間に線形な関係がある）となる。
- t値
$$t = \frac{\hat{\beta}_1 – 0}{SE(\hat{\beta}_1)}$$
このt値は自由度$n-2$のt分布に従う。このt値からp値を計算し、有意水準（通常は0.05）と比較することで、帰無仮説を棄却するかどうかを判断する。
回帰直線の傾きの95%信頼区間
推定された回帰係数$\hat{\beta}_1$が、真の回帰係数$\beta_1$をどの程度の範囲で含むかを95%信頼区間として示す。これは、もし同じような分析を100個のサンプルで、100回繰り返した場合、そのうち95回は真の回帰係数がこの区間内に収まるという可能性がある区間を意味する（手持ちのデータは、100 個のサンプルのうちの一つということになる。真の回帰係数は、不明だが一つなので、いくつも存在しうるわけではないことに注意）
$$\hat{\beta}_1 \pm t_{\alpha/2, n-2} \times SE(\hat{\beta}_1)$$
ここで、$t_{\alpha/2, n-2}$は自由度$n-2$のt分布における上側$\alpha/2$点（95%信頼区間では$\alpha=0.05$なので上側0.025点）である。
Rスクリプトでの実装と lm との比較
それでは、上記で解説した計算をRスクリプトで実装し、Rの標準関数であるlmの出力と比較してみよう。
# データの準備
set.seed(123)
X <- 1:100 + rnorm(100, 0, 10)
Y <- 2 * X + 5 + rnorm(100, 0, 20)
data_df <- data.frame(X = X, Y = Y)
# 1. 最小二乗法による回帰係数の推定
n <- length(X)
mean_X <- mean(X)
mean_Y <- mean(Y)
print(c(n, mean_X, mean_Y))
# β1 (傾き) の推定
beta1_hat <- sum((X - mean_X) * (Y - mean_Y)) / sum((X - mean_X)^2)
print(beta1_hat)
# β0 (切片) の推定
beta0_hat <- mean_Y - beta1_hat * mean_X
print(beta0_hat)
cat("--- 手動計算による回帰直線の推定値 ---\n")
cat(sprintf("切片 (beta0_hat): %.4f\n", beta0_hat))
cat(sprintf("傾き (beta1_hat): %.4f\n\n", beta1_hat))
# 2. 残差分散の推定と標準誤差の計算
# 予測値 (回帰直線上のYの値)
Y_pred <- beta0_hat + beta1_hat * X
# 残差平方和 (RSS)
RSS <- sum((Y - Y_pred)^2)
# 残差分散の推定 (MSE)
sigma_squared_hat <- RSS / (n - 2)
# 回帰係数 (beta1) の標準誤差
se_beta1 <- sqrt(sigma_squared_hat / sum((X - mean_X)^2))
cat("--- 手動計算による回帰直線の傾きの標準誤差 ---\n")
cat(sprintf("傾き (beta1) の標準誤差: %.4f\n\n", se_beta1))
# 3. 回帰直線の傾きに関する検定 (t検定)
# t値
t_value_beta1 <- beta1_hat / se_beta1
# 自由度
df <- n - 2
# p値 (両側検定)
p_value_beta1 <- 2 * pt(abs(t_value_beta1), df, lower.tail = FALSE)
cat("--- 手動計算による回帰直線の傾きの検定結果 ---\n")
cat(sprintf("傾き (beta1) のt値: %.4f\n", t_value_beta1))
cat(sprintf("自由度: %d\n", df))
ifelse(p_value_beta1 < 0.001, "p値: <0.001", cat(sprintf("p値: %.3f\n\n", p_value_beta1)))
# 4. 95%信頼区間の計算
# 有意水準 α = 0.05
alpha <- 0.05
# t分布の上側 α/2 点
t_critical <- qt(1 - alpha/2, df)
# 信頼区間の下限
conf_int_lower <- beta1_hat - t_critical * se_beta1
# 信頼区間の上限
conf_int_upper <- beta1_hat + t_critical * se_beta1
cat("--- 手動計算による95%信頼区間 ---\n")
cat(sprintf("95%%信頼区間: [%.4f, %.4f]\n\n", conf_int_lower, conf_int_upper))
# lm関数による結果との比較
cat("--- lm関数による回帰分析結果 ---\n")
lm_model <- lm(Y ~ X, data = data_df)
print(summary(lm_model))
confint(lm_model)
実行結果：
> cat("--- 手動計算による回帰直線の推定値 ---\n")
--- 手動計算による回帰直線の推定値 ---
> cat(sprintf("切片 (beta0_hat): %.4f\n", beta0_hat))
切片 (beta0_hat): 0.4487
> cat(sprintf("傾き (beta1_hat): %.4f\n\n", beta1_hat))
傾き (beta1_hat): 2.0467
> cat("--- 手動計算による回帰直線の傾きの標準誤差 ---\n")
--- 手動計算による回帰直線の傾きの標準誤差 ---
> cat(sprintf("傾き (beta1) の標準誤差: %.4f\n\n", se_beta1))
傾き (beta1) の標準誤差: 0.0626
> cat("--- 手動計算による回帰直線の傾きの検定結果 ---\n")
--- 手動計算による回帰直線の傾きの検定結果 ---
> cat(sprintf("傾き (beta1) のt値: %.4f\n", t_value_beta1))
傾き (beta1) のt値: 32.6750
> cat(sprintf("自由度: %d\n", df))
自由度: 98
> if (p_value_beta1 < 0.001) {
+ cat("p値: <0.001\n\n")
+ } else {
+ cat(sprintf("p値: %.3f\n\n", p_value_beta1))
+ }
p値: <0.001
> cat("--- 手動計算による95%信頼区間 ---\n")
--- 手動計算による95%信頼区間 ---
> cat(sprintf("95%%信頼区間: [%.4f, %.4f]\n\n", conf_int_lower, conf_int_upper))
95%信頼区間: [1.9224, 2.1710]
> cat("--- lm関数による回帰分析結果 ---\n")
--- lm関数による回帰分析結果 ---
> lm_model <- lm(Y ~ X, data = data_df)
> summary(lm_model)
Call:
lm(formula = Y ~ X, data = data_df)
Residuals:
Min 1Q Median 3Q Max
-38.532 -13.507 -1.954 11.070 66.859
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.44871 3.75826 0.119 0.905
X 2.04670 0.06264 32.675 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 19.38 on 98 degrees of freedom
Multiple R-squared: 0.9159, Adjusted R-squared: 0.9151
F-statistic: 1068 on 1 and 98 DF, p-value: < 2.2e-16
> confint(lm_model)
2.5 % 97.5 %
(Intercept) -7.009437 7.906858
X 1.922393 2.170999
上記の実行結果を見ると、手動で計算した推定値、標準誤差、t値、p値、95%信頼区間が、lm関数のsummary()およびconfint()関数が出力する結果とほぼ完全に一致していることがわかる。わずかな小数点以下の違いは、浮動小数点計算の丸め誤差によるものだ。
まとめ
本記事では、Rのlm関数を使わずに、単回帰分析における回帰直線の検定（具体的には傾きである回帰係数の検定）と信頼区間の計算をステップバイステップで実装した。これにより、lm関数が内部で行っている統計的計算のロジックを深く理解することができたのではないだろうか。
統計ソフトウェアは非常に便利なツールだが、その裏側の原理を理解することで、より適切に分析結果を解釈し、データ分析に自信を持って取り組むことができるようになるだろう。今回の解説が、回帰直線の理解を深める一助となれば幸いだ。
コメント