How to balance/interpret increasing "good fit" statistics, but seeing more departing from linear regression assumptions?

Question

I noticed that I can substantially improve my model's R$^2$ and Residual standard error values by adding some interaction terms.

My model's statistics go from:

Residual standard error: 52.06 on 790 degrees of freedom
Multiple R-squared:  0.5862,    Adjusted R-squared:  0.4773 
F-statistic: 5.381 on 208 and 790 DF,  p-value: < 0.00000000000000022

to

Residual standard error: 45.22 on 106 degrees of freedom
Multiple R-squared:  0.9581,    Adjusted R-squared:  0.6057 
F-statistic: 2.718 on 892 and 106 DF,  p-value: 0.0000000009271

but look what happens to the plots:

from:

to:

tells me that adding the interaction terms is causing something very bad to the model. Or could this signify that actually my data set is e.g. highly non-linear, which is why I would get a "better fit" (according to the summaries above), eventhough the plots start to show more non-linear behaviour?

The reason for believing in "non-linearity" is that:

In linear models I've heard there to be an assumption of the errors having to be normally distributed, that would mean that the QQ-plot should be as close to the straight line as possible. Right? So my thinking of "non-linearity" came from thinking that "linearity assumptions" are violated.

Also, does violating the linearity assumptions reflect to the predictions made on "new data", but not on the data (observations) that's used for the fit? So the reason to take linearity assumptions seriously would be so that the model wouldn't be dependent on the original data set and can work in predicting "future observations"?

How should I balance and interpret these results and what model should I continue to develop, the first one or the second one?

Answer 1

A non-straight Q-Q plot indicates non-normal residuals, not "nonlinearity in the model". That may be due to a variety of model issues. In this case the reason for the flat sections in the Q-Q plot seem clear enough.

The hundreds of interaction terms in the model (684 df worth) result in you getting exact or nearly exact fits to a substantial fraction of your data. The resulting zero residuals lead to the "line" across the residuals vs fitted and to the large flat part in the normal Q-Q plot.

This could occur (for example) when you have many factor combinations corresponding to a single observation (which would then be fitted exactly). Smaller flat sections might correspond to factor combinations with two observations (or even 3 or 4 if the response is discrete).

It doesn't necessarily suggest anything is nonlinear; it might be purely to the enormous number of interaction terms.

If you want to add all of those terms I'd suggest considering some form of regularization.