Calculating confidence intervals for the variance of the residuals in R

I have three variables:

• Number of house sales
• Month (in couples)
• Region of a city (N-W-E-S)

and I want to calculate confidence intervals for the residual of the errors. So, given the data:

month <- c("1", "1", "1", "1", "2", "2", "2", "2", "3", "3", "3",
"3", "4", "4", "4", "4", "5", "5", "5", "5", "6", "6", "6",
"6")
region <-c("1", "2", "3", "4", "1", "2", "3", "4", "1", "2", "3",
"4", "1", "2", "3", "4", "1", "2", "3", "4", "1", "2", "3",
"4")
sales <-c(85, 107, 61, 22, 40, 65, 58, 51, 60, 41, 45, 27, 15,
30, 68, 63, 28, 3, 57, 12, 36, 21, 10, 16)

data <- cbind(sales, month, region)
data <- as.data.frame(data)

salesmod <- lm(sales ~ month + region, data=data)
summary(salesmod)
anova(salesmod)


We can check the degrees of freedom and the sum of the residuals from the summary(salesmod) and anova(salesmod).


Call:
lm(formula = sales ~ month + region, data = data)

Residual standard error: 22.88 on 15 degrees of freedom
Multiple R-squared:  0.4855,    Adjusted R-squared:  0.2111
F-statistic: 1.769 on 8 and 15 DF,  p-value: 0.1623

> anova(salesmod)
Analysis of Variance Table

Response: sales
Df Sum Sq Mean Sq F value Pr(>F)
month      5 6368.7 1273.74  2.4325 0.0835 .
region     3 1042.8  347.60  0.6638 0.5871
Residuals 15 7854.5  523.63
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


So the variance of the residuals would be $$523.63/15=34.90867$$, but how do I compute a confidence intervals for this value (of given 95% confidence).

You seem to have made a small mistake. The variance of residuals is $$7854.5/15=523.63$$ (you have divided twice).

Now, what you are looking for is distribution of the estimate of the variance of true errors ($$\varepsilon$$) so that you can construct a confidence interval for it.

First let $$\boldsymbol{\varepsilon} \sim N(\mathbf{0},\sigma^2I)$$. Now the estimate for $$\sigma^2$$ is (where $$K$$ is the number of paramters in the model):

\begin{align} s^2 \equiv \frac{\mathbf{e'e}}{n-K} = \frac{\boldsymbol{\varepsilon'}M\boldsymbol{\varepsilon}}{n-K} \end{align}

So for this you can use the following property that (taken directly from this wiki page)

If $$Y$$ is a vector of $$k$$ i.i.d. standard normal random variables and $$A$$ is a $$k\times k$$ symmetric, idempotent matrix with rank $$k-n$$, then the quadratic form $$Y^{T}AY$$ is chi-square distributed with $$k-n$$ degrees of freedom.

Since $$M$$ is idempotent and symmetric with rank $$n-K$$, it is clear that:

$$\Big(\frac{\boldsymbol{\varepsilon'}}{\sigma}\Big)M\Big(\frac{\boldsymbol{\varepsilon}}{\sigma}\Big) = \frac{n-K}{\sigma^2}s^2 \sim \chi^2_{(n-K)}$$

This can be used now to construct the confidence interval:

$$Pr\bigg(q_{\alpha/2}<\frac{n-K}{\sigma^2}s^2 < q_{1-\alpha/2}\bigg)=1-\alpha$$

$$Pr\bigg(\frac{n-K}{q_{1-\alpha/2}}s^2 < \sigma^2 < \frac{n-K}{q_{\alpha/2}}s^2 \bigg) = 1-\alpha$$

where $$q_p$$ is the $$p^{th}$$ quantile of $$\chi^2_{(n-K)}$$

Based on your data, we can use qchisq in R to get these quantilies.

# lower bound:
sum(salesmod$residuals^2)/qchisq(0.975,df=15) [1] 285.7373 # upper_bound sum(salesmod$residuals^2)/qchisq(0.025,df=15)
[1] 1254.277

• This is a quite non-robust interval, it would be interesting to see the answer augmented with a more robust alternative! Commented May 8, 2021 at 14:41
• @kjetilbhalvorsen: that would be great! I am certainly interested in a better answer as this was just an attempt from me. Commented May 8, 2021 at 16:27