# predict(lm) with two interactions possible or incorrect?

in order to understand the output of one of my lme models I produced a little simpler example using lm (so no random factor). I noticed that my fitted model does not seem to fit the data correctly, as the predicted y-values deviate strongly from the given values when using predict(lm).

My data set is:

a = c(rep(1:10, 4))
b = c(10,20,30,40,50,60,70,80,90,100, 5,8,10,14,17,22,27,35,42,50, 90,82,73,64,56,48,40,33,25,18, 5,6,8,10,12,14,17,20,23,26)
c = c(rep("male", 10), rep("female", 10), rep("male", 10), rep("female", 10))
d = c(rep("low", 20), rep("high", 20))
e = data.frame(yval = b, xval = a, sex = c, education = d)


Graphically it looks like this:

library(car)
scatterplot(yval~xval | education, smooth = F, grid = T, spread = F, reg.line = T, data = e, xlab = "x", ylab = "y")
scatterplot(yval~xval | sex, smooth = F, grid = T, spread = F, reg.line = T, data = e, xlab = "x", ylab = "y")


and the linear model is:

lm2 = lm(yval~xval+sex+xval:sex+education+xval:education, data = e)
summary(lm2)


Using

e_pred = e
e_pred$pred = predict(lm2)  I get the predicted values which do not match the real data at all: e = cbind(e,e_pred$pred)


Is this due to the fact, that there is more than one significant interaction?

There is something very wrong with that model in terms of the fit to the actual data. It is always a good idea to look at plots of aspects of the model fit, particularly in this case for patterns in residuals and influential observations.

To look at several diagnostics plots, do

layout(matrix(1:4, ncol = 2))
plot(lm2)
layout(1)


This produces the plot below

The clear patterns in the residuals suggests a missing covariate or term from the model. As the residuals aren't independent, that calls into question the inferences in the table produced by summary(). Note also the lower right plot, where most of your observations lie outside the the red (dashed) lines indicating strong influence and leverage on the model coefficients.

Now we need to understand what might be causing this. If we reuse the scatterplot code but with resid(lm2) as the y-axis variable we can examine how the residuals vary with xval conditional upon the two factors

scatterplot(resid(lm2) ~ xval | education, smooth = F, grid = T, spread = F,
reg.line = T, data = e)
scatterplot(resid(lm2) ~ xval | sex, smooth = F, grid = T, spread = F,
reg.line = T, data = e)


The two plots are shown below, first conditioned by education

and next by sex

This shows that the effect of xval on yval varies considerably (change of slope) depending which sets of observations we look at. If you were to plot the residuals vs sex or vs education (as boxplots, not shown here) they don't look too bad, so the issue looks to be with xval.

As you mention you are doing this to look at the results of a mixed model, one assumes the observations are clustered by subject. This could account for the apparent differences in the slope of the relationship for xval on yval, where different individuals had different intercepts and slopes.

To that end you should draw plots of the residuals conditioned on Subject.

If there is no Subject effect in the data here, then the pattern of residuals suggests that there is another, unmeasured (or not provided here) factor that is modifying the relationship between yval and xval.

But to answer the specific question, there are effects in the data that are not accounted for by the simple linear model fitted you fitted. One option to follow up upon is whether there are missing higher order interactions. For example, if we allow for 3-way interactions of the variables in the model, a strong 3-way interaction is observed and the xval term is now significant.

lm4 <- lm(yval ~ (xval + sex + education)^3, data = e)

> summary(lm4)

Call:
lm(formula = yval ~ (xval + sex + education)^3, data = e)

Residuals:
Min      1Q  Median      3Q     Max
-3.5394 -0.7318  0.0000  0.3379  4.8545

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)                 1.0667     1.2218   0.873  0.38915
xval                        2.3697     0.1969  12.034 2.03e-13 ***
sexmale                    96.0667     1.7279  55.598  < 2e-16 ***
educationlow               -5.1333     1.7279  -2.971  0.00559 **
xval:sexmale              -10.4121     0.2785 -37.390  < 2e-16 ***
xval:educationlow           2.5515     0.2785   9.162 1.84e-10 ***
sexmale:educationlow      -92.0000     2.4436 -37.649  < 2e-16 ***
xval:sexmale:educationlow  15.4909     0.3938  39.335  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.789 on 32 degrees of freedom
Multiple R-squared:  0.9965,    Adjusted R-squared:  0.9957
F-statistic:  1295 on 7 and 32 DF,  p-value: < 2.2e-16


The data are almost perfectly fitted now, all terms in the model are significant and the model diagnostics plots look a lot better, thought there are still some influential observations and the distribution of the residuals is alittle heavier in the tails:

layout(matrix(1:4, ncol = 2))
plot(lm4)
layout(1)


which produces

So not perfect, but perhaps illustrates the issue of missing predictors or terms from the original model quite nicely.