Hoe robuuste standaardfouten in r te berekenen

Eén van de aannames van lineaire regressie is dat de modelresiduen gelijkelijk verspreid zijn op elk niveau van de voorspellende variabele.

Wanneer niet aan deze veronderstelling wordt voldaan, is er sprake van heteroskedasticiteit in een regressiemodel.

Wanneer dit gebeurt, worden de standaardfouten van de regressiecoëfficiënten van het model onbetrouwbaar.

Om hiermee rekening te houden, kunnen we robuuste standaardfouten berekenen, die „robuust“ zijn tegen heteroskedasticiteit en ons een beter idee kunnen geven van de werkelijke standaardfoutwaarden voor de regressiecoëfficiënten.

Het volgende voorbeeld laat zien hoe u robuuste standaardfouten voor een regressiemodel in R kunt berekenen.

Voorbeeld: het berekenen van robuuste standaardfouten in R

Stel dat we het volgende dataframe in R hebben dat informatie bevat over de gestudeerde uren en examenscores behaald door 20 studenten in een klas:

 #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

We kunnen de functie lm() gebruiken om een regressiemodel in R te fitten dat uren als voorspellende variabele en score als responsvariabele gebruikt:

 #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

De eenvoudigste manier om visueel te controleren of heteroskedasticiteit een probleem is in het regressiemodel, is door een residuele plot te maken:

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

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

Op de x-as worden de gefitte waarden van de responsvariabele weergegeven en op de y-as de bijbehorende residuen.

Uit de grafiek kunnen we zien dat de variantie van de residuen toeneemt naarmate de aangepaste waarden toenemen.

Dit geeft aan dat heteroscedasticiteit waarschijnlijk een probleem is in het regressiemodel en dat de standaardfouten van de modelsamenvatting onbetrouwbaar zijn.

Om robuuste standaardfouten te berekenen, kunnen we de functie coeftest() uit het lmtest- pakket en de vcovHC()- functie uit het sandwichpakket als volgt gebruiken:

 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

Merk op dat de standaardfout voor de urenvoorspellingsvariabele is toegenomen van 1,075 in de vorige modelsamenvatting naar 1,2072 in deze modelsamenvatting.

Omdat heteroscedasticiteit aanwezig is in het oorspronkelijke regressiemodel, is deze standaardfoutschatting betrouwbaarder en moet deze worden gebruikt bij het berekenen van een betrouwbaarheidsinterval voor de urenvoorspellingsvariabele .

Opmerking : het meest voorkomende type schatting dat in de functie vcovHC() moet worden berekend, is ‚HC0‘, maar u kunt de documentatie raadplegen om andere typen schattingen te vinden.

Aanvullende bronnen

In de volgende tutorials wordt uitgelegd hoe u andere veelvoorkomende taken in R kunt uitvoeren:

Hoe White’s test voor heteroscedasticiteit in R uit te voeren
Hoe lineaire regressie-uitvoer in R te interpreteren
Hoe maak je een restplot in R