R でロバストな標準誤差を計算する方法

線形回帰の仮定の 1 つは、モデルの残差が予測変数の各レベルで均等に分散しているということです。

この仮定が満たされない場合、回帰モデルには不均一分散性が存在すると言えます。

これが起こると、モデルの回帰係数の標準誤差は信頼できなくなります。

これを説明するために、不均一分散性に対して「堅牢」であるロバストな標準誤差を計算することができ、回帰係数の真の標準誤差値についてより適切なアイデアを得ることができます。

次の例は、R の回帰モデルの堅牢な標準誤差を計算する方法を示しています。

例: R でのロバストな標準誤差の計算

R に、クラスの 20 人の生徒の学習時間と取得した試験の得点に関する情報を含む次のデータ フレームがあるとします。

 #create data frame
df <- data. frame (hours=c(1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4,
                         4, 5, 5, 5, 6, 6, 7, 7, 8),
                 score=c(67, 68, 74, 70, 71, 75, 80, 70, 84, 72,
                         88, 75, 95, 75, 99, 78, 99, 65, 96, 70))

#view head of data frame
head(df)

  hours score
1 1 67
2 1 68
3 1 74
4 1 70
5 2 71
6 2 75

lm()関数を使用して、時間を予測変数として使用し、スコアを応答変数として使用する回帰モデルを R で近似できます。

 #fit regression model
fit <- lm(score ~ hours, data=df)

#view summary of model
summary(fit)

Call:
lm(formula = score ~ hours, data = df)

Residuals:
    Min 1Q Median 3Q Max 
-19,775 -5,298 -3,521 7,520 18,116 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 71.158 4.708 15.11 1.14e-11 ***
hours 1.945 1.075 1.81 0.087 .  
---
Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.48 on 18 degrees of freedom
Multiple R-squared: 0.154, Adjusted R-squared: 0.107 
F-statistic: 3.278 on 1 and 18 DF, p-value: 0.08696

不均一分散性が回帰モデルの問題であるかどうかを視覚的に確認する最も簡単な方法は、残差プロットを作成することです。

 #create residual vs. fitted plot
plot(fitted(fit), reside(fit))

#add a horizontal line at y=0 
abline(0,0) 

X 軸は応答変数の近似値を示し、Y 軸は対応する残差を示します。

グラフから、近似値が増加するにつれて残差の分散が増加することがわかります。

これは、不均一分散性が回帰モデルの問題である可能性が高く、モデルの要約の標準誤差が信頼できないことを示しています。

堅牢な標準誤差を計算するには、次のようにlmtestパッケージのcoeftest()関数とSandwichパッケージのvcovHC()関数を使用できます。

 library (lmtest)
library (sandwich)

#calculate robust standard errors for model coefficients
coeftest(fit, vcov = vcovHC(fit, type = ' HC0 '))

t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 71.1576 3.3072 21.5160 2.719e-14 ***
hours 1.9454 1.2072 1.6115 0.1245    
---
Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

時間予測変数の標準誤差が、前のモデルの概要の 1.075 から、このモデルの概要では 1.2072 に増加していることに注意してください。

元の回帰モデルには不均一分散性が存在するため、この標準誤差推定値はより信頼性が高く、時間予測変数の信頼区間を計算するときに使用する必要があります。

注: vcovHC()関数で計算する最も一般的な推定タイプは「HC0」ですが、 ドキュメントを参照して他のタイプの推定を見つけることができます。

追加リソース

次のチュートリアルでは、R で他の一般的なタスクを実行する方法について説明します。

R で不均一分散性のホワイト検定を実行する方法
R で線形回帰出力を解釈する方法
R で残差プロットを作成する方法