重回帰分析における外れ値の処理と信頼区間のシミュレーション検証

Sep 27, 2025 Biostatistics

#biostatistics #linear-regression #outlier #simulation #r

背景 (Background)

観察データに基づく統計モデリング、とりわけ重回帰分析（Multiple Linear Regression）は、交絡因子を調整しながら複数の説明変数が目的変数に与える影響を評価するための基盤となる手法です。しかし、現実の臨床データや疫学データには入力ミスや特異な性質を持つ「外れ値（Outlier）」が混入していることが多く、これらは最小二乗推定量（OLS）に過度な影響（てこ比）を与え、モデルの予測精度や統計的推論を歪める危険性があります。

本プロジェクトでは、18歳時の体重（WT18）を目的変数とし、過去の身体指標（9歳時の体重 WT9 と18歳時の身長 HT18）を説明変数とする重回帰モデルを構築しました。その過程で、回帰診断プロット（Diagnostic Plots）を用いて外れ値を特定・除外することによるモデル適合度の改善を定量的に評価しました。さらに、 $M = 10, 000$ 回のモンテカルロシミュレーションを通じて、回帰係数の「95% 信頼区間」が持つ頻度論的な意味（被覆率が理論通り 95% になること）を実証的に検証しました。

回帰モデリングと外れ値の診断

まず、R の alr4 パッケージに含まれる BGSgirls データセットを用いて、初期の重回帰モデルを構築しました。

library(alr4)
data(BGSgirls)
BGSg <- subset(x = BGSgirls, select = c(HT9, WT9, HT18, WT18))

# 初期の重回帰モデル
mreg2 <- lm(WT18 ~ WT9 + HT18, data = BGSg)
summary(mreg2)

初期モデルでは、残差の標準誤差（Residual Standard Error）が $6.032$ 、調整済み決定係数（Adjusted $R^{2}$ ）が $0.5153$ でした。ここで、回帰の前提条件（線形性、正規性、等分散性）を確認するために回帰診断プロットを出力します。

外れ値削除前の回帰診断プロット

診断プロット（特に右下の「Residuals vs Leverage」プロット）から、データ番号 "134" の観測値が極端な外れ値であり、回帰直線の推定に大きな影響（高い Cook の距離）を与えていることが視覚的に確認できました。

外れ値の除外とモデルの改善

特定された外れ値 "134" をデータセットから除外し、再度モデルを適合させました。

# 外れ値の削除
BGSg2 <- BGSg[rownames(BGSg) != "134", ]

# 修正後の重回帰モデル
mreg2_new <- lm(WT18 ~ WT9 + HT18, data = BGSg2)
summary(mreg2_new)

外れ値削除後の回帰診断プロット

結果の比較

外れ値の除外により、モデルの適合度は改善しました：

残差の正規性と等分散性の改善: 診断プロット（図後）を見ると、Q-Q プロットが理論直線（ $y = x$ ）に近づき、残差の分布が 0 を中心によりランダムに散布するようになりました。
適合度指標の向上:
- 残差標準誤差（RSE）が $6.032$ から $4.614$ に減少しました。
- 調整済み決定係数（Adjusted $R^{2}$ ）が $0.5153$ から $0.611$ へ上昇しました。モデルが目的変数のばらつきをより多く説明できるようになりました。
モデル全体の有意性の向上: モデル全体の有意性を表す F 検定の $p$ 値が $1.083 \times 1 0^{- 11}$ から $1.102 \times 1 0^{- 14}$ へと小さくなりました（個々の説明変数の $t$ 検定の $p$ 値とは区別される）。

この結果は、データ解析における前処理（外れ値の適切な処理）が、モデルの予測精度と推論の安定性にとっていかに重要であるかを示しています。

信頼区間のシミュレーション検証

統計学において「95% 信頼区間」とは、「同じ母集団から無作為抽出を繰り返し、そのたびに信頼区間を計算した場合、そのうちの 95% の区間が真のパラメータ（母数）を含む」という頻度論的な性質を指します。これを実感として理解するために、R を用いて $M = 10, 000$ 回のシミュレーションを行いました。

既知の真のパラメータ $β = (1, 2, - 0.5)$ を設定し、ノイズを含む仮想データを $10, 000$ セット生成して毎回回帰分析を行い、 $β_{2}$ の 95% 信頼区間が真の値 $2$ を含むかどうかを判定するプログラムを実装しました。

M <- 10000
n <- 100
beta <- c(1, 2, -0.5)
set.seed(44)

count_in <- 0   # 真の値を含む回数
count_out <- 0  # 真の値を含まない回数

for (m in 1:M) {
  x1 <- runif(n, -3, 3)
  x2 <- runif(n, -3, 3)
  error <- rnorm(n, 0, sqrt(5))
  y <- beta[1] + beta[2]*x1 + beta[3]*x2 + error
  
  reg <- lm(y ~ x1 + x2)
  ci.beta1 <- confint(reg)[2,] # beta[2] (x1) の信頼区間
  
  if (ci.beta1[1] > beta[2] | ci.beta1[2] < beta[2]) {
    count_out <- count_out + 1
  } else {
    count_in <- count_in + 1
  }
}

cat("真の値を含む回数:", count_in, "\n")
cat("真の値を含まない回数:", count_out, "\n")
cat("被覆率:", count_in / M, "\n")

シミュレーション結果

真の値を含む回数: 9513 
真の値を含まない回数: 487 
被覆率: 0.9513

$10, 000$ 回の繰り返しの結果、真のパラメータを含む区間の割合（被覆率: Coverage Probability）は $0.9513$ となり、理論上の確率 $0.95$ と極めてよく一致しました。これにより、信頼区間の数理的な定義を計算機実験を通して数値的に確認することができました。

まとめ

本プロジェクトでは、重回帰モデルの構築を通じて、外れ値の検出とその除外がモデルの適合度と予測精度にもたらす影響を実践的に評価しました。さらに、モンテカルロシミュレーションにより、推測統計の根幹である「信頼区間の被覆率」を数値的に確認しました。

これらの分析手法とシミュレーション技術は、観測データのノイズやバイアスに対処しつつ、統計モデルの信頼性を確保するための基盤となります。これは、のちに因果推論や臨床試験データの解析において交絡因子を適切にモデル化する際にも不可欠な素養です。

ソースコードとデータ

本解析で使用した R コードは、以下の GitHub リポジトリにて公開している：

GitHub

multiple-linear-regression-analysis

重回帰分析における外れ値の診断と、信頼区間被覆率のモンテカルロシミュレーションの R コード。