How to calculate the interaction standard error of a linear regression model in R?

I have the following model and want to make a table with the interpretation of the interaction effects as suggested by Bambor and Clark in this paper. However, I have no idea on how to calculate the $cov(\hat{\beta_1}, \hat{\beta_5})$ in the formula.

The model:

reg <- lm( log(Y) ~ as.factor(X1) +  as.factor(X2) + log(as.numeric(X3)) + log(as.numeric(X4)) + as.factor(X1)*as.factor(X2), data=DB)


The results:

Residuals:
Min      1Q  Median      3Q     Max
-5.1091 -0.3036  0.0294  0.3396  3.6537

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)                                     -0.08848    0.09523  -0.929   0.3531
as.factor(X1)1                                  -0.12795    0.06227  -2.055   0.0402 *
as.factor(X2)1                                   0.05666    0.06694   0.846   0.3976
log(as.numeric(X3))                              0.03602    0.02121   1.699   0.0898 .
log(as.numeric(X4))                              0.97546    0.02671  36.514   <2e-16 ***
as.factor(X1)1:as.factor(X2)1                    0.10733    0.11790   0.910   0.3629
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6729 on 795 degrees of freedom
Multiple R-squared: 0.9161, Adjusted R-squared: 0.9155
F-statistic:  1735 on 5 and 795 DF,  p-value: < 2.2e-16


The formula:

$\sqrt{var(\hat{\beta_1}) + var(\hat{\beta_5}) + 2cov(\hat{\beta_1}\hat{\beta_5})}$

-
Ah, sorry, didn't fully answer your question. $cov(\hat{\beta}_1,\hat{\beta}_2)$ would be the entry in the matrix with column name, in your case, as.factor(X1)1 row name as.factor(X2)1, assuming what you wished to find the covariance for were the coefficients with values -0.12795 and 0.05666, and similarly for the other covariances. –  jbowman Jul 28 '12 at 15:14
Thanks @jbowman, I think now I understand! I must use 0.00340945299? Ok?I just added the covariance matrix to the question. –  Davi Moreira Jul 28 '12 at 15:31
Yes, you have got it! –  jbowman Jul 28 '12 at 15:32
Thanks a lot!!!!!!! –  Davi Moreira Jul 28 '12 at 15:34
I also edited the question to avoid confusion. Thanks @jbowman. –  Davi Moreira Jul 28 '12 at 17:57

You will need a little more information than summary(reg) provides, namely, the covariance matrix of the estimates. vcov(reg) will give that to you:

x1 <- rnorm(100)
x2 <- 0.7*rnorm(100) + 0.7*x1
y <- x1 + x2 + rnorm(100)

reg <- lm(y~x1+x2)
vcov(reg)
(Intercept)           x1           x2
(Intercept)  0.009780556 -0.002229766  0.001652152
x1          -0.002229766  0.016996594 -0.012423096
x2           0.001652152 -0.012423096  0.018662900


The covariance between the coefficient of x1 and x2 is in the cell with row label x1 and column label x2 (or vice versa), and the variance terms are on the diagonal.

The difference between vcov(reg) and summary(reg)$cov is that the latter is not scaled by$\hat{\sigma}^2\$.

-
Thanks again @jbowman! –  Davi Moreira Jul 28 '12 at 17:29