Jak obliczyć solidne błędy standardowe w r

Jednym z założeń regresji liniowej jest to, że reszty modelu są równomiernie rozproszone na każdym poziomie zmiennej predykcyjnej.

Kiedy to założenie nie jest spełnione, w modelu regresji mówi się, że heteroskedastyczność jest obecna.

Kiedy tak się dzieje, błędy standardowe współczynników regresji modelu stają się niewiarygodne.

Aby to uwzględnić, możemy obliczyć solidne błędy standardowe , które są „odporne” na heteroskedastyczność i mogą dać nam lepsze pojęcie o prawdziwych wartościach błędów standardowych dla współczynników regresji.

Poniższy przykład pokazuje, jak obliczyć solidne błędy standardowe dla modelu regresji w R.

Przykład: obliczanie solidnych błędów standardowych w R

Załóżmy, że mamy następującą ramkę danych w języku R, która zawiera informacje o przestudiowanych godzinach i wynikach egzaminów uzyskanych przez 20 uczniów w klasie:

 #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

Możemy użyć funkcji lm() , aby dopasować model regresji w R, który wykorzystuje godziny jako zmienną predykcyjną i wynik jako zmienną odpowiedzi:

 #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

Najłatwiejszym sposobem wizualnego sprawdzenia, czy heteroskedastyczność jest problemem w modelu regresji, jest utworzenie wykresu reszt:

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

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

Oś x pokazuje dopasowane wartości zmiennej odpowiedzi, a oś y pokazuje odpowiednie reszty.

Z wykresu widać, że wariancja reszt wzrasta wraz ze wzrostem dopasowanych wartości.

Wskazuje to, że heteroskedastyczność jest prawdopodobnie problemem w modelu regresji i że błędy standardowe podsumowania modelu są niewiarygodne.

Aby obliczyć solidne błędy standardowe, możemy użyć funkcji coeftest() z pakietu lmtest i funkcji vcovHC() z pakietu kanapkowego w następujący sposób:

 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

Należy zauważyć, że błąd standardowy zmiennej predykcyjnej godzin wzrósł z 1,075 w poprzednim podsumowaniu modelu do 1,2072 w tym podsumowaniu modelu.

Ponieważ w pierwotnym modelu regresji występuje heteroskedastyczność, to oszacowanie błędu standardowego jest bardziej wiarygodne i powinno być stosowane przy obliczaniu przedziału ufności dla zmiennej predykcyjnej godzin .

Uwaga : Najpopularniejszym typem szacunków obliczanych w funkcji vcovHC() jest „HC0”, ale możesz zapoznać się z dokumentacją , aby znaleźć inne typy szacunków.

Dodatkowe zasoby

Poniższe samouczki wyjaśniają, jak wykonywać inne typowe zadania w języku R:

Jak wykonać test White’a na heteroskedastyczność w R
Jak interpretować wyniki regresji liniowej w R
Jak utworzyć wykres rezydualny w R