Standard Errors for 2X2 Regression If we have two dichotomous variables variables and run an interaction constructing the means for each cell seems straight forward, but not the standard errors. I would like to use these standard errors to compute if the differences between any two given cells is significantly different.
Y = Variable of interest (continuous and normally distributed) 
T = Treatment (dichotomous) 
A = Variable splitting sample by Age of 21 and older (dichotomous)

With the following (fake) output from an OLS regression :
         b    se   N
T       0.5   0.1  90
A       2.0   0.4  100
T*A     0.9   0.3  75
Cons    0.8   0.1  110

I interpret the above for the mean Y. The Cons representing mean Y for those under 21 and untreated (A = 0, T = 0); T representing mean Y for the treated individuals under 21 (A = 0, T = 1); A representing the mean Y for those 21 and older and untreated (A = 1, T = 0); T*A representing the mean Y for those for those 21 and older and treated (A = 1, T = 1).
1) Is this interpretation of means correct?
2) How are the corresponding standard errors of the means computed
3) Is it statistically acceptable to then compare the means across this 2X2?
I can check the above by hands and calculating means, but it becomes unclear when you include controls in the regression.
 A: *

*Yes, your interpretation of the coefficients is correct.

*This gets a bit more tricky / abstract. 

*Yes, you can compare estimated coefficients.


Question 2: where do standard errors come from?
Refresher:
Recall some standard OLS formulas. Estimate of coefficients $\mathbf{b}$:
$$ \mathbf{b} = (X'X)^{-1}(X'\mathbf{y}) $$
Regular OLS estimate of the variance of the estimate (i.e. estimate of $Var(\mathbf{b})$):
$$ V_{\mathrm{ols}} = (X'X)^{-1} s^2 \quad \quad s^2 = \frac{\boldsymbol{e}'\boldsymbol{e}}{n-k}$$
Regular OLS estimate of the variance of the estimate (i.e. estimate of $Var(\mathbf{b})$):
$$ V_{\mathrm{robust}} = (X'X)^{-1} \left(X'\Omega X \right)  (X'X)^{-1}\frac{n}{n-k} \quad \quad \Omega = \mathrm{diag}\left(u^2_1, u^2_2, \ldots, u^2_n \right)$$
Special case of indicators: example
Imagine we have:


*

*5 observations

*2 categories


Also for simplicity, let's construct a data matrix with exclusive indicators rather than a constant + a single indicator. (For simplicity, I want the columns of X orthogonal.) Eg. our data matrix $X$ might be:
$$ X = \left[ \begin{array}{cc} 1 & 0 \\ 1 & 0 \\ 0 & 1 \\ 0 & 1 \\ 0 & 1 \end{array}\right]$$
That is, the first two observations are category A and 2nd three observations are category B. Hence:
$$ X'X = \left[ \begin{array}{cc} 2 & 0 \\ 0 & 3 \end{array} \right] \quad \quad (X'X)^{-1} = \left[ \begin{array}{cc} \tfrac{1}{2} & 0 \\ 0 & \tfrac{1}{3} \end{array} \right]$$
$X'X$ is simply the counts of variables in each category! $X'\mathbf{y}$ sums up observations by category and multiplying by $(X'X)^{-1}$ divides by the counts. $\mathbf{b} = (X'X)^{-1}X'\mathbf{y}$ in this case is simply the category by category mean (i.e. this is your Q1).
The regular OLS, homoskedastic standard errors would be
$$V_{\mathrm{ols}} = (X'X)^{-1} s^2 = \left[ \begin{array}{cc} \tfrac{s^2}{2} & 0 \\ 0 & \tfrac{s^2}{3} \end{array} \right]$$
That is, it's the estimate of the variance of the residuals divided by the number of observations (we're more certain for categories where we have more observations).
The robust covariance matrix is a bit trickier, but in some sense, easier to interpret! Let $u = [u_1, u_2, u_3, u_4, u_5]'$ be a vector of residuals.
The robust estimator turns into (dropping the finite sample degrees of freedom correction $\frac{n}{n-k}$ to keep it simpler):
$$\begin{align*} V_{\mathrm{robust}} &= \left[ \begin{array}{cc} \tfrac{1}{2} & 0 \\ 0 & \tfrac{1}{3} \end{array} \right] \left[ \begin{array}{cc} u_1^2 + u_2^2 & 0 \\ 0 & u_3^2 + u_4^2 + u_5^2 \end{array} \right] \left[ \begin{array}{cc} \tfrac{1}{2} & 0 \\ 0 & \tfrac{1}{3} \end{array} \right] \\
&= \left[ \begin{array}{cc} \tfrac{1}{2} \left( \frac{u_1^2 + u_2^2 }{2} \right)  & 0 \\ 0 & \tfrac{1}{3}\left( \frac{u_4^2 + u_5^2 + u_6^2}{3} \right) \end{array} \right]
\end{align*}$$
That is, it's an estimate of the variance for each category estimated separately, divided by the number of observations of that category!
Basically, (1) OLS with indicators and robust standard errors in a single regression is equivalent to (2) calculating the mean and standard error for each category separately!
Question 3:
Most any statistical package will have a way to output matrix $V$ described earlier.
Let's say you're interested in some $\mathbf{a}' \mathbf{b}$. Asymptotically, we'll have:
$$ \mathbf{a}' \mathbf{b} \sim \mathcal{N}\left(\mathbf{a}' \boldsymbol{\beta}, \mathbf{a}'V\mathbf{a}  \right) $$
Example: let's say you're interested in $b_1$ - $b_3$. The standard error would be given by:
$$ \sqrt{\left[\begin{array}{ccc}1 & 0 & -1 \end{array} \right] V \left[\begin{array}{c}1 \\ 0 \\ -1 \end{array} \right]} $$
A: In most statistical packages I believe Cons (mean for A = 0, T = 0), Cons+T (mean for A = 0, T = 1), Cons+A (mean for A = 1, T = 0) and Cons+A+T+A*T (mean for A = 1, T = 1) would have the interpretations you describe. (2) is complicated, because the estimates are not independent, but most statistical packages will have an option for giving you either the contrasts you want or outputting the estimated covariance matrix for the estimates that let's you calculate it by hand.
I must be missing some crucial point regarding the third question - of course one can compare the means across the 2x2 table. Or is the question how to do this, or whether some multiplicity adjustment is needed?
