Heteroskedasticity Question

Question

I have a model that's affected by Heteroskedasticity:

bptest(m1)

studentized Breusch-Pagan test

data: m1 BP = 65.055, df = 6, p-value = 4.205e-12

In this model I had 2 variables (yrs.since.phd and yrs.service) that have significant coefficients:

summary(m1)

Call:
lm(formula = salary ~ ., data = data)

    Residuals:
       Min     1Q Median     3Q    Max 
    -65248 -13211  -1775  10384  99592 

    Coefficients:
                  Estimate        Std.Error       tvalue        Pr(>|t|)    

    (Intercept)    78862.8     4990.3  15.803  < 2e-16 ***

    rankAsstProf  -12907.6     4145.3  -3.114  0.00198 ** 

    rankProf       32158.4     3540.6   9.083  < 2e-16 ***

    disciplineB    14417.6     2342.9   6.154 1.88e-09 ***

    yrs.since.phd    535.1      241.0   2.220  0.02698 *  

    yrs.service     -489.5      211.9  -2.310  0.02143 *  

    sexMale         4783.5     3858.7   1.240  0.21584    
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 22540 on 390 degrees of freedom
Multiple R-squared:  0.4547,    Adjusted R-squared:  0.4463 
F-statistic:  54.2 on 6 and 390 DF,  p-value: < 2.2e-16

I tried to fix this model by using robust standard errors, and now the two variables mentioned above (yrs.since.phd and yrs.service) are not significant (or barely significant):

coeftest(m1, df = Inf, vcovHC(m1, omega = NULL, type = "HC4"))

z test of coefficients:

     Estimate Std. Error z value  Pr(>|z|)    

(Intercept)    78862.82    3936.67 20.0329 < 2.2e-16 ***

rankAsstProf  -12907.59    2229.92 -5.7884 7.108e-09 ***

rankProf       32158.41    2343.04 13.7251 < 2.2e-16 ***

disciplineB    14417.63    2320.44  6.2133 5.188e-10 ***

yrs.since.phd    535.06     320.96  1.6670   0.09551 .  

yrs.service     -489.52     315.23 -1.5529   0.12045    

sexMale         4783.49    2457.27  1.9467   0.05157 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Could anyone please explain me what's happened? All other variables Std.errors decreased, while happened the opposite for these two...

Thanks!

James

Answer 1

Likely, the heteroscedasticity led to high influence points in which the regression model had anticonservative estimates of residual variance.

I would believe the robust standard error estimates, since it is robust and gives approximately correct standard errors regardless of the actual distribution of the error term.

Other considerations would be whether you neglected to use an appropriate change of variable. I rarely think it's important to justify such a change based on the findings in the analyses themselves, it's evident "p-hacking" to do so. However, sometimes a log transform or square root transform makes scientific sense. A biproduct of such an approach is that one may be able to stabilize variance.

One sensitivity analysis you should consider is histograms of the dfbetas, called using dfbeta. We would be concerned if the trends were being driven by a single observation, leading to irregular DF beta plots. Then we would consider presenting analyses which do and do not include the high influence observation.