Как рассчитать надежные стандартные ошибки в r

Одно из предположений линейной регрессии заключается в том, что остатки модели одинаково разбросаны на каждом уровне переменной-предиктора.

Если это предположение не выполняется, говорят, что в регрессионной модели присутствует гетероскедастичность .

Когда это происходит, стандартные ошибки коэффициентов регрессии модели становятся ненадежными.

Чтобы учесть это, мы можем рассчитать надежные стандартные ошибки , которые «устойчивы» к гетероскедастичности и могут дать нам лучшее представление об истинных значениях стандартных ошибок для коэффициентов регрессии.

В следующем примере показано, как вычислить надежные стандартные ошибки для регрессионной модели в 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 показаны соответствующие остатки.

На графике мы видим, что дисперсия остатков увеличивается по мере увеличения подобранных значений.

Это указывает на то, что гетероскедастичность, вероятно, является проблемой в регрессионной модели и что стандартные ошибки сводки модели ненадежны.

Для расчета надежных стандартных ошибок мы можем использовать функцию coeftest() из пакета lmtest и функцию 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