Why do my (coefficients, standard errors & CIs, p-values & significance) change when I add a term to my regression model? Lots of people seem to be asking this. They often seem to get shallow answers that merely assert what is true, instead of drawing or explaining the mechanism. They also seem to not find each other -- if you don't know the answer, it's not obvious that these are all the same question.
Has someone given a high-quality answer that we can refer all these askers to? If not, give it your best shot below.
Here are some related questions. This list is incomplete. Feel free to edit and add to it.
Linear Regression Coefficient changes with additional variables
GAM Interactions : Individual and Combined Interactions are different
Nonsignificant interaction still causes main effect to flip?
Interpreting main effect coefficient in different models
Interactions make terms significant in regression when they should not be
Related but not bang on:
What if interaction wipes out my direct effects in regression?
Understanding Simpson's paradox: Andrew Gelman's example with regressing income on sex and height
 A: The ordinary least squares solution is simply given by:
$$\beta = (X'X)^{-1}X'y$$
Let's imagine we augment $X_{n\times p}$ with one or more variables $\tilde X_{n\times \tilde p}$, appending its corresponding values as columns, and call the resulting matrix ${X^*}_{n\times p^*}$, $p^* = p + \tilde p$.
Now, given enough degrees of freedom, coefficients will be given by:
$$\beta^* = ({X^*}'{X^*})^{-1}{X^*}'y$$
If $\left\{A\right\}_{k,:}$ represents the $k$-th row of a matrix $A$, the $k$-th coefficient in each vector corresponds to:
$$\beta_k = \left\{(X'X)^{-1}\right\}_{k,:}X'y$$
$$\beta^*_k = \left\{({X^*}'{X^*})^{-1}\right\}_{k,:}{X^*}'y$$
Since $X^* \neq X$, then there is no reason for $\beta_k$ and $\beta^*_k$ to be equal.
Notice, however, that it can still be true in a specific case.
If the added variables are orthogonal the original independent variables, then $(\beta_k = \beta^*_k | k\leq p)$.
This result can be achieved through:
$${X^*}'{X^*}=\left[
\matrix{
{X}'{X} & X' \tilde X \\
\tilde X' X & \tilde X' \tilde X
}
\right]$$
Rows from the inverse of this covariance matrix left multiply ${X^*}'$.
If $X' \tilde X = \tilde X' X = 0$ the covariance become a block matrix, and it can thus be inverted blockwise:
$$({X^*}'{X^*})^{-1}=\left[
\matrix{
{X}'{X} & \color{red}{0} \\
\color{red}{0} & \tilde X' \tilde X
}
\right]^{-1}=
\left[
\matrix{
({X}'{X})^{-1} & 0 \\
0 & (\tilde X' \tilde X) ^{-1}
}
\right]$$
The full solution becomes:
$$
\begin{align}
\beta^* &= ({X^*}'{X^*})^{-1}{X^*}'y=\\
&=
\left[
\matrix{
({X}'{X})^{-1} & 0 \\
0 & (\tilde X' \tilde X) ^{-1}
}
\right]
\left[\matrix{X'y \\ {\tilde X'y}}\right]=\\
&=\left[\matrix{({X}'{X})^{-1}X'y \\ {(\tilde X' \tilde X) ^{-1}\tilde X'y}}\right]=\left[\matrix{\beta \\ \tilde \beta}\right]
\end{align}
$$
Thus keeping the identity between both results, as the entries pertaining to $\tilde X$ do not affect the coefficients pertaining to the original $X$.
Since coefficients change, so do CIs and p-values. A more in-depth look into how express things in terms of the hat matrix will lead to it all as well.
A: Coefficient change
Let some there be some data distributed according to a quadratic curve:
$$y \sim \mathcal{N}(\mu = a+bx+cx^2, \sigma^2 = 10^{-3})$$
For instance with $x \sim \mathcal{U}(0,1)$ and $a=0.2$, $b=0$ and $c=1$. Then a linear curve and a polynomial curve will have very different coefficients for the linear term.

set.seed(1)
x <- runif(100, 0,1)
y <- rnorm(100, mean = 0.2+0*x+1*x^2,
                sd = 10^-1.5)
plot(x,y, ylim = c(0,1.5),
     pch = 21, col = 1 , bg = 1, cex = 0.7)

mod1 <- lm(y~x)
mod2 <- lm(y~poly(x,2, raw =TRUE))

xs <- seq(0,10,0.01)
lines(xs,predict(mod1,newdata = list(x = xs)), lty = 2)
lines(xs,predict(mod2,newdata = list(x = xs)),lty =1)

legend(0,1.5,c("y = 0.009 + 1.023 x", "y = 0.193 + 0.016 x + 0.994 x^2"), lty = c(2,1))


Correlation
The reason is that the variables/regressors $x$ and $x^2$ correlate.
The coefficient estimates computed with a linear regression are not a simple correlation (perpendicular projection onto each regressor seperately):
$$\hat{\beta} \neq \alpha = \mathbf{X^t} y$$
(this would give coefficients $\alpha_1$ and $\alpha_2$ in the image below, and these coordinates/coefficients/correlations do not change when you add or remove other regressors)
Using the correlation/projection $\mathbf{X^t}y$ is wrong, because if there is a correlation between the vectors in $\mathbf{X}$, then there will be an overlap between some vectors. This part that overlaps will be redundant and added too much. The predicted value $\hat{y} = \alpha \mathbf{X}$ would be too large.
For this reason there is a correction with a term $(\mathbf{X^t}\mathbf{X})^{-1}$ that accounts for the overlap/correlation between the regressors. This might be clear in the image below which stems from this question: Intuition behind $(X^TX)^{-1}$ in closed form of w in Linear Regression

Intuitive view
So the regressors $x$ and $x^2$ both correlate with the data $y$ and they both will be able to express the variation in the dependent data. But when we use them together then we are not gonna add them according to their single independent effects (according to correlation with $y$) because that would be too much.
If we use both $x$ and $x^2$ in the regression then obviously the coefficient for the linear term $x$ should be very small since this is the same in the true relation.
However, when we are not using the quadratic term $x^2$ in the regression (or otherwise add a bias to the coefficient for the quadratic term), then the coefficient for $x$ which correlates somewhat with $x^2$ will partly take correct this (take over) and... the value of the estimate for the coefficient of the linear term will change.
Standard error change (and confidence intervals and p-values)
The errors of the variables may be correlated leading to very large errors in some coefficient when they strongly correlate with others. The matrix $(X^TX)^{−1}$
describes this correlation.
Error in the regression line
The image below shows intuitively how this changes when adding other regressors.
The intercept is the point where a regression line crosses $x=0$.

*

*On the left the error of the intercept is the error of the mean of the population.

*On the right the error of the intercept is the error of the regression line intercept.


Confidence regions for correlated parameters
The next image displays the confidence regions (contrasting with confidence intervals) of the above regression in a 2-D plot. Here it takes into account the correlation between the parameters.
The ellipse shows the confidence region which is a based on a multivariate distribution of the slope and intercept which may be related via a correlation matrix. For illustration an alternative type of region is also show. This is depicted by the box which is based on two single variate distributions assuming independence (now the confidence for the single variables is $\sqrt{0.95}$).
By changing the model from $y = a + bx$ to a shifted model $y = a + b(x-35.5)$ we see that the correlation between the slope and intercept changes. Now the "intercept" coincides with the standard error of the line around the point $x=35.5$ which you see in the image above is smaller.

#used model and data
set.seed(1)

xt <- seq(0,40,0.1)
x <- c(1:10)+30
y <- 10+0.5*x+rnorm(10,0,3)

See also:

*

*why does the same variable have a different slope when incorporated into a linear model with multiple x variables


*regression with multiple independent variables vs multiple regressions with one independent variable


*Why is the intercept in multiple regression changing when including/excluding regressors?


*Why and how does adding an interaction term affects the confidence interval of a main effect?


*Why is the intercept changing in a logistic regression when all predictors are standardized?


*Intuition behind $(X^TX)^{-1}$ in closed form of w in Linear Regression


*Does adding more variables into a multivariable regression change coefficients of existing variables?


*Estimating $b_1 x_1+b_2 x_2$ instead of $b_1 x_1+b_2 x_2+b_3x_3$


*Why do regression coefficients change when excluding variables?


*Does the order of explanatory variables matter when calculating their regression coefficients?


*Intercept changing after adding an interaction
