Як обчислити стійкі стандартні помилки в 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) 

На осі абсцис відображаються підігнані значення змінної відповіді, а на осі у – відповідні залишки.

З графіка ми бачимо, що дисперсія залишків збільшується зі збільшенням підігнаних значень.

Це вказує на те, що гетероскедастичність, ймовірно, є проблемою в регресійній моделі і що стандартні помилки підсумкової моделі є ненадійними.

Щоб обчислити надійні стандартні помилки, ми можемо використати функцію 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