回帰係数の検定
参考サイト.
矢内勇生「回帰分析のシミュレーション」
津田裕之「Rによる統計入門 8.回帰分析」
1 回帰の有意性の検定.
単回帰式 y=a+bx において、b=0 を検定します。
- x に対して y は有意に変化しているか.
- それとも変化は誤差の範囲内なのか.
判断するのが有意性の検定です。
dat<- data.frame(X=c(1,3,4,5,6,7,9,10,11,12),Y=c(4,4,5,6,5,7,6,8,9,7))
mean(dat$Y)## [1] 6.1回帰分析では、x は固定された定数です。
y の平均は 6.1 なので、帰無仮説は 「b=0 , a=6.1」 です。
＜帰無仮説＞ b=0 , a=mean(Y).
回帰モデル \(y=a+bx+\varepsilon\) を作るには、誤差分散 \(V(\varepsilon)=\sigma^2\) が必要ですが、残差分散で代用します。
fit = lm( formula = Y ~ X, data = dat )
summary(fit)##
## Call:
## lm(formula = Y ~ X, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.15217 -0.73829 0.09699 0.63043 1.24247
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.41639 0.61769 5.531 0.000553 ***
## X 0.39465 0.08097 4.874 0.001233 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8855 on 8 degrees of freedom
## Multiple R-squared: 0.7481, Adjusted R-squared: 0.7166
## F-statistic: 23.76 on 1 and 8 DF, p-value: 0.001233Residual standard error: 0.8855
回帰残差 s=0.8855 としてモデルを構成します。
＜帰無仮説＞ 「回帰変動はない（b=0）」 のもとで .
xによる変動がなく誤差だけの変動ならば、誤差変動 \(\varepsilon_i \sim
N(0,~0.8855^2)\)として
\(y_i=6.1+0\cdot
x_i+\varepsilon_i\) という単回帰モデルになります。
1.1 １回シミュレーション.
X=c(1,3,4,5,6,7,9,10,11,12)
e <- rnorm(10,0,0.8855)
Y <- 6.1+0*X+e
data <- data.frame(X=X,Y=Y);data| X | Y |
|---|---|
| 1 | 5.000929 |
| 3 | 6.608828 |
| 4 | 6.345525 |
| 5 | 6.399468 |
| 6 | 6.439493 |
| 7 | 6.700170 |
| 9 | 5.986318 |
| 10 | 5.666699 |
| 11 | 6.105494 |
| 12 | 5.898244 |
fit = lm( formula = Y ~ X, data = data )
b <- fit$ coefficients[2];b#回帰係数## X
## 0.003346726a <- fit$ coefficients[1];a#回帰係数## (Intercept)
## 6.092359plot(data)
abline(fit)
回帰係数 b=0 のモデルから生成しても、回帰分析で出力される b は0とは限りません。なので、回帰が有意かどうかの検定が必要になるわけです。
1.2 千回シミュレーション.
X=c(1,3,4,5,6,7,9,10,11,12)
B=c()
A=c()
for (i in 1:1000) {
e <- rnorm(10,0,0.8855)
Y <- 6.1+0*X+e
data <- data.frame(X=X,Y=Y)
fit = lm( formula = Y ~ X, data = data )
b <- fit$ coefficients[2]#回帰係数
a <- fit$ coefficients[1]#回帰係数
B=c(B,b)
A=c(A,a)
}
hist(B,prob=T,breaks="sturges",xlim=c(-0.4,0.4))
abline(v=0.39,col="red")
hist(A,prob=T,breaks="sturges",xlim=c(3,9))
abline(v=3.41,col="red")
1.3 回帰係数の標準偏差とt値.
帰無仮説 b=0 のもとでシミュレーションを繰り返したとき、Rが出力するbの期待値は0と考えられますが、bの標準偏差はどのくらいでしょうか？
sd(B)#標準偏差## [1] 0.08140881\(t=\dfrac{b-0}{sd(B)}\sim t(n-2)\) が有意差を超えれば b=0 は棄却されることになります。
0.39/sd(B)#t値## [1] 4.790636p208で求めた b=0.39 に対する t値
はヒストグラムの右端に位置しており有意です。
summaryを見ると、
- 回帰係数bの標準偏差 Std. Error=0.08097.
- t値 t value=4.874.
1000回のシミュレーションの平均にほぼ一致しています.
a=3.4 の有意性についてはどうでしょうか？
sd(A)#標準偏差## [1] 0.640259(3.4-6.1)/sd(A)#t値## [1] -4.217044(3.4-0)/sd(A)#t値## [1] 5.310351summaryは a=0 を帰無仮説として 3.4 のt値を求めているようです。
ここではシミュレーションによってsd(B)を求めましたが、理論による推定値は.
\(sd=\sqrt{\frac{\sigma^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2}}\) です.
X=c(1,3,4,5,6,7,9,10,11,12)
sqrt(0.8855^2/(var(X)*9))## [1] 0.08096978一回のデータから推測したbの標準偏差ですが、シミュレーションから求めた値と一致しています。
2 回帰係数の差の検定.
２変量 x,y
のデータが与えられたとき、想定される母集団は2変量正規分布であり楕円状の分布を持ちます。（前節「回帰と相関」参照）.
回帰分析は2変量正規分布から調査によって抽出したデータをもとに
xからyを予測する式を求めようとするものです。ここで、2変量正規分布と回帰分析との関係を整理しておきましょう。
- 2変量正規分布でxを固定した条件付き期待値は yのxによる回帰直線に等しい.
- yの条件付き期待値は縦の弦の中点.
- 回帰直線は縦の弦の中点の軌跡.
（回帰直線は2変量正規分布で図形的な意味を持ちます）.
\(\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)\sim N\left(\mu_y+\rho \sigma_y\left(\dfrac{x-\mu_x}{\sigma_x}\right),(1-\rho^2)\sigma_y^2\right)\).
回帰の残差分散と2変量正規分布の分散との関係.
- 回帰の残差分散は、\(s^2=(1-\rho^2)\sigma_y^2\).
残差分散はxに依存しません。
具体例で見ていきましょう。
3 2変量正規乱数.
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="grey")
abline(model)
abline(v=1.5,col="blue")
3.1 回帰モデル.
https://selfboast.hatenablog.jp/entry/2022/08/18/114604
回帰の残差分散は、\(s^2=(1-\rho^2)\sigma_y^2=0.75\) なので.
回帰モデルは \(y=0.5x+\varepsilon~~,~~\varepsilon\sim N(0,0.75)\).
回帰の Residual standard error =\(\sqrt{0.75}=0.866\).
X <-rep(seq(-3,3,by= 1),10)
e <- rnorm(70,0,0.866)
Y <- 0.5*X+e
dat <- data.frame(X=X,Y=Y)
plot(x,asp=1,col="grey",xlim=c(-4,4),ylim=c(-4,4))
abline(model)
abline(v=1.5,col="blue")
par(new=T)
plot(dat,asp=1,col="blue",xlim=c(-4,4),ylim=c(-4,4))
点の配置は回帰直線に平行な帯の内側になり、Yの推定値の信頼区間をイメージすることができます。
3.2 標本抽出.
実際の観測で得られた2変量データの回帰係数と母集団である2変量正規分布から得られる母回帰直線とはどれくらい違うのでしょうか。シミュレーションで確かめます。
S1 <- sample_n(tbl = data, size = 100)
head(S1)| X | Y |
|---|---|
| 0.1654005 | 1.0773641 |
| 0.2180630 | -0.8015362 |
| -0.2217801 | -0.3830044 |
| -0.5026906 | -1.6139025 |
| -0.1289488 | 0.8854419 |
| 1.3133249 | -0.2228231 |
model1 <- lm(Y ~ X, data = S1)
summary(model1)##
## Call:
## lm(formula = Y ~ X, data = S1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.09100 -0.60664 -0.01742 0.48548 2.12630
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.09804 0.08613 -1.138 0.258
## X 0.54221 0.09428 5.751 1.01e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8596 on 98 degrees of freedom
## Multiple R-squared: 0.2523, Adjusted R-squared: 0.2447
## F-statistic: 33.08 on 1 and 98 DF, p-value: 1.008e-07cov(S1)## X Y
## X 0.8397549 0.4553221
## Y 0.4553221 0.9783348 \(\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}\).
\(V(\beta)=\dfrac{\sigma^2}{T_{xx}}=\dfrac{\sigma^2}{ns_{x}^2}\)
cov(S1)[1,2]/cov(S1)[1,1]#回帰係数## [1] 0.5422083sqrt((sum((model1$residuals)^2)/98)/(var(S1$X)*99))#標準誤差## [1] 0.09427675回帰係数とその標準誤差がsummaryの値と一致していることがわかります。
3.3 信頼区間.
母集団の回帰直線の信頼区間は回帰直線に平行な帯ですが、標本から推定した回帰直線では係数に誤差が含まれるので信頼区間は左右にいくと広がります。
plot(x,asp=1,col="grey",xlim=c(-4,4),ylim=c(-4,4))
abline(model)
for (i in 1:100){
S <- sample_n(tbl = data, size = 100)
model <- lm(Y ~ X, data = S)
par(new=T)
abline(model,col="lightblue")
}
信頼区間は.
\(y_i\sim N\left(\alpha+\beta x_i\:,\:\dfrac{\sigma^2}{n}+\dfrac{\sigma^2}{T_{xx}}(x_i-\overline{x})^2\right)\) から得られる。.
式から.
- 標本数が大きければ予測の幅が小さくなること.
- xの中心から離れると予測の幅が小さくなること.
が分かります。
3.5 同一母集団からの標本
S1 <- sample_n(tbl = data, size = 100)
model1 <- lm(Y ~ X, data = S1)
S2 <- sample_n(tbl = data, size = 100)
model2 <- lm(Y ~ X, data = S2)
plot(x,asp=1,col="grey",xlim=c(-4,4),ylim=c(-4,4),xlab="",ylab="")
abline(model)
par(new=T)
plot(S1,asp=1,col="blue",xlim=c(-4,4),ylim=c(-4,4))
abline(model1,col="blue")
par(new=T)
plot(S2,asp=1,col="red",xlim=c(-4,4),ylim=c(-4,4))
abline(model2,col="red")
model1$coef## (Intercept) X
## 0.03102363 0.43485116model2$coef## (Intercept) X
## 0.1491536 0.5452817異なる標本であっても、母集団が同じならば回帰直線は前項のlightblueの帯に収まります。
3.6 異なる母集団からの標本
回帰係数は正規分布に従う。
回帰係数に差があるかどうかは、２種類の正規集団に差があるかを検定する。
例7_3(2) については、２つの標本の回帰係数について２標本t検定を行う。
数理的計算はテキストに詳しく説明されている。
ここではRを用い、標本の所属を表す変数 y とTBGの値 x
を説明変数とし、甲状腺ホルモンT4の値 z
を目的変数とする重回帰分析を行う。
標本によって回帰係数に差があるかどうかは、交互作用の検定に出力される。
参考サイト 初学者の統計学習ノート様「回帰係数の差の検定」.
3.7 例7_3_2(p220)
x <- c(54,50,39,44,35,41,32,46,38,30,43,39,49,42,35,31,36,34,27)
y <- c(rep(1,10),rep(2,9))
z <- c(17,15,15,14,13,12,11,10,10,9,22,21,19,18,16,15,13,11,9)
data <- data.frame(x=x,y=y,z=z)S <- data[,-2]
model <- lm(z ~ x, data = S)
S1 <- data[1:10,-2]
model1 <- lm(z ~ x, data = S1)
S2 <- data[11:19,-2]
model2 <- lm(z ~ x, data = S2)
plot(S,asp=1,col="grey",xlim=c(20,60),ylim=c(5,25),xlab="",ylab="")
abline(model)
par(new=T)
plot(S1,asp=1,col="blue",xlim=c(20,60),ylim=c(5,25))
abline(model1,col="blue")
par(new=T)
plot(S2,asp=1,col="red",xlim=c(20,60),ylim=c(5,25))
abline(model2,col="red")
黒が全体での回帰直線、赤と青は標本の回帰直線。
3.7.1 重回帰分析.
res <- lm(z ~ x*y, data = data)
summary(res)##
## Call:
## lm(formula = z ~ x * y, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7880 -1.7779 0.2802 1.5614 4.1341
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.54210 10.33563 0.923 0.371
## x -0.05367 0.25470 -0.211 0.836
## y -6.46937 6.75830 -0.957 0.354
## x:y 0.28661 0.17161 1.670 0.116
##
## Residual standard error: 2.513 on 15 degrees of freedom
## Multiple R-squared: 0.6557, Adjusted R-squared: 0.5868
## F-statistic: 9.522 on 3 and 15 DF, p-value: 0.0009101交互作用 x:y は、 t=1.670 となっており、p222のt値と一致している。.