logo

多重回帰分析における残差の分散の推定量と回帰係数の標準誤差 📂統計的分析

多重回帰分析における残差の分散の推定量と回帰係数の標準誤差

定理

[y1y2⋮yn]=[1x11⋯xp11x12⋯xp2⋮⋮⋱⋮1x1n⋯xpn][β0β1⋮βp]+[ε1ε2⋮εn] \begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{bmatrix} = \begin{bmatrix} 1 & x_{11} & \cdots & x_{p1} \\ 1 & x_{12} & \cdots & x_{p2} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{1n} & \cdots & x_{pn} \end{bmatrix} \begin{bmatrix} \beta_{0} \\ \beta_{1} \\ \vdots \\ \beta_{p} \end{bmatrix} + \begin{bmatrix} \varepsilon_{1} \\ \varepsilon_{2} \\ \vdots \\ \varepsilon_{n} \end{bmatrix} pp 個の独立変数とnn 個のデータが与えられた時、線形重回帰モデルを計画行列で表すと上のようになり、簡単にY=Xβ+εY = X \beta + \varepsilon と表そう。回帰係数の推定値はβ^=(XTX)−1XTY\hat{\beta} = \left( X^{T} X \right)^{-1} X^{T} Y であるから、適合値のベクトルY^\hat{Y} は Y^=Xβ^=X(XTX)−1XTY \hat{Y} = X \hat{\beta} = X \left( X^{T} X \right)^{-1} X^{T} Y とわかっている。便宜上P:=X(XTX)−1XTP := X \left( X^{T} X \right)^{-1} X^{T} としよう。一方、残差は線形性を持つこと、つまりε1,⋯ ,εn\varepsilon_{1} , \cdots , \varepsilon_{n} の母平均は00 と仮定する。

残差平方和の期待値

  • [1]: もし残差が等分散性も持っているなら、つまりある定数σ>0\sigma > 0 に対してε1,⋯ ,εn∼(0,σ2)\varepsilon_{1} , \cdots , \varepsilon_{n} \sim \left( 0, \sigma^{2} \right) が成立するなら、SSESSEsum of squared errorの期待値は次のようになる。 E(SSE)=E[∑i=1n(yi−y^i)2]=nσ2−∑i,jE(yiyj)Pij E \left( SSE \right) = E \left[ \sum_{i=1}^{n} \left( y_{i} - \hat{y}_{i} \right)^{2} \right] = n \sigma^{2} - \sum_{i,j} E \left( y_{i} y_{j} \right) P_{ij}

残差平方和の分散に対する不偏推定量

  • [2]: もし残差が独立性も持っているなら、つまりε1,⋯ ,εn∼iid(0,σ2)\varepsilon_{1} , \cdots , \varepsilon_{n} \overset{\text{iid}}{\sim} \left( 0, \sigma^{2} \right) が成立するなら、SSESSE の分散に対する不偏推定量σ2^\widehat{\sigma^{2}}は次のようになる。 Eσ2^=E[1n−p−1∑i=1n(yi−y^i)2]=σ2 E \widehat{\sigma^{2}} = E \left[ {{ 1 } \over { n-p-1 }} \sum_{i=1}^{n} \left( y_{i} - \hat{y}_{i} \right)^{2} \right] = \sigma^{2}

回帰係数の標準誤差

  • [3] もし残差が正規性も持っているなら、つまりε1,⋯ ,εn∼iidN(0,σ2)\varepsilon_{1} , \cdots , \varepsilon_{n} \overset{\text{iid}}{\sim} N \left( 0, \sigma^{2} \right)が成立するなら、回帰係数の標準誤差は次のようになる。 s.e.(β^k)=σ^[(XTX)−1]kk \text{s.e.} \left( \hat{\beta}_{k} \right) = \hat{\sigma} \sqrt{ \left[ \left( X^{T} X \right)^{-1} \right]_{kk} }

説明

ほとんどの統計学部の学生は、学校で回帰分析を初めて学ぶ時には、プロジェクトや他の科目に追われて、こうした数理統計学的な理論展開を大雑把に見過ごす場合が多い。意欲やモチベーションとは無関係に、内容が理解するには難しすぎるので、無理に勉強するのも効率が良くないと思う。証明が最初から理解できないと感じたら、失望せずに立ち去っても大丈夫。

しかし、修士以上で学業を続け、学部科目を復習するなら、ここに上手く整理された内容を見ることを強くお勧めする。重回帰分析のモデル診断で最も重要なのは線形性であり、次に等分散性、その次は独立性、その後が正規性であるが、回帰分析でのt-検定、F-検定を導くためには、正確にその順序で仮定が追加されなければならない。直感や経験からその序列を理解できないかもしれないが、理論学習だけで納得できるというのは、幸いである。

証明 1

戦略: そう簡単ではないかもしれない。数理統計学はともかくとして、最低でも行列代数について十分に学んでいなければならない。定理のステートメントで簡単にP:=X(XTX)−1XTP := X \left( X^{T} X \right)^{-1} X^{T} と示したPP が冪等idempotent、つまり射影作用素であること、つまり P2=X(XTX)−1XT⋅X(XTX)−1XT=X(XTX)−1(XTX)(XTX)−1XT=X(XTX)−1XT=P=PT \begin{align*} P^{2} =& X \left( X^{T} X \right)^{-1} X^{T} \cdot X \left( X^{T} X \right)^{-1} X^{T} \\ =& X \left( X^{T} X \right)^{-1} \left( X^{T} X \right) \left( X^{T} X \right)^{-1} X^{T} \\ =& X \left( X^{T} X \right)^{-1} X^{T} \\ =& P \\ =& P^{T} \end{align*} したがってP2=P=PTP^{2} = P = P^{T} であり、その直交射影作用素 (I−P)(I-P) も射影作用素であるため、(I−P)2=(I−P)\left( I - P \right) ^{2} = \left( I - P \right) が成り立つという事実を補助定理として使う。ここから難しいと感じるなら、無理に今この証明を見ようとせず、数年勉強してから再び戻ることをお勧めする。

[1] 2

クロネッカーデルタδij={1,if i=j0,if i≠j\delta_{ij} = \begin{cases} 1 & , \text{if } i = j \\ 0 & , \text{if } i \ne j \end{cases} について、次のことが成立する： E[∑i=1n(yi−y^i)2]=E[(Y−PY)T(Y−PY)]=E[[(I1+p−P)Y]T[(I1+p−P)Y]]=E[YT(I1+p−P)T(I1+p−P)Y]=E[YT(I1+p−P)(I1+p−P)Y]=E[YT(I1+p−P)2Y]=E[YT(I1+p−P)Y]=E[∑i,jyiyj(δij−Pij)]=∑i,jE[yiyjδij]−∑i,jE[yiyjPij]=∑iE[yi2]−∑i,jE[yiyj]Pij=nσ2−∑i,jE[yiyj]Pij \begin{align*} E \left[ \sum_{i=1}^{n} \left( y_{i} - \hat{y}_{i} \right)^{2} \right] =& E \left[ \left( Y - P Y \right)^{T} \left( Y - P Y \right) \right] \\ =& E \left[ \left[ \left( I_{1+p} - P \right) Y \right] ^{T} \left[ \left( I_{1+p} - P \right) Y \right] \right] \\ =& E \left[ Y^{T} \left( I_{1+p} - P \right)^{T} \left( I_{1+p} - P \right) Y \right] \\ =& E \left[ Y^{T} \left( I_{1+p} - P \right) \left( I_{1+p} - P \right) Y \right] \\ =& E \left[ Y^{T} \left( I_{1+p} - P \right)^{2} Y \right] \\ =& E \left[ Y^{T} \left( I_{1+p} - P \right) Y \right] \\ =& E \left[ \sum_{i,j} y_{i} y_{j} \left( \delta_{ij} - P_{ij} \right) \right] \\ =& \sum_{i,j} E \left[ y_{i} y_{j} \delta_{ij} \right] - \sum_{i,j} E \left[ y_{i} y_{j} P_{ij} \right] \\ =& \sum_{i} E \left[ y_{i}^{2} \right] - \sum_{i,j} E \left[ y_{i} y_{j} \right] P_{ij} \\ =& n \sigma^{2} - \sum_{i,j} E \left[ y_{i} y_{j} \right] P_{ij} \end{align*}

■

[2]

残差が独立ということは、i≠ji \ne jの場合、yiy_{i} とyjy_{j} も相関関係がないということで、i≠ji \ne j の時E[yiyj]=0E \left[ y_{i} y_{j} \right] = 0 であり、i=ji = j の時、残差の線形性と等分散性によりE[yiyj]=σ2E \left[ y_{i} y_{j} \right] = \sigma^{2} だから、次を得る： E[∑i=1n(yi−y^i)2]=nσ2−∑i,jE[yiyj]Pij=nσ2−∑iσ2Pii \begin{align*} E \left[ \sum_{i=1}^{n} \left( y_{i} - \hat{y}_{i} \right)^{2} \right] =& n \sigma^{2} - \sum_{i,j} E \left[ y_{i} y_{j} \right] P_{ij} \\ =& n \sigma^{2} - \sum_{i} \sigma^{2} P_{ii} \end{align*}

トレースの巡回性： Tr(ABC)=Tr(BCA)=Tr(CAB) \text{Tr}(ABC) = \text{Tr}(BCA) = \text{Tr}(CAB)

∑iPii\sum_{i} P_{ii} はPPのトレースtrP\text{tr} P だから、 E[∑i=1n(yi−y^i)2]=nσ2−σ2∑iPii=σ2(n−trP)=σ2(n−trX(XTX)−1XT)=σ2(n−trXTX(XTX)−1)=σ2(n−trI1+p)=σ2(n−(1+p)) \begin{align*} E \left[ \sum_{i=1}^{n} \left( y_{i} - \hat{y}_{i} \right)^{2} \right] =& n \sigma^{2} - \sigma^{2} \sum_{i} P_{ii} \\ =& \sigma^{2} \left( n - \text{tr} P \right) \\ =& \sigma^{2} \left( n - \text{tr} X \left( X^{T} X \right)^{-1} X^{T} \right) \\ =& \sigma^{2} \left( n - \text{tr} X^{T} X \left( X^{T} X \right)^{-1} \right) \\ =& \sigma^{2} \left( n - \text{tr} I_{1+p} \right) \\ =& \sigma^{2} \left( n - (1+p) \right) \end{align*} を得る。両辺を(n−p−1)(n-p-1) で割ると、 1n−p−1E[∑i=1n(yi−y^i)2]=σ2 {{ 1 } \over { n-p-1 }} E \left[ \sum_{i=1}^{n} \left( y_{i} - \hat{y}_{i} \right)^{2} \right] = \sigma^{2} だから、σ2\sigma^{2} の不偏推定量σ2^=∑(yi−y^i)2/(n−p−1)\widehat{\sigma^{2}} = \sum \left( y_{i} - \hat{y}_{i} \right)^{2} / (n-p-1)を得る。

■

[3]

回帰係数ベクトルの多変量正規性： β^∼N1+p(β,σ2(XTX)−1) \hat{\beta} \sim N_{1+p} \left( \beta , \sigma^{2} \left( X^{T} X \right)^{-1} \right)

残差がiidで正規分布に従う場合、β^=(β^0,⋯ ,β^p)\hat{\beta} = \left( \hat{\beta}_{0} , \cdots , \hat{\beta}_{p} \right)のkk番目の成分β^k\hat{\beta}_{k}のマージナル確率分布も以下のような単変量正規分布に従う： β^k∼N(βk,σ2[(XTX)−1]kk) \hat{\beta}_{k} \sim N \left( \beta_{k} , \sigma^{2} \left[ \left( X^{T} X \right)^{-1} \right]_{kk} \right)

標準誤差の一般的な定義：ある推定量estimator TTに対し、TT の標準偏差の推定値estimateを標準誤差standard errorという。 s.e.(T):=Var⁡(T)^ \text{s.e.} \left( T \right) := \sqrt{ \widehat{ \operatorname{Var} \left( T \right) } }

Var⁡β^k=σ2[(XTX)−1]kk\operatorname{Var} \hat{\beta}_{k} = \sigma^{2} \left[ \left( X^{T} X \right)^{-1} \right]_{kk} だから、次を得る。 s.e.(β^k)=Var⁡(β^k)^=σ2[(XTX)−1]kk^=1n−p−1∑i=1n(yi−y^i)2[(XTX)−1]kk=σ^[(XTX)−1]kk \begin{align*} \text{s.e.} \left( \hat{\beta}_{k} \right) =& \sqrt{ \widehat{ \operatorname{Var} \left( \hat{\beta}_{k} \right) } } \\ =& \sqrt{ \widehat{\sigma^{2} \left[ \left( X^{T} X \right)^{-1} \right]_{kk} } } \\ =& \sqrt{ {{ 1 } \over { n-p-1 }} \sum_{i=1}^{n} \left( y_{i} - \hat{y}_{i} \right)^{2} \left[ \left( X^{T} X \right)^{-1} \right]_{kk} } \\ =& \hat{\sigma} \sqrt{ \left[ \left( X^{T} X \right)^{-1} \right]_{kk} } \end{align*}

■

  1. Hadi. (2006). Regression Analysis by Example(4th Edition): p82~84. ↩︎

  2. https://math.stackexchange.com/questions/2277739/why-is-the-denominator-n-p-1-in-estimation-of-variance ↩︎