Confidence intervals of linear model with several factors

I am trying to interpret the coefficients and confidence intervals of a model that has several categorical variables. I cannot share this model but I managed to narrow down my situation to a simple example.

Let's say that I want to model the weight of animals from three categorical features of the animals: cat/dog, house/wild and friendly/unfriendly. The following toy dataset creates such

> require(dplyr)
> set.seed(1234)
> toy = expand.grid(animal=c('cat', 'dog'),
habitat=c('house', 'wild'),
friend=c(T, F),
observation=seq(100)) %>%
mutate(weight=ifelse(animal=='cat', 15, 20) +
ifelse(habitat=='house', 0, 5) +
rnorm(n=nrow(.))) %>%
select(-observation)  # drop observation, not used

animal habitat friend   weight
1    cat   house   TRUE 13.79293
2    dog   house   TRUE 20.27743
3    cat    wild   TRUE 21.08444
4    dog    wild   TRUE 22.65430
5    cat   house  FALSE 15.42912
6    dog   house  FALSE 20.50606
7    cat    wild  FALSE 19.42526
8    dog    wild  FALSE 24.45337

So I create a model weight ~ animal + habitat + friend. Since there are three categories, coding them as 2^3=8 columns would give a rank-deficient matrix, but R knows this and it creates the correct matrix design by grouping these binary factors. We can verify this with the model.matrix function:

> lm1 = lm(weight ~ animal + habitat + friend, data=toy)
> model.matrix(lm1) %>% head(8)
(Intercept) animaldog habitatwild friendTRUE
1           1         0           0          1
2           1         1           0          1
3           1         0           1          1
4           1         1           1          1
5           1         0           0          0
6           1         1           0          0
7           1         0           1          0
8           1         1           1          0

And the model results:

> summary(lm1)

Call:
lm(formula = weight ~ animal + habitat + friend, data = toy)

Residuals:
Min      1Q  Median      3Q     Max
-3.3347 -0.6526 -0.0069  0.6284  3.2479

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.96788    0.07094 211.003   <2e-16 ***
animaldog    4.98014    0.07094  70.205   <2e-16 ***
habitatwild  4.99066    0.07094  70.354   <2e-16 ***
friendTRUE   0.03934    0.07094   0.555    0.579
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.003 on 796 degrees of freedom
Multiple R-squared:  0.9254,    Adjusted R-squared:  0.9252
F-statistic:  3293 on 3 and 796 DF,  p-value: < 2.2e-16

So this makes sense to me, but the intercept absorbs all cases where the binary variables are zero. In other words, the intercepts represents the mean weight for unfriendly house cats. This is not really a problem, I just need to be careful when interpreting the model coefficients.

For example, unfriendly house cats weigh about 14.96, but if I want to know the weight of friendly wild dogs it would be 14.96 + 4.98 + 4.99 + 0.03 = 24.96.

Now, I turn to confidence intervals, I have them calculated with a bootstrap method, but here let's assume that they have the same results as confint:

> confint(lm1)
2.5 %     97.5 %
(Intercept) 14.82863826 15.1071290
animaldog    4.84089540  5.1193861
habitatwild  4.85141476  5.1299055
friendTRUE  -0.09990898  0.1785817

The question:

What would be the confidence interval for the weight of friendly wild dogs?

Naturally, I thought that this could be estimated as the sum of all confidence intervals above. For unfriendly house dogs, it would be the confidence interval of (Intercept) + the interval of animaldog.

While this seems logic to me (it would be like the calculation of the coefficients in the first part), I am not entirely sure and it feels like I am missing something.

Method 1: Get the estimated covariance matrix of $$\hat \beta$$. (maybe the software will not give it to you by default, but you can force software by adding some options,...)
Use $$Var(\sum \hat \beta_i)=\sum(Var(\hat \beta_i) + 2\sum_{i>j} Cov(\hat \beta_i,\hat \beta_j)$$.
Get the CI by $$\sum \hat \beta_i \pm Z \sqrt{Var(\sum \hat \beta_i)}$$.
• Interesting. Method 2 makes a lot of sense, but it would cost me some computation time. I could use Method 1 since I do have the covariance matrix of $\hat\beta$, but I wonder if this has some underlying normal assumption? Do you have a link/reference that one could read concerning this approach? – YuppieNetworking Oct 17 '18 at 18:30
• Method 1 is based on the assumption that $\hat\beta$ has multivariate normal distribution. If $\epsilon$ follows normal, of course $\hat\beta$ will be normal. But even $\epsilon$ not normal, $\hat\beta$ will be asymptotic normal if your sample size is large enough. From your results on se and CI on $\hat\beta$ it seems assuming $\hat\beta$ be normal is acceptable. I have no link/ref because this kind of thing is very old and common. – user158565 Oct 17 '18 at 18:43