Como calcular erros padrão robustos em r

Uma das suposições da regressão linear é que os resíduos do modelo estão igualmente dispersos em cada nível da variável preditora.

Quando esta suposição não é atendida, diz-se que a heterocedasticidade está presente em um modelo de regressão.

Quando isso acontece, os erros padrão dos coeficientes de regressão do modelo tornam-se não confiáveis.

Para dar conta disso, podemos calcular erros padrão robustos , que são “robustos” contra a heterocedasticidade e podem nos dar uma ideia melhor dos verdadeiros valores de erro padrão para os coeficientes de regressão.

O exemplo a seguir mostra como calcular erros padrão robustos para um modelo de regressão em R.

Exemplo: cálculo de erros padrão robustos em R

Suponha que temos o seguinte quadro de dados em R que contém informações sobre as horas estudadas e as notas dos exames obtidas por 20 alunos em uma turma:

 #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

Podemos usar a função lm() para ajustar um modelo de regressão em R que usa horas como variável preditora e pontuação como variável de resposta:

 #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

A maneira mais fácil de verificar visualmente se a heterocedasticidade é um problema no modelo de regressão é criar um gráfico residual:

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

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

O eixo x mostra os valores ajustados da variável resposta e o eixo y mostra os resíduos correspondentes.

No gráfico podemos ver que a variância dos resíduos aumenta à medida que os valores ajustados aumentam.

Isto indica que a heterocedasticidade é provavelmente um problema no modelo de regressão e que os erros padrão do resumo do modelo não são confiáveis.

Para calcular erros padrão robustos, podemos usar a função coeftest() do pacote lmtest e a função vcovHC() do pacote sanduíche da seguinte forma:

 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

Observe que o erro padrão para a variável preditora de horas aumentou de 1,075 no resumo do modelo anterior para 1,2072 neste resumo do modelo.

Como a heterocedasticidade está presente no modelo de regressão original, esta estimativa do erro padrão é mais confiável e deve ser utilizada no cálculo de um intervalo de confiança para a variável preditora de horas .

Nota : O tipo de estimativa mais comum para calcular na função vcovHC() é ‘HC0’, mas você pode consultar a documentação para encontrar outros tipos de estimativa.

Recursos adicionais

Os tutoriais a seguir explicam como realizar outras tarefas comuns em R:

Como realizar o teste de White para heterocedasticidade em R
Como interpretar a saída da regressão linear em R
Como criar um gráfico residual em R