13章_1

本書は「続2・高校生の統計学」の続編です。
高校の探究的学習から大学の統計学へ繋ぐ教材としてご利用ください。

1 はじめに

＜回帰分析について＞.
　単回帰分析は、高校数学B、情報Ⅰで扱われています。しかし、回帰の数理モデルの理解には行列や多変量正規分布の知識が必要なため、高校での扱いは最小二乗法の紹介にとどまっています。
　本書のねらいは、シミュレーションによる実証的な学習により、回帰分析の本質を理解しようというものです。高校生の方は、数理の難しい箇所は公式の確認にとどめ、シミュレーション結果が公式に合致しているかどうかに注目してください。公式が成り立つことをシミュレーションで確かめておけば、統計ソフトで実際の回帰分析を行う際、自信を持って分析を行うことができるでしょう。

＜ひょっこり方式＞.
　従来の学習は、「数理的な証明のあと活用を学ぶ」というものでした。本書が提案する学習法（ひょっこり方式）は、シミュレーションを活用して体験的に理解し、数理をあとにしようというものです。
　将来の再学習を前提とするため、数理は行列を用い簡潔に説明しています。大学で線形代数と微積分を学んだあとで本書に立ち戻れば、更に理解が深まり重回帰分析へと進むことができるでしょう。行列による説明は重回帰分析への拡張が容易です。

皆さんが統計学を深く理解し、様々な場面で活用されることを願います。

2 回帰と相関

回帰と相関は2変量 x,y の関係ですが、その違いは何なのでしょうか。
市原清志先生は「バイオサイエンスの統計学」のなかで、その違いについて説明しておられます。

	回帰	相関
誤　差	y方向のみ	x,y両方が確率変数
信頼域	yの信頼区間	x,yの分布の中心からの等確率距離
目　的	どのように直線に関係づけられるか	どの程度直線的か
単　位	単位に依存	単位に依存しない

本書では、確率モデルによって生成したデータを分析します。
確率モデルは、目的とする調査・研究によって変ります。
そこで、回帰分析と相関分析の目的を明確にし、適用する確率モデルを考えます。

2.1 相関分析のモデル.

相関分析は、2つのデータの間に関係性があるかどうかを相関係数を用いて調べます。
簡単のため正規乱数を用いることにすると、相関分析の確率モデルは２変量正規分布です。
平均を mean、分散共分散行列を sigma とすると、

　　　rmvnorm(n,mean,sigma)　によって n 個の２変量データを生成することができます。

library(mvtnorm)
sigma <- matrix(c(1,0.5,0.5,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)  
data <- data.frame(X=x[,1],Y=x[,2])
model <- lm(Y ~ X, data = data)
plot(x,asp=1,col="lightblue")
abline(model,col="red")
abline(v=1.5,col="blue")

相関分析のモデルとした２変量正規分布では、(x,y) は2変量の確率変数であり、x,y 各々も確率変数です。

一方、回帰分析は 2変量データを直線に当てはめ、xからyを予測しようというものです。
xを青の縦線に固定したとき、予測される y は楕円の弦の中点であるので、回帰直線は楕円の縦の弦の中点の軌跡（赤の直線）になります。

2.2 回帰分析のモデル

回帰分析では、説明変数 xから目的変数 yを予測します。
関数でいうと、xは独立変数、yは従属変数であり、yはxの一次式で表されます。
関数との違いは、誤差変動が加わることです。

回帰分析では、説明変数は定数として固定して考えます。
真の回帰式（母回帰式）として、
　　　\(y=\alpha+\beta x\)　を考え、

固定した説明変数 \(x_i:i=1,2,\cdots,n\) のもとで、
誤差変動 \(\varepsilon_i \sim N(0,~\sigma^2) ~:~ i.i.d.\) を加えた

　　　\(y_i=\alpha+\beta x_i+\varepsilon_i\)　　を単回帰モデルとします。

＜具体例＞

　　説明変数：　\(X=(x_i)=(0,1,2,3,4)\)
　　母回帰式：　\(y=1+0.5 x\)
　　誤差変動：　\(\varepsilon_i \sim N(0,~0.5^2)\)

として y を出力します。

X <- c(0,1,2,3,4)
e <- rnorm(5,0,0.5)
Y <- 0.5*X+1+e
data <- data.frame(X=X,Y=Y);data

X	Y
0	1.282598
1	3.020207
2	2.419595
3	2.367939
4	2.675535

２変量データ (x,y) をプロットし、母回帰直線を書き加えると、

plot(data,ylim=c(0,4))
curve(0.5*x+1,add=T)

誤差項があるため プロットした点は母回帰直線から少しずれています。

今は、確率モデルにより真の回帰直線が \(y=1+0.5 x\)　 と分かっています。
しかし、分析の場面では、この５つの点をもとに真の回帰直線を推定しなければなりません。
推定方法は最小二乗法といい、Rでは関数 lm( ) によって回帰直線を推定することができます。
（最小二乗法の数理については後述します）
実際に、lm( )を用いて回帰直線を推定し、赤で書き加えます。
（矢内勇生先生「統計学２」によれば、lm( )は linear methodの略です）.

fit = lm( formula = Y ~ X, data = data )
plot(data,ylim=c(0,4))
curve(0.5*x+1,add=T)
abline(a =fit$ coefficients[1], b =fit$ coefficients[2], col = "red")

推定した回帰直線（赤）は、真の回帰直線（黒）からずれています。
理由は、５つの点はたまたま得られた標本であり、推定した回帰直線の傾き及び切片は誤差を含んでいるからです。
推定した回帰直線を 標本回帰直線 といいます。

2.3 標本回帰直線.

私たちは「平均には標本平均と母平均があり、標本平均の期待値が母平均になる」ことを学んでいます。
手元にあるデータは、母集団からたまたま得られた標本であり、標本平均は変動する確率変数です。

回帰分析のために得られたデータも、母集団からたまたま得られたものであり、その背後に無限の母集団があることを想定しなければなりません。

母回帰式：\(y=1+0.5 x\) に誤差変動を加えた確率モデル（母回帰モデル）から選られたデータは１つの標本に過ぎず、そこから得られる回帰直線は母回帰直線に一致している訳ではありません。標本から選られたのは標本回帰直線なのです。

しかし、無限回データ収集を行えば、得られる回帰直線は母回帰直線に近いはずです。

　　「母回帰モデルのもとでデータを生成し、回帰直線を推定する」.

無限回とはいきませんが、100回のシミュレーションを行います。
即ち、前述の確率モデルによるデータ生成を100回繰り返します。データは rbind( )関数でつないでいきます。

DATA=data.frame(X=c(),Y=c())
X <- c(0,1,2,3,4)
for (i in 1:1000){
  e <- rnorm(5,0,0.5)
  Y <- 0.5*X+1+e
  data <- data.frame(X=X,Y=Y)
  DATA <- rbind(DATA, data)
}

100回の調査を繋いだ DATA に対し lm( ) によって当てはまる model を求め、回帰直線を abline( ) によって赤で表示します。母回帰直線は黒で表示します。

model = lm( formula = Y ~ X, data = DATA )
plot(DATA,asp=1,col="grey")
abline(model, col = "red")
abline(a =1, b =0.5)

表示した標本回帰直線（赤）は、母回帰直線 \(y=1+0.5 x\) （黒）とほぼ一致していることが分かります。

3 正規方程式

標本回帰直線の係数を、標本回帰係数と言います。

回帰モデル \(y=a+bx\) による y の予測値と、標本による y の観測値の、差の平方和が最小となるように係数 a, b を定めます。この方法を最小二乗法といいます。係数 a, b を求める連立方程式が得られ、得られた連立方程式を正規方程式といいます。

3.1 最小二乗法

簡単のため、２点に当てはまる回帰直線を最小二乗法で計算します。

X <- c(0,1)
Y <- c(1,2)
data <- data.frame(X=X,Y=Y);data

X	Y
0	1
1	2

plot(data,ylim=c(0,3))

２点 (0,1) (1,2) に最も当てはまる直線は、その２点を通る直線 \(y=1+x\) そのものだと想像できます。
この自明な例で、最小二乗法が有効に働くことを確かめます。

最小二乗法で２点\((x_1,y_1),~(x_2,y_2)\) に最もよく当てはまる直線を推定するには、

\(y=a+bx\) による予測値 \(a+bx_1,~a+bx_2\) と実測値 \(y_1,~y_2\) との差の平方和

　　　\(Q=(y_1-a-bx_1)^2+(y_2-a-bx_2)^2\).

を最小とする a,b の値を求めます。
微分を学んでいれば、是非、微分を用いて a,b を求める式を求めてください。
文字で計算すると構造が分かり、一般化が容易になります。
aで微分するときは bは定数として扱い、bで微分するときは aは定数として扱います。（偏微分といいます）.

　　　\(2a+(x_1+x_2)b=y_1+y_2\)
　　　\((x_1+x_2)a+(x_1^2+x_2^2)b=x_1y_1+x_2y_2\)
　　正規方程式（タイプ１）と呼ぶことにします。

各式の両辺を 2 で割ると更にシンプルな式になります。

　　　\(1a+\bar{x}b=\bar{y}\)
　　　\(\bar{x}a+\bar{x^2}b=\bar{xy}\)
　　正規方程式（タイプ２）と呼ぶことにします。

正規方程式は標本回帰係数を求める連立方程式です。

今、\((x_1,y_1)=(0,1),~(x_2,y_2)=(1,2)\) であるので、

　　　\(Q=(1-a)^2+(2-a-b)^2\)　　\(\therefore a=1,~b=1\)　

微分を用いなくても、標本回帰直線は \(y=1+x\) であることが容易にわかります。
また、正規方程式（タイプ１）は、

　　　\(2a+b=3\)
　　　\(~~a+b=2\).
これを解くと \(a=1,~b=1\) であり、標本回帰直線　\(y=1+x\)　が求まります。

＜探求課題＞
３点 \((x_1,y_1),~(x_2,y_2),~(x_3,y_3)\) について正規方程式を求め、２点の場合と同じになることを確かめましょう。

3.2 行列表現

現行の高校教育課程では行列は扱われていませんが、旧の教育課程（1985年～1996年）では「行列と一次変換」が「代数・幾何」という文理共通科目で扱われていました。
行列は多変量のデータ分析に欠かせないツールだからです。
少し背伸びして、行列による表現を見ておきましょう。高校生の方も行列の美しさに触れておくと、学習するときのモチベーションになることでしょう。

正規方程式は、
　　　\(\begin{pmatrix} 1 & \overline{x} \\ \overline{x} & \overline{x^2} \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \overline{y} \\ \overline{xy} \end{pmatrix}\)

　　　\(\begin{pmatrix} 1 & \overline{x} \\ \overline{x} & \overline{x^2} \end{pmatrix}^{-1} = \dfrac{1}{s_{x}^2} \begin{pmatrix} \overline{x^2} & \ \ -\overline{x} \\ -\overline{x} & \ \ \ \ 1 \end{pmatrix}\)　より、　\(\begin{pmatrix} a \\ b \end{pmatrix} =\dfrac{1}{s_{x}^2} \begin{pmatrix} \overline{x^2} & \ \ -\overline{x} \\ -\overline{x} & \ \ \ \ 1 \end{pmatrix} \begin{pmatrix} \overline{y} \\ \overline{xy} \end{pmatrix}\)

　\(\begin{pmatrix} a \\ b \end{pmatrix} = \dfrac{1}{s_{x}^2} \begin{pmatrix} \overline{x^2}\overline{y}-\overline{x}\cdot \overline{xy} \\ -\overline{x}\overline{y}+\overline{xy} \end{pmatrix} = \dfrac{1}{s_{x}^2} \begin{pmatrix} (s_x^2+\overline{x}^2)\overline{y}-\overline{x}\cdot (s_{xy}+\overline{x}\cdot\overline{y})\\ -\overline{x}\overline{y}+\overline{xy} \end{pmatrix}\)

\(\therefore \begin{pmatrix} a \\ b \end{pmatrix} = \dfrac{1}{s_{x}^2} \begin{pmatrix} s_x^2\overline{y}-s_{xy}\cdot \overline{x}\\ s_{xy} \end{pmatrix} = \begin{pmatrix} \overline{y}-\dfrac{s_{xy}}{s_{x}^{2}}\overline{x} \\ \dfrac{s_{xy}}{s_{x}^2} \end{pmatrix} = \begin{pmatrix} \overline{y}-r_{xy}\dfrac{s_{y}}{s_{x}}\overline{x} \\ r_{xy}\dfrac{s_{y}}{s_{x}} \end{pmatrix}\)

\(s_x^2,~s_y^2\) は x,y の分散、\(s_{xy}\) は x,y の共分散、\(r_{xy}\) は x,y の相関係数です。

行列を用いると、更に、簡潔な表現が可能です。

　\(Y=\begin{pmatrix} y_{1} \\ y_{2} \end{pmatrix}\:\:,\:\:X=\begin{pmatrix} 1 & x_{1} \\ 1 & x_{2} \end{pmatrix} \:\:,\:\: \beta=\begin{pmatrix} a\\b \end{pmatrix}\)　とおくとき、

最小二乗法の目的は、\(Y= X\beta\) に最もよく当てはまるよう係数 \(\beta\)を定めることです。

\(X\beta=Y\)の両辺に Xの転置行列 X’ を乗ずると、\(X'X\beta=X'Y\).
（\(X'X\) は正方行列で、逆行列を持ちます）

　　　\(X'X=\begin{pmatrix} 1 & 1 \\ x_{1} & x_{2} \end{pmatrix} \begin{pmatrix} 1 & x_{1} \\ 1 & x_{2} \end{pmatrix}= \begin{pmatrix} 2 & x_1+x_2 \\ x_1+x_2 & x_1^2+x_2^2 \end{pmatrix}\)

　　　\(X'Y=\begin{pmatrix} 1 & 1 \\ x_{1} & x_{2} \end{pmatrix} \begin{pmatrix} y_{1} \\ y_{2} \end{pmatrix}= \begin{pmatrix} y_1+y_2 \\ x_1y_1+x_2y_2 \end{pmatrix}\)

正規方程式（タイプ１）は　\(X'X\beta=X'Y\) と表され、

　　　\(\beta=(X'X)^{-1}X'Y\) 　が回帰係数の推定値です。

3.3 回帰モデル

前節では、２点に当てはまる回帰直線について正規方程式を構成しました。
実は、データが n個の点になっても、同じ形の正規方程式が導かれます。
これが数学の面白いところです。

\(X=\begin{pmatrix} 1 & x_{1}\\ \vdots & \vdots \\ 1 & x_{n} \end{pmatrix}\) 、\(x= \begin{pmatrix} x_{1}\\ \vdots\\ x_{n} \end{pmatrix}\)、\(Y= \begin{pmatrix} y_{1}\\ \vdots\\ y_{n} \end{pmatrix}\)、\(1= \begin{pmatrix} 1\\ \vdots\\ 1 \end{pmatrix}\) とします。

母回帰係数：\(\beta=\begin{pmatrix} \beta_1\\ \beta_2 \end{pmatrix}\)　　誤差変動：\(u=\begin{pmatrix} \varepsilon_{1} \\ \vdots\\ \varepsilon_{n} \end{pmatrix}\) とすると、

回帰モデルは　\(Y=X\beta+u\)　で表されます。

　　　但し、\(\varepsilon_1,\cdots,\varepsilon_n \sim N(0,~\sigma^2) ~:~ i.i.d.\) （互いに独立）
　

＜正規方程式を構成＞.
正規方程式は回帰モデルから容易に構成できます。
ざっくりとした理解として、

　　「\(Y\fallingdotseq X\beta\) だから、回帰係数は \(\beta\fallingdotseq X^{-1}Y\)でよいのではないか」

と思われるかもしれません。
しかし、Xは正方行列ではないので逆行列はありません。
そこで、\(X\beta \fallingdotseq Y\) の両辺に転置行列 \(X'\)を乗じ、

　　　\(X'X\beta\fallingdotseq X'Y\) （正規方程式）とすることで.

\(\beta\) の係数 \(X'X\) を正方行列にします。
（逆行列については線形代数で学びます）

　　　\(X'X=\begin{pmatrix} 1' \\ x' \end{pmatrix} \begin{pmatrix} 1 & x \end{pmatrix}= \begin{pmatrix} 1'\cdot 1 & 1'\cdot x \\ x'\cdot1 & x'\cdot x \end{pmatrix} = \begin{pmatrix} n & n\overline x \\ n\overline x & n\overline {x^2} \end{pmatrix}\).

　　　\(X'Y=\begin{pmatrix} 1' \\ x'\end{pmatrix}Y =\begin{pmatrix} 1'\cdot Y \\ x'\cdot Y \end{pmatrix} = \begin{pmatrix} n\overline{y} \\ n\overline{xy} \end{pmatrix}\).

\(X'X\beta= X'Y\) は最小二乗法から導いた正規方程式と一致します。

注意.
詳しくは、
\(X\beta+u=Y\) より　\(X'X\beta+X'u=X'Y\).
誤差の性質から　\(X'u=0\) となります。
興味のある方は回帰分析の流れを参照してください。

4 回帰係数の変動

\(\beta=(X'X)^{-1}X'Y\) 　は回帰係数の推定値です。

Yは誤差変動のある確率変数なので、\((X'X)^{-1}X'Y\) も確率変数です。

Y の変動による回帰係数の推定値 \((X'X)^{-1}X'Y\) の変動を、回帰モデルを用いて考えます。

誤差がなく \(Y=X\beta\) であれば、
　　　\((X'X)^{-1}X'Y=(X'X)^{-1}X'X\beta\)　は母回帰係数 \(\beta\) となり、

誤差があり \(Y=X\beta+u\) であれば、
　　　\((X'X)^{-1}X'Y=(X'X)^{-1}X'(X\beta+u)\)　は標本回帰係数 \(\widehat{\beta}\) になります。

標本回帰係数の変動は、

　　　\(\begin{split} \widehat{\beta }-\beta &=\left( X'X\right) ^{-1}X'u \end{split}\)

回帰モデルの誤差\(u\)が回帰係数の変動に影響していることが分かります。

\(E\left(\widehat{\beta }-\beta\right) =\left( X'X\right) ^{-1}X'E(u)=0\)　であるので、
　　　\(E\left(\widehat{\beta}\right)-\beta=0\)　　\(\therefore E\left(\widehat{\beta }\right)=\beta\)

標本回帰係数 \(\widehat{\beta}\) の期待値は 母回帰係数となります。

4.1 分散共分散行列

2変量確率変数 \(X=\begin{pmatrix}X_1\\ X_2\end{pmatrix}\)　の期待値を \(m=\begin{pmatrix}m_1\\ m_2\end{pmatrix}\)　としたとき、

　　\(Cov\left(X\right)=E\left((X-m)(X-m)'\right)\) を分散共分散行列といます。

成分で計算すると、

　\(\begin{split}(X-m)(X-m)'&=\begin{pmatrix} X_1-m_1 \\ X_2-m_2 \end{pmatrix} \begin{pmatrix} X_1-m_1 & X_2-m_2 \end{pmatrix}\\&= \begin{pmatrix} (X_1-m_1)^2 & (X_1-m_1)(X_2-m_2) \\ (X_1-m_1)(X_2-m_2) & (X_2-m_2)^2 \end{pmatrix} \end{split}\)

\(\therefore E\left((X-m)(X-m)'\right)=\begin{pmatrix} V(X_1) & Cov(X_1,X_2) \\ Cov(X_1,X_2) & V(X_2) \end{pmatrix}\)

対角成分は分散、対角でない成分は共分散になります。
２変量の確率変数では分散のほかに変量相互の関連を示す共分散を求める必要があり、 分散共分散行列は、これらを同時に求める有効な手段です。

標本回帰係数 \(\widehat\beta =\begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix}\) は2変量の確率変数なので、行列を用いて分散共分散行列を計算できます。

\(\begin{split} \begin{aligned}Cov\left( \widehat{\beta }\right)=E\left( \left( \widehat{\beta }-\beta \right) \ \left( \widehat{\beta }-\beta \right) '\right) &=\left( X'X\right) ^{-1}X'E\left( uu'\right) X\left( X'X\right) ^{-1}\\ &=\left( X'X\right) ^{-1}X'\sigma^2 IX\left( X'X\right) ^{-1}=\sigma ^{2}\left( X'X\right) ^{-1} \ \end{aligned}\\ \end{split}\)

回帰係数の変動は、回帰モデルの誤差変動 \(E(uu')=\sigma^2I\) に集約されます。
高校生の皆さんは、この美しさを感じてください。疑問の残るところは行列を学べば自然と解消します。

＜まとめ＞.

　標本回帰係数 \(\widehat\beta =\begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix}\)　は標本によって変動する2変量確率変数であり、

　　　\(Cov\left( \widehat{\beta }\right)=\sigma ^{2}\left( X'X\right) ^{-1}\)　は標本回帰係数の分散共分散行列です。

　　　\((X'X)^{-1}=\begin{pmatrix} n & \ \ n\overline{x} \\ n\overline{x} &\ \ n\overline{x^2} \end{pmatrix}^{-1} = \dfrac{1}{n} \begin{pmatrix} 1 & \ \ \overline{x} \\ \overline{x} & \ \ \overline{x^2} \end{pmatrix}^{-1} = \dfrac{1}{ns_{x}^2} \begin{pmatrix} \overline{x^2} & \ \ -\overline{x} \\ -\overline{x} & \ \ \ \ 1 \end{pmatrix}\)

\({ns_{x}^2}\)は x の偏差平方和であるので、
　　　\(\sigma ^{2}\left( X'X\right) ^{-1}= \dfrac{\sigma^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2} \begin{pmatrix} \overline {x^2} & -\overline x \\ -\overline x & 1 \end{pmatrix}\)

簡単のため、\(T_{xx}=\sum_{i=1}^{n}( x_{i}-\overline{x})^2\) とおくと、

　　　\(Cov\left( \widehat{\beta }\right)=\sigma ^{2}\left( X'X\right) ^{-1}= \dfrac{\sigma^2}{T_{xx}} \begin{pmatrix} \overline {x^2} & -\overline x \\ -\overline x & 1 \end{pmatrix}\)

とまとめることができます。

4.2 標本回帰直線の変動

［出典］　統計検定１級 統計応用 2012年 共通①［単回帰の基本］

＜回帰モデル＞
　\(y_i=\alpha+\beta x_i+\varepsilon_i\:\:(i=1,2,\cdots,n)\) 　
　　　　 \(x_i\:\:\:(i=1,2,\cdots,n)\) は定数
　　　　 \(\varepsilon_i\sim N(0,\:\sigma^2)\) は互いに独立

今、モデルのパラメーターは全て既知であるとし、モデルからの標本による標本回帰直線の変動範囲を考えます。

• \(\begin{pmatrix}\alpha\\ \beta\end{pmatrix}\)は真の回帰係数であり定数です。

• \(\begin{pmatrix}a\\ b\end{pmatrix}\)はモデルからの標本による標本回帰係数で誤差による変動を伴います。

最小二乗法では、
残差　\(e_i=y_i-(a+bx_i)\) の平方和が最小となるように \(a,b\) を定めます。（OLS推定量）.

\(y_i=a+bx_i+e_i\:\:(i=1,2,\cdots,n)\)より、残差 \(e_i\) は誤差 \(\varepsilon_i\) の実現値と考えられます。

　　　\(\displaystyle s^2=\dfrac{1}{n-2}\sum_{i=1}^n e_i^2\) 　は \(\sigma^2\) の不偏推定量かつ一致推定量

　　　\(Cov\left[\begin{pmatrix}a\\ b\end{pmatrix}\right]=\sigma ^{2}\left( X'X\right) ^{-1}= \dfrac{\sigma^2}{T_{xx}} \begin{pmatrix} \overline {x^2} & -\overline x \\ -\overline x & 1 \end{pmatrix}\).

予測値 \(a+bx_i\:\:(i=1,2,\cdots,n)\) について、

　　　\(V(a+bx_i)=V(a)+x_i^2V(b)+2x_i\cdot Cov(a,b)\)　であり、

　　　\(V(a)=\dfrac{\sigma^2}{T_{xx}}\overline{x^2} =\dfrac{\sigma^2}{T_{xx}}\left(\dfrac{T_{xx}}{n}+\overline{x}^2\right) =\dfrac{\sigma^2}{n}+\dfrac{\sigma^2}{T_{xx}}\overline{x}^2\).
　　　\(V(b)=\dfrac{\sigma^2}{T_{xx}}\)　　\(Cov(a,b)=-\dfrac{\sigma^2\overline{x}}{T_{xx}}\)　であるので、

　\(V(a+bx_i)=\dfrac{\sigma^2}{n}+\dfrac{\sigma^2}{T_{xx}}\overline{x}^2+\dfrac{\sigma^2}{T_{xx}}x_i^2 -2\dfrac{\sigma^2}{T_{xx}}x_i\overline{x}=\dfrac{\sigma^2}{n}+\dfrac{\sigma^2}{T_{xx}}(x_i-\overline{x})^2\)

　　　\(\therefore~a+bx_i\sim N\left(\alpha+\beta x_i\:,\:\dfrac{\sigma^2}{n}+\dfrac{\sigma^2}{T_{xx}}(x_i-\overline{x})^2\right)\)

\(\sigma^2\) が既知であるので、Z変換により95%の確率で、
　　　\(|(a+bx_i)-(\alpha+\beta x_i)|<1.96\sqrt{\dfrac{\sigma^2}{n}+\dfrac{\sigma^2}{T_{xx}}(x_i-\overline{x})^2}\)

標本回帰直線による予測値の上限と下限は、

　　　　　　\(a+bx_i=\alpha+\beta x_i\pm1.96\sqrt{\dfrac{\sigma^2}{n}+ \dfrac{\sigma^2}{T_{xx}}(x_i-\overline{x})^2}\)

標本回帰直線による予測値は、母回帰直線の値を中心として、標準偏差 1.96倍の範囲内にあることが分かります。
具体的な例でシミュレーションしましょう。

＜例＞
「2.2 回帰分析のモデル」を用います。

　　説明変数：　\(X=(x_i)=(0,1,2,3,4)\)
　　母回帰式：　\(y=1+0.5 x\)
　　誤差変動：　\(\varepsilon_i \sim N(0,~0.5^2)\)

モデルをもとに100回シミュレーションを行い、1回ごとに lm( )関数で標本回帰直線を描きます。
また、各回ごとの標本回帰係数 a, b を求めておくことで、標本回帰係数の分布を調べます。
母回帰直線を赤い直線で、信頼区間の上限・下限を赤丸で表示します。

plot(2,2,asp = 1,xlim=c(-2,6),ylim=c(-1,5),xlab="",ylab="")
S=c()
A=c()
B=c()
DATA=data.frame(X=c(),Y=c())
X <- c(0,1,2,3,4)
for (i in 1:100){
  e <- rnorm(5,0,0.5)
  Y <- 0.5*X+1+e
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ X, data = data )
  f <- fit$ fitted.values#予測値
  a <- fit$ coefficients[1]#回帰係数
  b <- fit$ coefficients[2]#回帰係数  
  e <- fit$ residuals#残差
  se <- sum(e^2)#残差平方和
  s <- se/3#残差分散
  S <- c(S,s)
  A <- c(A,a)
  B <- c(B,b)
  DATA <- rbind(DATA, data)
  abline(fit,col="lightblue")
}
abline(a=1,b=0.5,col="red")
#信頼区間
vy <- 0.5^2/5+0.5^2/10*(X-2)^2 #分散
yu <- 1+0.5*X+sqrt(vy)*1.96 #上限
yl <- 1+0.5*X-sqrt(vy)*1.96 #下限
x <- c(X,X)
y <- c(yu,yl)
par(new=T)
plot(x,y,asp = 1,xlim=c(-2,6),ylim=c(-1,5),col="red")
segments(X, yu, X, yl,lty ="dotted",col="red")

赤の点を縦に結んだ区間が \(a+bx\) の95%信頼区間です。
中心から離れるほど信頼区間が広くなっていることが分かります。
青の直線は、標本からlm( )関数を用いて求めた標本回帰直線です。

＜理論値と実現値の比較＞

標本回帰係数の分布について、理論値とシミュレーション結果を比較します。

標本回帰係数 \(\widehat\beta =\begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix}\) の各分散と共分散は、

• \(V\left(\widehat\beta_{1}\right)=\dfrac{\sigma^2\overline{x^2}}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)　　　\(V\left(\widehat\beta_{2}\right)=\dfrac{\sigma^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)、

• \(Cov(\widehat\beta_{1},\widehat\beta_{2})=\dfrac{-\sigma^2\overline{x}}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)

公式の右辺は、回帰モデルから計算される理論値です。
公式の左辺は、シミュレーション結果から計算されます。
理論値とシミュレーションによる実現値を比較します。

• 標本回帰係数.

mean <- c(mean(A),mean(B));mean

## [1] 0.9737466 0.5230003

　　　母回帰係数 1 と 0.5 に一致しています。

• 回帰係数の分散・共分散.

varB <- c(0.5^2/sum((X-mean(X))^2),var(B));varB

## [1] 0.02500000 0.02941424

varA <- c(0.5^2*mean(X^2)/sum((X-mean(X))^2),var(A));varA

## [1] 0.1500000 0.1706244

covAB <- c(0.5^2*mean(-X)/sum((X-mean(X))^2),cov(A,B));covAB

## [1] -0.05000000 -0.05865562

　　　分散・共分散の理論値は実現値とほぼ一致しています。

標本回帰係数の分散共分散行列を数理的に求めるのは行列と線形代数を学んでからですが、シミュレーションで確かめておくことは重要です。

次に、理論分布に基づく標本回帰係数を、２次元正規乱数 rmvnorm( ) を用いて出力し、標本回帰直線をピンクで描きます。

sigma <- matrix(c(0.15,-0.05,-0.05,0.025), ncol=2)
plot(2,2,asp = 1,xlim=c(-2,6),ylim=c(-1,5),xlab="",ylab="")
for (i in 1:100){
  X <- rmvnorm(n=1, mean=c(1,0.5), sigma=sigma)  
  abline(a=X[1],b=X[2],col="pink")
}
abline(a=1,b=0.5,col="red") #母回帰直線
par(new=T)
plot(x,y,asp = 1,xlim=c(-2,6),ylim=c(-1,5),col="blue")
segments(x, yu, x, yl,lty ="dotted",col="blue")

ピンクの直線群は青の直線群とほぼ一致しており、分散共分散行列による理論値が妥当なものであることが分かります。

4.3 母回帰式の予測

実際の問題では、母回帰式は分かっていません。標本回帰式をもとに母回帰式の値（真の値）を区間推定する必要があります。
即ち、標本回帰式の値から真値の信頼区間を求めなければなりません。

前出の式を、
　　　\(\alpha+\beta x_i=a+bx_i\pm1.96\sqrt{\dfrac{\sigma^2}{n}+ \dfrac{\sigma^2}{T_{xx}}(x_i-\overline{x})^2}\)

と考えれば、標本回帰式の値 \(a+bx_i\) から真値 \(\alpha+\beta x_i\) の信頼区間を求めることができます。

4.3.1 誤差分散が既知のとき

まず、\(\sigma^2\) が既知の場合から考えます。

標本回帰式を青で示し、標本回帰式を中心とした信頼区間を青の点線で示します。
赤の母回帰直線は 95% の確率で、青の点線の信頼区間を通るはずです。

plot(2,2,asp = 1,xlim=c(-2,6),ylim=c(-1,5),col="red",xlab="",ylab="",pch=16)
abline(a=1,b=0.5,col="red") #母回帰直線
X <- c(0,1,2,3,4)
e <- rnorm(5,0,0.5)
Y <- 0.5*X+1+e
data <- data.frame(X=X,Y=Y)
fit = lm( formula = Y ~ X, data = data )
a <- fit$ coefficients[1]#回帰係数
b <- fit$ coefficients[2]#回帰係数  
e <- fit$ residuals#残差
se <- sum(e^2)#残差平方和
s <- se/3#残差分散
abline(fit,col="blue")
par(new=T)
plot(2,mean(Y),asp = 1,xlim=c(-2,6),ylim=c(-1,5),col="blue",xlab="",ylab="",pch=16)
vy <- s/5+s/10*(X-2)^2 #分散
yu <- a+b*X+sqrt(vy)*qnorm(0.975) #上限
yl <- a+b*X-sqrt(vy)*qnorm(0.975)  #下限
x <- c(X,X)
y <- c(yu,yl)
par(new=T)
plot(x,y,asp = 1,xlim=c(-2,6),ylim=c(-1,5),col="blue")
segments(X, yu, X, yl,lty ="dotted",col="blue")

母回帰直線は (2,2) を通り、標本回帰直線は (2,mean(Y)) を通ります。（図では赤点、青点で表示）
これらは一致するとは限りません。

課題
　上記のプログラムをRで何度も実行すると、赤の母回帰直線は 95% の確率で青の信頼区間内を通過するはずです。
このことをRstudioやGoogleColab で確かめましょう。

4.4 誤差分散が未知のとき

誤差分散 \(\sigma^2\) が未知の場合、標本の残差分散 \(s^2\) で代用しなければなりません。
信頼区間は 標準正規分布の qnorm(0.975)=1.96 の代替として、自由度３のt分布 t(0.975, 3)=3.18 を用いなければなりません。

前出の式は、
　　　\(\alpha+\beta x_i=a+bx_i\pm dt(0.975,3)\sqrt{\dfrac{s^2}{n}+ \dfrac{s^2}{T_{xx}}(x_i-\overline{x})^2}\)

と修正し、標本回帰式の値 \(a+bx_i\) から真値 \(\alpha+\beta x_i\) の信頼区間を求めることができます。

シミュレーションで確かめます。
標本回帰式を青で示し、標本回帰式を中心とした信頼区間を青の点線で示します。

plot(2,2,asp = 1,xlim=c(-2,6),ylim=c(-1,5),col="red",xlab="",ylab="",pch=16)
abline(a=1,b=0.5,col="red") #母回帰直線
X <- c(0,1,2,3,4)
e <- rnorm(5,0,0.5)
Y <- 0.5*X+1+e
data <- data.frame(X=X,Y=Y)
fit = lm( formula = Y ~ X, data = data )
a <- fit$ coefficients[1]#回帰係数
b <- fit$ coefficients[2]#回帰係数  
e <- fit$ residuals#残差
se <- sum(e^2)#残差平方和
s <- se/3#残差分散
abline(fit,col="blue")
par(new=T)
plot(2,mean(Y),asp = 1,xlim=c(-2,6),ylim=c(-1,5),col="blue",xlab="",ylab="",pch=16)
vy <- s/5+s/10*(X-2)^2 #分散
yu <- a+b*X+sqrt(vy)*qt(0.975, 3)  #上限
yl <- a+b*X-sqrt(vy)*qt(0.975, 3) #下限
x <- c(X,X)
y <- c(yu,yl)
par(new=T)
plot(x,y,asp = 1,xlim=c(-2,6),ylim=c(-1,5),col="blue")
segments(X, yu, X, yl,lty ="dotted",col="blue")

信頼区間の幅が大きくなっていることが読み取れます。

• 母回帰直線は、データの中心 (2,2) を通る直線です。

• 標本回帰直線は、標本の中心 (2,mean(Y)) を通る直線です。

5 回帰統計量

回帰分析における基本的な統計量について整理しておきます。

計算は lm関数で行われ summary( )で出力されます。

しかし、summaryで出力される数値の意味は理解できているでしょうか。
数値の意味を解釈するには、基本的な数理の理解が必要です。

また、lm関数から統計量を抽出する方法も知っておく必要があります。

• \(fit = lm( formula = Y \sim X, data = data )\)

fit とおいた lm( )の出力には様々な統計量が含まれ、summary(fit) によって要約が出力されます。
統計量が知りたい場合は「fit $ 統計量」によって抽出します。

5.1 残差

残差は「観測値と予測値の差」です。
回帰分析の summary では “Residuals” として表示されています。

最小二乗法では、残差平方和が最小となる回帰係数を、微分を用いて計算します。
そのプロセスを見ておくことは、回帰統計量への理解を深めるのに役立つでしょう。
行列で簡潔に表現します。高校生の方は線形代数を学んでから再学習してください。

\(X=\begin{pmatrix} 1 & x_{1}\\ \vdots & \vdots \\ 1 & x_{n} \end{pmatrix}\) 、\(x= \begin{pmatrix} x_{1}\\ \vdots\\ x_{n} \end{pmatrix}\)、\(Y= \begin{pmatrix} y_{1}\\ \vdots\\ y_{n} \end{pmatrix}\)、\(1= \begin{pmatrix} 1\\ \vdots\\ 1 \end{pmatrix}\)

回帰係数：\(\beta=\begin{pmatrix} \beta_1\\ \beta_2 \end{pmatrix}\)　 残差：\(e=\begin{pmatrix} e_{1} \\ \vdots\\ e_{n} \end{pmatrix}\)　とし、

残差平方和を行列を用いて表現すると、

　　　\(h(\beta)=(Y-X\beta)'(Y-X\beta)\)　

最小二乗法は、残差平方和が最小となるよう回帰係数 \(\widehat{\beta}\) を定めるので、

　　　\(h'(\widehat{\beta})=2(-X')(Y-X\widehat{\beta})=0\) 　

残差は　\(e=Y-X\widehat{\beta}\)　であるので、

　　　\(X'e=0\)　　（よって \(X,~e\) は無相関）

予測値は \(\widehat{Y}=X\widehat{\beta}\) であり

　　　\(e'\widehat{Y}=e'X\widehat{\beta}=0~\widehat{\beta}=0\)　　（よって \(\widehat{Y},~e\) は無相関） 　　　

無相関となる理由について、もう少し詳しく説明します。
相関係数の定義式.
　　　\(Cor(x,y)=\dfrac{Cov(x,y)}{s_x~s_y}\) 　の分子は共分散です。

相関係数が 0 であることと、共分散が 0 であることは同値です。

\(Cov(x,y)=E(xy)-E(x)E(y)\)　であるので、

　　　\(Cov(X,e)=E(Xe)-E(X)E(e)=0-0=0\)

　　　\(Cov(\widehat{Y},e)=E(\widehat{Y}e)-E(\widehat{Y})E(e)=0-0=0\)

よって、\((X,~e)\) は無相関、\((\widehat{Y},~e)\) も無相関となります。

共分散について、もう少し詳しく見ておきます。

5.2 共分散

共分散は線形性をもつ良い統計量です。

• \(Cov(X+Y,Z)=Cov(X,Z)+Cov(Y,Z)\)
  
• \(Cov(aX,Y)=aCov(X,Y)\)

偏差積和は \((n-1)Cov(X,Y)\) で計算されます。
偏差積和も線形性を持ちます。

残差と予測値の共分散は０なので、残差と予測値の偏差積和も 0 になります。

残差と観測値の共分散は、

　　　\(Cov(e,Y)=Cov(e,\widehat{Y}+e)=Cov(e,\widehat{Y})+Cov(e,e)=V(e)\)

予測値と観測値の共分散は、
　　　\(Cov(\widehat{Y},Y)=Cov(\widehat{Y},\widehat{Y}+e)=Cov(\widehat{Y},\widehat{Y})+Cov(\widehat{Y},e)=V(\widehat{Y})\)

観測値の分散は、
　　　\(Cov(Y,Y)=Cov(\widehat{Y}+e,\widehat{Y}+e)=Cov(\widehat{Y},\widehat{Y})+Cov(e,e)\)

共分散は偏差積和をデータ数で割ったものなので、共分散にデータ数をかければ偏差積和になります。

従って、

• 残差と観測値の偏差積和は、残差平方和に等しく、
  
• 予測値と観測値の偏差積和は、予測値の偏差平方和（回帰変動）に等しい。
  
• 観測値の偏差平方和は、予測値の偏差平方和と残差平方和の合計に等しい。
  （\(S_{T}=S_R+S_e\) で表され、平方和の分解公式といいます）

これらの事項は、シミュレーションで確かめることができます。

5.3 標本回帰係数

n個の指定された説明変数 \(x_i\) に対し \(y_i\) が観測されたとき、最小二乗法によって標本回帰係数 \(\widehat\beta\) が計算されます。

行列を用いると　\(\widehat\beta=(X'X)^{-1}X'Y\)

標本回帰係数の従う分布は、
　　　\(\widehat{\beta }\sim N\left(\beta\ ,\ \sigma ^{2}\left( X'X\right) ^{-1}\right)\) 　（\(\beta\) は母回帰係数）

　　　\(\sigma^2\left( X'X\right) ^{-1}= \dfrac{\sigma^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2} \begin{pmatrix} \overline {x^2} & -\overline x \\ -\overline x & 1 \end{pmatrix}\)

\(\widehat\beta =\begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix}\) の分散は次のようになります。

• \(V\left(\widehat\beta_{1}\right)=\dfrac{\sigma^2\overline{x^2}}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)　　　\(V\left(\widehat\beta_{2}\right)=\dfrac{\sigma^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)

\(\sigma^2\) を \(s^2=\dfrac{1}{n-2} e'e\) （eは残差）で代用すると、

• \(V\left(\widehat\beta_{1}\right)=\dfrac{s^2\overline{x^2}}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)　　　\(V\left(\widehat\beta_{2}\right)=\dfrac{s^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)
  　

\(\sqrt{V\left(\widehat\beta_{1}\right)},~\sqrt{V\left(\widehat\beta_{2}\right)}\) は標本回帰係数の標準誤差です。

6 統計量の抽出

　　\(fit = lm( formula = Y \sim X, data = data )\)

\(fit\) から回帰係数、残差などを抽出し必要な回帰統計量を構成します。

• 標本回帰係数は　“fit $ coefficients”
  fit $ coefficients[1] は定数項を、fit $ coefficients[2] は一次の係数を出力します。
  

• 残差は　e = ” fit$ residuals “
  残差の自由度は (n-2) なので、残差分散は \(s^2=var(e)(n-1)/(n-2)\)
  残差標準誤差 \(s\) は Residual standard errorといいます。

• 残差分散は
  \(s^2=sum(e^2)/(n-2)\) でも計算されます。
  残差分散 \(s^2\) は誤差分散 \(\sigma^2\) の推定値です。

• 回帰の予測値は　f=“fit $ fitted.values” で抽出されます。

実際に出力してみましょう。
引き続き「2.2 回帰分析のモデル」を用います。

　　説明変数：　\(X=(x_i)=(0,1,2,3,4)\)
　　母回帰式：　\(y=1+0.5 x\)
　　誤差変動：　\(\varepsilon_i \sim N(0,~0.5^2)\)

\(\overline{x}=2,~~~~\overline{x^2}=6,~~~~\sum_{i=1}^{n}( x_{i}-\overline{x})^2=10\) ですが、この値は \(var(X)*4\) で計算できます。

X <- c(0,1,2,3,4)
e <- rnorm(5,0,0.5)
Y <- 0.5*X+1+e
data <- data.frame(X=X,Y=Y)
fit = lm( formula = Y ~ X, data = data )
f <- fit$ fitted.values#予測値
a <- fit$ coefficients[1]#回帰係数
b <- fit$ coefficients[2]#回帰係数  
e <- fit$ residuals#残差

lm( )関数から求めた統計量とsummaryの内容が符合することを確かめます。

Residuals（残差）

e <- fit$ residuals;e #残差

##          1          2          3          4          5 
##  0.5487052 -0.4771162 -0.5314481  0.2994238  0.1604352

Coefficients（標本回帰係数）

a <- fit$ coefficients[1]#回帰係数（切片）
b <- fit$ coefficients[2]#回帰係数
c(a,b)

## (Intercept)           X 
##   0.8280505   0.4798930

Residual standard error（残差標準誤差）

se <- sum(e^2)#残差平方和
s2 <- se/3#残差分散
sqrt(s2)#残差標準誤差

## [1] 0.5557423

残差標準誤差は summary の下部　“Residual standard error” に表示される数値で、回帰モデルの誤差の推定値です。
残差は平均が0なので、var( )関数を用いた計算も可能です。自由度が \(n-2=3\) であることに注意しましょう。

sqrt(var(e)*4/3)

## [1] 0.5557423

Std. Error（標準誤差）
“Coefficients”に表示される標本回帰係数の標準誤差です。

　　　\(V\left(\widehat\beta_{1}\right)=\dfrac{s^2\overline{x^2}}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)　　　\(V\left(\widehat\beta_{2}\right)=\dfrac{s^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)

に基づいて計算します。

vb <- s2/(var(X)*4)#標本回帰係数の分散
va <- mean(X^2)*s2/(var(X)*4)#標本回帰係数の分散
c(sqrt(va),sqrt(vb))#標本回帰係数の標準誤差

## [1] 0.4304761 0.1757411

6.1 summary（概要）

統計量が summary に表示されていることを確認しましょう。

summary(fit)

## 
## Call:
## lm(formula = Y ~ X, data = data)
## 
## Residuals:
##       1       2       3       4       5 
##  0.5487 -0.4771 -0.5314  0.2994  0.1604 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.8281     0.4305   1.924   0.1501  
## X             0.4799     0.1757   2.731   0.0719 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5557 on 3 degrees of freedom
## Multiple R-squared:  0.7131, Adjusted R-squared:  0.6175 
## F-statistic: 7.457 on 1 and 3 DF,  p-value: 0.0719

7 理論値との整合性

データサイエンスの学習は公式を使うだけのものになりがちです。
しかし、自ら発見すれば、公式は違った輝きを帯びてきます。
予測値と観測値の偏差積和が回帰変動に等しいことを発見したときは感動しました。
そして、重相関係数と決定係数の関係がすんなり理解できました。

抽出した回帰統計量と理論値との整合性を確かめることは、探究的活動として意味のあることではないかと思います。

7.1 残差の無相関性

相関係数は単位に影響されません。
分散は nで割るか、(n-1)で割るかという問題がありますが、相関係数はどちらでも値は変わりません。
相関係数は優れた統計量と言えるでしょう。
ここでは「残差の無相関性」に注目すると、

• 残差は、説明変数と無相関
  
• 残差は、予測値と無相関

であることが確かめられます。

cor(e,X)#残差と説明変数の相関

## [1] 5.471009e-17

cor(e,f)#残差と予測値の相関

## [1] 2.090087e-16

残差と説明変数の相関係数は 0 に近く、無相関であることが分かります。
残差と予測値の相関係数も 0 に近く、無相関です。

Note.
残差と観測値との相関係数は 0 にはなりません。

cor(e,Y)#残差と観測値の相関係数

## [1] 0.5356301

cov(e,Y)#残差と観測値の共分散

## [1] 0.2316371

共分散も 0 になりません。
次のセクションで、残差と観測値の共分散・偏差積和について考えます。

7.2 観測値との共分散

観測値との共分散については、

• 観測値と残差の共分散は、残差の分散に等しく、
  
• 観測値と予測値の共分散は、予測値の分散に等しく、
  
• 観測値の分散は、上記２つの合計になります。

これは、「5.2 共分散」でも説明しましたが、

• \(Cov(Y,e)=V(e)\)
  
• \(Cov(Y,\widehat{Y})=V(\widehat{Y})\)
  
• \(Cov(Y,\widehat{Y}+e)=V(\widehat{Y})+V(e)\)

第３式の左辺は、
　　　\(Cov(Y,\widehat{Y}+e)=Cov(Y,Y)=V(Y)\) なので、

• \(V(Y)=V(\widehat{Y})+V(e)\)

と簡潔に説明されます。

出力された値で確かめると、

＜Cov(観測値, 残差) と V(残差)＞

c(cov(e,Y),var(e))

## [1] 0.2316371 0.2316371

＜Cov(観測値,予測値) と V(予測値)＞.

c(cov(f,Y),var(f))

## [1] 0.5757433 0.5757433

＜V(観測値) と V(予測値)＋V(残差)＞.

c(var(Y),var(f)+var(e))

## [1] 0.8073804 0.8073804

このように、公式を出力データで確かめることができます。

Rには偏差積和の関数はありませんが、共分散を(n-1)倍すれば偏差積和になるので、

• 観測値と残差の偏差積和は、残差の偏差平方和に等しく、
  
• 観測値と予測値の偏差積和は、予測値の偏差平方和に等しく、
  
• 観測値の偏差平方和は、上記２つの合計.

が確かめられたことになります。
１つ目の値を残差変動 \(S_e\)、２つ目の値を回帰変動 \(S_R\)、３つ目の値を全変動 \(S_T\) と呼ぶことにすれば、

　　　\(S_{T}=S_R+S_e\)　　（平方和の分解公式）.

これらの事項は、変動という概念を用いて次のようにまとめられます。

7.3 平方和の分解公式

全変動 (Total Sum of Squares, TSS):
　　\(\displaystyle S_{T}=\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}\)

回帰変動 (Explained Sum of Squares, ESS):
　　\(\displaystyle S_R=\sum ^{n}_{i=1}\left( (\alpha+ \beta x_i)-\overline{y}\right)^2\)

残差変動 (Residual Sum of Squares, RSS):
　　\(\displaystyle S_e=\sum ^{n}_{i=1}\left( y_{i}-(\alpha+ \beta x_i)\right)^2\)

分解公式
　　\(S_{T}=S_R+S_e\)

（導出はひょっこり統計「回帰分析の流れ」 参照）

導出は行列を用いると簡明であり多変量への拡張も容易ですが、ここではシミュレーションの出力から確かめておくことにしましょう。

st <- var(Y)*4
sr <- var(f)*4#回帰変動
se <- var(e)*4#残差平方和
c(st,sr+se)

## [1] 3.229521 3.229521

\(S_{T}\) と \(S_R+S_e\) は一致していることが分かります。

7.4 重相関係数

観測値と推測値の相関係数を重相関係数（Multiple R）といいます。

\(Cov(Y,\widehat{Y})=V(\widehat{Y})\) であることから、

　　　\(R=Cor(Y,\widehat{Y})=\dfrac{Cov(Y,\widehat{Y})}{\sqrt{V(Y)}\sqrt{V(\widehat{Y})}}=\dfrac{V(\widehat{Y})}{\sqrt{V(Y)}\sqrt{V(\widehat{Y})}}=\sqrt{\dfrac{V(\widehat{Y})}{V(Y)}}\)
　　　\(R^2=\dfrac{V(\widehat{Y})}{V(Y)}=\dfrac{S_R}{S_{T}}\).

これは全変動に対する回帰変動の寄与率である決定係数と一致しています。

cor(f,Y)^2

## [1] 0.7131004

var(f)/var(Y)

## [1] 0.7131004

Note.
　重相関係数の二乗 \(R^2\) “Multiple R-squared” は、決定係数と一致します。

8 有意性と適合度

回帰が有意とは、説明変数が有効に働き、説明変数に従って目的変数が有意に変化することです。単回帰では、回帰係数のt検定で判断します。

回帰の適合度が高いとは、説明変数による目的変数の変動が大きいことです。誤差変動の割合が小さいとも言えます。これは、回帰変動の全変動に対する割合（決定係数）によって評価します。

予測値と実測値との相関係数である重相関係数も適合度の指標として用いられます。
単回帰では、重相関係数は相関係数の絶対値に等しく、二乗は決定係数に一致するので決定係数と同等の指標と考えられます。

8.1 回帰係数のt値

回帰分析のモデルでは説明変数が固定されており、変動するのは y のみです。
単回帰分析では、「標本回帰係数＝標本相関係数」であり、 回帰の有意性は相関の有意性と同等になります。
散布図に傾きがあるというイメージです。
回帰が有意であることを示すには、
　　　「帰無仮説：回帰係数＝０」
を検定すればよいことになります。検定には t検定を用います。

t value
「帰無仮説：回帰係数＝０」のもとでの、標本回帰係数の t値です。
標本回帰係数を \(a,b\) 、標準誤差を\(va, vb\) とすると

　　　\(\dfrac{a-0}{\sqrt{va}}\sim t(n-2)\)　　　\(\dfrac{b-0}{\sqrt{vb}}\sim t(n-2)\).

これが各回帰係数の t値になります。

c(a/sqrt(va),b/sqrt(vb))

## (Intercept)           X 
##    1.923569    2.730681

次に t 値に対応するP値を求めます。

Pr(>|t|)

• pt(t値, 自由度) は t分布の分布関数の値であるので、
  P値は 1-pt(t値, 自由度)で計算されます。

c(2*(1-pt(a/sqrt(va),3)),2*(1-pt(b/sqrt(vb),3)))

## (Intercept)           X 
##  0.15011057  0.07189947

t値に対応するP値が有意水準より小さければ回帰係数は０ではないと判断します。

標本回帰係数 b の t値は、２乗すると F 値に等しくなります。
数理的説明は「11. t値とF値の関係」で行います。ここでは、シミュレーションの値で確かめておきましょう。

Ft <- var(f)*4/s2 #F値 
t <- b/sqrt(vb) #t値
c(Ft,t^2)

##               X 
## 7.45662 7.45662

8.2 決定係数

回帰分析の目的は、説明変数から目的変数を予測することです。
x に応じて y が変化し、y の変化が誤差より有意に大きいとき、回帰は適合しているといいます。
変動を表す統計量は偏差平方和であり、平方和の分解公式が成立します。

全変動に対する回帰変動の寄与率 \(\dfrac{S_R}{S_{T}}\) を決定係数といい、
これが大きければ回帰は適合していると判断します。

sr/st

## [1] 0.7131004

Note.
　単回帰の場合、決定係数は相関係数の２乗に等しくなります。実際、

(cor(X,Y))^2

## [1] 0.7131004

数理的説明は以下のようになります。

\(\displaystyle \beta=\frac{\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) \left( y_{i}-\overline{y}\right) }{ \sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) ^{2}}=\dfrac{s_{xy}}{s_x^2}=\dfrac{s_{xy}}{s_xs_y}\cdot\dfrac{s_y}{s_x}=r_{xy}\cdot\dfrac{s_y}{s_x}\)　であるので、

\(\begin{split} 　　\dfrac{S_R}{S_T}&=\frac {\sum ^{n}_{i=1}\left( (\alpha+ \beta x_i)-\overline{y}\right)^2 } {\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}} =\frac {\sum ^{n}_{i=1}\left( (\alpha+ \beta x_i)-(\alpha+\beta\overline{x})\right)^2 } {\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}}\\ &=\beta^2\frac {\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) ^{2}} {\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}} =r_{xy}^2\cdot\dfrac{s_y^2}{s_x^2}\cdot\dfrac{s_x^2}{s_y^2} =r_{xy}^2 \end{split}\).

さて、6.1（概要） summary をみてみましょう。

• 決定係数の英訳 “coefficient of determination” の記載はありません。
• 寄与率の英訳 “contribution rate” も記載はありません。
• 統計量は summary の下から２行目に”Multiple R-squared” として記載されています。

“Multiple R-squared” とはどういう意味なのでしょう？

• “Multiple R” は重相関係数といい、予測値と観測値の相関係数です。
  
• “Multiple R-squared” は「重相関係数の２乗」という意味になります。

「7.4 重相関係数」で示したように、重相関係数の２乗は決定係数と一致するので、決定係数でなく”Multiple R-squared” としてsummaryに記載されています。

重相関係数について、詳しく見ていきます。

8.3 重相関係数

予測値と観測値の相関係数を重相関係数といいます。

　　　R=cor（予測値, 観測値）

重相関係数は（予測値, 観測値） を散布図で表示したとき、どの程度直線的であるかを示す指標です。
重相関係数が１に近ければ、予測値は観測値に近いと考えられます。

この重相関係数は “Multiple R” と呼ばれ、その平方（squared）である \(R^2\)が決定係数と等しいため、Multiple R-squared と表示されているのです。

予測値を \(\widehat{y}=a+bx\) 、観測値を \(y\) とおくと、

　　\(Cor(y,~\widehat{y})=\dfrac{Cov(y,~\widehat{y})}{\sigma_y \cdot \sigma_{\widehat{y}}} = \dfrac{Cov(y,~bx)}{\sigma_y \cdot |b|\sigma_{x}} =\dfrac{b}{|b|}Cor(x,~y)\)
\(\dfrac{b}{|b|}=\pm 1\) であるので、 重相関係数と相関係数は符合を除いて一致します。

相関係数は「どの程度直線的かの指標」ですから、相関係数が１だからといって x と y の値が一致するわけではありません。
\(y=0.5x+1\) でも相関係数は１になります。
しかし、予測値は y と値が近く、\((y,~\widehat{y})\) は直線 y=x の近くに分布しているので重相関係数の符号は正で値は１に近く、予測値と観測値が近いことの指標になります。

シミュレーションで確かめましょう。

\(y=-0.5x+1\) の場合、

X <- c(0,1,2,3,4)
e <- rnorm(5,0,0.5)
Y <- -0.5*X+1+e
data <- data.frame(X=X,Y=Y)
fit = lm( formula = Y ~ X, data = data )
f <- fit$ fitted.values#予測値
a <- fit$ coefficients[1]#回帰係数
b <- fit$ coefficients[2]#回帰係数  
e <- fit$ residuals#残差
cor(Y,f)#重相関係数

## [1] 0.9667597

cor(X,Y)#相関係数

## [1] -0.9667597

x,と y の相関が負でも、重相関係数は正であることがわかります。
誤差変動を 0.5 から 0.1 に小さくすると、

X <- c(0,1,2,3,4)
e <- rnorm(5,0,0.1)
Y <- -0.5*X+1+e
data <- data.frame(X=X,Y=Y)
fit = lm( formula = Y ~ X, data = data )
f <- fit$ fitted.values#予測値
a <- fit$ coefficients[1]#回帰係数
b <- fit$ coefficients[2]#回帰係数  
e <- fit$ residuals#残差
cor(Y,f)#重相関係数

## [1] 0.9898794

cor(X,Y)#相関係数

## [1] -0.9898794

重相関係数は１に近くなります。

Note

• 重相関係数は単回帰係数ではありません。
  
• 単回帰では、重相関係数と相関係数は符合を除いて一致します。
  （重相関係数の二乗）＝（相関係数の二乗）
• 単回帰では、相関係数の二乗は決定係数に等しいので、
  （重相関係数の二乗）＝（決定係数）
  決定係数は “Multiple R-squared” としてsummaryに記載されます。

重相関係数の２乗が決定係数に等しいことは、出力データでも確かめられます。

R <- (cor(Y,f))^2;R#重相関係数の二乗

## [1] 0.9798613

R <- (var(f)*4)/(var(Y)*4);R#決定係数

## [1] 0.9798613

9 F統計量

回帰の有意性を示すもう一つの統計量として、F統計量があります。
説明が長くなるので、取り出して説明します。

決定係数は寄与率ともいい、説明変数が説明する割合を表します。
決定係数にはP値のような確率的評価がありません。

一方、F統計量は回帰モデル全体の有意性を検定する統計量であり、P値による確率評価を行います。ざっくり言えば、

• F統計量= 回帰分散／残差分散

残差変動に対して回帰変動が大きいかどうかを、分散比を用いて確率評価するものです。

分散の計算では自由度に注意する必要があります。

＜単回帰分析の自由度＞
　データ数 n のとき、帰無仮説（回帰係数０）のもとでの自由度は、

　　　回帰変動：自由度１　　　残差変動：自由度 n-2.

となります。以下で、この理由について考えます。

9.1 偏差平方和の自由度.

自由度は様々な文脈で用いられます。
ここでは、偏差平方和の自由度に絞って説明します。

ある確率変数とその母分散が与えられているとき、

• 偏差平方和が母分散の何倍であるか、その倍数が自由度です。

偏差平方和が「何個分の母分散を変動としてを持っているか」に注目していることになります。
当然、自由度で割れば母分散の推定量になる訳です。

• n個の偏差平方和の各項が独立ならば、分散の加法性により自由度は n. 

しかし、偏差平方和の計算に標本平均を用いた場合、各項は独立ではありません。
この場合、偏差平方和は \((n-1)\) 個の独立な偏差平方和に書き直されるので、値は母分散の \((n-1)\) 倍になり、自由度は \((n-1)\) になります。つまり、数理としては「何個の独立な偏差平方和になるか、その個数が自由度」ということになります。

＜まとめ＞.
偏差平方和の自由度は、

• 数理的には、何個の独立な確率変数の平方和で表されるか、その個数ということができます。
  
• シミュレーションでは、偏差平方和の期待値が母分散の何倍であるか、その倍数ということができます。

偏差平方和／母分散　という統計量を考えれば、この値が自由度ということになります。
この統計量は、偏差平方和の自由度を自由度とするカイ二乗分布に従います。
この統計量は、F統計量の構成に用いられるので重要です。（ひょっこり統計「F検定」参照）

自由度は、統計学を学ぶ人が頭を悩ませるところです。
詳しい説明は、ひょっこり統計「高校生の統計学第２章」 を参照してください。

ここでは、シミュレーションによる説明をします。
数理は難しいのですが、シミュレーションで確かめることは容易です。.

＜帰無仮説の回帰モデル＞.

　　説明変数：　\(X=(x_i)=(0,1,2,3,4)\)
　　母回帰式：　\(y=1+0\cdot x\)
　　誤差変動：　\(\varepsilon_i \sim N(0,~0.5^2)\)

• 回帰変動（予測値の偏差平方和）を変数 SR に、
  
• 残差変動（残差の偏差平方和）を変数 SE に格納し、

1000回のシミュレーション平均を求めます。

SR <- c()
SE <- c()
X <- c(0,1,2,3,4)
for (i in 1:1000){
  e <- rnorm(5,0,0.5)
  Y <- 0*X+1+e
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ X, data = data )
  f <- fit$ fitted.values#予測値
  a <- fit$ coefficients[1]#回帰係数
  b <- fit$ coefficients[2]#回帰係数  
  e <- fit$ residuals#残差
  se <- sum(e^2)#残差平方和
  SR <- c(SR,var(f)*4)
  SE <- c(SE,se)
}
mean(SR)

## [1] 0.2326054

mean(SE)

## [1] 0.7634947

• 回帰変動は母分散 0.25 の１倍なので、自由度は１
  
• 残差変動は母分散 0.25 の３倍なので、自由度は３

と判断できます。

9.2 自由度の数理

単回帰の分解公式　\(S_{T}=S_R+S_e\)　において、\(S_{T}\) の自由度は(n-1) です。

\(S_R\) の自由度が 1 であることは、どうやって確かめたら良いのでしょうか？

「4.1 共分散」のまとめから、
標本回帰係数　\(\widehat\beta =\begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix}\)　は標本によって変動する2変量確率変数であり、

　　　\(V\left(\widehat\beta_{2}\right)=\dfrac{\sigma^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)　　\(\therefore \displaystyle \sum_{i=1}^{n}( x_{i}-\overline{x})^2=\dfrac{\sigma^2}{V\left(\widehat\beta_{2}\right)}\)

　　　\(\displaystyle S_{R}={\sum ^{n}_{i=1}\left( (\widehat\beta_{1}+ \widehat\beta_{2} x_i)-\overline{y}\right)^2 } =\widehat\beta_{2}^2{\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) ^{2}}=\left(\widehat\beta_{2}\right)^2 \dfrac{\sigma^2}{V\left(\widehat\beta_{2}\right)}\)

帰無仮説（回帰係数０）のもとでは、\(\dfrac{\widehat\beta_{2}}{\sqrt{V(\widehat\beta_{2})}} \sim N(0,1)\)
よって、　\(\left(\dfrac{\widehat\beta_{2}}{\sqrt{V(\widehat\beta_{2})}}\right)^2 \sim \chi^2(1)\).
自由度１のカイ二乗分布の期待値は１であり、\(S_R\) は \(\sigma^2\) の約１倍となるので、\(S_R\) の自由度は１と考えられます。

偏差平方和／母分散　という統計量を考えれば、

　　　\(\dfrac{S_{R}}{\sigma^2}=\left(\dfrac{\widehat\beta_{2}}{\sqrt{V(\widehat\beta_{2})}}\right)^2 \sim \chi^2(1)\)　　　

回帰変動／母分散 は、標本回帰係数を確率変数と見た標準化変数の平方となり、自由度１のカイ二乗分布に従います。

このように、\(\dfrac{S_{R}}{\sigma^2}\) が何個の独立な標準正規分布の平方和で表されるかを考えれば、従うカイ二乗分布の自由度が分かり、それが回帰の自由度になっています。

自由度が１であることは、\(\dfrac{S_{R}}{\sigma^2}\) が１個の標準正規分布に従う確率変数の二乗で表されることを意味します。
\(\dfrac{S_{R}}{\sigma^2}\) は 単回帰係数を確率変数と見た標準化式の二乗で表されたことになります。
　

＜残差変動の自由度＞.

「残差平方和／母分散」は、(n-2)個の独立な標準正規分布の平方和で表されることが導かれます。
これは、9.1 のシミュレーションから分かりますが、数理的導出には冪等行列や固有値など線形代数の知識が必要です。

　　　\(\dfrac{S_e}{\sigma ^{2}}=v_1^2+v_2^2+\cdots +v_{n-2}^2\)　　　\(v_i\sim N(0,1)\) は互いに独立.
　　\(\therefore \dfrac{S_e}{\sigma ^{2}}\sim \chi ^{2}\left( n-2\right)\).

これは線形代数を学んでから学習するとよいでしょう。
（ ひょっこり統計「回帰分析の流れ」参照）.

9.3 回帰変動の期待値.

\(t=\dfrac{\widehat\beta_{2}}{\sqrt{V(\widehat\beta_{2})}}\) とおくと、\(V(t)=1\).

　　　\(V\left(\widehat\beta_{2}\right)=\dfrac{\sigma^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)　は定数です。

母回帰係数を \(\beta_2\) とすると、\(E\left(\widehat\beta_{2}\right)= \beta_2\) であるので、
　　　\(E(t)^2=\beta_2^2~\dfrac{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}{\sigma^2}\)

　　　\(E(t^2)=V(t)+E(t)^2=1+\beta_2^2~\dfrac{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}{\sigma^2}\).

回帰変動は、\(\displaystyle S_{R}=t^2 \sigma^2\) であるので、

　　　\(E(S_R)=\sigma^2+\beta_2^2~{\displaystyle \sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)

• 帰無仮説：\(\beta_2=0\) のもとでは、
  　　　回帰変動は \(\sigma^2\) に近い値になりますが、

• 対立仮説：\(\beta_2\neq 0\) のもとでは、
  　　　\(\sigma^2\) より \(\beta_2^2~{\displaystyle \sum_{i=1}^{n}( x_{i}-\overline{x})^2}\) だけ大きくなります。

シミュレーションで確かめます。

＜対立仮説の回帰モデル＞.

　　説明変数：　\(X=(x_i)=(0,1,2,3,4)\)
　　母回帰式：　\(y=1+0.5\cdot x\)
　　誤差変動：　\(\varepsilon_i \sim N(0,~0.5^2)\)

SR <- c()
SE <- c()
X <- c(0,1,2,3,4)
for (i in 1:1000){
  e <- rnorm(5,0,0.5)
  Y <- 0.5*X+1+e
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ X, data = data )
  f <- fit$ fitted.values#予測値
  a <- fit$ coefficients[1]#回帰係数
  b <- fit$ coefficients[2]#回帰係数  
  e <- fit$ residuals#残差
  se <- sum(e^2)#残差平方和
  SR <- c(SR,var(f)*4)
  SE <- c(SE,se)
}
mean(SR)

## [1] 2.743447

mean(SE)

## [1] 0.7795318

• \(E(S_{R})=E(t^2)\sigma^2\).
  
• \(E(t^2)=1+\beta_2^2~\dfrac{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}{\sigma^2}\)

今、\(E(t^2)=1+10=11\)

対立仮説のもとでは、

• 回帰変動 \(S_R\) は、\(\sigma^2=0.25\) の約11倍になりますが、
  
• 残差変動 \(S_e\) は、\(\sigma^2=0.25\) の約3倍で変化ありません。

自由度で割って分散にすると、
残差分散は \(\sigma^2\)と変わらないのに、回帰分散は \(\sigma^2\)の11倍になっています。

回帰分散／残差分散 は、11 に近い大きな値を取ることになり、このことを利用して回帰の有意性を検出しようというのがF検定です。つまり、対立仮説のもとでは回帰変動が大きくなることを用いて、回帰の有意性を調べようという訳です。

しかし、帰無仮説のもとでのP値を求めるには、F値の理論分布を調べておく必要があります。

9.4 F値の分布.

F統計量は分散比です。

帰無仮説のもとでは、
　回帰変動を自由度１で割った回帰分散 \(S_R/1\) は誤差分散に近く、
　残差変動を自由度３で割った残差分散 \(S_e/3\) は誤差分散に近く、
その商は1に近い統計量になる。

対立仮説のもとでは、
　回帰分散の増加量は、回帰係数に比例して大きくなり、
　残差分散は、影響を受けず不変であるので、
回帰分散／残差分散　が大きくなることを用いて回帰を検定する。

このような考え方になりますが、回帰の有意性を確率評価するには、帰無仮説のもとでどのような確率分布に従うのかを明確にしておく必要があります。

＜F値の定義＞.

　\(U_1\sim \chi^2(k_1)~,~U_2\sim \chi^2(k_2):i.i.d.\)　であるとき、

　　　\(F=\dfrac{U_1/k_1}{U_2/k_2}\sim F(k_1,k_2)\).

F値は ２つの独立なカイ二乗統計量の商によって構成されます。
F分布の確率密関数は分かっており、自由度に従ってP値を正確に求めることができます。
つまり、有意性の確率評価が可能になります。

＜回帰分析のF値＞.

単回帰分析の 回帰変動 \(S_R\) 、残差変動 \(S_e\) に適用すると、

\(E(S_R)=1\sigma^2\)より、
　　　\(S_R\)の自由度は 1 　\(\therefore \dfrac{S_R}{\sigma^2}\sim \chi^2(1)\).

\(E(S_e)=(n-2)\sigma^2\)より、
　　　\(S_e\)の自由度は \((n-2)\) 　\(\therefore \dfrac{S_e}{\sigma^2}\sim \chi^2(n-2)\).

辺々割ったとき\(\sigma^2\)が約分されることに注意して、

　　　\(F=\dfrac{S_R/1}{S_e/(n-2)}\sim F(1,n-2)\).

これが単回帰分析のF値です。

9.1「帰無仮説の回帰モデル」では、
　　　分散比 \(F=\dfrac{S_R}{S_e/3} \sim F(1,3)\)
　　　自由度(1,3)のF分布に従います。

Rでは、自由度 \((k_1, k_2)\) のF分布の確率密度関数は、

　　　関数　\(df(x,~k_1, k_2)\)　　で与えられます。

帰無仮説のもとでは、F分布の理論曲線と分散比のヒストグラムがほぼ一致していることを確認しましょう。 .

SR <- c()
SE <- c()
X <- c(0,1,2,3,4)
for (i in 1:500){
  u <- rnorm(5,0,0.5)
  Y <- 0*X+1+u
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ X, data = data )
  f <- fit$ fitted.values#予測値
  a <- fit$ coefficients[1]#回帰係数
  b <- fit$ coefficients[2]#回帰係数  
  e <- fit$ residuals#残差
  se <- sum(e^2)#残差平方和
  SR <- c(SR,var(f)*4)
  SE <- c(SE,se) 
}
hist(SR/(SE/3),breaks="scott",prob=T)
curve(df(x,1,3),add=T)

mean(SR/(SE/3))

## [1] 3.260003

Note.
F分布 \(F(n,m)\) の理論上の期待値は \(m/m-2\) になります。
上記の例では、\(m/m-2=3/1=3\).
回帰分散と残差分散がほぼ等いのに、その比が1にならないのは奇妙な感じがするのですが、これはデータ数が極端に少ないためで、十分なデータ数を分析すれば分散比は1に近くなります。

説明変数の値を増やし、残差の自由度が 5 になるようなモデルでは、
F値の期待値は　\(m/m-2=5/3\fallingdotseq 1.67\)　になるはずです。
確かめてみましょう。

SR <- c()
SE <- c()
X <- c(0,1,2,3,4,5,6)
for (i in 1:500){
  u <- rnorm(7,0,0.5)
  Y <- 0*X+1+u
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ X, data = data )
  f <- fit$ fitted.values#予測値
  a <- fit$ coefficients[1]#回帰係数
  b <- fit$ coefficients[2]#回帰係数  
  e <- fit$ residuals#残差
  se <- sum(e^2)#残差平方和
  SR <- c(SR,var(f)*6)
  SE <- c(SE,se) 
}
hist(SR/(SE/5),breaks="scott",prob=T)
curve(df(x,1,5),add=T)

mean(SR/(SE/5))

## [1] 1.478439

ほぼ理論値に近い値になっています。

10 回帰の検出力.

対立仮説「\(\beta_2\neq 0\)」のモデルによってデータを得たとき、F値がどのように反応するかシミュレーションで確かめましょう。

• 残差分散は誤差分散 \(\sigma^2\) に近い値ですが、
  
• 回帰分散は　\(E(S_R)=\sigma^2+\beta_2^2~{\displaystyle \sum_{i=1}^{n}( x_{i}-\overline{x})^2}\).

F値＝回帰分散／残差分散　が大きくなれば、回帰は有意であると判断されます。

9.3の対立仮説では \(\beta_2=0.5\) であるので、

• \(E(S_R)=0.5^2+0.5^2\cdot 10=0.5^2\cdot 11\).
• \(E(S_e/3)=0.5^2\).

\(\therefore~ F=\dfrac{S_R}{S_e/3}\fallingdotseq 11\)　と予想されます。

シミュレーションで確かめましょう。

SR <- c()
SE <- c()
X <- c(0,1,2,3,4)
for (i in 1:50){
  e <- rnorm(5,0,0.5)
  Y <- 0.5*X+1+e
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ X, data = data )
  f <- fit$ fitted.values#予測値
  a <- fit$ coefficients[1]#回帰係数
  b <- fit$ coefficients[2]#回帰係数  
  e <- fit$ residuals#残差
  se <- sum(e^2)#残差平方和
  s <- se/3#残差分散
  SR <- c(SR,var(f)*4)
  SE <- c(SE,se/3)
}
c(mean(SR),mean(SE))

## [1] 3.2402631 0.2113171

mean(SR)/mean(SE)

## [1] 15.33366

平均の比は 11 に近いですが、

mean(SR/SE)

## [1] 38.19264

比の平均は、かなり大きくなっています。
状況を把握するため、(SR,SE) の散布図とSR/SEの値をプロットしてみます。

par(mfrow=c(1,2))
plot(SR,SE)
plot(SR/SE)

分散比にはかなり大きい値があり、この極端に大きい値が平均値に影響しているようです。

10.1 検出力の算出

「対立仮説のもとでの分散比のヒストグラム」と「帰無仮説のもとで従うF分布曲線」を描き、検出力を算出します。

• 自由度(1,3)のF分布の有意水準5％棄却境界は “qf(0.95,1,3)” で出力され、

• 分散比が棄却境界を超えるデータ数は “sum(SR/SE>a)” で出力されます。

a <- qf(0.95,1,3)
hist(SR/SE,prob=T,xlim=c(0,100),breaks=seq(0, 3000, by = 5))
curve(df(x,1,3),add=T)
abline(v=a,col="red")

sum(SR/SE>a)/50#検出力

## [1] 0.68

赤の縦線の右側が棄却域であり、ヒストグラム全体に対するこの部分の割合が検出力です。
この場合、検出力は高くありません。

説明変数は \(X=0,1,2,3,4~\)で、各々に対して1回のデータがあるのみ、データ数は５個です。
検出力が低い理由として、データ数が少ないことも考えられます。

10.2 検出力を上げる.

各説明変数に対する出力を2回行い、データ数を2倍の10個にしてみます。

SR <- c()
SE <- c()
X <- c(0,1,2,3,4)
for (i in 1:50){
  DATA=data.frame(X=c(),Y=c())
  for (j in 1:2){
    e <- rnorm(5,0,0.5)
    Y <- 0.5*X+1+e
    data <- data.frame(X=X,Y=Y)
    DATA <- rbind(DATA, data)
  }
  fit = lm( formula = Y ~ X, data = DATA )
  f <- fit$ fitted.values#予測値
  a <- fit$ coefficients[1]#回帰係数
  b <- fit$ coefficients[2]#回帰係数  
  e <- fit$ residuals#残差
  se <- sum(e^2)#残差平方和
  SR <- c(SR,var(f)*9)
  SE <- c(SE,se/8)
}
c(mean(SR),mean(SE))

## [1] 5.618201 0.282105

c(mean(SR)/mean(SE),mean(SR/SE))

## [1] 19.91528 24.51279

par(mfrow=c(1,2))
plot(SR,SE)
plot(SR/SE)

a <- qf(0.95,1,8)
hist(SR/SE,prob=T,xlim=c(0,100),breaks=seq(0, 3000, by = 5))
curve(df(x,1,8),add=T)
abline(v=a,col="red")

sum(SR/SE>a)/50#検出力

## [1] 1

説明変数に対し2回データを与えただけで、検出力は跳ね上がっています。
検出力を上げるには、データ数も重要です。

Rには検出力を分析するパッケージ “pwr” があり、回帰分析に必要なサンプルサイズを算出できるということです。
将来、Rに慣れてきたら挑戦してください。

参考サイト.
初学者の統計学習ノート.

11 t値とF値の関係.

単回帰分析の場合、回帰係数のt値の平方はF値に一致します。

簡潔に説明します。

帰無仮説のもとでの回帰係数の標準化変数は、

　　　\(z=\dfrac{\widehat\beta_{2}}{\sqrt{V(\widehat\beta_{2})}} \sim N(0,1)\)
\(V\left(\widehat\beta_{2}\right)=\dfrac{\sigma^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}\)　ですが、誤差分散 \(\sigma^2\) を推定値である残差分散 \(s^2=\dfrac{S_e}{n-2}\) で代用すると、自由度(n-2)のt分布に従います。

即ち、
　　　\(z^2=\dfrac{\widehat\beta_{2}^2}{V(\widehat\beta_{2})}=\dfrac{\widehat\beta_{2}^2~\sum_{i=1}^{n}( x_{i}-\overline{x})^2}{\sigma^2}\)　において、\(\sigma^2\) を残差分散 \(s^2=\dfrac{S_e}{n-2}\) で代用すると、\(z^2\)は \(t^2\)になります。

　　　\(t^2=\dfrac{\widehat\beta_{2}^2~\sum_{i=1}^{n}( x_{i}-\overline{x})^2}{S_e/(n-2)}\).

\(\displaystyle S_{R} =\widehat\beta_{2}^2~{\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) ^{2}}\)　であるので、.
　　　\(t^2=\dfrac{S_R}{S_e/(n-2)}=F\).

＜t分布の定義＞.
　\(Z\sim N(0,1)\)：標準正規分布 、
　\(W\sim \chi^2(m)\)：自由度 \(m\) のカイ二乗分布、
　\(Z,W\) は互いに独立のとき、
　　　\(t=\dfrac{Z}{\sqrt{\left. W \right/ m}} \sim t(m)\)：自由度 \(m\) のt分布
この定義に従った説明をすると、

• \(z=\dfrac{\widehat\beta_{2}}{\sqrt{V(\widehat\beta_{2})}} \sim N(0,1)\)
  
• \(\dfrac{S_{e}}{\sigma^2} \sim \chi^2(n-2)\)より、
  \(t=\dfrac{\dfrac{\widehat\beta_{2}}{\sqrt{V(\widehat\beta_{2})}}}{\sqrt{\left. \dfrac{S_{e}}{\sigma^2} \right/ (n-2)}} \sim t(n-2)\).

\(\dfrac{S_{R}}{\sigma^2}=\left(\dfrac{\widehat\beta_{2}}{\sqrt{V(\widehat\beta_{2})}}\right)^2\) であるので、
　　　\(t^2=\dfrac{\dfrac{S_{R}}{\sigma^2}}{\left. \dfrac{S_{e}}{\sigma^2} \right/ (n-2)}=\dfrac{S_R}{S_e/(n-2)}\).
一方、

• \(\dfrac{S_{R}}{\sigma^2}=\left(\dfrac{\widehat\beta_{2}}{\sqrt{V(\widehat\beta_{2})}}\right)^2 \sim \chi^2(1)\).
  
• \(\dfrac{S_{e}}{\sigma^2} \sim \chi^2(n-2)\)より、
  \(F=\dfrac{S_R/\sigma^2/1}{S_e/\sigma^2/(n-2)}=\dfrac{S_R}{S_e/(n-2)}\)

従って、\(t^2=F\) が成立します。

11.1 決定係数との関係.

単回帰分析の t 値と F 値は回帰の有意性の指標であり、\(t^2=F\) によって繋がっています。

適合度の指標である決定係数 \(R^2\)との関係はどうでしょうか。

回帰係数が有意ならば相関係数も有意であり、相関係数が1に近ければ \(x,y\) は直線的に分布することになります。回帰の適合度も高くなると想像できます。

\(t^2=\dfrac{S_R}{S_e/(n-2)}\)　及び　\(R^2=\dfrac{S_R}{S_T}\)　より、

　　　\(t^2=\dfrac{S_R}{(S_T-S_R)/(n-2)}=\dfrac{S_R/S_T}{(1-S_R/S_T)/(n-2)}\).

即ち、\(t^2=\dfrac{R^2}{(1-R^2)/(n-2)}\)　　　\(\therefore~R^2=\dfrac{t^2}{t^2+(n-2)}\).

決定係数は、t値から算出されます。（nは標本数）
適合度の指標が有意性の指標と繋がったことになります。

しかし、この式では nが大きくなると決定係数は小さくなるように見えます。
標本数 n が大きければ t 値も大きくなり、\(R^2\)も１に近くなると考えられますが、標本数 n が t値と\(R^2\)に与える影響を見ておく必要がありそうです。

11.2 標本数との関係.

引き続き「2.2 回帰分析のモデル」を用います。

　　説明変数：　\(x_i=(0,1,2,3,4)\)
　　回帰モデル：\(y_i=1+0.5 x_i+\varepsilon_i\)
　　誤差変動：　\(\varepsilon_i \sim N(0,~0.5^2)\)

• このモデルから５個の標本を生成し、data とします。
  data にlm( )関数を適用すると、t値と \(R^2\)が得られます。

• 再度、生成した標本をつけ加えた10個の標本を DATA とします。
  このDATAに lm( )関数を適用すると、t値と \(R^2\)が得られます。

この操作を 繰り返すことで標本数を順次増やし、t値と \(R^2\)がどのように変化するかをシミュレーションします。

繰り返して得られた t値の列を TT、P値の列を PP、\(R^2\)の列をRR とします。
次のコードでは、繰り返しを10回行い、変化をplot( )関数で描画します。
t値、P値、決定係数は summary から取り出します。

PP=c()
RR=c()
TT=c()
DATA=data.frame(X=c(),Y=c())
X <- c(0,1,2,3,4)
for (i in 1:10){
  e <- rnorm(5,0,0.5)
  Y <- 0.5*X+1+e
  data <- data.frame(X=X,Y=Y)
  DATA <- rbind(DATA, data)
  fit = lm( formula = Y ~ X, data = DATA )
  t_value <- summary(fit)$coefficients[, "t value"]
  p_values <- summary(fit)$coefficients[, "Pr(>|t|)"]
  r_squared <- summary(fit)$r.squared
  p <- as.numeric(p_values[2])
  t <- as.numeric(t_value[2])
  R2 <- as.numeric(r_squared)
  TT <- c(TT,t)
  RR <- c(RR,R2)
  PP <- c(PP,p_values)
}

＜t値とP値の変化＞.

par(mfrow=c(1,2))
plot(TT,type="l",ylab="t-value")
plot(PP,type="l",ylab="p-value")

標本数が増えると t値は大きくなり、P値は急速に0に近くなります。
標本数が少ないと t値が安定しないことも、シミュレーションを繰り返すと観察されます。
安定した結果を得るには、標本数を大きくすることが必要なようです。

＜決定係数とF値の変化＞.

par(mfrow=c(1,2))
plot(RR,type="l",ylim=c(0,1),ylab="R-squared")
plot(TT^2,type="l",ylab="F-value")

シミュレーションを繰り返すと、決定係数が減少することもあります。
標本数が増えれば、t値とF値は大きくなりますが、決定係数は大きくなるとは限らないようです。

標本数が少ないのに決定係数が大きくなることも、シミュレーションを繰り返すと観察されます。
有意性と適合度は分析の２つの側面であり、別の指標と考える必要があるようです。

11.3 回帰係数との関係.

適合度は回帰係数の大きさと関係がある可能性もあります。
回帰係数が大きければ目的変数の変動が大きくなり、回帰直線からのズレが大きくなる可能性があります。

回帰係数を 0.2 としてシミュレーションします。

PP=c()
RR=c()
TT=c()
DATA=data.frame(X=c(),Y=c())
X <- c(0,1,2,3,4)
for (i in 1:10){
  e <- rnorm(5,0,0.5)
  Y <- 0.2*X+1+e
  data <- data.frame(X=X,Y=Y)
  DATA <- rbind(DATA, data)
  fit = lm( formula = Y ~ X, data = DATA )
  t_value <- summary(fit)$coefficients[, "t value"]
  p_values <- summary(fit)$coefficients[, "Pr(>|t|)"]
  p <- as.numeric(p_values[2])
  t <- as.numeric(t_value[2])
  R2 <- t^2/(t^2+3)
  TT <- c(TT,t)
  RR <- c(RR,R2)
  PP <- c(PP,p_values)
}

＜t値とP値の変化＞.

par(mfrow=c(1,2))
plot(TT,type="l",ylab="t-value")
plot(PP,type="l",ylab="p-value")

t値とP値の変化は、係数が 0.5 のときとは異なっています。

＜決定係数とF値の変化＞.

par(mfrow=c(1,2))
plot(RR,type="l",ylim=c(0,1),ylab="R-squared")
plot(TT^2,type="l",ylab="F-value")

標本数とともに決定係数は緩やかに上昇します。 F値は t値と同じ形の変動となります。

このほか、モデルの誤差変動の大きさも影響すると考えられます。
これらは探究課題として取り組んでください。

シミュレーションをしておくことで、分析結果の判断は慎重に行わなければならないことが分かります。

［note］.
　このシミュレーションから、いくつかの教訓が得られます。

• 標本数が増えると回帰直線から離れたデータが増えるため、決定係数が低下する場合がある。
  しかし、モデルが悪いとは限らない。
• 標本数が少ないと決定係数が高く出ることがある。
  しかし、結果は安定せず信頼性に欠ける。

モデルが適合しているとは、
　　「標本数が大きいときでも決定係数が高く安定していることである」
と言えそうです。

12 ２変量正規分布との関係

回帰分析では、得られた2変量データ \((x,y)\) に最もよく当てはまる回帰直線を求めます。
2変量の標本データの背後には、2変量正規分布に従う母集団が存在します。

• 2変量正規分布に従う母集団からの標本データに、回帰分析を適用して得られるのが標本回帰直線

• 2変量正規分布に従う母集団全体に、回帰分析を適用して得られるのが母回帰直線

と考えられます。

回帰分析の目的は y の予測であり、予測値は期待値であると考えれば、

• 2変量正規分布で x を指定したときの y の条件付き期待値は、母回帰直線の値.

となるはずです。

これらの数理は少しハードルが高いです。しかし、２次元で捉えることで理解は飛躍的に深まります。高い次元から俯瞰すると関連が見えてくるからです。

Rには rmvnorm( ) 関数という多変量正規乱数を生成する関数があります。
この関数をツールとして、2変量正規分布を体験的に理解することに挑戦しましょう。

12.1 2変量正規乱数.

平均を mean、分散共分散行列を sigma とすると、

　　　rmvnorm(n,mean,sigma)　によって n 個の２変量正規乱数を生成することができます。

ここで、分散共分散行列について再度確認しておきます。（4.1 参照）

2変量確率変数 \(X=\begin{pmatrix}X_1\\ X_2\end{pmatrix}\)　の期待値を \(m=\begin{pmatrix}m_1\\ m_2\end{pmatrix}\)　としたとき、

　　\(Cov\left(X\right)=E\left((X-m)(X-m)'\right)\) を分散共分散行列といます。

成分で計算すると、

　\(\begin{split}(X-m)(X-m)'&=\begin{pmatrix} X_1-m_1 \\ X_2-m_2 \end{pmatrix} \begin{pmatrix} X_1-m_1 & X_2-m_2 \end{pmatrix}\\&= \begin{pmatrix} (X_1-m_1)^2 & (X_1-m_1)(X_2-m_2) \\ (X_1-m_1)(X_2-m_2) & (X_2-m_2)^2 \end{pmatrix} \end{split}\)

\(\therefore Cov\left(X\right)=\begin{pmatrix} V(X_1) & Cov(X_1,X_2) \\ Cov(X_1,X_2) & V(X_2) \end{pmatrix}\)

対角成分は分散、対角でない成分は共分散になります。
２変量の確率変数では分散のほかに変量相互の関連を示す共分散を求める必要があります。

＜Note＞
分散共分散行列は、一次元の分散の自然な拡張になっています。実際、

１変量確率変数 \(X=\begin{pmatrix}X_1\end{pmatrix}\)　の期待値を \(m=\begin{pmatrix}m_1\end{pmatrix}\)　としたとき、

　\(\begin{split}(X-m)(X-m)'&=\begin{pmatrix} X_1-m_1 \end{pmatrix} \begin{pmatrix} X_1-m_1 \end{pmatrix}\\&= \begin{pmatrix} (X_1-m_1)^2 \end{pmatrix} \end{split}\)

\(\therefore Cov\left(X\right)=\begin{pmatrix} V(X_1) \end{pmatrix}\)

• １変量確率変数の確率分布は、平均と分散が与えられれば決まり、

• ２変量確率変数の確率分布は、平均ベクトルと分散共分散行列が与えられれば決まります。

12.1.1 共分散と相関係数

それでは、実際に２変量正規乱数を出力し、散布図でイメージをつかみましょう。

rmvnorm( ) 関数を使用するには、mvtnorm パッケージを導入します。
“library(mvtnorm)” がプログラムコードです。

指定するパラメータは \((x,y)\) の平均ベクトル mean と分散共分散行列 sigma ですが、x,y の関連の度合いである共分散 \(C(x,y)\) が必要になります。
共分散は x,y の連動性を表す指標で、高校では相関係数のセクションで学習しています。

• 共分散：\(\displaystyle Cov(x,y)=\dfrac{1}{n}\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})\).

• 相関係数：\(r_{xy}=\dfrac{Cov(x,y)}{s_x~s_y}\).

標準偏差 \(s_x,~s_y\) が固定されていると、共分散が決まれば相関係数が決まります。
この状況では、共分散の指定と相関係数の指定は同等の意味を持ちます。
（詳しくは ひょっこり統計「相関係数」を参照）

特に、 \(s_x=1,s_y=1\) という状況では、共分散と相関係数は等しくなります。
分散共分散行列を指定することによって、分散と相関係数を指定しているとイメージすると良いでしょう。

12.1.2 r=0 及び r=0.7

相関係数が r=0 の場合と r=0.7 の場合に、rmvnorm( ) 関数を用いて２変量正規乱数を出力します。
r=0 の場合の分散共分散行列を sigma.zero、
r=0.7 の場合の分散共分散行列を sigma とすると、

sigma.zero <- matrix(c(1,0,0,1), ncol=2);sigma.zero

##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

sigma <- matrix(c(1,0.7,0.7,1), ncol=2);sigma

##      [,1] [,2]
## [1,]  1.0  0.7
## [2,]  0.7  1.0

rmvnorm( )を用いて乱数を各々10000個出力し、散布図を並列します。

sigma.zero <- matrix(c(1,0,0,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma.zero) 
par(mfrow = c(1,2))
plot(x,xlim=c(-4,4),ylim=c(-4,4),asp=1,col="grey") 
sigma <- matrix(c(1,0.7,0.7,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)  
plot(x,asp=1,col="grey")

r=0.7 の場合は、xとyに関係性があることが分かります。

12.2 母回帰直線.

このデータを母集団とした母回帰直線を lm( )関数を用いて求めます。

r=0 のとき.

sigma.zero <- matrix(c(1,0,0,1), ncol=2)
x0 <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma.zero) 
data <- data.frame(X=x0[,1],Y=x0[,2])
model <- lm(Y ~ X, data = data)
plot(x0,asp=1,col="grey")
abline(model,col="red")
abline(v=1.5,col="blue")

model$coef

##   (Intercept)             X 
## -0.0048748692  0.0006233494

母回帰直線は、\(y=0\cdot x\)
回帰係数は０に近く、回帰が有意でない状態です。

r=0.7 のとき

sigma <- matrix(c(1,0.7,0.7,1), ncol=2)
x0.7 <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma) 
data <- data.frame(X=x0.7[,1],Y=x0.7[,2])
model <- lm(Y ~ X, data = data)
plot(x,asp=1,col="grey")
abline(model,col="red")
abline(v=1.5,col="blue")

model$coef

## (Intercept)           X 
## 0.001370394 0.702840334

母回帰直線は、\(y=0.7x\)
標準化されているので、母回帰係数と相関係数は一致します。
母回帰直線は、x を固定したときの y 推定値を与えます。

例えていうと、相関係数はサツマイモの太さのようなもので、右斜めに置いてサツマイモを切ったときに、切り口の中点を繋いだ線が母回帰直線です。

12.3 回帰モデル

実際の回帰分析は、この母集団において x を指定して抽出した標本に対して行われ、結果として標本回帰直線を得ることになります。

x を指定して抽出した標本を \((x,y)\) とするとき、 \(y\) はどのような確率分布に従うでしょうか。

この分布は x を指定したときの y の条件付き分布と呼ばれ、
　　条件付き期待値：\(E(y|x)\)、条件付き分散：\(V(y|x)\) を持ちます。

二変量正規分布で x を指定したときの y の条件付き期待値は母回帰直線の値に一致。

　\(\begin{pmatrix} x \\ y \end{pmatrix} \sim N_2\left( \begin{pmatrix} \mu_x \\ \mu_y \end{pmatrix},\begin{pmatrix} \sigma_x^2 & \sigma_{xy}\\ \sigma_{xy} & \sigma_y^2 \end{pmatrix}\right) = N_2\left( \begin{pmatrix} \mu_x \\ \mu_y \end{pmatrix},\begin{pmatrix} \sigma_x^2 & \rho\sigma_x \sigma_y \\ \rho\sigma_x \sigma_y & \sigma_y^2 \end{pmatrix}\right)\) のとき、

条件付き確率分布：\(f(y|x)\) は、

　　　\(f(y|x)\sim N\left(\mu_y+\rho \sigma_y\left(\dfrac{x-\mu_x}{\sigma_x}\right),(1-\rho^2)\sigma_y^2\right)\)

• \(E(y|x)=\mu_y+\rho \sigma_y\left(\dfrac{x-\mu_x}{\sigma_x}\right)\)
  
• \(V(y|x)=(1-\rho^2)\sigma_y^2\)

導出については、ひょっこり統計「2変量正規分布と線形回帰分析」を参照してください。
（高校生の方は微積・線形代数の学習が進んでから再学習することを勧めます）

従って、
　x を指定して抽出した標本 \((x,y)\) が従う回帰モデルは、

• \(y_i=\mu_y+\rho \sigma_y\left(\dfrac{x_i-\mu_x}{\sigma_x}\right)+\varepsilon_i\)
  　　 \(\varepsilon_i \sim N\left(0,~(1-\rho^2)\sigma_y^2\right) ~:~ i.i.d.\).

　　　条件付き分散が回帰モデルの誤差変動になります。

「12.2　r=0.7 のとき」 を例とすると、

• 回帰係数は相関係数と同じ 0.7

• 回帰モデルの誤差分散は、
  \(V(y|x)=(1-0.7^2)\cdot 1^2=0.51\)、　標準偏差は \(\sqrt{0.51}=0.714\)

この２変量正規分布に対して定まる回帰モデルは、

　　　\(y_i=0.7 x_i+\varepsilon_i\)　　\(\varepsilon_i \sim N(0,~0.714^2) ~:~ i.i.d.\)

２変量正規分布が与えられると母回帰直線が決まり、回帰モデルが一意に定まることが分かります。

標本回帰直線がどのような動きをするか、シミュレーションでみてみましょう。

12.4 標本回帰直線

２変量正規分布から得られた回帰モデル

　　　\(y_i=0.7 x_i+\varepsilon_i\) ： \(\varepsilon_i \sim N(0,~0.51)\) 　に従い、

　\(X=(x_i)\) に対して \(y_i\) を出力します。

ここで注目するのは、

• 説明変数 X の取り方によって標本回帰直線がどのように変わるか.
  
• 標本の個数によって標本回帰直線の変動はどのようになるか.

これらは回帰分析の精度に関わることで、調査を設計する上で重要です。

12.4.1 X ={ -2,-1,0,1,2 } のとき.

抽出した標本を青点で、標本回帰直線を青線で表します。
母回帰直線（赤）とのずれに注目します。

sigma <- matrix(c(1,0.7,0.7,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)  
data <- data.frame(X=x[,1],Y=x[,2])
plot(x,asp=1,col="grey",xlim=c(-4,4),ylim=c(-4,4),xlab="",ylab="")
X <-rep(seq(-2,2,by= 1),1)
e <- rnorm(5,0,0.714)
Y <- 0.7*X+e
dat <- data.frame(X=X,Y=Y)
par(new=T)
plot(X,Y,asp=1,col="blue",xlim=c(-4,4),ylim=c(-4,4))
model <- lm(Y ~ X, data = dat)
abline(model,col="blue")
abline(a=0,b=0.7,col="red")

model$coef

## (Intercept)           X 
## -0.07541866  0.37717420

抽出された点に対して標本回帰直線が求まり、母回帰直線から若干変動していることが分かります。
この誤差変動を小さくするように回帰分析を設計することが重要になります。

方法としては、x ごとの調査回数を増やすことが考えられます。

12.4.2 X ={ -2,-1,0,1,2 } を n 回.

調査回数を n とします。
X の値を n 回繰り返し、同じ x に対し y を n 回出力します。

プログラムコードではリピート関数 rep( ) を用い、
　　　X <- rep(seq(-2,2,by= 1),n) 　としています。

こうすることで、n を変えてシミュレーションすることが可能になります。

まず、n=1 として実行します。
n=1 のもとでの試行を100回繰り返し、標本回帰直線の変動を観察します。

plot(x,asp=1,col="grey",xlim=c(-4,4),ylim=c(-4,4))
n <- 1
X <-rep(seq(-2,2,by= 1),n)
for (i in 1:100){
  e <- rnorm(5*n,0,0.714)
  Y <- 0.7*X+e
  dat <- data.frame(X=X,Y=Y)
  par(new=T)
  model <- lm(Y ~ X, data = dat)
  par(new=T)
  abline(model,col="lightblue")
}
abline(a=0,b=0.7,col="red")

標本回帰直線の変動はかなり大きく、y の予測精度は良くありません。

次に、n=4 として実行します。
n=4 のもとでの試行を100回繰り返し、標本回帰直線の変動を観察します。

plot(x,asp=1,col="grey",xlim=c(-4,4),ylim=c(-4,4))
n <- 4
X <-rep(seq(-2,2,by= 1),n)
for (i in 1:100){
  e <- rnorm(5*n,0,0.714)
  Y <- 0.7*X+e
  dat <- data.frame(X=X,Y=Y)
  par(new=T)
  model <- lm(Y ~ X, data = dat)
  par(new=T)
  abline(model,col="lightblue")
}
abline(a=0,b=0.7,col="red")

標本回帰直線の変動は小さくなります。

予測の精度を上げるには、標本数を多くすることが有効と考えられます。

このほか、xの数を増やすことも有効です。

12.4.3 ＜探究課題＞.

上の例では、X ={ -2,-1,0,1,2 } は5個ですが、
X ={ -4,-3,-2,-1,0,1,2,3,4 } の9個に増やすとどうでしょうか。
n=1 のとき、n=4 のとき、それぞれシミュレーションしましょう。

12.4.3.1 X ={ -4,-3,-2,-1,0,1,2,3,4 } を n 回.

＜n=1 とすると＞、

plot(x,asp=1,col="grey",xlim=c(-4,4),ylim=c(-4,4))
n <- 1
X <-rep(seq(-4,4,by= 1),n)
for (i in 1:100){
  e <- rnorm(9*n,0,0.714)
  Y <- 0.7*X+e
  dat <- data.frame(X=X,Y=Y)
  par(new=T)
  model <- lm(Y ~ X, data = dat)
  par(new=T)
  abline(model,col="lightblue")
}
abline(a=0,b=0.7,col="red")

＜n=4 とすると＞、

plot(x,asp=1,col="grey",xlim=c(-4,4),ylim=c(-4,4))
n <- 4
X <-rep(seq(-4,4,by= 1),n)
for (i in 1:100){
  e <- rnorm(9*n,0,0.714)
  Y <- 0.7*X+e
  dat <- data.frame(X=X,Y=Y)
  par(new=T)
  model <- lm(Y ~ X, data = dat)
  par(new=T)
  abline(model,col="lightblue")
}
abline(a=0,b=0.7,col="red")

標本回帰直線の変動は小さく、予測精度は高くなります。

［Note］

• 調査する x を幅広く取り、調査回数を増やせば、推定精度は上がる。.

13 平均への回帰.

同じ対象に検査を2回行うとします。

• 1回目の検査結果 \(X\) の母集団平均からの偏差が \(|X-\mu|\) であるとき、
  
• 2回目の検査結果 \(Y\) の偏差 \(|Y-\mu|\) は、1回目の偏差より小さい傾向があります。

この現象を「平均への回帰」といいます。
イギリスのゴルトン（1822-1911）は「平均への回帰」の数値化を試み、その結果生まれたのが「相関係数」と言われています。
（参照：統計Web「平均への回帰、相関係数―統計学史(2)」）.

どのようにして、この現象は起こるのでしょうか。

13.1 2変量正規分布による説明.

もし、2変量正規分布や相関係数の理論が確立されていれば、説明は簡潔です。

\(X,Y,\) が同じ平均 \(\mu\)と分散 \(\sigma^2\)を持つ正規分布に従い、相関係数が \(\rho\) であるとします。
\(Cov(X,Y)=\rho\sigma^2\)　であるので、.
　　　\(\begin{pmatrix} X \\ Y \end{pmatrix} \sim N_2\left( \begin{pmatrix} \mu \\ \mu \end{pmatrix},\begin{pmatrix} \sigma^2 & \rho\sigma^2\\ \rho\sigma^2 & \sigma^2 \end{pmatrix}\right)\).

\(X\) を指定したときの \(Y\) の予測値は、\(X\)に対する条件付き期待値であり、母回帰直線の値に一致するので、

　　予測値：\(Y=\mu+\rho(X-\mu)\)　　\(\therefore Y-\mu=\rho(X-\mu)\).

この式から予測値 \(Y\) の偏差は \(X\) の偏差より小さく( \(\rho\)倍)、平均へ回帰することが分かります。

しかし、相関係数と２変量正規分布の理論を確立したのは、ゴルトンのアイディアを継承したピアソンでした。
また、\(X,Y\) が同じ平均と分散を持つ正規分布に従っていても、\((X,Y)\)が2変量正規分布に従うかという問題があります。
（ひょっこり統計「多変量正規分布」参照）.
ゴルトンが注目したのは「測定誤差」であると言われています。 この誤差に注目したアプローチを考えます。

統計検定１級 2022年統計応用の問5に「平均への回帰」についての出題があります。
「平均への回帰」を誤差の確率モデルで説明しようとするものです。
（問題については統計検定公式問題集を参照してください）.

13.2 血圧測定の平均回帰.

［出典］統計検定１級 統計応用 2022年⑤.［平均への回帰］
　血圧値は測定ごとに変動するが、各人は個人ごとに血圧の真値（その人では定数）を持つとし、各人の測定値を、

　1回目：\(X=\theta+\varepsilon_1\)、　２回目：\(Y=\theta+\varepsilon_2\)　とする。

• \(\theta\sim N(\mu,~\tau^2)\)：各人の真値の母集団分布.
  
• \(\varepsilon_i\sim N(0,~\psi^2)\)：誤差変動（真値と独立）.

今、\(\mu=120\)　\(\tau^2=108,~~\psi^2=36\) として1000個の \(X,Y\) を出力し、散布図を描きます、

u <- 36 #誤差分散
theta <- rnorm(1000,mean=120,sd=sqrt(108))
e1 <- rnorm(1000,mean=0,sd=sqrt(u))
e2 <- rnorm(1000,mean=0,sd=sqrt(u))
x <- theta+e1
y <- theta+e2
plot(x,y,col="grey")
abline(a=0,b=1,col="red")

rmvnorm( )関数を使わないで散布図出力されたことに驚きました。
赤の直線は \(y=x\) です。散布図はこの直線について対称になっているように見えます。

分散、共分散、相関係数を計算します。

　　　\(E(X)=120,~~V(X)=108+36=144\).
　　　\(E(Y)=120,~~V(Y)=108+36=144\).
　　　（全分散は、真値分散と誤差分散の合計）

\(Cov(X,Y)=Cov(\theta+\varepsilon_1,\theta+\varepsilon_2)=Cov(\theta,\theta)=\tau^2\)から、
　相関係数　\(Cor(X,Y)=\dfrac{\tau^2}{\tau^2+\psi^2}=\dfrac{108}{144}=0.75\)
（相関係数は、全分散に対する真値の分散の割合になっています）

これで、問題文と同じパラメーターを持つ 2変量確率変数 \((X,Y)\) が得られたことになります。
\((X,Y)\) は 2変量正規分布に従うとして話を進めます。

13.2.1 誤差分散 小.

誤差変動が微小ならば、1回目と2回目の計測値はほぼ一致するはずです。
誤差分散を u <- 0.01 とします。

u <- 0.01 #誤差分散
theta <- rnorm(1000,mean=120,sd=sqrt(108))
e1 <- rnorm(1000,mean=0,sd=sqrt(u))
e2 <- rnorm(1000,mean=0,sd=sqrt(u))
x <- theta+e1
y <- theta+e2
plot(x,y,col="grey")
abline(a=0,b=1,col="red")

誤差がなければ相関係数は1に近く「平均への回帰」が起こる余地はありません。

ゴルトンは「誤差の存在が平均への回帰に関係している」と考えたのかもしれません。

ゴルトンのデータは、
　HistData “Galton’s data on the heights of parents and their children, by child” で見ることができます。
パソコンのない時代、分析は大変だったと思います。

13.2.2 誤差分散 大.

誤差分散を「u <- 36」と大きくし、統計検定の問題の数値に戻します。

ゴルトンの目的は「親の身長から、子の身長を予測する」ことでした。
x が観測されたときに y を予測することは、条件付き期待値 \(E(y|x)\) を求めることに相当します。

1回目の測定値が \(X=132\) である \((X,Y)\) について、 \(Y\) の平均を求めます。

• 出力した \((X,Y)\) をデータフレームに格納、
  
• \(132-0.1<X<132+0.1\) を満たす \((X,Y)\) を抽出し.
• 抽出した点の \(Y\) 平均を求めます。

条件付き期待値をシミュレーションで求めます。

u <- 36 #誤差分散
theta <- rnorm(10000,mean=120,sd=sqrt(108))
e1 <- rnorm(10000,mean=0,sd=sqrt(u))
e2 <- rnorm(10000,mean=0,sd=sqrt(u))
x <- theta+e1
y <- theta+e2
plot(x,y,col="grey",xlim=c(80,160),ylim=c(80,160))
abline(a=0,b=1,col="red")
#データフレームに格納
data <- data.frame(X=x,Y=y)
#X=132に近い点を抽出
DATA <- data[data$X <132+0.1 & 132-0.1 < data$X,]
par(new=T)
plot(DATA,col="red",xlim=c(80,160),ylim=c(80,160))

mean(DATA[,2])#抽出点のY平均

## [1] 127.9206

\(X=132\) のもとでの \(Y\) の平均は 132より小さくなっており、「平均への回帰」が起きているようです。

＜真値の条件付き分布＞.

測定値は誤差を含むので、 \(X=132\) であったとしても、その人の真値 \(\theta\) は132とは限りません。
真値の予測値は \(E(\theta\mid X=132)\) で計算されます。

\(Y=\theta+\varepsilon_2\) であるので、

　　　\(E(Y\mid X=132)=E(\theta\mid X=132)+E(\varepsilon_2\mid X=132)\).

\(E(\varepsilon_2\mid X=132)=0\)　であるので、

　　　\(E(\theta\mid X=132)＝E(Y\mid X=132)\)

\(X\) に対する条件付き期待値が母回帰直線の値に一致することを用いれば、

　　　\(E(Y\mid X=132)=120+0.75(132-120)=129\)

従って、\(E(\theta\mid X=132)=129\).

かくして「平均への回帰」は次のように説明されます。

• 観測値が \(X=132\) であったとき、真値は\(\theta=129\)と予想され、
  
• 2回目の計測値は、真値 129 に誤差変動を加えたものになる.
• よって、2回目 \(Y=129+\varepsilon_2\) は1回目の132より小さい。

2変量正規分布の理論を使うことは循環論の心配がありますが、ご容赦ください。
因みに、統計検定公式問題集の解説は、ベイズの定理を用いて \((X,\theta)\) が2変量正規分布に従うことにまで言及しています。

14 重回帰分析.

説明変数が複数ある回帰分析を重回帰分析といいます。
目的変数に影響する要因（説明変数）が複数ある場合に用いられます。

簡単のため、説明変数は２個とします。（データは n 個）
n 個のデータをベクトル成分で

\(1=\begin{pmatrix} 1\\ \vdots\\ 1 \end{pmatrix}, X_{2}=\begin{pmatrix} x_{21} \\ \vdots \\ x_{2n} \end{pmatrix}, X_{3}=\begin{pmatrix} x_{31} \\ \vdots \\ x_{3n} \end{pmatrix}, Y=\begin{pmatrix} Y_1\\ \vdots\\ Y_n \end{pmatrix}\)　と表し、

重回帰モデル.
　　　\(Y=\beta_1 1+\beta_2X_2+\beta_3X_3+u\) 　を成分で表現すると、

　　　\(\begin{pmatrix} Y_{1} \\ \vdots \\ Y_{n} \end{pmatrix}=\begin{pmatrix} 1 & x_{21} & x_{31} \\ \vdots & \vdots & \vdots \\ 1 & x_{2n} & x_{3n} \end{pmatrix}\ \begin{pmatrix} \beta_{1} \\ \beta_{2} \\ \beta_{3} \end{pmatrix}+\begin{pmatrix} u_{1} \\ \vdots \\ u_{n} \end{pmatrix}\)

　　　\(X=\begin{pmatrix} 1 & x_{21} & x_{31} \\ \vdots & \vdots & \vdots \\ 1 & x_{2n} & x_{3n} \end{pmatrix}\)　 \(\beta=\begin{pmatrix} \beta_{1} \\ \beta_{2} \\ \beta_{3} \end{pmatrix}\)　とおけば、

簡潔に　\(Y=X\beta+u\)　と表すことができます。

これは 単回帰と同じ式になります。

つまり、単回帰で成り立つ公式は、ほぼそのまま重回帰でも成り立ちます。
行列を用いた説明は重回帰にも当てはまると考えてよいでしょう。
以下で確認する重要事項は、単回帰と照らし合わせると理解がスムーズです。

14.1 正規方程式.

予測値との差である残差は　\(u=Y-X\beta\) です。

残差の大きさが最小となるよう最小二乗法で回帰係数を定めます。

＜残差平方和の行列表示＞

　\(Y'\) は \(Y\) の転置行列を表すとします。

　　　\(Y'Y=Y_1^2+Y_2^2+\cdots +Y_n^2\)　　です。

このように、平方和は転置行列との積で表すことができます。

残差平方和は　\(Q=(Y-X\beta)'(Y-X\beta)\)

　　　\(h\left( \beta\right) =Q=\left( X\beta-Y\right) '\left( X\beta-Y\right)\) とおいて \(\beta\) で微分します。

ベクトルの微分では、
　　　\(\dfrac{\partial }{\partial x}\ \left\| Ax-b\right\| ^{2}=\dfrac{\partial }{\partial x}\ \left( Ax-b\right)'\times \left( Ax-b\right) =2A'\left( Ax-b\right)\)

よって、　\(h'\left( \beta\right) =2X'\left( X\beta-Y\right)\)

\(h'( \widehat\beta ) =0\)とすると、\(X'\left( X\widehat\beta-Y\right)=0\)
　　\(\therefore~X'X\widehat\beta=X'Y\)　これを正規方程式といいます。

　　標本回帰係数は　\(\widehat\beta=(X'X)^{-1}X'Y\)　となります。

\(\begin{split} \widehat{\beta }-\beta &=\left( X'X\right) ^{-1}X'u \end{split}\)　であるので、

　　\(Cov\left( \widehat{\beta }\right)=E\left( \left( \widehat{\beta }-\beta \right) \ \left( \widehat{\beta }-\beta \right) '\right) =\sigma ^{2}\left( X'X\right) ^{-1}\)

説明は「4 回帰係数の変動」を参照してください。
1変量での説明が、そのまま多変量でも通用する。これが数学の醍醐味です。

＜残差の性質（「5.1 残差」参照）＞

　残差は　\(e=Y-X\widehat{\beta}\)　であるので、

　　　\(X'e=0\)　　（よって \(X,~e\) は無相関）

予測値は \(\widehat{Y}=X\widehat{\beta}\) であり

　　　\(e'\widehat{Y}=e'X\widehat{\beta}=0~\widehat{\beta}=0\)　　（よって \(\widehat{Y},~e\) は無相関） 　　　

共分散を用いると、
　　　\(Cov(X,e)=0\)　　　\(Cov(\widehat{Y},e)=0\)

＜平方和の分解公式＞

\(Cov(\widehat{Y},e)=0\) から、

• \(Cov(Y,e)=Cov(\widehat{Y}+e,e)=V(e)\).

• \(Cov(Y,\widehat{Y})=Cov(\widehat{Y}+e,\widehat{Y})=V(\widehat{Y})\).

• \(Cov(Y,Y)=Cov(\widehat{Y}+e,\widehat{Y}+e)=V(\widehat{Y})+V(e)\).

　第３式は、 \(V(Y)=V(\widehat{Y})+V(e)\)　であり、データ数をかければ偏差平方和になるので、

• \(S_{T}=S_R+S_e\)

平方和の分解公式が成立します。（「5.2 共分散」参照）

14.2 回帰の適合度.

単回帰だけなら、相関係数だけで事足ります。
決定係数は \((x,y)\) の相関係数の二乗、重相関係数も相関係数の二乗であり、相関係数を見れば回帰の適合度が判断できます。

しかし、重回帰分析では複数の説明変数があります。
説明変数と目的変数の相関係数は複数個あり、どう評価するかが問題になります。

回帰の適合度を１つの統計量で評価するには、決定係数・重相関係数が必要になります。
また、p値を用いて回帰の有意性を判断するにはF統計量が必要です。

＜重相関係数＞

観測値と推測値の相関係数を重相関係数（Multiple R）といいます。
観測値が大きくなれば推測値も大きくなるので、

　　　重相関係数　\(R>0\)　と考えられます。

また　\(Cov(Y,\widehat{Y})=V(\widehat{Y})\) であることから、

　　　\(R=Cor(Y,\widehat{Y})=\dfrac{Cov(Y,\widehat{Y})}{\sqrt{V(Y)}\sqrt{V(\widehat{Y})}}=\dfrac{V(\widehat{Y})}{\sqrt{V(Y)}\sqrt{V(\widehat{Y})}}=\sqrt{\dfrac{V(\widehat{Y})}{V(Y)}}\)

この式からも、重相関係数は常に正であることが分かります。

　　　\(R^2=\dfrac{V(\widehat{Y})}{V(Y)}=\dfrac{S_R}{S_{T}}\).

これは全変動に対する回帰変動の寄与率である決定係数と一致しています。

＜決定係数＞

　（全変動）=（回帰変動）+（残差変動）　において、

全変動に対する回帰変動の割合を決定係数といいます。

分解公式　\(S_{T}=S_R+S_e\)　では、

　　　\(R^2=S_R/S_T\)　が決定係数になります。

\(R^2=\dfrac{S_T-S_e}{S_T}=1-\dfrac{S_e}{S_T}\)　であるので、

残差平方和 \(S_e\) が小さくなれば、決定係数は大きくなります。

重回帰モデルの残差平方和 \(S_e\) は、

　　　\(\displaystyle h(\beta)=\sum _{i}\left\{ Y_{i}-\left( \beta _{1}\cdot 1+\beta _{2}\cdot x_{2i}+\beta _{3}x_{3i}\right) \right\} ^{2}\)　の最小値です。
　　　これは \(\beta=(\beta_1,\beta_2,\beta_3)\) という３変数関数の最小値です。

単回帰モデルの残差平方和 \(S_e\) は、
　　　\(\beta_3=0\) という条件のもとでの条件付き最小値で、

　　　\(\displaystyle h(\beta)=\sum _{i}\left\{ Y_{i}-\left( \beta _{1}\cdot 1+\beta _{2}\cdot x_{2i}\right) \right\} ^{2}\)　の最小値です。
　　　これは \(\beta=(\beta_1,\beta_2,0)\) という２変数関数の最小値です。

重回帰の\(\beta\)の範囲は、単回帰の条件付き\(\beta\)の範囲を含むので、

　　　（重回帰の\(S_e\)） \(\leq\)（単回帰の\(S_e\)）.
従って、
　　　（重回帰の決定係数）\(\geq\)（単回帰の決定係数）.

• 説明変数が増えれば決定係数は大きくなります。
  
• 決定係数の平方根である重相関係数も大きくなります。

［注意］.
　関数 \(y=x\) について、

　　　\([-1,1]\) における最小値よりも、\([-2,2]\)における最小値の方が小さいです。

一般に、同じ関数で定義域が広がれば最小値は小さくなります。

＜自由度調整決定係数＞

説明変数の数を増やせば決定係数は大きくなりますが、統計的に有意でない説明変数を増やす必要はありません。
予測においては、説明変数は独立であることが望ましく、説明変数相互の相関が高くなると回帰係数の推定が不安定になり予測精度が下がると言われています。この現象を多重共線性といいます。.

この欠点を補うため、説明変数の数を考慮して決定係数を補正するのが自由度調整決定係数です。

重回帰分析の自由度については次のようになります。

　データは n 個、説明変数は \(k\) 個とします。

　帰無仮説（偏回帰係数はすべて０）のもとでの自由度は、

　　　回帰変動：自由度 \(k\) 　　　残差変動：自由度 \(n-k-1\)

\(R^2=1-\dfrac{S_e}{S_T}\)　の第２項を　\(\dfrac{S_e/(n-k-1)}{S_T(n-1)}\)　とした、

　　　\(R^{*2}=1-\dfrac{S_e/(n-k-1)}{S_T(n-1)}\)

を自由度調整済み決定係数（Adjusted R-squared）と呼びます。

自由度調整済み決定係数が大きくなることは、重回帰モデルの改善を意味し、
説明変数の数が異なるモデル同士を比較する際に、より公平な判断基準となります。
より少ない変数で、より高い説明力を持つモデルを評価できます。

14.3 回帰の有意性.

単回帰では、「回帰係数の有意性」は「相関係数の有意性」と同等です。

「3.1 行列表現」から、回帰モデル \(y=a+bx\) の x の標本回帰係数は、

　　　\(b=r_{xy}\dfrac{s_{y}}{s_{x}}\)
　　　（\(s_x^2,~s_y^2\) は x,y の分散、\(r_{xy}\) は相関係数）

　\(b=0\) を検定することは、\(r_{xy}=0\) を検定することと同等です。

しかし、重回帰では回帰係数が複数あります。
回帰の有意性を総合的に判断するには、F統計量が必要になります。

＜F統計量＞

F統計量は回帰モデル全体の有意性を検定する統計量であり、P値による確率評価を行います。

• F統計量= 回帰分散／残差分散

単回帰と同じ定義です。（「9 F統計量」を参照）

重回帰分析の自由度は次のようになります。

　データは n 個、説明変数は２個とします。

　帰無仮説（偏回帰係数はすべて０）のもとでの自由度は、

　　　回帰変動：自由度 2　　　残差変動：自由度 n-3.

単回帰の場合の類比と考えてください。
シミュレーションで確かめることは、各自の探究課題としましょう。

単回帰分析の 回帰変動 \(S_R\) 、残差変動 \(S_e\) とすると、

\(E(S_R)=2\sigma^2\)より、
　　　\(S_R\)の自由度は 2 　\(\therefore \dfrac{S_R}{\sigma^2}\sim \chi^2(2)\).

\(E(S_e)=(n-3)\sigma^2\)より、
　　　\(S_e\)の自由度は \((n-3)\) 　\(\therefore \dfrac{S_e}{\sigma^2}\sim \chi^2(n-3)\).

辺々割ったとき\(\sigma^2\)が約分されることに注意して、

　　　\(F=\dfrac{S_R/2}{S_e/(n-3)}\sim F(2,n-3)\).

これが説明変数が２個の重回帰分析のF値です。

F統計量では自由度が考慮されています。
自由度調整済み決定係数との違いは、F分布を用いた確率評価ができることです。
即ち、p値を用いた有意性の判断が可能になります。

［参考］
説明変数が k 個のときの F統計量と決定係数との関係は、

　　　\(F=\dfrac{S_R/k}{S_e/(n-k-1)}=\dfrac{R^2/k}{(1-R^2)/(n-k-1)}\).

　　　F 統計量と決定係数との関係　参照.

F統計量と決定係数のどちらを使うべきかは、評価の目的によって異なります。

• 自由度調整済み決定係数は、変数の数による偏りを補正してモデルの適合度を評価します。
  しかし、モデルが統計的に有意かどうかは分かりません。
  
• F統計量はモデル全体の統計的有意性を検定したいときに使用します。
  しかし、個々の変数の重要度やモデルの適合度は分かりません。

15 多重共線性.

多重共線性は、どのようにして起こるのでしょうか。
多重共線性は、説明変数相互の相関が強くなると回帰の予測精度が下がる現象です。

相関の強い２変量正規乱数を生成し、これを説明変数とする重回帰のシミュレーションを行います。
単回帰では予測値は回帰直線として視覚化されますが、重回帰では回帰平面になります。

＜回帰モデル＞

　\(z_i=\beta_1+\beta_2x_{i}+\beta_3y_{i}+\varepsilon_i\:\:(i=1,2,\cdots,n)\) 　

　　　　 \(\begin{pmatrix} x \\ y \end{pmatrix} \sim N_2\left( \begin{pmatrix} \mu_x \\ \mu_y \end{pmatrix},\begin{pmatrix} \sigma_x^2 & \rho\sigma_x \sigma_y \\ \rho\sigma_x \sigma_y & \sigma_y^2 \end{pmatrix}\right)\)

　　　　 \(\varepsilon_i\sim N(0,\:\sigma^2)\) は互いに独立

誤差がなければ、\(z=\beta_1+\beta_2x+\beta_3y\)　は空間における平面の方程式を表します。

平面の方程式に誤差を加えて、重回帰モデルを構成しています。

15.1 相関が強いとき.

説明変数 \((x,y)\) の相関が強い場合を考えます。パラメーターは、

• x の標準偏差　sx=2 、　y の標準偏差　sy=2.
  
• \((x,y)\) の相関係数　r=0.95.
  
• x の平均 mx=0 、y の平均 my=0.

このパラメーターを持つ２変量正規乱数を n=20 個出力し、説明変数とします。

library(mvtnorm)
sx <- 2
sy <- 2
r <- 0.98
n <- 20
sigma <- matrix(c(sx^2,r*sx*sy,r*sx*sy,sy^2), ncol=2)
x <- rmvnorm(n, mean=c(0,0), sigma=sigma) 
plot(x)
abline(a=0,b=0.98,col="red")

説明変数 \((x,y)\) は \(y=0.95x\) を中心に、幅の狭い帯状の範囲に出力されています。
これは \((x,y)\) の相関が強いことを意味します。

次に、回帰モデルとして、

　　　\(z_i=1+0.5x_{i}+0.3y_{i}+\varepsilon_i~~~~(~\varepsilon_i \sim N\left(0,~0.2^2\right)~)\).

を用いて目的変数 z を出力し、得られた \((x,y,z)\) に対し lm( )関数を用いて回帰分析を行います.

e <- rnorm(n,0,0.2)
z <- 0.5*x[,1]+0.3*x[,2]+1+e
data <- data.frame(X=x[,1],Y=x[,2],Z=z)
fit = lm( formula = Z ~ X+Y, data = data )
summary(fit)

## 
## Call:
## lm(formula = Z ~ X + Y, data = data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.232757 -0.069390  0.007144  0.094500  0.234676 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.08287    0.03286  32.950  < 2e-16 ***
## X            0.66293    0.08578   7.729 5.82e-07 ***
## Y            0.12867    0.08644   1.489    0.155    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1362 on 17 degrees of freedom
## Multiple R-squared:  0.9923, Adjusted R-squared:  0.9914 
## F-statistic:  1100 on 2 and 17 DF,  p-value: < 2.2e-16

回帰平面を描画するため、パッケージ “rgl” を導入します。

• rgl パッケージは3次元グラフィックスを生成するためのパッケージです。
  
• options (rgl.printRglwidget = TRUE).
  rglで生成された3次元グラフィックスの動的操作（回転）を可能にします。

library(rgl)
options(rgl.printRglwidget = TRUE) 
x <- data$X
y <- data$Y
z <- data$Z
plot3d(x = x, y = y, z = z,type = "s", col = "blue", size = 1)
#標本回帰平面の描画
coefs <- coef(fit)
#coefs
a <- coefs["X"]
b <- coefs["Y"]
c <- -1
d <- coefs["(Intercept)"]
#ax+by＋cz+d=0で表される平面
planes3d(a, b, c, d, alpha = 0.3,col="lightblue")
planes3d(0.5, 0.3, c, 1, alpha = 0.5,col="red")

薄赤は母回帰平面、薄青は標本回帰平面です。
若干のズレがあります。

標本\((x,y,z)\)が細長い棒状に並んでいるので、棒の上に乗る回帰平面は安定しないことが予想されます。
これは、説明変数がほとんど１次元の幅しか持たないためです。
x と y がほぼ同じ値なら、片方の変数のみで回帰できる可能性があります。

＜AIC：赤池情報量基準＞.
日本の統計学者である赤池弘次（1927-2002）が提唱した統計量です。
AIC（Akaike’s Information Criterion）を用いると「予測精度が高く、説明変数の少ないモデル」を選ぶことができます。
AIC最小のモデルが、バランスの取れた良いモデルであると判断されます。
ここでは、AIC の数理には立ち入りません。Rのstep( )関数で AICを利用してみましょう。

step( )関数はステップワイズ法を実行する関数です。説明変数を段階的に減少させ AICの変化を見ます。

• は説明変数を減らさないモデル.
• 「-X」は X を除外したモデルです。

step_model <- step(fit)

## Start:  AIC=-76.98
## Z ~ X + Y
## 
##        Df Sum of Sq     RSS     AIC
## <none>              0.31558 -76.982
## - Y     1   0.04113 0.35671 -76.531
## - X     1   1.10885 1.42443 -48.839

• AIC最小のモデルが、最適モデルとして summary に表示されます。

summary(step_model)

## 
## Call:
## lm(formula = Z ~ X + Y, data = data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.232757 -0.069390  0.007144  0.094500  0.234676 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.08287    0.03286  32.950  < 2e-16 ***
## X            0.66293    0.08578   7.729 5.82e-07 ***
## Y            0.12867    0.08644   1.489    0.155    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1362 on 17 degrees of freedom
## Multiple R-squared:  0.9923, Adjusted R-squared:  0.9914 
## F-statistic:  1100 on 2 and 17 DF,  p-value: < 2.2e-16

　

＜標本回帰平面の変動＞.

上記の操作を5回繰り返して、標本回帰平面の変動を見ます。

for (i in 1:5){
  x <- rmvnorm(n, mean=c(0,0), sigma=sigma)
  e <- rnorm(n,0,0.2)
  z <- 0.5*x[,1]+0.3*x[,2]+1+e
  data <- data.frame(X=x[,1],Y=x[,2],Z=z)
  fit = lm( formula = Z ~ X+Y, data = data )
  coefs <- coef(fit)
  a <- coefs["X"]#回帰係数
  b <- coefs["Y"]#回帰係数
  c <- -1
  d <- coefs["(Intercept)"]
  planes3d(a, b, c, d, alpha = 0.3,col="lightblue")
}
planes3d(0.5, 0.3, c, 1, alpha = 0.5,col="red")

15.2 相関が弱いとき.

x,y が無相関に近い r=0.1 の場合を考えます。

library(mvtnorm)
sx <- 2
sy <- 2
r <- 0.1
n <- 20
sigma <- matrix(c(sx^2,r*sx*sy,r*sx*sy,sy^2), ncol=2)
x <- rmvnorm(n, mean=c(0,0), sigma=sigma) 
e <- rnorm(n,0,0.2)
z <- 0.5*x[,1]+0.3*x[,2]+1+e
data <- data.frame(X=x[,1],Y=x[,2],Z=z)
x <- data$X
y <- data$Y
fit = lm( formula = Z ~ X+Y, data = data )
coefs <- coef(fit)
a <- coefs["X"]#回帰係数
b <- coefs["Y"]#回帰係数
c <- -1
d <- coefs["(Intercept)"]
plot3d(x = x, y = y, z = z,type = "s", col = "blue", size = 1)
planes3d(a, b, c, d, alpha = 0.3,col="lightblue")
planes3d(0.5, 0.3, c, 1, alpha = 0.5,col="red")

説明変数の相関が低ければ、標本回帰平面と母回帰平面のズレはほとんどなく、予測精度が高くなっていることがわかります。

step_model <- step(fit)

## Start:  AIC=-66.06
## Z ~ X + Y
## 
##        Df Sum of Sq     RSS     AIC
## <none>               0.5448 -66.060
## - Y     1    5.4943  6.0391 -19.949
## - X     1   16.6312 17.1761   0.956

グラフでは差が感じられないかもしれませんが、AIC の数値は.

　説明変数間の相関が強のAIC ＞ 説明変数間の相関が弱のAIC.

となっており、説明変数間の相関は弱いほうが良いことが分かります。

16 終わりに.

乱数シミュレーションを活用した学習方法を提案し、その方法を用いて本書を書いてきました。
様々な疑問が解消し、学ぶことの奥深さを感じました。

• 単回帰で回帰変動の自由度が１である理由.
  
• 平均へ回帰する理由.
  
• 単回帰のF値とt値の関係.
  
• 重相関係数の意味.
  
• 多重共線性の意味.

これらは、統計学を学ぶ人が引っ掛かりを覚えるところではないかと思います。
教材をシミュレーションで確かめることで理解が深まりました。

• 共分散を利用した分解公式の説明.
  
• 正規方程式の帰納的導入.

これらも高校生向けに考えたもので、自分ではよかったと思います。

実際の分析で最も重要で大変なのは、データを揃えることです。
実際のデータで分析する機会は限られてきます。
「自分でモデルを設定し、生成したデータ分析する」.
そこでは予測や検証といったクリティカルシンキングが必要になります。
この方法は、探究的な学習においても有用ではないかと思います。

16.1 高校生の統計学シリーズ.

「高校生の統計学」全４巻が完成しました。
これらのシリーズは「数学や統計が分からないが理解したい」という方に向いています。
是非、GoogleColabやRstudioでコードを動かしながら学んでください。

この4冊の作成は、自分の目と手で確かめながら学ぶ「看て学ぶ」の連続でした。
分かっていなかったことが見えてきて興奮する場面もありました。
「書き手が面白いんだから、読者だって面白いはず」と思うのですが、果たしてどうでしょうか。
しかし、学びというのは自分楽しむことが大切なのは確かなようです。
興味が興味を呼び、疑問が疑問を呼ぶ。分からないことは苦痛ではない。
学んでいって視界が開けた時の喜びはどんなものでしょう。
これだけでも、やってきてよかったと思います。

次のシリーズとして「情報Ⅱの統計学」を予定しています。
情報Ⅱについては、まだ具体的な教材が不足しているという現状があります。
情報Ⅱを高校生が学ぶには、様々な角度からの教材へのアプローチが必要です。
確率モデルとシミュレーションという手法を用いて教材化を進めたいと思います。