Formula used in confidence intervalle on R's lm function

Question

Does someone know how confidence interval on factor values are computed on R (lm function), here is a simple example :

data

df <- data.frame(c("Male", "Female", "Male","Female"),c(900,600,1200,800)) names(df) <- c("gender","salary") df$gender <- as.factor(df$gender)

model

model <- lm(salary ~ gender, data = df)

Male confidence intervale

predict(model, list(gender="Male"), interval = "confidence")

fit lwr upr 1050 501.5172 1598.483

Female confidence intervale

predict(model, list(gender="Female"), interval = "confidence")

fit lwr upr 700 151.5172 1248.483

How these lower and upper bounds have been computed ?

nb : this question is linked to Confidence interval in R lm function with factors values

Answer 1

Note that you are in fact estimating an ANOVA model with "treatment contrasts" (the baseline is set to "Female"). This simple ANOVA is equivalent to a simple linear regression model with some dummy-variable encoding: By running lm() you are in fact running a regression of $y$ on $X$, where $y=(900,600,1200,800)^T$ and $X$ is given by: $$X=\begin{pmatrix}1 & 1\\ 1 & 0 \\ 1 & 1 \\ 1 & 0 \end{pmatrix}.$$ (see model.matrix(model) to get $X$).

Hence, you are estimating a regression line in a simple linear regression model of the form $y= \alpha + \beta x + \epsilon_i$ and are interested in a predicition interval of the a mean response at $x=1$ (this is NOT at a new observation!).

It is well known that under the usual assumptions (including the normality assumption on $\epsilon_i$), the prediction interval for the mean response at a in-sample $x$ is then given by $$\widehat{y}_x \pm t_{1-\frac{\alpha}{2},n-2}SE(\widehat{y}_x),$$ where $t_{\alpha,n}$ is the $\alpha$ quantile of a $t$-distribution with $n$-degrees of freedem, $s= \sqrt{\frac{1}{n-2}\sum_{i=1}^n e_i^2}$, $\widehat{y}_x = \widehat{\alpha} + \widehat{\beta} x$ and $SE(\widehat{y})$ is given by:: $$SE(\widehat{y}_x) = s \sqrt{\frac{1}{n}+\frac{(x-\overline{x})^2}{\sum(x_i-\overline{x})^2}}.$$

In your case, for $x=1$, you have $\widehat{\alpha} = 700$, $\widehat{\beta}=350$, $s\approx 180.2776$, $SE(\widehat{y}_{x=1}) \approx 127.4755$, $t_{0.975,2}\approx 4.302653 $. Hence your $0.95$-prediction interval is given by: $$\widehat{y}_{x=1} \pm t_{1-\frac{\alpha}{2},n-2}SE(\widehat{y}_{x=1}) = 1050 \pm 548.4828.$$