Come calcolare gli errori standard robusti in r

Uno dei presupposti della regressione lineare è che i residui del modello siano equamente dispersi a ciascun livello della variabile predittrice.

Quando questa ipotesi non è soddisfatta, si dice che l’eteroschedasticità è presente in un modello di regressione.

Quando ciò accade, gli errori standard dei coefficienti di regressione del modello diventano inaffidabili.

Per tenere conto di ciò, possiamo calcolare gli errori standard robusti , che sono “robusti” contro l’eteroschedasticità e possono darci un’idea migliore dei veri valori dell’errore standard per i coefficienti di regressione.

L’esempio seguente mostra come calcolare gli errori standard robusti per un modello di regressione in R.

Esempio: calcolo degli errori standard robusti in R

Supponiamo di avere il seguente frame di dati in R che contiene informazioni sulle ore studiate e sui punteggi degli esami ottenuti da 20 studenti in una classe:

 #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

Possiamo utilizzare la funzione lm() per adattare un modello di regressione in R che utilizza le ore come variabile predittiva e il punteggio come variabile di risposta:

 #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

Il modo più semplice per verificare visivamente se l’eteroschedasticità è un problema nel modello di regressione è creare un grafico dei residui:

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

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

L’asse x mostra i valori adattati della variabile di risposta e l’asse y mostra i residui corrispondenti.

Dal grafico possiamo vedere che la varianza dei residui aumenta all’aumentare dei valori adattati.

Ciò indica che l’eteroschedasticità è probabilmente un problema nel modello di regressione e che gli errori standard del riepilogo del modello non sono affidabili.

Per calcolare errori standard robusti, possiamo utilizzare la funzione coeftest() del pacchetto lmtest e la funzione vcovHC() del pacchetto sandwich come segue:

 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

Si noti che l’errore standard per la variabile predittore delle ore è aumentato da 1,075 nel riepilogo del modello precedente a 1,2072 in questo riepilogo del modello.

Poiché l’eteroschedasticità è presente nel modello di regressione originale, questa stima dell’errore standard è più affidabile e dovrebbe essere utilizzata quando si calcola un intervallo di confidenza per la variabile predittrice ore .

Nota : il tipo di stima più comune da calcolare nella funzione vcovHC() è “HC0”, ma puoi fare riferimento alla documentazione per trovare altri tipi di stima.

Risorse addizionali

I seguenti tutorial spiegano come eseguire altre attività comuni in R:

Come eseguire il test di White per l’eteroschedasticità in R
Come interpretare l’output della regressione lineare in R
Come creare una trama residua in R