【発展】単回帰分析と標本回帰係数の母分散

 単回帰分析の標本回帰係数の母分散を、Σを用いて求めます。数式が多いです。

要約

(1)問い

 「単回帰分析の標本回帰係数の母分散」を数式で求めると、どうなるでしょうか?

(2)要約

 結論としては、下の気持ち悪い数式になります。

$$ 標本回帰係数の母分散Var (\widehat{\beta_1})$$

$$= \frac{ \sum\limits_{i=1}^n (x_i-\overline{x})^2 Var(U_i) }{ \left(\sum\limits_{i=1}^n (x_i-\overline{x})^2 \right) ^2} +2 \frac{ \sum\limits_{i=1}^n\sum\limits_{j=1}^n (x_i-\overline{x})(x_j-\overline{x})Cov(U_i,U_j) }{ \left(\sum\limits_{i=1}^n (x_i-\overline{x})^2 \right) ^2 }$$

$$ただし i≠j$$

(3)着眼点

 繰り返し期待値の法則が突破口です。

前提

(1)仮定する条件

 この記事では

  • 仮定1:データ(X,Y)は独立で同一の分布に従い、単回帰モデルに従う

$$単回帰モデル:Y_i=\beta_0 +\beta_1 X_i +誤差項U_i$$

$$i:データ番号(1〜n)、n:サンプル・サイズ$$

  • 仮定2:標本が無作為抽出
  • 仮定3:説明変数Xの分散が0でない

に加えて

$$仮定4:外生性:E(U_i|X_i)=0$$

(2)仮定しない条件

 以下の条件は仮定しません

  • 系列相関なし
  • 均一分散
  • 誤差項の正規性

(3)事実

最小二乗推定量

 仮定1〜3が満たされ、説明変数Xが確定しているとき、単回帰分析の最小二乗推定量は

$$\widehat{\beta_1}=\beta_1+ \frac{\sum\limits_{i=1}^n(x_i-\overline{x})U_i } {\sum\limits_{i=1}^n(x_i-\overline{x})^2 }・・・式(1)$$

となります。詳しくは「単回帰分析/ 最小二乗推定量とパラメーターの関係」をご覧ください。

不偏性

 外生性が仮定されるとき、不偏性

$$E( \widehat{\beta_1}) =\beta_1 ・・・式(2)$$

が成り立ちます。詳しくは単回帰分析と不偏性をご覧ください。

最小二乗推定量には分散がある

 標本抽出できたデータによって、様々な推定値が計算されます。したがって、推定量は確率変数で、分散(Variance、 Var)が存在します。

 例えば、下図はパラメーターは1で生成した標本2つです。それぞれの推定値は0.898や1.130と、ばらつきがあります。

 標本を1000個生成して、それぞれで推定値を計算すると、次のような分布になりました。最小二乗推定量には分散があるのです。

分散と共分散の定義

 確率変数Xの分散とは

$$分散Var(X)=E\left[ (X-E(X))^2 \right]$$

です。確率変数X、Yの共分散とは

$$共分散Cov(X,Y)=E \left[(X-E(X)) (Y-E(Y)) \right]$$

です。

結果

最小二乗推定量の母分散

 最小二乗推定量の母分散は

$$ Var (\widehat{\beta_1}) = E \left[(\widehat{\beta_1}- E(\beta_1))^2 \right]$$

ですが、式(2)の不偏性をもつので

$$= E \left[(\widehat{\beta_1}- \beta_1)^2 \right]・・・式(3)$$

です。

①E( )の中身の計算

 式(1)より

$$\widehat{\beta_1}-\beta_1=\frac{\sum\limits_{i=1}^n(x_i-\overline{x})U_i } {\sum\limits_{i=1}^n(x_i-\overline{x})^2 }$$

になります。したがって、

$$(\widehat{\beta_1}-\beta_1)^2=\left( \frac{\sum\limits_{i=1}^n(x_i-\overline{x})U_i } {\sum\limits_{i=1}^n(x_i-\overline{x})^2 } \right)^2・・・式(4)$$

 次をCとおくと

$$C =\sum\limits_{i=1}^n(x_i-\overline{x})^2 $$

 式(4)は

$$(\widehat{\beta_1}-\beta_1)^2=\left( \frac {\sum\limits_{i=1}^n (x_i-\overline{x}) U_i }{C} \right)^2$$

Cはiによって変化しないので

$$=\left( \sum\limits_{i=1}^n\frac{x_i-\overline{x}}{C}U_i \right)^2・・・式(5)$$

となります。

 式(5)の内部を展開すると

$$\left( \sum\limits_{i=1}^n\frac{x_i-\overline{x}}{C}U_i \right)^2$$

$$=(\frac{x_1-\overline{x}}{C})^2U_1^2 +(\frac{x_2-\overline{X}}{C})^2U_2^2+\cdots+(\frac{x_n-\overline{x}}{C})^2U_n^2$$

$$+2(\frac{x_1-\overline{x}}{C})(\frac{x_2-\overline{x}}{C})U_1U_2$$

$$+2(\frac{x_1-\overline{x}}{C})(\frac{x_3-\overline{x}}{C})U_1U_3$$

$$+\cdots$$

$$+2(\frac{x_{n-1}-\overline{x}}{C})(\frac{x_n-\overline{x}}{C})U_{n-1}U_n・・・式(6)$$

になります。これはuの二乗部分の前半と、異なるuの積の後半に分けられます。

②中身をE( )に入れる

 式(6)を式(3)に代入すると

$$ 標本回帰係数の母分散Var (\widehat{\beta_1})$$

$$=E \Bigg[ (\frac{x_1-\overline{x}}{C})^2U_1^2 +\cdots+(\frac{x_n-\overline{x}}{C})^2 U_n^2$$

$$+2(\frac{x_1-\overline{x}}{C})(\frac{x_2-\overline{x}}{C})U_1U_2+\cdots+2(\frac{x_{n-1}-\overline{x}}{C})(\frac{x_n-\overline{x}}{C})U_{n-1}U_n  \Bigg]$$

 xは確定しているのでEの外に出せて

$$=(\frac{x_1-\overline{x}}{C})^2 E (U_1^2) +\cdots+(\frac{x_n-\overline{x}}{C})^2 E(U_n^2)$$

$$+2(\frac{x_1-\overline{x}}{C})(\frac{x_2-\overline{x}}{C})E(U_1U_2)+\cdots+2(\frac{x_{n-1}-\overline{x}}{C})(\frac{x_n-\overline{x}}{C})E(U_{n-1}U_n) $$

シグマによる表現

 上式は次のように前半や後半に分けられます。

$$ 標本回帰係数の母分散Var (\widehat{\beta_1})$$

$$= \frac{ \sum\limits_{i=1}^n (x_i-\overline{x})^2 E (U_i^2) }{ \left(\sum\limits_{i=1}^n (x_i-\overline{x})^2 \right) ^2} +2 \frac{ \sum\limits_{i=1}^n\sum\limits_{j=1}^n (x_i-\overline{x})(x_j-\overline{x})E(U_iU_j) }{ \left(\sum\limits_{i=1}^n (x_i-\overline{x})^2 \right) ^2 }$$

$$ただし、i≠j $$

 ΣΣについては次のように考えてください。

$$\sum\limits_{i=1}^n \sum\limits_{j=1}^n A_i B_j$$

$$=A_1B_1+A_2 B_1+\cdots+A_nB_1$$

$$+A_1B_2+A_2 B_2+\cdots+A_nB_2$$

$$\cdots$$

$$+A_1B_n+A_2 B_n+\cdots+A_nB_n$$

分散、共分散による表現

 外生性を仮定しているので

$$E (U_i)=0$$

です。(←証明は繰り返し期待値の法則を使う。E(U)=E(E(U|X))=E(0)=0)したがって、

$$E (U_i^2)=E \left[(U_i-E(U_i)^2 \right]=Var(U_i)・・・誤差項の分散$$

$$E(U_iU_j)=E \left[(U_i-E(U_i)) (U_j-E(U_j)) \right]=Cov(U_i,U_j)・・・誤差項の共分散$$

と表せます。したがって、結論としては

$$ 標本回帰係数の母分散Var (\widehat{\beta_1})$$

$$= \frac{ \sum\limits_{i=1}^n (x_i-\overline{x})^2 Var(U_i) }{ \left(\sum\limits_{i=1}^n (x_i-\overline{x})^2 \right) ^2} +2 \frac{ \sum\limits_{i=1}^n\sum\limits_{j=1}^n (x_i-\overline{x})(x_j-\overline{x})Cov(U_i,U_j) }{ \left(\sum\limits_{i=1}^n (x_i-\overline{x})^2 \right) ^2 }$$

$$ただし、i≠j $$

です。

考察

(1)長い数式

 結果は次の式でした。

$$ 標本回帰係数の母分散Var (\widehat{\beta_1})$$

$$= \frac{ \sum\limits_{i=1}^n (x_i-\overline{x})^2 Var(U_i) }{ \left(\sum\limits_{i=1}^n (x_i-\overline{x})^2 \right) ^2} +2 \frac{ \sum\limits_{i=1}^n\sum\limits_{j=1}^n (x_i-\overline{x})(x_j-\overline{x})Cov(U_i,U_j) }{ \left(\sum\limits_{i=1}^n (x_i-\overline{x})^2 \right) ^2 }$$

$$ただし、i≠j $$

です。

(2)外的な意義

 標準誤差について理解する出発点となります。

付録:Rコード

 シミュレーションは、R言語を用いて行いました。R言語については「しまうまのRでデータ分析入門」をご覧ください。

 標本を作成します。↓

trial_number <- 1000 #試行回数
sample_size <- 300 #サンプル・サイズ
beta_1 <- rep(NA, trial_number) #推定値

se_classical <-rep(NA, trial_number) #標準誤差
x <- rnorm(n=sample_size, mean =15, sd=3) #xの分布
u <- rep(NA, sample_size)

#{}を繰り返す。標本を試行回数作る。
for(j in 1:trial_number){ 
  
#{}を繰り返す。誤差項をサンプル・サイズ作る。冗長なコードだが、不均一分散に拡張できるようにしてある
  for(i in 1:sample_size){ 
    xx <- x[i]
    uu <- rnorm(n=1, mean =0, sd=2) #均一分散
    u[i] <- uu
  }
  
  y <- 2 + 1*x+ u
  answer <- lm(y~x)
  
  beta_1[j] <- summary(answer)$coef["x","Estimate"]
}

 回帰分析と推定量の標準偏差を計算します。なお、回帰分析は、最後に生成されたサンプルについて行います。↓

lm(y ~ x) #回帰分析
sd(beta_1) #推定量の標準偏差

 ggplot2用にデータフレームを作成します。↓

category <- 
  ifelse(x <= 10,"~10",
         ifelse(x <= 15 ,"10~15",
                ifelse(x <= 20, "15~20","20~")))

beta <- data.frame(beta_1)
data0 <- data.frame(x,y,u,category) 

 ggplot2で作図します。なお、散布図で作図されるのは、最後に生成されたサンプルです。↓

#パッケージの呼び出し。未インストールならインストール!
library(ggplot2) 

#散布図
ggplot(data0, aes(x = x,y=y)) + 
  geom_point(aes(colour=category) ) +
  scale_x_continuous(limits = c(0, 30),breaks = seq(0,30,5))+
  scale_y_continuous(limits = c(0, 30),breaks = seq(0,30,5))+
  stat_smooth(method = "lm", formula='y~x', se = FALSE)

#ヒストグラム
ggplot(beta, aes(x = beta_1)) +
  geom_histogram(colour="white", fill = "skyblue2",bins=50)+
  scale_x_continuous(limits = c(0.75, 1.25))

コメント欄 お気軽にコメントをお寄せください!

タイトルとURLをコピーしました