# Calculating combined confidence intervals for categorical (dummy) interactions

I have an OLS regression with interacted categorical variables. One of these variables encodes $$n$$ cross-sectional units, the other indicates time along $$m$$ years. I am interested in estimating unit-time specific effects (conditional on other variables) to get a sense of relative levels and changes over time for each unit.

The model I am fitting has the following form:

$$Y = \alpha + \sum\mathbf{\beta_{n-1}X} + \sum\mathbf{\gamma_{m-1}T} + \sum\mathbf{\delta_{(n-1)(m-1)}XT} + ... + \epsilon$$,

where $$\mathbf{X}$$ is a matrix of $$n-1$$ dummy variables that capture the different units (relative to a natural baseline category), and $$\mathbf{T}$$ is a matrix of $$m-1$$ dummy variables for the years (relative to the first year). The data actually has more than two dimensions, so I am not estimating a single fixed effect for each observation but I ignore this here for readability.

The sum of the coefficient estimates of $$\beta_2 + \gamma_3 + \delta_6$$ (for instance) tells me how unit 2 in year 3 differs from baseline unit 0 in baseline year 0, and so forth.

What I would like to know is: How can I calculate the combined confidence intervals for the sum of these three (and other triples of) coefficients?

There are two related posts here and here, which address this question for combining confidence intervals for two coefficient estimates from an interaction. Based on these posts, I have tried to adapt the procedure to the case with three coefficient estimates, but I am uncertain as to whether my approach is correct.

Here is what I did, illustrated on some toy data:

> set.seed(1)
> unit <- sample(0:2, 100, replace = T)
> yr <- sample(0:3, 100, replace = T)
>
> # for simplicity, create Y based on a linear DGP
> Y <- 1 + .5*unit + 1.5*yr - .9*unit*yr + rnorm(100, sd = 4)
>
> summary(reg <- lm(Y ~ factor(unit) * factor(yr)))

Call:
lm(formula = Y ~ factor(unit) * factor(yr))

Residuals:
Min     1Q Median     3Q    Max
-6.093 -2.373 -0.201  2.433  8.717

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)                  1.374      1.135    1.21   0.2291
factor(unit)1               -2.694      1.702   -1.58   0.1171
factor(unit)2                3.787      2.123    1.78   0.0779 .
factor(yr)1                 -0.429      1.768   -0.24   0.8088
factor(yr)2                  0.534      1.605    0.33   0.7399
factor(yr)3                  4.789      1.853    2.58   0.0114 *
factor(unit)1:factor(yr)1    4.714      2.454    1.92   0.0580 .
factor(unit)2:factor(yr)1   -1.394      2.789   -0.50   0.6183
factor(unit)1:factor(yr)2    3.314      2.293    1.45   0.1519
factor(unit)2:factor(yr)2   -6.800      2.763   -2.46   0.0158 *
factor(unit)1:factor(yr)3    3.121      2.623    1.19   0.2374
factor(unit)2:factor(yr)3   -8.864      2.818   -3.15   0.0023 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.59 on 88 degrees of freedom
Multiple R-squared:  0.302, Adjusted R-squared:  0.214
F-statistic: 3.45 on 11 and 88 DF,  p-value: 0.0005


Thus, for unit 2 in year 3 I get:

coef(reg)[c(3,6,12)]
factor(unit)2               factor(yr)3 factor(unit)2:factor(yr)3
3.78706                   4.78857                  -8.86370


So, $$\beta_2 + \gamma_3 + \delta_6$$ equals:

> sum(coef(reg)[c(3,6,12)])
[1] -0.288069


My attempt to calculate combined confidence intervals for the three coefficients looks like this:

# get the degrees of freedom
mf <- model.frame(reg)
DoF <- nrow(mf) - ncol(mf)
ca_1.96 <- qt(0.975, DoF)

# get the relevant part of the variance-covariance matrix
sigma2mat  <- vcov(reg)[c(3,6,12), c(3,6,12)]
sigma2mat
factor(unit)2 factor(yr)3 factor(unit)2:factor(yr)3
factor(unit)2                   4.50659     1.28760                  -4.50659
factor(yr)3                     1.28760     3.43359                  -3.43359
factor(unit)2:factor(yr)3      -4.50659    -3.43359                   7.94018

# calculate standard errors
SE <- sqrt(sum(sigma2mat))
SE
[1] 1.60474

# calculate confidence intervals
CIs <- sum(coef(reg)[c(3,6,12)]) + c(-ca_1.96, ca_1.96) * SE
CIs
[1] -3.47303  2.89690


Edit:

sum(coef(reg)[c(1,3,6,12)])
[1] 1.08606

CIs <- sum(coef(reg)[c(1,3,6,12)]) + c(-ca_1.96, ca_1.96) * SE
CIs
[1] -1.16605  3.33817


I would be grateful for any pointers as to whether this is the right way of going about this - and if not - how it could be done. Thank you!

• Just a note: when simulating your data, your call to rnorm() should return a number of samples that is equal to the number of rows. You've got 100 rows, but only pull 3 samples. So your "random error" will be repeated to fill up the 100 slots (and you get a warning like longer object length is not a multiple of shorter object length). You should have Y <- 1 + .5*unit + 1.5*yr - .9*unit*yr + rnorm(100, sd = 4) Sep 28, 2021 at 16:59
• Thanks for the note, my mistake. I updated the post accordingly. Sep 28, 2021 at 19:07
• Your simple sum of the covariance matrix entries to get the variance of this combination of predictor values is just a special case of the general formula for the variance of a linear combination of correlated variables. Your special case is with all weights equal to 1. You can thus generalize to any logically consistent number of predictors, predictions based on coefficients for continuous variables, etc. Keep in mind the distinction between confidence intervals and prediction intervals.
– EdM
Sep 28, 2021 at 20:18

You've left out the intercept and its part of the variance-covariance matrix from your predictions - otherwise this looks fine, but you've probably worked harder than is necessary. The predict() function can do the hard part for you:

set.seed(1)
unit <- sample(0:2, 100, replace = T)
yr <- sample(0:3, 100, replace = T)

# for simplicity, create Y based on a linear DGP
Y <- 1 + .5*unit + 1.5*yr - .9*unit*yr + rnorm(100, sd = 4)

reg <- lm(Y ~ factor(unit) * factor(yr))

# now create a row of data to be used for prediction:
pred_data = data.frame( unit=factor(2, levels=0:2), yr=factor(3, levels=0:3) )

# if you want a confidence interval (on the fitted mean line):
predict(reg, pred_data, level=0.95, interval="confidence")

# if you want a predictive interval (adds the estimated residual variance)
predict(reg, pred_data, level=0.95, interval="prediction")


Results are:

> predict(reg, pred_data, level=0.95, interval="confidence")
fit       lwr      upr
1 1.086058 -1.168965 3.341081
> predict(reg, pred_data, level=0.95, interval="prediction")
fit       lwr      upr
1 1.086058 -6.393008 8.565124

• Thanks a lot for your answer and the hint about the intercept. If I include it, I get (almost) the same results as the predict() approach you suggest yields (see: Edit). One thing that I'm still unsure about, however, is how this would travel to a setting where there are additional variables in the model (which is the case in my actual application). Specifically, what would I do a) with the intercept (which then reflect not only the dummy baseline category anymore) and b) with the other covariates when I would want to use predict() to obtain CIs in such a setting? Sep 28, 2021 at 19:21
• I'm sorry, but I'm a bit confused about what you want. You've demonstrated that you can generate estimates and confidence intervals for an arbitrary combination of factor levels. There's nothing particularly different if you have a third factor variable. In that case, your manual calculations or the predict() function would work as well. Can you please explain where you have encountered problems? Sep 28, 2021 at 19:33
• Are you saying that you also have some non-factor variables? If so, your approach would still work but would become more slightly complicated. I can demonstrate, if that is the case. Sep 28, 2021 at 19:38
• Sorry for the confusion. Yes exactly, I have continuous variables in the model, too. I just played around with my toy data (to which I added a continuous variable). It seems that I need to pass a newdata argument to predict() that contains a cell with a 0 for each continuous variable. So if I have one additional variable z, this would like like this: pred_data = data.frame( unit=factor(2, levels=0:2), yr=factor(3, levels=0:3), z = 0 ). Do I understand this correctly? Sep 28, 2021 at 19:47
• Yes, that would give you the predictions where z is equal to zero. You could also set z to other specific values (such as its mean value). The overall variance of the fitted value will be lowest where z is equal to its mean value. Sep 28, 2021 at 19:52