I am curious about how interaction terms are calculated within a multiple linear regression model. I understand that in a simple linear regression, the following holds:
$$y_i=\beta_0+\beta_1x_i+\epsilon_i $$ $$\hat{\beta_1}=\frac{\Sigma(x_i-\bar{x})y_i}{\Sigma(x_i-\bar{x})^2}$$ $$var(\hat{\beta_1})=\frac{\sigma^2}{\Sigma(x_i-\bar{x})^2}$$
Now, if we extend this to a multiple linear regression model, I understand that estimates corresponding to the above are:
$$y_i=\beta_0+\beta_1x_i+\beta_2z_i+\epsilon_i $$
$$\hat{\beta_1}=\frac{\Sigma\hat{r}_{i1}y_i}{\Sigma\hat{r}^2_{i1}}$$ and $$var(\hat{\beta_1})=\frac{\sigma^2}{\Sigma(x_i-\bar{x})^2(1-R^2_1)}$$
where $r^2_1$ are the residuals from a simple regression of $x$ on $z$, and $R^2_1$ is the R-squared from regressing $x$ on $z$. This is following a partialling out interpretation of multiple linear regression.
My question is how the above can be expressed for an interaction term. Say we have the model:
$$y_i=\beta_0+\beta_1x_i+\beta_2z_i+\beta_3xz_i+\epsilon_i $$
Do I substitute the product of $x$ and $z$ into the above equations? This would seem strange, as it would rely upon regressing the product of $x$ and $z$ upon $x$ and $z$ to calculate $R^2_1$. Otherwise is there some aspect of the estimation procedure of which I am unaware? It would be great to know from where these estimates magically appear.
