「続3・高校生の統計学」のプログラムコードを掲載しています。

  • RstudioやPosit Cloud.
  • GoogleColabのR環境.

においてコピー＆ペーストで利用できます。
初めは指導者のもとで行うのが良いでしょう。
（電子書籍はデジタル著作権管理（DRM）仕様のためコピーできません）

インストールせずにブラウザでプログラムを実行するには、次の方法があります。

以下、電子書籍との章番号を対応させるため、項目名は残してあります。

  • 必要なパッケージを用いるのに必要なコードで書籍に表示されていないものもありますが、以下では表示してあります。
    例えば、library(mvtnorm)、library(scatterplot3d).

1 はじめに

2 回帰と相関

2.1 相関分析のモデル

平均を 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")

2.2 回帰分析のモデル

＜具体例＞

　　説明変数：　\(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.4165223
1 0.7168174
2 2.3176022
3 2.6494690
4 2.5660058

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

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

実際に、lm( )を用いて回帰直線を推定し、赤で書き加えます。

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 標本回帰直線

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

無限回とはいきませんが、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)

3 正規方程式

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\) そのものだと想像できます。

3.2 行列表現

3.3 回帰モデル

4 回帰係数の変動

4.1 分散共分散行列

4.2 標本回帰直線の変動

＜例＞
「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.9942717 0.4965107

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

  • 回帰係数の分散・共分散
varB <- c(0.5^2/sum((X-mean(X))^2),var(B));varB
## [1] 0.02500000 0.02622608
varA <- c(0.5^2*mean(X^2)/sum((X-mean(X))^2),var(A));varA
## [1] 0.1500000 0.1455142
covAB <- c(0.5^2*mean(-X)/sum((X-mean(X))^2),cov(A,B));covAB
## [1] -0.05000000 -0.04966014

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

次に、理論分布に基づく標本回帰係数を、２次元正規乱数 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)) を通ります。（図では赤点、青点で表示）
これらは一致するとは限りません。

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 回帰統計量

5.1 残差

5.2 共分散

5.3 標本回帰係数

5.4 統計量の抽出

引き続き「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.37121143 -0.13549639  0.95754208 -0.02374928 -0.42708499

Coefficients（標本回帰係数）

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

Residual standard error（残差標準誤差）

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

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

sqrt(var(e)*4/3)
## [1] 0.647047

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.5012005 0.2046142

5.5 summary（概要）

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

summary(fit)
## 
## Call:
## lm(formula = Y ~ X, data = data)
## 
## Residuals:
##        1        2        3        4        5 
## -0.37121 -0.13550  0.95754 -0.02375 -0.42708 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   1.1681     0.5012   2.331    0.102
## X             0.2200     0.2046   1.075    0.361
## 
## Residual standard error: 0.647 on 3 degrees of freedom
## Multiple R-squared:  0.2781, Adjusted R-squared:  0.03746 
## F-statistic: 1.156 on 1 and 3 DF,  p-value: 0.3611

6 理論値との整合性

6.1 残差の無相関性

ここでは「残差の無相関性」に注目すると、

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

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

cor(e,X)#残差と説明変数の相関
## [1] -6.265327e-17
cor(e,f)#残差と予測値の相関
## [1] 1.068128e-16

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

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

cor(e,Y)#残差と観測値の相関係数
## [1] 0.8496507
cov(e,Y)#残差と観測値の共分散
## [1] 0.3140024

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

6.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.3140024 0.3140024

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

c(cov(f,Y),var(f))
## [1] 0.1209604 0.1209604

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

c(var(Y),var(f)+var(e))
## [1] 0.4349628 0.4349628

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

6.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] 1.739851 1.739851

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

6.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.2780936
var(f)/var(Y)
## [1] 0.2780936

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

7 有意性と適合度

7.1 回帰係数の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 
##    2.330611    1.075018

次に 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.1020909   0.3611308

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

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

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

7.2 決定係数

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

sr/st
## [1] 0.2780936

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

(cor(X,Y))^2
## [1] 0.2780936

7.3 重相関係数

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

\(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.9306195
cor(X,Y)#相関係数
## [1] -0.9306195

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.9878814
cor(X,Y)#相関係数
## [1] -0.9878814

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

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

R <- (cor(Y,f))^2;R#重相関係数の二乗
## [1] 0.9759097
R <- (var(f)*4)/(var(Y)*4);R#決定係数
## [1] 0.9759097

8 F統計量

8.1 偏差平方和の自由度

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

　　説明変数：　\(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.2383108
mean(SE)
## [1] 0.7673836
  • 回帰変動は母分散 0.25 の１倍なので、自由度は１
  • 残差変動は母分散 0.25 の３倍なので、自由度は３

と判断できます。

8.2 自由度の数理

8.3 回帰変動の期待値

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

　　説明変数：　\(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.737044
mean(SE)
## [1] 0.7540288
  • \(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倍で変化ありません。

8.4 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] 2.318399

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.543441

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

9 回帰の検出力

対立仮説「\(\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] 2.5822503 0.2170489
mean(SR)/mean(SE)
## [1] 11.89709

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

mean(SR/SE)
## [1] 36.15701

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

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

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

9.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.54

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

9.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] 4.9508999 0.2710212
c(mean(SR)/mean(SE),mean(SR/SE))
## [1] 18.26757 23.62334
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] 0.96

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

10 t値とF値の関係

10.1 決定係数との関係

10.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値は大きくなりますが、決定係数は大きくなるとは限らないようです。

10.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値と同じ形の変動となります。

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

11.1 2変量正規乱数

11.1.1 共分散と相関係数

11.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に関係性があることが分かります。

11.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.01144738  0.01816098

母回帰直線は、\(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.002035705 0.689183398

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

11.3 回帰モデル

11.4 標本回帰直線

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

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

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

ここで注目するのは、

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

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

11.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.2195162   0.5716489

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

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

11.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")