回帰係数の標準誤差 / 均一分散vs不均一分散

 回帰係数の標準誤差は「不均一分散に頑健な標準誤差(ロバスト標準誤差)」を使うべきと言われます。

 では、通常の標準誤差を使うとどうなるのでしょうか?

要約

(1)問い

  • 通常の標準誤差
  • 不均一分散に頑健な標準誤差

を使った場合に、何が異なるのでしょうか?

(2)結論

 誤差項が不均一分散のときは、「不均一分散に頑健な標準誤差」の方がバイアスが小さいです。

 シミュレーションの例は下です。黒が推定量の母標準偏差、青が「通常の標準誤差」の分布、赤が「不均一分散に頑健な標準誤差」の分布です。

黒が推定量の母標準偏差青が均一分散の標準誤差赤が不均一分散の標準誤差

(3)意義

 「不均一分散に頑健な標準誤差を使うべき」という主張が腑に落ちます。

前提

(1)定義

 標準誤差

  • 推定量の母標準偏差SDの推定値

です。詳しくは「標準誤差」をご覧ください。

(2)仮定する条件

推定量の望ましさはバイアスの小ささ

$$バイアス=E(SE)-SD$$

が小さいとき、その標準誤差は望ましいと考えることにします。なお、標準誤差SE、標準誤差の期待値E(SE)、推定量の母標準偏差SDです。

共通の仮定する条件

 今回は

  1. データ(X,Y)は独立で同一の分布に従い、単回帰モデルに従う
  2. 標本が無作為抽出
  3. 説明変数Xの分散が0でない
  4. 外生性がある
  5. 誤差項の系列相関がない

を仮定します。

場合によって仮定する条件

 そして、

  • 均一分散
  • 不均一分散

の両方の場合について考えます。

(3)事実

標準誤差のRによる計算

 通常の標準誤差(Std. Error)は、Rでは次のlm関数で求められます。

#均一分散が仮定された標準誤差を用いた回帰分析
lm(y ~ x)
 
#回帰
lm0 <- lm(y ~ x)  #回帰分析をlm0と名づける
summary(lm0) #lm0を詳しく表示して「Std. Error」が標準誤差

不均一分散に頑健な標準誤差の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と同じこともできます

方法

(1)場合分け

 誤差項の分散が

  1. 均一分散
  2. xが大きいと誤差項の分散が大きくなる不均一分散
  3. xが平均に近いと誤差項の分散が小さくなる不均一分散
  4. xが平均に近いと誤差項の分散が大きくなる不均一分散

について場合分けして考えます。

(2)シミュレーション

 シミュレーションは次のように行います。

  1. 誤差項の分散について仮定
  2. サンプル・サイズ1000の標本を生成
  3. 最小二乗推定値を計算
  4. 標準誤差を計算
  5. 「2」〜「4」を2000回繰り返す

 この結果、

  • 2000個の最小二乗推定値
  • 2000個の標準誤差

が得られます。

  • 推定値の標準偏差=母標準偏差SD (黒線
  • 標準誤差の平均=標準誤差の期待値E(SE) (均一分散は青い線不均一分散は赤い線

とみなします。

結果

(1)均一分散①

 誤差項の分散の仮定を「均一分散」とします。その結果、生成された標本の一つが↓です。

 この標本を2000回作った結果、最小二乗推定量の標準偏差は0.0623でした。そこで、最小二乗推定量の母標準偏差は0.0623とみなします。

$$SD=0.0623$$

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

黒が推定量の母標準偏差青が均一分散の標準誤差赤が不均一分散の標準誤差

 それぞれの標準誤差SEの平均を、標準誤差の期待値E(SE)とみなします。結果は、以下の通りでした。

$$E(SE_{均一分散})=0.0624$$

$$E(SE_{不均一分散HC0})=0.0622$$

$$E(SE_{不均一分散stata})=0.0623$$

$$E(SE_{不均一分散HC2})=0.0623$$

$$E(SE_{不均一分散HC3})=0.0625$$

 どれもバイアスは小さいですが、均一分散なら「均一分散の標準誤差」を用いた方がよいですね。

(2)不均一分散②:xが大きいと分散大

 誤差項の分散の仮定を「不均一分散②:xが大きいと誤差項の分散が大きくなる」とします。その結果、生成された標本の一つが↓です。

 この標本を2000回作った結果、最小二乗推定量の標準偏差は0.0503でした。そこで、最小二乗推定量の母標準偏差は0.0503とみなします。

$$SD=0.0503$$

 このときの標準誤差の分布は、下図でした。青が均一分散の標準誤差で、赤が不均一分散の標準誤差(stata, HC1)です。縦線は、黒が推定量の標準偏差、青が均一分散の標準誤差の平均、赤が不均一分散の標準誤差の平均です。

黒が推定量の母標準偏差青が均一分散の標準誤差赤が不均一分散の標準誤差

 それぞれの標準誤差SEの平均を、標準誤差の期待値E(SE)とみなします。結果は、以下の通りでした。

$$E(SE_{均一分散})=0.0484$$

$$E(SE_{不均一分散HC0})=0.0501$$

$$E(SE_{不均一分散stata})=0.0501$$

$$E(SE_{不均一分散HC2})=0.0502$$

$$E(SE_{不均一分散HC3})=0.0503$$

 不均一分散なら「不均一分散に頑健な標準誤差」を用いた方がよいですね。ただ、不均一分散に頑健な標準誤差を用いても、0.048が0.050になるだけですし、あんまり大きな影響はありませんね。それに不均一分散に頑健な標準誤差の分散の大きさも気になります。

(3)不均一分散③:xが平均近傍で分散小

 誤差項の分散の仮定を「不均一分散③:xが平均に近いと誤差項の分散が小さくなる」とします。その結果、生成された標本の一つが↓です。

 この標本を2000回作った結果、最小二乗推定量の標準偏差は0.146でした。そこで、最小二乗推定量の母標準偏差は0.146とみなします。

$$SD=0.146$$

 このときの標準誤差の分布は、下図でした。青が均一分散の標準誤差で、赤が不均一分散の標準誤差(stata, HC1)です。縦線は、黒が推定量の標準偏差、青が均一分散の標準誤差の平均、赤が不均一分散の標準誤差の平均です。

黒が推定量の母標準偏差青が均一分散の標準誤差赤が不均一分散の標準誤差

 それぞれの標準誤差SEの平均を、標準誤差の期待値E(SE)とみなします。結果は、以下の通りでした。

$$E(SE_{均一分散})=0.0771$$

$$E(SE_{不均一分散HC0})=0.145$$

$$E(SE_{不均一分散stata})=0.145$$

$$E(SE_{不均一分散HC2})=0.146$$

$$E(SE_{不均一分散HC3})=0.146$$

 今回は、明らかに「不均一分散に頑健な標準誤差」を用いた方がよいですね。不均一分散に頑健な標準誤差を用いると、0.077が0.146に修正されます。

(4)不均一分散④:xが平均近傍で分散大

 誤差項の分散の仮定を「不均一分散④:xが平均に近いと誤差項の分散が大きくなる」とします。その結果、生成された標本の一つが↓です。

 この標本を2000回作った結果、最小二乗推定量の標準偏差は0.359でした。そこで、最小二乗推定量の母標準偏差は0.359とみなします。

$$SD=0.359$$

 このときの標準誤差の分布は、下図でした。青が均一分散の標準誤差で、赤が不均一分散の標準誤差(stata, HC1)です。縦線は、黒が推定量の標準偏差、青が均一分散の標準誤差の平均、赤が不均一分散の標準誤差の平均です。

黒が推定量の母標準偏差青が均一分散の標準誤差赤が不均一分散の標準誤差

 それぞれの標準誤差SEの平均を、標準誤差の期待値E(SE)とみなします。結果は、以下の通りでした。

$$E(SE_{均一分散})=0.414$$

$$E(SE_{不均一分散HC0})=0.357$$

$$E(SE_{不均一分散stata})=0.357$$

$$E(SE_{不均一分散HC2})=0.357$$

$$E(SE_{不均一分散HC3})=0.358$$

 今回も「不均一分散に頑健な標準誤差」を用いた方がよいですね。不均一分散に頑健な標準誤差を用いると、0.41が0.35に修正されます。

考察

(1)結論

 不均一分散に頑健な標準誤差は、不均一分散の状況下でのバイアスが小さいです。

(2)前提評価

 この記事では、バイアスの大小だけで、標準誤差の計算方法の望ましさを判断しました。これには一定の合理性があります。

 ただし、バイアスの大小以外にも、標準誤差の計算方法の望ましさを判断する方法はあるようにみえます。

 例えば、下図を見ればわかるように、不均一分散に頑健な標準誤差の分散は大きく、推定量の母標準誤差とかけ離れた推定値になる確率も高いです。ならバイアスがあったとしても、分散の小さい通常の標準誤差を使用する方がいい場合もありそうです。

黒が推定量の母標準偏差青が均一分散の標準誤差赤が不均一分散の標準誤差

(3)意義

  バイアスの大小に注目するならば、「不均一分散に頑健な標準誤差を使うべき」という主張が腑に落ちたのではないでしょうか?

付録1:頑健な標準誤差と有意性

結論

「不均一分散に頑健な標準誤差を使うと、常に、仮説検定にて有意な結果が出にくくなる」は誤りです。

反例

 なぜなら、次の反例があるからです。

 不均一分散④の場合、不均一分散に頑健な標準誤差を使うと、仮説検定にて有意な結果が出やすくなります。↓

黒が推定量の母標準偏差青が均一分散の標準誤差赤が不均一分散の標準誤差

 不均一分散④は

  • Xの平均に近いところの誤差項の分散が大きく
  • Xの平均から遠いところの誤差項の分散が小さい

場合です。↓

↑④

数式による理解

「不均一分散に頑健な標準誤差を使うと、仮説検定にて有意な結果が出にくくなる、とは限らない」理由は、数式で理解できます。

 単回帰分析と「均一分散の標準誤差」によると均一分散の標準誤差は

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

$$= \sqrt {\frac{ \sum\limits_{i=1}^n (x_i-\overline{x})^2 \widehat{\sigma^2} }{ \left(\sum\limits_{i=1}^n (x_i-\overline{x})^2 \right) ^2} }$$

です。単回帰分析と「不均一分散に頑健な標準誤差」によると不均一分散に頑健な標準誤差SEは

$$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} }$$

となります。標準誤差の数式の分子に注目すると、

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

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

となっています。これは④のように平均値から外れた誤差項の分散が小さいタイプの不均一分散では、

$$\sum\limits_{i=1}^n(x_i-\overline{x})^2 \widehat{\sigma^2}  >  \sum\limits_{i=1}^n(x_i-\overline{x})^2 \widehat{\sigma_i^2}$$

となり、

  • 不均一分散に頑健な標準誤差 < 均一分散の標準誤差

がおきます。逆に、③のように平均値から外れた誤差項の分散が大きいタイプの不均一分散では、

  • 不均一分散に頑健な標準誤差 > 均一分散の標準誤差

となるのです。

付録2: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)

コメント欄 お気軽にコメントをお寄せください!

タイトルとURLをコピーしました