回帰分析や分散分析の手法は本当に妥当なのか？　どのくらいの有意差を検知できるのか、確率モデルに基づいた乱数からデータを生成し、妥当性と検出力を確認します。

1 帰無仮説のもとで

1.1 回帰分析

帰無仮説のもとで、回帰分散／残差分散 は F分布に従います。

回帰分散、残差分散はともに誤差の母分散の推定値であり、その自由度倍は独立にカイ二乗分布に従うので、帰無仮説のもとで 回帰分散／残差分散 はF分布に従うことが保障されています。（少し難しいですが、当サイト回帰分析の流れ などを参照）
母分散の最尤推定量については、 ウサギさんの統計サロン　【統計学】最尤推定・最尤推定量 を参照してください。

ここでは実際に帰無仮説（回帰係数０）のデータを作成してみましょう。
続いて、シミュレーションで1000回データを作成し、回帰分散と誤差分散の期待値が母分散になっていること、回帰分散と誤差分散が独立になっていることを確かめます。（独立の確認は難しいですが、相関係数が０に近いことを確かめます）

1.1.1 データの生成

x=0,1,2,3 とし、１つの x に対してデータを５個、誤差（標準偏差）0.5で出力します。

X <- c(rep(0,5),rep(1,5),rep(2,5),rep(3,5))
e <- rnorm(20,0,0.5)
Y <- 0*X+e
data <- data.frame(X=X,Y=Y)
#data
plot(data)

1.1.2 aov による回帰分析

aov は分散分析の関数です。R は因子型（factor）の説明変数に対し、特性関数（indicator function）を用いて回帰分析を行い、分散分析として出力しているそうです。このため、水準数が４ならば、使う特性関数は３つ、もう一つの水準は x=0 のときの切片の値として使うので、自由度は 4-1=3 、分散分析の場合の自由度が 水準数−１であるのは、こういう理由によります。
説明変数 x を数値型（numeric）として aov を適応すると、回帰分析が行われ、説明変数の自由度は１になります。
（factor(x)として因子型にすると、自由度は因子数−１になります）

res <- aov(Y ~ X, data = data)
summary(res)
##             Df Sum Sq Mean Sq F value Pr(>F)
## X            1 0.3431  0.3431    2.07  0.167
## Residuals   18 2.9832  0.1657

１変数回帰変動の自由度は１です。
aov の代わりに線形回帰の関数 lm （linear model）でも同じ結果が得られます。 summary を使えば必要な情報が整理して出力されますが、予測値、回帰係数、残差など個々の値を抽出し、統計量を自分で計算することも可能です。
予測値の分散に自由度19をかければ回帰平方和になり、回帰の自由度１で割れば回帰分散になります。残差平方和を残差の自由度18で割れば残差分散になります。

fit = lm( formula = Y ~ X, data = data )
f <- fit$ fitted.values;#予測値
a <- fit$ coefficients;#回帰係数
e <- fit$ residuals;#e#残差
se <- sum(e^2);se#残差平方和
## [1] 2.983236
s <- sqrt(se/18);s#残差標準偏差
## [1] 0.4071061
ff<- var(f)*19/(se/18);ff#F値
## [1] 2.070152
1-pf(ff, 1, 18)#p値
## [1] 0.1673689
summary(fit)
## 
## Call:
## lm(formula = Y ~ X, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.97372 -0.26231  0.05935  0.20141  0.62889 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.07431    0.15233   0.488    0.632
## X           -0.11715    0.08142  -1.439    0.167
## 
## Residual standard error: 0.4071 on 18 degrees of freedom
## Multiple R-squared:  0.1031, Adjusted R-squared:  0.05332 
## F-statistic:  2.07 on 1 and 18 DF,  p-value: 0.1674

1.2 誤差分散と回帰分散

残差分散と回帰分散の期待値は母分散になるはずですが、１回の試行では分かりません。試行を1000回行って1000回の平均を計算し、母分散に近い値になっているか確かめましょう。

S2=c();Vf=c();FF=c()
for (i in 1:1000){
  X <- c(rep(0,5),rep(1,5),rep(2,5),rep(3,5))
  e <- rnorm(20,0,0.5)
  Y <- 0*X+e
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ X, data = data )
  f <- fit$ fitted.values;#f#予測値
  e <- fit$ residuals;#e#残差
  se <- sum(e^2);se#残差平方和
  s2 <- se/18;#残差分散
  S2=c(S2,s2)
  vf <- var(f)*19#回帰分散
  Vf=c(Vf,vf)
  FF=c(FF,vf/s2)
}
hist(S2,prob=T,xlim=c(0,0.8),breaks = "scott")

mean(S2)
## [1] 0.2502946
hist(Vf,prob=T,xlim=c(0,0.8),breaks = "scott")

mean(Vf)
## [1] 0.2391788
cor(S2,Vf)
## [1] 0.04890166

誤差の母分散は \(0.5^2=0.25\) であり、回帰分散Vf、誤差分散S2とも、1000回の平均を求めると母分散 0.25に近い値になっています。相関係数も０に近く、回帰変動と残差変動は独立であると予想されます。回帰分散のヒストグラムが偏っているのは、回帰の自由度が１で小さいからと考えられます。

F値のヒストグラムと理論上の自由度(1,18)のF分布のグラフを描き、ほぼ一致していることを確かめます。

hist(FF,prob=T,xlim=c(0,10),breaks = "scott")
curve(df(x, 1, 18, ),add=T)

1.3 分散分析

同じモデルのデータを分散分析すると、結果はどうなるでしょうか。
説明変xを factor(x)として因子型で計算すると、aovもlm も分散分析を行ってくれ、説明変数の自由度は（水準数）−１＝３ になります。

X <- c(rep(0,5),rep(1,5),rep(2,5),rep(3,5))
e <- rnorm(20,0,0.5)
Y <- 0*X+e
data <- data.frame(X=X,Y=Y)
#data
plot(data)

res <- aov(Y ~ factor(X), data = data)
summary(res)
##             Df Sum Sq Mean Sq F value Pr(>F)
## factor(X)    3  0.801  0.2668   1.322  0.302
## Residuals   16  3.231  0.2019
fit = lm( formula = Y ~ factor(X), data = data )
f <- fit$ fitted.values;#予測値
a <- fit$ coefficients;#回帰係数
e <- fit$ residuals;#e#残差
se <- sum(e^2);se#残差平方和
## [1] 3.230722
s <- sqrt(se/16);s#残差標準偏差
## [1] 0.4493552
vs <- sum(e^2)/16;vs#残差分散
## [1] 0.2019201
vf <- var(f)*19/3;vf#回帰分散
## [1] 0.2668434
ff <- vf/vs;ff#F値
## [1] 1.321529
summary(fit)
## 
## Call:
## lm(formula = Y ~ factor(X), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.59480 -0.29010 -0.06436  0.33595  0.67688 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.1521     0.2010   0.757   0.4600  
## factor(X)1   -0.5298     0.2842  -1.864   0.0807 .
## factor(X)2   -0.1469     0.2842  -0.517   0.6124  
## factor(X)3   -0.3436     0.2842  -1.209   0.2442  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4494 on 16 degrees of freedom
## Multiple R-squared:  0.1986, Adjusted R-squared:  0.04831 
## F-statistic: 1.322 on 3 and 16 DF,  p-value: 0.3021

分散分析では、群内平均が予測値であること、説明変数factorをダミー変数として扱うので自由度が３（x=0 を除くので 4-1=3）となっていることが分かります。

1.4 誤差分散と分散分析

S2=c();Vf=c();FF=c()
for (i in 1:1000){
  X <- c(rep(0,5),rep(1,5),rep(2,5),rep(3,5))
  e <- rnorm(20,0,0.5)
  Y <- 0*X+e
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ factor(X), data = data )
  f <- fit$ fitted.values;f#予測値
  e <- fit$ residuals;e#残差
  se <- sum(e^2);se#残差平方和
  s2 <- se/16;#残差分散
  S2=c(S2,s2)
  vf <- var(f)*19/3#回帰分散
  Vf=c(Vf,vf)
  FF=c(FF,vf/s2)
}
hist(S2,prob=T,xlim=c(0,0.8),breaks = "scott")

mean(S2)
## [1] 0.2543187
hist(Vf,prob=T,xlim=c(0,0.8),breaks = "scott")

mean(Vf)
## [1] 0.2482396
cor(S2,Vf)
## [1] -0.04502646

誤差の母分散は \(0.5^2=0.25\) であり、回帰分散、誤差分散ともに近い値になっています。相関係数も０に近く、独立であると予想されます。回帰分散のヒストグラムが偏っているのは、回帰の自由度が１であるからと考えられます。

F値のヒストグラムと理論上の自由度(1,18)のF分布のグラフを描き、ほぼ一致していることを確かめます。

hist(FF,prob=T,xlim=c(0,10),breaks = "scott")
curve(df(x, 3, 16 ),add=T)

2 対立仮説のもとで

回帰係数が０でないモデルの場合、回帰分散や誤差分散はどのようになるでしょうか。

2.1 回帰分析

S2=c();Vf=c();FF=c()
for (i in 1:1000){
  X <- c(rep(0,5),rep(1,5),rep(2,5),rep(3,5))
  e <- rnorm(20,0,0.5)
  Y <- 0.3*X+e
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ X, data = data )
  f <- fit$ fitted.values;#f#予測値
  e <- fit$ residuals;#e#残差
  se <- sum(e^2);se#残差平方和
  s2 <- se/18;#残差分散
  S2=c(S2,s2)
  vf <- var(f)*19#回帰分散
  Vf=c(Vf,vf)
  FF=c(FF,vf/s2)
}
hist(S2,prob=T,xlim=c(0,0.8),breaks = "scott")

mean(S2)
## [1] 0.2516408
hist(Vf,prob=T,breaks = "scott")

mean(Vf)
## [1] 2.580031
cor(S2,Vf)
## [1] -0.05164554
mean(FF)
## [1] 11.52313
1-pf(mean(FF), 1, 18)
## [1] 0.003230094
hist(FF,prob=T,breaks = "scott")
curve(df(x, 1, 18, ),add=T)

誤差分散は変わりませんが、回帰分散が大きくなり、F値が大きくなることで異常が検知されます。

2.2 分散分析

S2=c();Vf=c();FF=c()
for (i in 1:1000){
  X <- c(rep(0,5),rep(1,5),rep(2,5),rep(3,5))
  e <- rnorm(20,0,0.5)
  Y <- 0.3*X+e
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ factor(X), data = data )
  f <- fit$ fitted.values;f#予測値
  e <- fit$ residuals;e#残差
  se <- sum(e^2);se#残差平方和
  s2 <- se/16;#残差分散
  S2=c(S2,s2)
  vf <- var(f)*19/3#回帰分散
  Vf=c(Vf,vf)
  FF=c(FF,vf/s2)
}
hist(S2,prob=T,xlim=c(0,0.8),breaks = "scott")

mean(S2)
## [1] 0.2465759
hist(Vf,prob=T,breaks = "scott")

mean(Vf)
## [1] 1.005865
cor(S2,Vf)
## [1] -0.01898453
mean(FF)
## [1] 4.6473
1-pf(mean(FF), 3, 16)
## [1] 0.01606317
hist(FF,prob=T,breaks = "scott")
curve(df(x, 3, 16 ),add=T)

やはり誤差分散は変わりませんが、回帰分散が大きくなり、F値が大きくなることで異常が検知されます。