6
$\begingroup$

In Applied Linear Statistical Models (Kutner et al) an example on weighted least squares is given.

Dataset: Age versus Diastolic Blood pressure

Regressing DBP on age results in the following model $$\hat Y = 56.157 + 0.58003 X$$

But a residual plots show the megaphone-alike shape which would imply heteroscedacity.

Residual plot

Since the absolute values of the residuals show a linear trend versus the age predictor OLS was performed on these residuals versus age, which resulted in $$\hat s = -1.54946 + 0.198172 X$$

Using these fitted values as weights $$w_i = \dfrac{1}{\hat s_i^2}$$ resulted in a new model: $$\hat Y = 55.566 + 0.59634 X$$.

Some questions remain, first denote the model as:

$$Y_i = \beta_0 + \beta_1 X_i + \epsilon_i$$

Why not use the residuals as estimator of the weights?

Since $w_i$ is defined as $\dfrac{1}{\sigma_i^2}$ where $\epsilon_i \stackrel{\text{d}}{=} N(0,\sigma^2_i)$ and $$\text{Var}(\epsilon_i) = \sigma^2_i = E(\epsilon_i^2) - E(\epsilon_i)^2 = E(\epsilon_i^2) $$ Now since $e_i = (\hat Y_i - Y_i)$ is an estimator of $\epsilon_i$ why not use $e_i^2$ as an estimator of $\sigma^2_i$? Which would imply $w_i = \dfrac{1}{e^2_i}$.

I don't really understand why we must regress the residuals once more to find the weights.

How much does Weighted Least Squares fix?

I performed the analysis above in R and plotted a new residual plot of the weighted least squares fit, which results in: Second residual plot

This doesn't seem much better though?

R-code

colnames(Blood_Pressure_Example) <- c("Age", "DBP");
attach(Blood_Pressure_Example)

plot(DBP~Age)
fit <- lm(DBP~Age)
abline(fit)
summary(fit)

plot(residuals.lm(fit)~Age)
plot(abs(residuals.lm(fit))~Age)

## Like in the book
fit.res <- lm(abs(residuals.lm(fit))~Age)
wii <- 1/predict.lm(fit.res)^2

fit.end <- lm(DBP~Age, weights = as.vector(wii))
summary(fit.end)
plot(residuals.lm(fit.end)~Age)

## why not use the residuals as weights?
residuals.lm(fit)
wii.test <- 1/(residuals.lm(fit)^2)
wii.test

fit.test <- lm(DBP~Age, weights=as.vector(wii.test))
summary(fit.test)
plot(residuals.lm(fit.test)~Age)  
$\endgroup$
  • 1
    $\begingroup$ The second plot looks identical to the first. $\endgroup$ – Alecos Papadopoulos Jan 12 '17 at 3:32
  • $\begingroup$ There is a tiny difference, but should it be a lot better? As I understand it the residuals where Age is high should be less in the weighted model, because it weighted them down? $\endgroup$ – dietervdf Jan 12 '17 at 3:34
  • $\begingroup$ You also implemented your idea. What was the residual plot here? You should upload it also $\endgroup$ – Alecos Papadopoulos Jan 12 '17 at 3:44
  • $\begingroup$ It also looked extremly similar, the fitted coefficients are a bit different though: $\hat Y = 57.03481 + 0.56127\cdot X$ $\endgroup$ – dietervdf Jan 12 '17 at 3:48
  • 3
    $\begingroup$ Normally the residual plots will look similar after running WLS unless high-leverage points have been severely downweighted. (The residual plot is not supposed to look homoscedastic: it's the plot of weighted residuals that should look like nice.) Often the interesting changes due to weighting occur in the standard errors, the p-values, and prediction and confidence limits. As far as using residuals for (inverse squared) weights goes, what do you suppose would happen if any of the original residuals was close to zero? $\endgroup$ – whuber Jan 12 '17 at 17:21
2
$\begingroup$

It would be too noisy to estimate weights as squared residuals. Consider this: you're estimating n weights using n observations. It's one observation per parameter. It's like estimating the variance of the population having a sample size one.

So, instead you observe that the variance seems to increase linearly with X, then you regress on X. You end up estimating just two coefficients from n observations. You get much more reliable estimate of the weight this way.

On the residuals looking the same after weighted least squares: that's not surprising given that the estimated weights are not that different.

$\endgroup$
1
$\begingroup$

I suspect the issue here is that the weighting $w_i = 1/e_i^2$ would make the regression estimates insufficiently sensitive to the response variable, since it would almost be tantamount to treating each deviation as having unit magnitude.

To see what I mean, suppose you let $e_i$ be the residuals from the first model fit, and suppose you then fit a new model using your specified weighting, and let $e_i^*$ be the residuals from this second model fit. (And more generally, I will put a star on all estimated quantities from the second model.) The coefficient vector in the second model fit is being estimated by minimising the objective:

$$\text{WRSS}(\boldsymbol{\hat{\beta}^*}) \equiv (\boldsymbol{y} - \boldsymbol{x} \boldsymbol{\hat{\beta}^*})^{\boldsymbol{\text{T}}} \boldsymbol{\text{W}} (\boldsymbol{y} - \boldsymbol{x} \boldsymbol{\hat{\beta}^*}) = \sum_{i=1}^n \Bigg( \frac {y_i - \sum_{k=0}^m \hat{\beta}_k^* x_{i,k}}{e_i} \Bigg)^2 = \sum_{i=1}^n \Bigg( \frac{e_i^*}{e_i} \Bigg)^2 .$$

Evaluating near the previous parameter estimate we have $\boldsymbol{\hat{\beta}^*} \approx \boldsymbol{\hat{\beta}}$ which gives $e_i^* \approx e_i$ so that:

$$\text{WRSS}(\boldsymbol{\beta}) \approx \sum_{i=1}^n \Bigg( \frac {e_i^*}{e_i} \Bigg) ^2 = \sum_{i=1}^n 1^2 =n.$$

So you can see that under this method, when you fit the second model, you are effectively treating each deviation from the regression line as having approximately unit magnitude, regardless of the actual response value in the data. That is, you are looking to minimise the weighted residuals, but these weighted residuals are (almost by definition) close to unit magnitude.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.