回帰分析において不均一分散に頑健な標準誤差を用いるべきか？

要約

　回帰分析では、不均一分散に頑健な標準誤差を用いるべきです。

　なぜなら、誤差項が不均一分散でも、不均一分散に頑健な標準誤差（赤）は、標本回帰係数の標準偏差(黒)を一致性を持って推定できます。しかし、通常の標準誤差（青）では一致性を持って推定できません。これは下図の例で確認できます。

全体像

（1）問題の構造

　「単回帰分析にて不均一分散に頑健な標準誤差を用いるべきか？」という問いを設定します。この問いを

①標準誤差とは

②均一分散の標準誤差とは

③不均一分散に頑健な標準誤差とは

④「用いるべき/べきでない」を判断する基準とは

⑤「均一分散の標準誤差」と「不均一分散に頑健な標準誤差」のどちらが④の判断基準に合致するか

に分解します。

（2）前提の選択

　①〜④については前提に含めてしまいます。

（3）論点の選択

　⑤をもう少し簡単にして「具体的な4つの場合でも、不均一分散に頑健な標準誤差が望ましさの基準を満たすか」を論点とします。

（4）付録一覧

　冗長さを避けて、可読性を上げるために、以下の内容は記事最後の付録に回します。

「不均一分散に頑健な標準誤差のRでの計算方法」「シミュレーションのRコード」

前提

①標準誤差とは

　標準誤差（standard error）とは「標本回帰係数の標準偏差」の推定値です。

$$SE=\widehat{ SD \left(\widehat{\beta} \right) } $$

$$標準誤差SE、標準偏差SD、推定値\widehat{　}$$

$$標本回帰係数\widehat{\beta}、母回帰係数\beta$$

②均一分散の標準誤差とは

　均一分散の標準誤差とは、Rやエクセルの初期設定で計算される標準誤差です。誤差項が均一に分散していると仮定しています。具体的な数式は、以下のとおりです。

$$SE_{均一分散}=\sqrt {\frac{ \widehat{\sigma^2 } }{ \sum\limits_{i=1}^n (x_i-\overline{x})^2 } }$$

　この記事では、Rにて次のコードで計算できる「Std. Error」を不均一分散に頑健な標準誤差と定義します。

lm(y ~ x) #単回帰分析。均一分散を仮定した標準誤差

③不均一分散に頑健な標準誤差とは

　不均一分散に頑健な標準誤差とは、誤差項が不均一に分散していても、標本回帰係数の標準偏差を一致性をもって計算できる標準誤差です。

$$SE_{不均一分散}= \sqrt {\frac{ \sum\limits_{i=1}^n (x_i-\overline{x})^2 \widehat{\sigma_i^2} }{ \left(\sum\limits_{i=1}^n (x_i-\overline{x})^2 \right) ^2} }$$

$$標準誤差SE,データ番号i、サンプル・サイズn$$

$$xの平均\overline{x}、誤差項の分散の推定値\widehat{\sigma_i^2}$$

　この記事では、Rにて次のコードで計算できる「Std. Error」を不均一分散に頑健な標準誤差と定義します。

library(estimatr) #パッケージestimatrの呼び出し
lm_robust(y ~ x, se_type = "stata") #単回帰分析。stataと同じ不均一分散に頑健な標準誤差を出力

④「用いるべき/べきでない」を判断する基準とは

・一致性をもつ　→　用いるべき

・一致性がない　→　用いないべき

という基準を設定します。この記事では

・大きなサンプルサイズ1000において

・得られる可能性のあった推定値が

・パラメーターの近傍に集まっている

ならば、一致性があるとみなします。なお、本来の一致性の定義とは異なるので注意してください（後述）。

方法

　論点は「『均一分散の標準誤差』と『不均一分散に頑健な標準誤差』のどちらが一致性をもつか」です。

　そこで

・大きなサンプルサイズ1000を想定し

・乱数の発生によって得られる可能性のあった標本回帰係数βと標準誤差SEを2000個作り出し

・「標本回帰係数βの標準偏差SD」の近傍に集まっているか否か

をヒストグラムをもちいてチェックします。

結果

均一分散の場合１

　次の単回帰モデルに従うサンプル・サイズ1000の標本を生成します。

$$単回帰モデル：Y=2+X+U$$

$$誤差項U　〜　N(0,2)$$

つまり、誤差項Uの分散は均一分散

　

　その結果、生成された標本の一つが↓です。

　この標本を2000回作った結果、標本回帰係数の標準偏差SDは0.0623でした。

　標準誤差の分布は、下図でした。青が均一分散の標準誤差で、赤が不均一分散に頑健な標準誤差（stata,HC1）です。

　均一分散の標準誤差も不均一分散に頑健な標準誤差の両方とも標本回帰係数の標準偏差の近傍に集まっており、一致性をもちそうです。

場合１	数値
標本回帰係数の標準偏差SD	0.0623
均一分散の標準誤差SEの平均	0.0624
不均一分散に頑健な標準誤差SEの平均	0.0623

不均一分散の場合２

　次の単回帰モデルに従うサンプル・サイズ1000の標本を生成します。

$$単回帰モデル：Y=2+X+U$$

$$誤差項U　〜　N(0,0.3X)$$

つまり、誤差項Uの分散は、xが大きいと誤差項の分散が大きくなる

　

　その結果、生成された標本の一つが↓です。

　この標本を2000回作った結果、標本回帰係数の標準偏差SDは0.0503でした。

　このときの標準誤差の分布は、下図でした。青が均一分散の標準誤差で、赤が不均一分散の標準誤差（stata, HC1）です。

　不均一分散に頑健な標準誤差は標本回帰係数の標準偏差の近傍に集まっており、一致性をもちそうです。一方で、均一分散の標準誤差は標本回帰係数の標準偏差の近傍に集まっているとは言えず、一致性をもたなそうです。

場合２	数値
標本回帰係数の標準偏差SD	0.0503
均一分散の標準誤差SEの平均	0.0484
不均一分散に頑健な標準誤差SEの平均	0.0501

不均一分散の場合３

　次の単回帰モデルに従うサンプル・サイズ1000の標本を生成します。

$$単回帰モデル：Y=2+X+U$$

$$誤差項U　〜　N(0,26-10X+X^2)$$

つまり、xが平均に近いと誤差項の分散が小さくなる

　

　その結果、生成された標本の一つが↓です。

　この標本を2000回作った結果、標本回帰係数の標準偏差SDは0.146でした。

　このときの標準誤差の分布は、下図でした。青が均一分散の標準誤差で、赤が不均一分散の標準誤差（stata, HC1）です。

　不均一分散に頑健な標準誤差は標本回帰係数の標準偏差の近傍に集まっており、一致性をもちそうです。一方で、均一分散の標準誤差は標本回帰係数の標準偏差の近傍に集まっているとは言えず、一致性をもたなそうです。

場合３	数値
標本回帰係数の標準偏差SD	0.146
均一分散の標準誤差SEの平均	0.0771
不均一分散に頑健な標準誤差SEの平均	0.145

不均一分散の場合４

　次の単回帰モデルに従うサンプル・サイズ1000の標本を生成します。

$$単回帰モデル：Y=2+X+U$$

$$誤差項U　〜　N(0,-11+10X-X^2)$$

$$ただし、-11+10X-X^2≦0ならば、標準偏差＝0.01$$

つまり、xが平均に近いと誤差項の分散が大きくなる

　

　その結果、生成された標本の一つが↓です。

　この標本を2000回作った結果、標本回帰係数の標準偏差SDは0.359でした。

　このときの標準誤差の分布は、下図でした。青が均一分散の標準誤差で、赤が不均一分散の標準誤差（stata, HC1）です。

　不均一分散に頑健な標準誤差は標本回帰係数の標準偏差の近傍に集まっており、一致性をもちそうです。一方で、均一分散の標準誤差は標本回帰係数の標準偏差の近傍に集まっているとは言えず、一致性をもたなそうです。

場合４	数値
標本回帰係数の標準偏差SD	0.359
均一分散の標準誤差SEの平均	0.414
不均一分散に頑健な標準誤差SEの平均	0.357

まとめ

　場合１〜４のすべてで、不均一分散に頑健な標準誤差は、一致性を持ちました。

　均一分散の場合１では、均一分散の標準誤差は一致性を持ちました。不均一分散の場合２〜４では、均一分散の標準誤差は一致性を持ちませんでした。

考察

（1）結論

　単回帰分析において、不均一分散に頑健な標準誤差を用いるべきです。なぜなら、不均一分散でも一致性をもつからです。

　データ分析前に誤差項の分布は未知ですので、誤差項が均一分散との仮定をおくのは不誠実です。

（2）妥当性評価

前提評価

　望ましさの指標を一致性と明確にしたのは、Goodです。

　一致性の定義を歪めたのがBadです。一致性の本来の定義は

$$標本の数nが無限大のとき$$

$$推定量\widehat{\theta}_nがパラメーター\thetaに確率収束する$$

です。これを

$$標本の数nが1000と大きいとき$$

$$推定量\widehat{\theta}_nがパラメーター\thetaの近傍に集まっている$$

と言い換えてしまいました。

方法評価

　視覚的にわかりやすく判断したのは、Goodです。

　具体例を4つあげて検証したにすぎず、一般性にかけるのがBadです。

結論評価

　不均一分散に頑健な標準誤差の名前の通り、不均一分散にも対応している標準誤差なので、こちらを使うべきという結論は正しいでしょう。なお、不均一分散に頑健な標準誤差の英訳は

・Heteroskedasticity-consistent standard errors

であり、直訳すると

・不均一（Hetero）な分散（skedasticity）で一致性（consistent）をもつ標準誤差（standard errors）

となります。HC標準誤差とも言えます。略称としてはロバスト標準誤差が多いです。

（3）意義

　回帰分析では不均一分散に頑健な標準誤差を使うべきということがわかりました。

付録：不均一分散に頑健な標準誤差のRでの計算方法

　不均一分散に頑健な標準誤差（Std. Error）は、Rでは次のlm_robust関数と、summary（）で得られます。↓

#パッケージの呼び出し。未インストールの場合はインストール！
library(estimatr)
　
#不均一分散に頑健な標準誤差を用いた回帰分析（4種）
lm_robust(y ~ x, se_type = "HC0")
lm_robust(y ~ x, se_type = "stata") #HC1でも可。stataと同じ
lm_robust(y ~ x, se_type = "HC2")
lm_robust(y ~ x, se_type = "HC3")
　
#均一分散が仮定された標準誤差を用いた回帰分析
lm_robust(y ~ x, se_type = "classical") #lmと同じこともできます

付録：Rコード

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

　以下のシミュレーションすべてにおいて、処理に数分の時間がかかります。少々お待ちください。

　使うパッケージは↓

library(estimatr)　#lm_robustを用いるため
library(ggplot2)  #ggplotを用いるため

①均一分散のコード↓

#均一分散
trial_number <- 2000
sample_size <- 1000
beta_1 <- rep(NA, trial_number)
　
se_classical <-rep(NA, trial_number)
se_HC0 <- rep(NA, trial_number)
se_stata <- rep(NA, trial_number)
se_HC2 <- rep(NA, trial_number)
se_HC3 <- rep(NA, trial_number)
se_HC4 <- rep(NA, trial_number)
　
u <- rep(NA, sample_size)
　　
#標本生成、推定値の計算、標準誤差の計算
for(j in 1:trial_number){
  x <- rnorm(n=sample_size, mean =5, sd=1)
　
  for(i in 1:sample_size){
    uu <- rnorm(n=1, mean =0, sd=2)  #均一分散
    u[i] <- uu
  }
  　
  y <- 2 + 1*x+ u
  answer <- lm(y~x)
  answer0 <- lm_robust(y ~ x,se_type = "HC0")
  answer1 <- lm_robust(y ~ x,se_type = "stata")
  answer2 <- lm_robust(y ~ x,se_type = "HC2")
  answer3 <- lm_robust(y ~ x,se_type = "HC3")
  　
  beta_1[j] <- summary(answer)$coef["x","Estimate"]
  　
  se_classical[j] <- summary(answer)$coef["x","Std. Error"]
  se_HC0[j] <- summary(answer0)$coef["x","Std. Error"]
  se_stata[j] <- summary(answer1)$coef["x","Std. Error"]
  se_HC2[j] <- summary(answer2)$coef["x","Std. Error"]
  se_HC3[j] <- summary(answer3)$coef["x","Std. Error"]
}
　
#データフレーム作成
data0 <- data.frame(x,y,u) 
beta <- data.frame(beta_1)
data_se <- data.frame(se_classical,se_HC0,se_stata,se_HC2,se_HC3) 
　
#(x,y)の散布図
ggplot(data0, aes(x = x,y= y)) +　
  geom_point() +
  scale_x_continuous(limits = c(0, 10),breaks = seq(0,10,2))+
  scale_y_continuous(limits = c(0, 15),breaks = seq(0,15,5))+
  stat_smooth(method = "lm", formula='y~x', se = FALSE)
　
#推定量の分布
ggplot(beta, aes(x = beta_1)) +
  geom_histogram(colour="white", fill = "skyblue2",bins=50)+
  scale_x_continuous(limits = c(0.6, 1.4))
　
#標準誤差の分布
ggplot(data_se) +
  geom_histogram(aes(x = se_classical),colour="white", fill = "skyblue2",bins=50)+
  geom_histogram(aes(x = se_stata),colour="white", fill = "red",alpha=0.2,bins=50)+
  scale_x_continuous(limits = c(0.055, 0.07))+
  geom_vline(xintercept =sd(beta_1) ) 
　
#推定量の標準偏差を計算
sd(beta_1) 
　
#標準誤差の平均を計算
mean(se_classical)
mean(se_HC0)
mean(se_stata)
mean(se_HC2)
mean(se_HC3)

②不均一分散のコード↓

#①不均一分散
trial_number <- 2000
sample_size <- 1000
beta_1 <- rep(NA, trial_number)
se_classical <-rep(NA, trial_number)
se_HC0 <- rep(NA, trial_number)
se_stata <- rep(NA, trial_number)
se_HC2 <- rep(NA, trial_number)
se_HC3 <- rep(NA, trial_number)
se_HC4 <- rep(NA, trial_number)
u <- rep(NA, sample_size)
#標本生成、推定値の計算、標準誤差の計算
for(j in 1:trial_number){
  x <- rnorm(n=sample_size, mean =5, sd=1)
  for(i in 1:sample_size){
    uu <- rnorm(n=1, mean =0, sd=0.3*x[i])  #不均一分散
    u[i] <- uu
  }
  
  y <- 2 + 1*x+ u
  answer <- lm(y~x)
  answer0 <- lm_robust(y ~ x,se_type = "HC0")
  answer1 <- lm_robust(y ~ x,se_type = "stata")
  answer2 <- lm_robust(y ~ x,se_type = "HC2")
  answer3 <- lm_robust(y ~ x,se_type = "HC3")
  
  beta_1[j] <- summary(answer)$coef["x","Estimate"]
  
  se_classical[j] <- summary(answer)$coef["x","Std. Error"]
  se_HC0[j] <- summary(answer0)$coef["x","Std. Error"]
  se_stata[j] <- summary(answer1)$coef["x","Std. Error"]
  se_HC2[j] <- summary(answer2)$coef["x","Std. Error"]
  se_HC3[j] <- summary(answer3)$coef["x","Std. Error"]
}
#データフレーム作成
data0 <- data.frame(x,y,u) 
beta <- data.frame(beta_1)
data_se <- data.frame(se_classical,se_HC0,se_stata,se_HC2,se_HC3) 
#(x,y)の散布図
ggplot(data0, aes(x = x,y= y)) +　
  geom_point() +
  scale_x_continuous(limits = c(0, 10),breaks = seq(0,10,2))+
  scale_y_continuous(limits = c(0, 15),breaks = seq(0,15,5))+
  stat_smooth(method = "lm", formula='y~x', se = FALSE)
#推定量の分布
ggplot(beta, aes(x = beta_1)) +
  geom_histogram(colour="white", fill = "skyblue2",bins=50)+
  scale_x_continuous(limits = c(0.6, 1.4))
#標準誤差の分布
ggplot(data_se) +
  geom_histogram(aes(x = se_classical),colour="white", fill = "skyblue2",bins=50)+
  geom_histogram(aes(x = se_stata),colour="white", fill = "red",alpha=0.2,bins=50)+
  scale_x_continuous(limits = c(0.04, 0.06))+
  geom_vline(xintercept =sd(beta_1) )+ #推定量の標準偏差
  geom_vline(xintercept =mean(se_classical),col="blue")+ #均一分散の標準誤差平均
  geom_vline(xintercept =mean(se_stata),col="red") #不均一分散の標準誤差平均
#推定量の標準偏差を計算
sd(beta_1) 
#標準誤差の平均を計算
mean(se_classical)
mean(se_HC0)
mean(se_stata)
mean(se_HC2)
mean(se_HC3)

③不均一分散のコード↓

#不均一分散②
trial_number <- 2000
sample_size <- 1000
beta_1 <- rep(NA, trial_number)
se_classical <-rep(NA, trial_number)
se_HC0 <- rep(NA, trial_number)
se_stata <- rep(NA, trial_number)
se_HC2 <- rep(NA, trial_number)
se_HC3 <- rep(NA, trial_number)
se_HC4 <- rep(NA, trial_number)
u <- rep(NA, sample_size)
min(23-10*x+x^2)
#標本生成、推定値の計算、標準誤差の計算
for(j in 1:trial_number){
  x <- rnorm(n=sample_size, mean =5, sd=1)
  for(i in 1:sample_size){
    uu <- rnorm(n=1, mean =0, sd=26-10*x[i]+x[i]^2)  #不均一分散
    u[i] <- uu
  }
  
  y <- 2 + 1*x+ u
  answer <- lm(y~x)
  answer0 <- lm_robust(y ~ x,se_type = "HC0")
  answer1 <- lm_robust(y ~ x,se_type = "stata")
  answer2 <- lm_robust(y ~ x,se_type = "HC2")
  answer3 <- lm_robust(y ~ x,se_type = "HC3")
  
  beta_1[j] <- summary(answer)$coef["x","Estimate"]
  
  se_classical[j] <- summary(answer)$coef["x","Std. Error"]
  se_HC0[j] <- summary(answer0)$coef["x","Std. Error"]
  se_stata[j] <- summary(answer1)$coef["x","Std. Error"]
  se_HC2[j] <- summary(answer2)$coef["x","Std. Error"]
  se_HC3[j] <- summary(answer3)$coef["x","Std. Error"]
}
#データフレーム作成
data0 <- data.frame(x,y,u) 
beta <- data.frame(beta_1)
data_se <- data.frame(se_classical,se_HC0,se_stata,se_HC2,se_HC3) 
#(x,y)の散布図
ggplot(data0, aes(x = x,y= y)) +　
  geom_point() +
  scale_x_continuous(limits = c(0, 10),breaks = seq(0,10,2))+
  scale_y_continuous(limits = c(-10, 20),breaks = seq(-10,20,5))+
  stat_smooth(method = "lm", formula='y~x', se = FALSE)
#推定量の分布
ggplot(beta, aes(x = beta_1)) +
  geom_histogram(colour="white", fill = "skyblue2",bins=50)+
  scale_x_continuous(limits = c(0.6, 1.4))
#標準誤差の分布
ggplot(data_se) +
  geom_histogram(aes(x = se_classical),colour="white", fill = "skyblue2",bins=50)+
  geom_histogram(aes(x = se_stata),colour="white", fill = "red",alpha=0.2,bins=50)+
  scale_x_continuous(limits = c(0.05, 0.2))+
  geom_vline(xintercept =sd(beta_1) )+ #推定量の標準偏差
  geom_vline(xintercept =mean(se_classical),col="blue")+ #均一分散の標準誤差平均
  geom_vline(xintercept =mean(se_stata),col="red") +#不均一分散の標準誤差平均
  xlab("標準誤差") + 
  ylab ("度数" )+
  theme_grey(base_family = "HiraKakuPro-W3")
#推定量の標準偏差を計算
sd(beta_1) 
#標準誤差の平均を計算
mean(se_classical)
mean(se_HC0)
mean(se_stata)
mean(se_HC2)
mean(se_HC3)

④不均一分散のコード↓

#均一分散
trial_number <- 2000
sample_size <- 1000
beta_1 <- rep(NA, trial_number)
se_classical <-rep(NA, trial_number)
se_HC0 <- rep(NA, trial_number)
se_stata <- rep(NA, trial_number)
se_HC2 <- rep(NA, trial_number)
se_HC3 <- rep(NA, trial_number)
se_HC4 <- rep(NA, trial_number)
u <- rep(NA, sample_size)
min(-11.9+10*x-x^2)
#標本生成、推定値の計算、標準誤差の計算
for(j in 1:trial_number){
  x <- rnorm(n=sample_size, mean =5, sd=1)
  for(i in 1:sample_size){
    sd0 <- -11+10*x[i]-x[i]^2
    uu <- rnorm(n=1, mean =0, sd=ifelse(sd0 > 0, sd0,0.01))  #不均一分散
    u[i] <- uu
  }
  
  y <- 2 + 1*x+ u
  answer <- lm(y~x)
  answer0 <- lm_robust(y ~ x,se_type = "HC0")
  answer1 <- lm_robust(y ~ x,se_type = "stata")
  answer2 <- lm_robust(y ~ x,se_type = "HC2")
  answer3 <- lm_robust(y ~ x,se_type = "HC3")
  
  beta_1[j] <- summary(answer)$coef["x","Estimate"]
  
  se_classical[j] <- summary(answer)$coef["x","Std. Error"]
  se_HC0[j] <- summary(answer0)$coef["x","Std. Error"]
  se_stata[j] <- summary(answer1)$coef["x","Std. Error"]
  se_HC2[j] <- summary(answer2)$coef["x","Std. Error"]
  se_HC3[j] <- summary(answer3)$coef["x","Std. Error"]
}
#データフレーム作成
data0 <- data.frame(x,y,u) 
beta <- data.frame(beta_1)
data_se <- data.frame(se_classical,se_HC0,se_stata,se_HC2,se_HC3) 
#(x,y)の散布図
ggplot(data0, aes(x = x,y= y)) +　
  geom_point() +
  scale_x_continuous(limits = c(0, 10),breaks = seq(0,10,2))+
  scale_y_continuous(limits = c(-40, 60),breaks = seq(-40,60,20))+
  stat_smooth(method = "lm", formula='y~x', se = FALSE)
#推定量の分布
ggplot(beta, aes(x = beta_1)) +
  geom_histogram(colour="white", fill = "skyblue2",bins=50)+
  scale_x_continuous(limits = c(0, 2))
#標準誤差の分布
ggplot(data_se) +
  geom_histogram(aes(x = se_classical),colour="white", fill = "skyblue2",bins=50)+
  geom_histogram(aes(x = se_stata),colour="white", fill = "red",alpha=0.2,bins=50)+
  scale_x_continuous(limits = c(0.3, 0.5))+
  geom_vline(xintercept =sd(beta_1) )+ #推定量の標準偏差
  geom_vline(xintercept =mean(se_classical),col="blue")+ #均一分散の標準誤差平均
  geom_vline(xintercept =mean(se_stata),col="red") +#不均一分散の標準誤差平均
  xlab("標準誤差") + 
  ylab ("度数" )+
  theme_grey(base_family = "HiraKakuPro-W3")
#推定量の標準偏差を計算
sd(beta_1) 
#標準誤差の平均を計算
mean(se_classical)
mean(se_HC0)
mean(se_stata)
mean(se_HC2)
mean(se_HC3)