# Using White's Robust Co-variance Matrix vs Weighted Least Squares to correct for heteroscedasticity

I've been trying to figure this out for a bit. How does using White's robust co-variance matrix in OLS vs weighted least squares affect mean response confidence intervals?

I've experimented with both of these and the mean responses are very similar, but my confidence intervals for the mean responses vary quite a bit.

Here is some sample code. I am using python statsmodels library. My features consist of 36 qualitative inputs and two quantitative inputs with 516 observations. I estimated the standard deviation function by regressing the residuals against the fitted values using OLS, then used those as the weights (as laid out by Kutner et al.) I use the get_prediction method of my 'reg' model instance below to get the mean and confidence interval.

Using White's Robust Covariance:

import statsmodels.api as sm
mod = sm.OLS(y, X)
reg = mod.fit(cov_type='HC0')
reg.get_prediction(exog=data).summary_frame()


Using Weighted Least Squares:

mod = sm.WLS(y, X, weights=weights)
reg = mod.fit()
reg.get_prediction(exog=data).summary_frame()

• What program are you using to perform your calculations? Can you paste in your code? Different programs have different ways of calculating these. Jan 29, 2019 at 21:39
• I updated with some additional details as well as some sample code. Jan 31, 2019 at 15:19

Under the usual Gauss-Markov Model where $$\boldsymbol{Y}=\boldsymbol{X\beta}+\boldsymbol{\epsilon}$$, it is assummed $$\boldsymbol{\epsilon}\sim N\left(0,\boldsymbol{I}\sigma^{2}\right)$$, where $$\boldsymbol{I}$$ is an $$n\times n$$ identity matrix, where $$n$$ is the number of observations in your dataset. This implies $$Var(\boldsymbol{Y})=Var\left(\boldsymbol{X\beta}+\boldsymbol{\epsilon}\right)=Var(\boldsymbol{\epsilon})=\boldsymbol{I}\sigma^{2}.$$ Then, under this model:

$$\begin{eqnarray*} Var\left(\hat{\boldsymbol{\beta}}\right) & = & Var\left[\boldsymbol{\left(X^{\prime}X\right)^{-}X^{\prime}Y}\right]\\ & = & Var\left[\boldsymbol{AY}\right]\\ & = & \boldsymbol{A}Var(\boldsymbol{Y})\boldsymbol{A^{\prime}}\\ & = & \boldsymbol{A}\left(\boldsymbol{I}\sigma^{2}\right)\boldsymbol{A}^{\prime}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(*)\\ & = & \left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\right]\left(\boldsymbol{I}\sigma^{2}\right)\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\right]^{\prime}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]\left[\boldsymbol{X}^{\prime}\boldsymbol{X}\right]\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]\\ & = & \sigma^{2}\boldsymbol{I}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]\\ & = & \sigma^{2}\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-} \end{eqnarray*}$$

And normally we don't know $$\sigma^{2}$$ so it is estimated from the residuals to give us: $$Var\left(\boldsymbol{\hat{\beta}}\right)=\hat{\sigma}^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\right]$$. Most statistical software will calculate the variance of the estimated coefficients by using this formula. But what happens now, if we use a model where we still have $$\boldsymbol{Y}=\boldsymbol{X\beta}+\boldsymbol{\epsilon}$$, but this time it is assummed $$\boldsymbol{\epsilon}\sim N\left(0,\boldsymbol{\Sigma}\right)$$? Well, in this case, the variance formula changes in (*) above to be:

$$\begin{eqnarray*} Var\left(\hat{\boldsymbol{\beta}}\right) & = & \boldsymbol{A\Sigma}\boldsymbol{A}^{\prime}\\ & = & \left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\right]\boldsymbol{\Sigma}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\right]^{\prime}\\ & = & \left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{\Sigma X}\left(\boldsymbol{X}^{\prime}\boldsymbol{X}\right)^{-} \end{eqnarray*}$$

So, now, you see that the formula has changed a bit. Since we no longer can assume homoscedastic variances, we have the covariance matrix, $$\boldsymbol{\Sigma}$$ sandwiched into the middle of the formula for the variance of the estimated coefficients (hence the nickname the 'Huber-White "Sandwich" Estimator'). Again, we typically do not know this, so we usually estimate $$\boldsymbol{\Sigma}$$ with $$\hat{\boldsymbol{\Sigma}}$$ from the residuals from an initial fitting by Ordinary Least Squares Regression and then once $$\hat{\boldsymbol{\Sigma}}$$ is estimated, the $$Var\left(\hat{\boldsymbol{\beta}}\right)$$ can be estimated. You may obtain the specifics of how each program goes about estimating $$\hat{\boldsymbol{\Sigma}}$$ exactly by referring to their documentation.

Now, under weighted least squares, the estimate of $$\boldsymbol{\beta}$$ is now weighted and takes on different form. It is of the form:

$$\begin{eqnarray*} \hat{\boldsymbol{\beta}}_{weighted} & = & \left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{WY} \end{eqnarray*}$$

where $$\boldsymbol{W}$$ is a diagonal matrix of weights (usually formed by taking the inverse of the fitted value from some variance function). The variance of this estimate is different. It is of the form:

$$\begin{eqnarray*} Var\left(\hat{\boldsymbol{\beta}}_{weighted}\right) & = & Var\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{WY}\right]\\ & = & Var\left[\boldsymbol{AY}\right]\\ & = & \boldsymbol{AY}\boldsymbol{A}^{\prime}\\ & = & \boldsymbol{A}\sigma^{2}\boldsymbol{W}\boldsymbol{^{-}A}^{\prime}\\ & = & \sigma^{2}\boldsymbol{A}\boldsymbol{W}^{-}\boldsymbol{A}^{\prime}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{W}\right]\boldsymbol{W}^{-}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{W}\right]^{\prime}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{W}\right]\boldsymbol{W}^{-}\boldsymbol{W}^{\prime}\boldsymbol{X}\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\boldsymbol{X}^{\prime}\boldsymbol{W}\right]\boldsymbol{X}\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\\ & = & \sigma^{2}\left[\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\right]\left[\boldsymbol{X}^{\prime}\boldsymbol{W}\boldsymbol{X}\right]\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\\ & = & \sigma^{2}\boldsymbol{I}\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-}\\ & = & \sigma^{2}\left(\boldsymbol{X}^{\prime}\boldsymbol{WX}\right)^{-} \end{eqnarray*}$$

and $$\sigma^{2}$$ is can be estimated here by a number of different methods, including use of replicates or near replicates, or by modelling the variance function. See Chapter 11 of Applied Linear Statistical Models by Kutner, Nachtsheim, Neter, and Li, 5th ed. for additional details.

• How is this possible? $Var(\boldsymbol{Y})=Var\left(\boldsymbol{X\beta}+\boldsymbol{\epsilon}\right)=Var(\boldsymbol{\epsilon})=\boldsymbol{I}\sigma^{2}$ Jan 30, 2019 at 0:45
• @HeteroskedasticJim, because $\boldsymbol{X}$ and $\boldsymbol{\beta}$ are considered constants in the Markov Model. Jan 30, 2019 at 21:56
• @StatsStudent. This is actually the exact book I am using that raised this question. What you have stated is fairly clear, but I'm still not understanding how this affects the confidence intervals for the mean response. The bottom of page 426 states that weighted least squares uses the modified variance for the mean response confidence intervals, but page 427 seems to imply that White's estimator is only used for coefficient std errors. I may not be connecting things through. I will post some code above which may help. Jan 31, 2019 at 14:53
• Perhaps this is an issue of different fields having different terminology/notation, but I've never heard a linear regression described as a "Markov model" before. I would use "Markov model" to describe an entirely different class of techniques. This is a bit of an aside from the main point, because the rest of your answer is wonderful, but I was curious as to whether you could elaborate or provide a reference for a linear regression being referred to as such. Could take it to the chat if mods deem it too extraneous. Jan 31, 2019 at 15:10
• @RyanSimmons, that should have read "Gauss-Markov model" instead of just "Markov-Model." This is the term you'll find in pure statistics due to the Gauss-Markov theorem. I'm travelling today, so don't my references in front of me, but I think you'll find this reference in "Applications of Linear and Nonlinear Models Fixed Effects, Random Effects, and Total Least Squares." by Grafarend and Awange and in "A Primer on Linear Models" by John F. Monahan. I'll also try to get to the the other Ryan's (the OP's) questions later today if time allows when I'm back to a computer. Jan 31, 2019 at 15:33