回帰分析 - モデルの評価

講義の予定

• 第1日: 回帰モデルの考え方と推定
• 第2日: モデルの評価
• 第3日: モデルによる予測と発展的なモデル

回帰分析の復習

線形回帰モデル

• 目的変数 を 説明変数 で説明する関係式を構成:
  • 説明変数: \(x_1,\dotsc,x_p\) (p次元)
  • 目的変数: \(y\) (1次元)

• 回帰係数 \(\beta_0,\beta_1,\dotsc,\beta_p\) を用いた一次式:

  \begin{equation} y=\beta_0+\beta_1x_1+\dotsb+\beta_px_p \end{equation}

• 誤差項 を含む確率モデルで観測データを表現:

  \begin{equation} y_i=\beta_0+\beta_1 x_{i1}+\cdots+\beta_px_{ip}+\epsilon_i \quad (i=1,\dotsc,n) \end{equation}

行列・ベクトルによる簡潔な表現

• デザイン行列:

  \begin{equation} \boldsymbol{X}= \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{pmatrix} \end{equation}

行列・ベクトルによる簡潔な表現

• ベクトル:

  \begin{equation} \boldsymbol{y}= \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix},\quad \boldsymbol{\epsilon}= \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{pmatrix},\quad \boldsymbol{\beta}= \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{pmatrix} \end{equation}

問題の記述

• 確率モデル:

  \begin{equation} \boldsymbol{y} =\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\epsilon} \end{equation}

• 回帰式の評価: 残差平方和 の最小化による推定

  \begin{equation} S(\boldsymbol{\beta}) =(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta})^\top (\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}) \end{equation}

解の表現

• 解の条件: 正規方程式

  \begin{equation} \boldsymbol{X}^\top\boldsymbol{X}\boldsymbol{\beta} =\boldsymbol{X}^\top\boldsymbol{y} \end{equation}

• 解の一意性: Gram 行列 \(\boldsymbol{X}^\top\boldsymbol{X}\) が正則

  \begin{equation} \boldsymbol{\hat{\beta}} = (\boldsymbol{X}^\top\boldsymbol{X})^{-1} \boldsymbol{X}^\top\boldsymbol{y} \end{equation}

最小二乗推定量の性質

• あてはめ値 \(\boldsymbol{\hat{y}}=\boldsymbol{X}\boldsymbol{\hat{\beta}}\) は \(X\) の列ベクトルの線形結合

• 残差 \(\boldsymbol{\hat{\epsilon}}=\boldsymbol{y}-\boldsymbol{\hat{y}}\) はあてはめ値 \(\boldsymbol{\hat{y}}\) と直交する

  \begin{equation} \boldsymbol{\hat{\epsilon}}\cdot\boldsymbol{\hat{y}} =0 \end{equation}

• 回帰式は説明変数と目的変数の 標本平均 を通る

  \begin{equation} \bar{y} = (1,\bar{\boldsymbol{x}}^\top)\boldsymbol{\hat{\beta}}, \quad \bar{\boldsymbol{x}} =\frac{1}{n}\sum_{i=1}^n\boldsymbol{x}_i, \quad \bar{y} =\frac{1}{n}\sum_{i=1}^ny_i, \end{equation}

寄与率

• 決定係数 (R-squared):

  \begin{equation} R^2 = 1-\frac{\sum_{i=1}^n\hat{\epsilon}_i^2}{\sum_{i=1}^n(y_i-\bar{y})^2} \end{equation}

• 自由度調整済み決定係数 (adjusted R-squared):

  \begin{equation} \bar{R}^2 = 1-\frac{\frac{1}{n{-}p{-}1}\sum_{i=1}^n\hat{\epsilon}_i^2} {\frac{1}{n{-}1}\sum_{i=1}^n(y_i-\bar{y})^2} \end{equation}

  (不偏分散で補正)

残差の統計的性質

あてはめ値と誤差の関係

• (A)あてはめ値は観測値の重み付けの和で表される
• (B)あてはめ値と観測値は誤差項の寄与のみ異なる

残差と誤差の関係

• (A)残差は誤差の重み付けの和で表される

ハット行列

• 定義:

  \begin{equation} H= \boldsymbol{X}(\boldsymbol{X}^{\top}\boldsymbol{X})^{-1}\boldsymbol{X}^{\top} \end{equation}
• ハット行列 \(H\) の性質:

  • あてはめ値や残差は \(H\) を用いて簡潔に表現される

    \begin{align} \boldsymbol{\hat{y}} &=H\boldsymbol{y}\\ \boldsymbol{\hat{\epsilon}} &=(I-H)\boldsymbol{\epsilon} \end{align}
  • 観測データの説明変数の関係を表す
  • 対角成分(テコ比 (leverage))は観測データが自身の予測に及ぼす影響の度合を表す

推定量の統計的性質

推定量の性質

• 推定量と誤差の関係

  \begin{align} \boldsymbol{\hat{\beta}} &=(\boldsymbol{X}^{\top}\boldsymbol{X})^{-1}\boldsymbol{X}^{\top}\boldsymbol{y}\\ &\qquad( \boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\epsilon} を代入) \\ &=(\boldsymbol{X}^{\top}\boldsymbol{X})^{-1}\boldsymbol{X}^{\top}\boldsymbol{X}\boldsymbol{\beta} +(\boldsymbol{X}^{\top}\boldsymbol{X})^{-1}\boldsymbol{X}^{\top}\boldsymbol{\epsilon}\\ &=\boldsymbol{\beta} +(\boldsymbol{X}^{\top}\boldsymbol{X})^{-1}\boldsymbol{X}^{\top}\boldsymbol{\epsilon} \end{align}

• 正規分布の重要な性質:

  正規分布に従う独立な確率変数の和は正規分布に従う

推定量の分布

• 誤差の仮定: 平均0，分散 \(\sigma^{2}\) の正規分布に従う

• 推定量は以下の多変量正規分布に従う:

  \begin{align} \mathbb{E}[\boldsymbol{\hat{\beta}}] &=\boldsymbol{\beta}\\ \mathrm{Cov}(\boldsymbol{\hat{\beta}}) &=\sigma^{2}(\boldsymbol{X}^{\top}\boldsymbol{X})^{-1} \end{align}

• 通常 \(\sigma^{2}\) は未知，必要な場合には不偏分散で代用

  \begin{equation} \hat{\sigma^{2}} =\frac{S}{n{-}p{-}1} =\frac{1}{n{-}p{-}1}\boldsymbol{\hat{\epsilon}}^{\top}\boldsymbol{\hat{\epsilon}} =\frac{1}{n{-}p{-}1}\sum_{i=1}^n\hat{\epsilon}_i^2 \end{equation}

分析の評価

寄与率 (再掲)

• 決定係数 (R-squared):
  (回帰式で説明できるばらつきの比率)

  \begin{equation} R^2 = 1-\frac{\sum_{i=1}^n\hat{\epsilon}_i^2}{\sum_{i=1}^n(y_i-\bar{y})^2} \end{equation}

• 自由度調整済み決定係数 (adjusted R-squared):
  (決定係数を不偏分散で補正)

  \begin{equation} \bar{R}^2 = 1-\frac{\frac{1}{n{-}p{-}1}\sum_{i=1}^n\hat{\epsilon}_i^2} {\frac{1}{n{-}1}\sum_{i=1}^n(y_i-\bar{y})^2} \end{equation}

標準誤差

• 推定されたパラメータの精度を評価:
  • 誤差の分布は平均0, 分散 \(\sigma^2\) の正規分布
  • \(\boldsymbol{\hat{\beta}}\) の分布は 平均 \(\boldsymbol{\beta}\) ， 共分散 \(\sigma^2(\boldsymbol{X}^{\top}\boldsymbol{X})^{-1}\) の \(p+1\) 変量正規分布
  • \(\hat{\beta}_j\) の分布は， 行列 \((\boldsymbol{X}^\top\boldsymbol{X})^{-1}\) の対角成分 を \(\xi_0,\xi_1,\dotsc,\xi_p\) とすると， 平均 \(\beta_j\), 分散 \(\sigma^2\xi_j\) の正規分布
  • 未知パラメータ \(\sigma^{2}\) は不偏分散 \(\hat{\sigma}^{2}\) で推定

• 標準誤差 (standard error): \(\hat{\beta}_j\) の精度の評価指標

  \begin{equation} \hat{\sigma}\sqrt{\xi}_j = \sqrt{\frac{1}{n{-}p{-}1}\sum_{i=1}^n\hat{\epsilon}_i^2} \cdot \sqrt{\xi}_j \end{equation}

演習: 標準誤差

• 05-se.r を確認してみよう

\(t\) -統計量

• 回帰係数の分布に関する定理:
  \(t\) -統計量 は自由度 \(n{-}p{-}1\) の \(t\) 分布に従う:

  \begin{equation} \text{(t-統計量)}\quad t=\frac{\hat{\beta}_j-\beta_j}{\hat{\sigma}\sqrt{\xi_j}} \end{equation}
• 証明には以下の性質を用いる:
  • \(\hat{\sigma}{}^2\) と \(\boldsymbol{\hat{\beta}}\) は独立となる
  • \((\hat{\beta}_j-\beta_j)/(\sigma\sqrt{\xi}_j)\) は標準正規分布に従う
  • \((n{-}p{-}1)\hat{\sigma}^2/\sigma^2=S/\sigma^2\) は自由度 \(n{-}p{-}1\) の \(\chi^{2}\) 分布に従う

\(t\) -統計量による検定

• 回帰係数 \(\beta_j\) が回帰式に寄与するか否かを検定:
  • 帰無仮説: \(\beta_j=0\) (\(t\) -統計量が計算できる)
  • 対立仮説: \(\beta_j\neq0\)

• \(p\) -値: 確率変数の絶対値が \(|t|\) を超える確率

  \begin{equation} \text{($p$-値)} = 2\int_{|t|}^\infty f(x)dx \quad\text{(両側検定)} \end{equation}
  • \(f(x)\) は自由度 \(n{-}p{-}1\) の \(t\) 分布の確率密度関数
  • 帰無仮説 \(\beta_j=0\) が正しければ \(p\) 値は小さくならない

演習: \(t\) -統計量

• 05-tpval.r を確認してみよう

\(F\) -統計量

• ばらつきの比に関する定理:
  \(\beta_1=\dotsb=\beta_p=0\) ならば， \(F\) -統計量 は自由度 \(p,n{-}p{-}1\) の \(F\) 分布に従う

  \begin{equation} \text{(F-統計量)}\quad F= \frac{\frac{1}{p}S_{r}}{\frac{1}{n{-}p{-}1}S} % =\frac{\frac{1}{p}\sum_{i=1}^n(\hat{y}_i-\bar{y})^2} % {\frac{1}{n{-}p{-}1}\sum_{i=1}^n(y_i-\hat{y}_i)^2} =\frac{n{-}p{-}1}{p}\frac{R^2}{1-R^2} \end{equation}
• 証明には以下の性質を用いる:
  • \(S_{r}\) と \(S\) は独立となる
  • \(S_{r}/\sigma^2\) は自由度 \(p\) の \(\chi^{2}\) 分布に従う
  • \(S/\sigma^2\) は自由度 \(n{-}p{-}1\) の \(\chi^{2}\) 分布に従う

\(F\) -統計量を用いた検定

• 説明変数のうち1つでも役に立つか否かを検定:
  • 帰無仮説: \(\beta_1=\dotsb=\beta_p=0\) (\(S_r\) が \(\chi^2\) 分布になる)
  • 対立仮説: \(\exists j\;\beta_j\neq0\)

• \(p\) -値: 確率変数の値が \(F\) を超える確率

  \begin{equation} \text{($p$-値)} = \int_{F}^\infty f(x)dx \quad\text{(片側検定)} \end{equation}
  • \(f(x)\) は自由度 \(p,n{-}p{-}1\) の \(F\) 分布の確率密度関数
  • 帰無仮説 \(\forall j\;\beta_j=0\) が正しければ \(p\) 値は小さくならない

演習: \(F\) -統計量

• 05-fstat.r を確認してみよう

演習

• 先週用いたデータの回帰分析の結果を， 寄与率・標準誤差・ \(t\) -統計量・ \(F\) -統計量の観点から 評価してみよう
  • datasets::airquality
  • datasets::LifeCycleSavings
• 人工データを用いたシミュレーションによって 推定量，\(t\) -統計量，\(F\) -統計量の分布を調べてみよう