Linear Regression: Why do the coefficients change on the original IVs when you interact them, and add that new interacted-variable to the model?

Basically I want to know how the 'constant' value differs in each of the following models:

Model 1: DV=income; IV1=gender (0=male, 1=female); IV2=location (0=east, 1=west) Here, I understand the constant is for 'Males in the East'

Model 2: DV=income; IV1=gender (0=male, 1=female); IV2=location (0=east, 1=west); IV3=gender*location (1=females in West, 0=everyone else)

What constitutes the constant term here? I don't understand why my coefficients for 'gender' and 'location', along with the 'constant' are drastically different from Model 1.

Please do help me out. Any insights would be greatly appreciated!

A better question might be why you would expect the terms to stay constant at all. Your uninteracted model represents a restriction on the form that the conditional expectation of income behaves. To see how, recall that the first model is $$Y = \beta_0 + \beta_1 female + \beta_2 west + e$$ where if we interpret $\beta_0,\beta_1,\beta_2$ as capturing correlations rather than any underlying causal structure, we let $\mathbb E[e | female, west] = 0$ by definition. Then $$\mathbb E\left[Y | \lnot female, \lnot west\right] = \beta_0$$ $$\mathbb E\left[Y | female, \lnot west\right] = \beta_0 + \beta_1$$ Therefore, $\beta_1$ is just the difference in the expected income of a female in the east and a male in the east. Similarly, $\beta_2$ is just the expected difference in income between a male in the west and a male in the east. But now, the first model also asserts that $$\mathbb E\left[Y | female, west\right] = \beta_0 + \beta_1 + \beta_2$$ In other words, the difference between a male in the east and a female in the west is $\beta_1 + \beta_2$. Therefore, the model asserts that $$\mathbb E\left[Y | female, west\right] = \mathbb E\left[Y | female, \lnot west\right] + \mathbb E\left[Y | \lnot female, west\right] - \mathbb E\left[Y | \lnot female, \lnot west\right]$$ It should be clear that not every distribution satisfies this. Consider as an example $$\mathbb E\left[Y | \lnot female, \lnot west\right] = 1$$ $$\mathbb E\left[Y | female, \lnot west\right] = 2$$ $$\mathbb E\left[Y | \lnot female, west\right] = 3$$ $$\mathbb E\left[Y | female, west\right] = 5$$ Under our model, we should have $\mathbb E\left[Y | female, west\right] = 4$ given the first 3 expectations, so our functional form is not flexible enough to allow that last expectation to be anything other than 4 fixing the first 3 expectations. OLS will therefore (in the plim) have to do its best to fit this distribution by getting it slightly wrong for each expectation (intuitively, the fact that OLS minimizes squared error means that it will prefer to be slightly wrong for each estimate rather than very wrong from some estimate and completely right for another). In particular, the fact that $\mathbb E\left[Y | female, west\right]$ is higher than our functional form would imply given the first three conditional expectations will tend to push the estimates of $\beta_1,\beta_2$ up and $\beta_0$ down to compensate (as an exercise, try to reason about why this would be the case). By contrast, if you fit model 2: $$Y = \beta_0 + \beta_1 female + \beta_2 west + \beta_3 female * west + e$$ you give OLS enough freedom to fit all of the features of the data, so the bias from the previous specification gets absorbed into the interaction term (so in this example, we can set $\beta_0 = 1,\beta_1 = 1,\beta_2 = 2,\beta_3 = 1$). If the violations of model 1 are more dramatic than my simple numerical example, the discrepancy can be even larger.