Diagnostics part 2: error covariance assumptions
10/21/2020
Visual diagnostics
Courtesy of: Faraway, J. Linear Models with R. 2nd Edition
- To the left is a plot where:
- The horizontal (x-axis) represents the fitted value \( \hat{\mathbf{Y}} \), as in some model \( \hat{\mathbf{Y}}=\mathbf{X}\hat{\boldsymbol{\beta}} \).
- The vertical axis (y-axis) represents the corresponding residual \( \hat{\boldsymbol{\epsilon}} \) value, equal to \( \mathbf{Y} - \hat{\mathbf{Y}} \).
- Suppose that the variances of the error \( \boldsymbol{\epsilon} \) are not fixed, i.e.,
- suppose that some observations have more variation and some observations have less variation around the signal, \( \mathbf{X}\boldsymbol{\beta} \).
- In this case, there is dependence of the variation (or error \( \boldsymbol{\epsilon} \)) on the observation.
- To the left, this is actuall a well behaved situation. In particular, the residuals show variation that doesn’t seem to depend on the value of the fitted value/ observation.
- The mean of the residuals appears to be zero as well, because there isn’t a clear preference for positive or negative residuals
- The situation to the left is denoted homoscedasticity.
- Later, we will develop tests to determine if the variance is “close-enough-to-constant”, so that we are satisfied.
- On the other hand, if there are non-constant variances (as discussed above), this will often show up as a pattern in the plot to the left.
Visual diagnostics
Courtesy of: Faraway, J. Linear Models with R. 2nd Edition
- The non-constant variance of \( \hat{\boldsymbol{\epsilon}} \) to the left is known as heteroscedasticity.
- In this case, there is a clear dependence on the variation of \( \hat{\boldsymbol{\epsilon}} \) on the fitted value/ the observation.
- Breaking the constant variance assumption, the Gauss-Markov theorem no longer applies
- Heteroscedasticity does not, in itself, cause ordinary least squares coefficient estimates to be biased;
- however, the theoretical estimates of the variance of the residuals (and therefore the standard errors) will become biased.
Visual diagnostics
Courtesy of: Faraway, J. Linear Models with R. 2nd Edition
- The bias in the standard errors complicates our ability to accurate quantify the uncertainty, and therefore to accurately:
- make hypothesis tests on the parameters for significance;
- provide confidence intervals for the parameters
- provide accurate prediction intervals and confidence intervals for the mean response;
- provide explanatory power in the relationship between the response and the explanatory variables.
Visual diagnostics
Courtesy of: Faraway, J. Linear Models with R. 2nd Edition
- Using the same plot as before, we can likewise determine if there is actually nonlinearity in the residuals.
- The plot at the left exhibits a nonlinear dependence of the residual on the fitted/ observed values
- The same technique can also be used replacing the fitted values \( \hat{\mathbf{Y}} \) in the horizontal axis with any other variable in the model \( \mathbf{X}_i \), to determine dependence of the residual on the explanatory variables
- Additionally, we may consider how the residual varies across variables \( \mathbf{X}_i \) that we have data for, but have not included as explanatory variables on the response.
- If there is a dependence structure for the residuals on the variable that was left out of the model, it suggests that we should consider its impact on the response and/or if it is a variable that is tightly correlated with our existing model variables.
Savings data example
- With this model, we plot the residuals versus the fitted values:
par(mai=c(1.5,1.5,.5,.5), mgp=c(3,0,0))
plot(fitted(lmod_savings),residuals(lmod_savings),xlab="Fitted",ylab="Residuals", cex=3, cex.lab=3, cex.axis=3)
abline(h=0)
- From the initial inspection, we do not notice anything particularly structured about the residuals – indeed they are roughly symmetric around zero and in the fitted value.
Half normal distributions
Courtesy of Nagelum CC BY-SA 4.0 via Wikimedia Commons
- Therefore, to investigate the constant variance assumption more closely, we will “zoom” in on the residuals.
- In particular, we will consider the plot of the square root of the absolute residuals, i.e., \[ \begin{align} \sqrt{\vert\hat{\boldsymbol{\epsilon}}\vert} \end{align} \]
- Under the assumption \( \boldsymbol{\epsilon}\sim N(0, \mathbf{I}\sigma^2) \), the values of \( \vert\hat{\boldsymbol{\epsilon}}\vert \) are distribted according to the half normal distribution
- the relationship between the normal distribution and the corresponding half-normal distribution is illustrated on the left.
- Having excluded the nonlinearity expressed in other plots, we will focus on the absolute values to increase the resolution of the residuals as a response of the fitted values.
Savings data example
- In the case below, we plot \( \sqrt{\vert\hat{\boldsymbol{\epsilon}}\vert} \) versus the fitted values
par(mai=c(1.5,1.5,.5,.5), mgp=c(3,0,0))
plot(fitted(lmod_savings),sqrt(abs(residuals(lmod_savings))), xlab="Fitted",ylab=expression(sqrt(hat(epsilon))), cex=3, cex.lab=3, cex.axis=1.5)
- We note in this case as well (irrespective of sign) the residuals have approximately even variation.
Plotting rediduals versus predictors
- We will now consider the plot of the residual versus the explanatory variables of the model:
par(mai=c(1.5,1.5,.5,.5), mgp=c(3,0,0))
plot(savings$pop15,residuals(lmod_savings), xlab="Population under 15", ylab="Residuals", cex=3, cex.lab=3, cex.axis=1.5)
abline(h=0)
We note that there are two distinct groups illustrated in the population below 15,
This reflects the multi-modality that was previously observed in the Q-Q plots.
Plotting rediduals versus predictors
- On the other hand, the residuals are of uniform variation with respect to the values for the population over 75.
par(mai=c(1.5,1.5,.5,.5), mgp=c(3,0,0))
plot(savings$pop75,residuals(lmod_savings),xlab="Population over 75",ylab="Residuals", cex=3, cex.lab=3, cex.axis=1.5)
abline(h=0)
- We are motivated thus to determine if there is non-constant variance in the residuals with respect to the pop15 predictor…
F-test continued…
Courtesy of IkamusumeFan CC BY-SA 4.0
- Consider the following null and alternative hypothesis: \[ \begin{align} H_0 : &\sigma^2_X = \sigma^2_Z \\ H_1 : & \sigma^2_X \neq \sigma^2_Z \end{align} \]
- Under the previous assumptions, and the null hypothesis \( H_0 \), the ratio of the two sample based variances is an F-statistic \[ \begin{align} \frac{S_X^2}{S_Z^2}\sim F_{(n-1,m-1)} \end{align} \] a distribution defined by the ratio of two \( \chi^2 \) variables, divided by their respective degrees of freedom.
- If the p-value of the above ratio is greater than the pre-specified significance (\( 5\% \) typically) then we fail to reject the null hypothesis, i.e., \( \sigma^2_X = \sigma^2_Z \).
- On the other hand, if we observe an extreme value for this ratio, this gives significant evidence to support inequality of the variances.
Returning to the regression analysis…
Courtesy of: Faraway, J. Linear Models with R. 2nd Edition
- We wish to determine if there is non-constant variance in the residuals with respect to the values of the explanatory variable “pop15”.
- This is of course, a violation of the Gauss-Markov theorem assumptions. We want to test this concretely with the F-test for variances.
- For the F-test discussed earlier we take
- \( \left\{X_i \right\}_{i=1}^n \) to be the residuals corresponding to a country with the percent of population under 15, less than \( 35\% \);
- \( \left\{Z_i \right\}_{i=1}^m \) to be the residuals corresponding to a country with the percent of population under 15, greater than \( 35\% \).
- Assuming that both sample populations are normally distributed with unknown variance, we will use an F-test to determine if the variances are equal.
- In our earlier analysis, we failed to reject the null hypothesis of Gaussian error distributions, so this assumption is justified.
An example of time-correlated observations
Ice core sample hole. Image courtesy of NOAA - National Centers for Environmental Information
- If we wish to model and understand the changing climate, we must bear in mind that there are only accurate measures of temperatures available since the early 1850’s.
- To estimate the average land surface temperature in a region for dates further back in time, we can use proxy data including, e.g,
- tree rings, which show the relative intensity of winters, summers and precipitation;
- ice core samples, which give stratified layers demonstrating the relative intensity of summer melting and precipitation events;
- corals, which by chemical properties of their skeletons can be used to estimate the ocean temperature in which they formed;
- and others…
- In this context, we can consider building a linear model for the average northern hemisphere surface temperature as the response, with the proxy data as the explanatory variables.
- Fitting the relationship where there are reliable observations of the temperature, we can predict the average surface temperature when we lack these futher back in time, based upon the proxies which extend much further.
An example of correlated observations
Q: What kind of plot would make sense in this instance to study the
time-correlation of the observations?A: we can once again plot the residuals, but in this context versus the time variable
par(mai=c(1.5,1.5,.5,.5), mgp=c(3,0,0))
plot(residuals(lmod) ~ year, na.omit(globwarm), ylab="Residuals", cex=3, cex.lab=3, cex.axis=1.5)
abline(h=0)
- Q: do we detect any structure to the residuals?
An example of correlated observations
Given the last diagnostic, it appears that the residuals are positively correlated in time,
- that is, the residual in the last time step appears to largely be similar to the outcome of the subsequent.
Therefore, we will plot the residual one step in time forward as a function of the residual one step in time backward:
par(mai=c(1.5,1.5,.5,.5), mgp=c(3,0,0))
n <- length(residuals(lmod))
plot(tail(residuals(lmod),n-1) ~ head(residuals(lmod),n-1), xlab=expression(hat(epsilon)[i]),ylab=expression(hat(epsilon)[i+1]), cex=3, cex.lab=3, cex.axis=1.5)
abline(h=0,v=0,col=grey(0.75))
Durbin-Watson test
Image courtesy of Geek3 CC BY 3.0 via Wikimedia Commons.
- Given the lag-one autocorrelation model in the residuals, \( \hat{\boldsymbol{\epsilon}}_t = \rho \hat{\boldsymbol{\epsilon}}_{t-1} + \nu_t \), we will consider the following hypotheses:
- \[ \begin{align} H_0: &\rho = 0 \\ H_1: &\rho \neq 0 \end{align} \] where the null corresponds to the absence of a time-correlation of the residuals.
- Under the null hypothesis, the following, \[ \begin{align} DW \triangleq \frac{\sum_{i=2}^n \left(\hat{\boldsymbol{\epsilon}}_{i} - \hat{\boldsymbol{\epsilon}}_{i-1}\right)^2 }{\sum_{i=1}^n \hat{\boldsymbol{\epsilon}}_i^2} \end{align} \] is distributed according to a linear combination of \( \chi^2 \) distributions (example Chi-squared distributions pictured left).