Intercept changing after adding an interaction Quite a basic multiple linear regression question.
I know that the intercept should describe the predicted value when all predictors equal 0. Then, why does adding a paired interaction term (X1X2) changes the intercept? Wouldn't the product term be zero when the predictors are zero?
Thanks for your help.

Consider a simple case with only a single numerical predictor and two groups. You are right that the intercept gives the fitted value when all predictors are zero, that is, when the numerical predictor has a value of zero and the group is the reference group (assuming standard dummy coding for factor variables).
If you plot your fits, the intercept then corresponds to the point where the fit for the reference group, plotted against the predictor, intersects the vertical axis. However, here is the key point: when you add an interaction to your model, you change the slopes and the intercepts of both fitting lines. Thus the $y$ axis intercepts of the two (now non-parallel) lines will change - and the "model" intercept, i.e., the constant term, will do so too, since it is just one of the two $y$ axis intercepts.

In the plot above, the intercept for the model plotted on the left is 0.3370 and for the one on the right, it is -0.0019. In both cases, it is the $y$ coordinate where the black dashed line intersects the $y$ axis.
R code:
set.seed(1)
nn <- 100
group <- sample(c("A","B"),size=nn,replace=TRUE)
predictor <- runif(nn)
dataset <- data.frame(group,predictor,
    outcome=predictor-1.4*predictor*(group=="B")+0.5*(group=="B")+rnorm(nn,0,0.2))
xx <- seq(0,1,0.01) # for plotting

opar <- par(mfrow=c(1,2),mai=c(1,1,.7,.1))
with(dataset,plot(predictor,outcome,col=(group=="B")+1,pch=19,las=1,main="Main effects only"))
legend("topleft",pch=19,col=1:2,legend=c("A","B"))
(model_1 <- lm(outcome~predictor+group,dataset))
lines(xx,predict(model_1,newdata=data.frame(group="A",predictor=xx)),col=1,lwd=2,lty=2)
lines(xx,predict(model_1,newdata=data.frame(group="B",predictor=xx)),col=2,lwd=2,lty=2)

with(dataset,plot(predictor,outcome,col=(group=="B")+1,pch=19,las=1,main="Model with interaction"))
legend("topleft",pch=19,col=1:2,legend=c("A","B"))
(model_2 <- lm(outcome~group*predictor,dataset))
lines(xx,predict(model_2,newdata=data.frame(group="A",predictor=xx)),col=1,lwd=2,lty=2)
lines(xx,predict(model_2,newdata=data.frame(group="B",predictor=xx)),col=2,lwd=2,lty=2)
par(opar)

