Violation of multiple linear regression assumptions after adding interactions

Question

I am working with a data set and properly fitted a model that satisfied the assumptions of lienar regression in a multi-variate setting. I have 8 predictors, and proceeded to test for significance of interaction using the anova function(comparing the original model to the one including the interaction term). After doing this I got results that stated, that some interactions are significant and should be included in the model. However now after updating the model I get that the model is no longer normal.

1. Is there a reason that this might be occurring?
2. Is there a method to zero in on what the cause of the problem is (perhaps an unwarranted interaction)?
3. Is there a separate way to test for interaction apart from anova comparison?

#my model
model23<-lm(strength2~blast2+flyash2+water2+superplast2+coarseagg2+fineagg2+age2)

I ran interactions between the predictors and included these interactions in the model and test against my initial model. The interaction models looked as follows:

int1<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+blast2*flyash)
int2<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+blast2*water2)
int3<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+blast2*superplast2)
int4<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+blast2*coarseagg2)
int5<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+blast2*age2)
int6<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+blast2*fineagg2)

int7<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+flyash2*water2)
int8<-lm(strength2 ~blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+flyash2*superplast2)
int9<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+flyash2*coarseagg2)
int10<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+flyash2*age2)
int11<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+flyash2*fineagg2)


int12<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+water2*superplast2)
int13<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+water2*coarseagg2)
int14<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+water2*age2)
int15<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+water2*fineagg2)



int16<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+superplast2*coarseagg2)
int17<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+superplast2*age2)
int18<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+superplast2*fineagg2)


int19<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+coarseagg2*age2)
int20<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+coarseagg2*fineagg2)

int21<- lm(strength2 ~blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+age2*fineagg2)

Then I ran the anova table against the original model letting #y indicate that we should include the interaction

anova(model23,int1) #N
anova(model23,int2) #N
anova(model23,int3) #Y
anova(model23,int4) #Y
anova(model23,int5) #Y
anova(model23,int6) #Y
anova(model23,int7) #y
anova(model23,int8) #Y
anova(model23,int9) #N
anova(model23,int10) #y
anova(model23,int11) #Y
anova(model23,int12) #y
anova(model23,int13) #y
anova(model23,int14) #Y
anova(model23,int15) #N
anova(model23,int16) #N
anova(model23,int17) #Y
anova(model23,int18) #N
anova(model23,int19) #Y
anova(model23,int20) #N
anova(model23,int21) #Y

Now I updated my model with the interaction terms:

model232<-lm(strength2 ~ blast2+flyash2+water2+superplast2+coarseagg2+age2+fineagg2+blast2*superplast2+blast2*coarseagg2+blast2*age2+blast2*fineagg2+flyash2*water2+flyash2*superplast2+flyash2*age2+flyash2*fineagg2+water2*superplast2+water2*coarseagg2+water2*age2+superplast2*age2+coarseagg2*age2+age2*fineagg2)

Before apply the interactions I ran a shapiro test, fitted vs residual plot,q-qplot,etc; and saw that the assumptions of linearity were satisfied. However if I run a shapiro test on the new model I get

shapriro.test(resid(model232)) #running this code



Shapiro-Wilk normality test

data:  resid(model232)
W = 0.99583, p-value = 0.00686

The plots of this model also showcase lack of satisfying assumption of linear regression.

Answer 1

There seem to be two things going on here.

First, with large data sets normality tests tend to fail. Real data seldom have exactly normal distributions. This page goes into extensive discussion. That doesn't mean normality tests are unimportant, but they have to be interpreted in terms of whether the deviations from normality are sufficiently large to invalidate your model for its intended uses. For example, normality of residuals isn't needed for least-squares regression to give the best linear unbiased estimates of regression coefficients. At most, normality of residuals is assumed by standard parametric tests of coefficient statistical significance. Even without strict normality of residuals the central limit theorem can still come to the rescue.

Second, you should rethink your approach toward evaluating which interactions to include. You have compared a large set of pairwise interactions individually against the model without interactions. With 21 separate tests for particular interaction pairs, you have (with a p < 0.05 cutoff) about a 66% chance of at least one false-positive finding. Then you included all of the individually identified interaction terms together into a model without regard to how much that combination of interaction terms is really helping. As your predictors are likely to be correlated I suspect that there is a lot of redundancy among the included interactions.

How to proceed depends on how you intend to use your model. The model without interactions might give adequate performance for some applications.

If the interactions themselves are of interest, you seem to have enough data to include all 28 possible pairwise interactions plus the 8 primary predictors in a single model. You would have nearly 30 times as many data points as the number of parameter values that model would estimate, typically enough to avoid overfitting. Then, depending on the purpose of your modeling, you could use that entire model, keep all interactions but with a penalized approach like ridge regression (to get a model potentially less likely to overfit and thus more likely to work for predicting new cases), or pull back carefully from the full model to remove the least useful interactions while all the remaining interactions are still taken into account.

A final warning: the names of your variables suggest that these data might be coming from a time series. If that's the case there are inherent problems with non-independence among observations that require different approaches than simple linear regression.