0
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ Since factor variables are dummy encoded automatically, in the same way as with continuous variables. Study the source code of predict.lm. $\endgroup$ – Roland Jan 8 at 14:02
1
$\begingroup$

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.$$

$\endgroup$
  • $\begingroup$ thanks for your answear, the input X seems strange, why is the input X 2 dimension vector ? With dummy coding on a 2 levels factor, i except to have only one column (arbitrarily isMale or isFemale) $\endgroup$ – philippe Jan 15 at 15:58
  • $\begingroup$ there are indeed different forms of modelling the influence of the factor variable and different corresponding $X$-matrices. however you want to measure two effects (the effect of eacht of the levels of the gender-factor) and hence you need two columns. Try for example: model <- lm(salary ~ gender, data = df,contrasts = list(gender = "contr.sum")) and take a look at the model.matrix(model). $\endgroup$ – chRrr Jan 15 at 16:34

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.