How can I calculate the variance of interaction term from an equivalent model without interaction?

Lets say that you have access to a model that estimates the mean of four independent groups like m2 below, but these groups have been formed from two factors (a & b) and you want to instead evaluate a model with two main effects (a & b) and an interaction (aXb) like m1.

n=10^3
a=rnorm(n, 0, 1)>0
b=rnorm(n, 0, 1)>0
e=rnorm(n, 0, 1)
y=1+a*.2 + b*-.2 + .3*a*b + e

aXb = a*b

ai=a==TRUE & b==FALSE
bi=b==TRUE & a==FALSE
abi=a==TRUE & b==TRUE

m1 = lm(y~a+b+aXb)
m2 = lm(y~ai+bi+abi)


Sample output below for m1 & m2:

(Intercept)  0.85847    0.06396  13.423  < 2e-16 ***
aTRUE        0.37947    0.09180   4.134 3.87e-05 ***
bTRUE        0.02101    0.08956   0.235    0.815
aXb          0.02870    0.12705   0.226    0.821

(Intercept)  0.85847    0.06396  13.423  < 2e-16 ***
aiTRUE       0.37947    0.09180   4.134 3.87e-05 ***
biTRUE       0.02101    0.08956   0.235    0.815
abiTRUE      0.42918    0.08873   4.837 1.53e-06 ***


The two models are equivalent but parameterized differently and most estimtates are identical except aXb from m1 & abi from m2. It is easy to calculate the point estimate of the interaction (aXb) in m1 from the independent groups modeled in m2 like this:

 summary(m2)$coefficients[4,1]-summary(m2)$coefficients[3,1]-summary(m2)$coefficients[2,1] [1] 0.0286971  And the other way around:  summary(m1)$coefficients[4,1]+summary(m1)$coefficients[3,1]+summary(m1)$coefficients[2,1]
[1] 0.4291785


But how can I calculate the variance/standard error of the interaction term in m1 using only information from m2 (and vice versa)?

• Aren't they the same model? – Aksakal Aug 13 '15 at 20:17
• Yeah, they are equivalent but I have not managed to get the variances from m2 to match that of m1 yet. The problem I am trying to solve is that I have access only to m2 but I need to calculate confidence interval for the interaction term in m1. I updated the question to make that more clear. – Poker Aug 13 '15 at 21:17
• @Aksakal: These are not the same model. a and ai differ; so do b and bi. That's why the var-covar matrices will be difficult to match--you need to compute the change-of-basis matrix from one model to the other. Incidentally, both models include (the same) "interaction" term, making it difficult to figure out what the title means--or even what the purpose of this exercise is. – whuber Aug 13 '15 at 22:22
• @whuber, the problem may seem odd but there is a good reason for the question. I have access to results from models fitted with four independent groups like m2 but the groups have been formed from two factors ( a & b) and thus should be evaluated like two main effects + interaction like m1. Interestingly, as you noticed, a and ai may be parameterized differently but they have identical estimates. in m1 you use a+b+aXb to predict the mean of abi from m2. But I cannot get the se/variance to match. – Poker Aug 13 '15 at 22:34
• Thank you for explaining the origin of the problem. I still wouldn't characterize it in terms of the presence vs the absence of an interaction term, though. If you look at it more generally, the change from one model to the other is simply a change of basis of the model space: including the (implicit) constant, it's a $4\times 4$ matrix. The relationships between the two models can all be written in terms of that change-of-basis matrix. – whuber Aug 13 '15 at 22:42

The details get in the way of seeing the simple relationships involved. This answer addresses the problem by generalizing it to the case of any two linear models that are related by a one-to-one linear reparameterization. The new coefficient estimates are related by the same reparametrization (equation $(1)$ below), while their covariance matrix is conjugated by that linear transformation (equation $(2)$).

The two models are parameterized differently. To understand how they are related, let $\beta$ be the parameters of the first model and $\gamma$ the parameters of the second. For the models to be equivalent, the parameters must be related by an invertible linear transformation $U$, so that

$$\gamma = U\beta;\ \beta = U^{-1}\gamma.$$

These models are

$$y = X\beta + \varepsilon = X(U^{-1}\gamma) + \varepsilon = (XU^{-1})\gamma + \varepsilon = Z\gamma + \varepsilon,$$

identifying $Z = XU^{-1}$ as the model matrix for the second model.

The least squares fits are solutions of the Normal equations,

$$\hat \beta = (X^\prime X)^{-1} X^\prime y$$

and

$$\hat \gamma = (Z^\prime Z)^{-1} Z^\prime y = ((U^\prime)^{-1} X^\prime X U^{-1})^{-1} (U^\prime)^{-1} X^\prime y = U(X^\prime X)^{-1} X^\prime y = U\hat \beta.\tag{1}$$

(Incidentally, once we find a formula for $U$, this will provide a straightforward way to convert between $\hat\beta$ and $\hat\gamma$.)

Letting $\sigma^2$ be the common variance of the $\varepsilon$, the familiar formula for the variance-covariance matrix of the estimates gives

\eqalign{ \text{Cov}(\hat \beta) &= \sigma^2 (X^\prime X)^{-1};\\ \text{Cov}(\hat \gamma) &= \sigma^2 (Z^\prime Z)^{-1} = \sigma^2 ((U^\prime)^{-1} X^\prime X U^{-1})^{-1} \\ &= \sigma^2 U (X^\prime X)^{-1} U^\prime \\ &= U\left(\text{Cov}(\hat \beta)\right)U^\prime. }\tag{2}

(The last equation employed the symmetry of $X^\prime X$.) In practice, $\sigma^2$ will be replaced by an estimate of it obtained from the residuals, $\hat\sigma^2$. That will not change this relationship between the two estimated covariances.

It all comes down to finding $U$. In principle, we already know it--it is determined by the two sets of variables we have created--but $U$, too, can be had by means of the Normal equations. The defining relationship $XU = Z$ is a (multivariate) regression problem that by construction of the variables in the two models has a perfect least squares fit given by

$$U = (X^\prime X)^{-1} X^\prime Z. \tag{3}$$

This, then, is the desired solution: from equations $(2)$ and $(3)$ we have obtained $\text{Cov}(\hat \gamma)$ in terms of $\text{Cov}(\hat \beta)$ and the two model matrices $X$ and $Z$.

To continue the example of the question, let's make the data generation reproducible by setting the seed at the outset:

set.seed(17)


After creating the variables a, b, ..., abi as before, here are the R commands to implement the foregoing calculations:

X <- cbind(1, a, b, aXb)                     # Model 1 matrix
Z <- cbind(1, ai, bi, abi)                   # Model 2 matrix
V <- zapsmall(solve(t(X) %*% X, t(X) %*% Z)) # This particular V has many exact zeros
U <- solve(V)                                # V^{-1}
colnames(U) <- rownames(U) <- colnames(Z)
fit.X <- lm(y ~ a + b + aXb)                 # Model 1
print(U %*% vcov(fit.X) %*% t(U), digits=13) # Model 2's covariance matrix


The last one outputs the covariance matrix of the second model in terms of the first one (which so far is the only one that has been calculated and stored in fit.X):

                                       ai                 bi                abi
0.003780312655181 -0.003780312655181 -0.003780312655181 -0.003780312655181
ai  -0.003780312655181  0.007430767242053  0.003780312655181  0.003780312655181
bi  -0.003780312655181  0.003780312655181  0.007832935967726  0.003780312655181
abi -0.003780312655181  0.003780312655181  0.003780312655181  0.007621353224501


We may check by fitting the second model and extracting its covariance matrix directly:

fit.Z <- lm(y ~ ai + bi + abi)     # Model 2
print(vcov(fit.Z), digits=13)      # Model 2's covariance matrix, again


Compare its output to the preceding:

                   (Intercept)             aiTRUE             biTRUE            abiTRUE
(Intercept)  0.003780312655181 -0.003780312655181 -0.003780312655181 -0.003780312655181
aiTRUE      -0.003780312655181  0.007430767242053  0.003780312655181  0.003780312655181
biTRUE      -0.003780312655181  0.003780312655181  0.007832935967726  0.003780312655181
abiTRUE     -0.003780312655181  0.003780312655181  0.003780312655181  0.007621353224501


They agree through the first 13 significant figures. Switch the roles of the models to convert covariances from the second model to the first.

• Thank you so much for your excellent answer. As you write, the details were confusing from a mathematical perspective because the interaction term was parameterized exactly the same in both models. It was the main effects, or groups groups a, ai & b bi, that were different. Interestingly, from the perspective of interpreting the result it was the opposite. – Poker Aug 14 '15 at 16:34

To get the standard error your are reporting, you can access the variance-covariance matrix of the lm object. To do that you simply need to do:

sqrt(diag(vcov(m1)))


The results of that will be in the order of the coefficients, so the fourth element should be the standard error for the interaction term.

Updated

Now that I better understand your question, the issue you are having derives from using TRUE/FALSE values in a regression in R. An interesting question is why do the logical values behave strangely, and that I don't have a good answer for, but if you do this:

aInt <- as.integer(a)
bInt <- as.integer(b)

m3 = lm(y ~ aInt + bInt + aInt*bInt)


You'll get the results to line up.

(Intercept)  0.85847    0.06396  13.423  < 2e-16 ***
aInt       0.37947    0.09180   4.134 3.87e-05 ***
bInt       0.02101    0.08956   0.235    0.815
aInt:bInt  0.42918    0.08873   4.837 1.53e-06 ***

• Hi there an thank you for your answer. Yes I know i can get the variances/se from a model like that but what I am trying to do is to use the equivalent model (m2) to reconstruct the model with the interaction (m1). And as I show above, the point estimate is easy but I haven't managed to get the variance/se from m2 to match that of m1. – Poker Aug 13 '15 at 21:13
• @Poker I've updated my response to answer your question now that I understand it better. – jknowles Aug 13 '15 at 21:53
• The issue is not integer/boolean types here. The issue is that while the two models are equivalent they are not the same. m2 does not have an interaction term but is parameterized differently as four independent groups, instead of two main effects and an interaction like m1. What I am trying to to is to calculate all effects from m1 using m2 and I can do it with the main effects (they are identical) and the point estimate of the interaction (as shown in my example) but I cannot get the variances/se. – Poker Aug 13 '15 at 22:09