How to deal with curvature in residuals plot

Question

I am trying to do a multiple linear regression in R but am having some problems. I have a set up where I am trying to develop a multiple linear regression model for one variable (y) using six other variables ($x_{1},...,x_{6}$), all of which are correlated to some degree.

Based on my understanding of the data I ran a multiple linear regression for y using $x_3$ and $x_5$. Here is an updated residual plot based on feedback:

$$y~x_3+x5+x_3*x_5$$

How can I fix this? Would feasible GLS work? Or have I selected a bad combination of independent variables?

By the way I have limited knowledge of regression and have never done a weighted regression before so this is new to me.

Thanks for your advice!

Edit: Here are the additional plots

Here are the residual plots for $x_1...x_6$ too:

Hope this helps!

Answer 1

It is sometimes helpful to plot the residuals against each $x_j$ variable individually (i.e., two plots in your initial situation, possibly up to 6 plots eventually). I am not sure if that will help you here, but it's worth a try.

Looking at your top lefthand plot ("Residuals vs. Fitted"), if you had a simple (i.e., 1 $x$ variable) regression, that picture would be a tell-tale sign that the relationship between $x_j$ and $y$ is curvilinear and you need to add a squared term ($x_j^2$). In this case, you have two variables, and that implies you may have an interaction. Try adding a new variable, $x_3\times x_5$, to your model. In R, the syntax would be: y~x3+x5+x3:x5.

To illustrate my point about a missing interaction, consider this simulation:

set.seed(9)
x1  = runif(200, min=0, max=100)
x2  = runif(200, min=0, max=100)
x12 = x1*x2
y   = 0 + 1*x1 + 1*x2 + .01*x1*x2 + rnorm(200, mean=0, sd=10)
mod = lm(y~x1+x2)

There is an interaction in the data generating process, but the model is mis-specified, it leaves the interaction out. Here is what the plots look like:

Note that a missing curvilinear term (x^2) could cause a similar picture. There is a little squirrelly curvature between y and x5, so a squared term might help there a little bit.

Answer 2

First of all get read of the outlier with too big cook distance and do the regression again.

Try to add x² in the regressor when y vs x does not look linear.

If the variance of your residuals is still increasing with y after that, perhaps an heteroscedastic regression would be better.

If it still does not work you can try a non-parametric regression.

Answer 3

If you don't want to abandon the lm setting because you're not familiar with regressions, my solution could help.

I've noticed you use R, so what I show you may be useful:

I suggest you car::residualPlots function to have a clearer idea of how the sample residuals behave against each predictor and against fitted values.
This function provides both the plots with the curve pattern and the tests for curvature of residuals. In your case, (some or all) tests will reject the NH of no curvature.
This does not imply that adding a quadratic term for each predictor will surely solve your problem, instead you'd rather think this tests as a proof that something is not fine. Then you can:

1) transform predictors (in this case car::powerTransform may provide useful hints);

2) add quadratic, cubic,... or interaction terms;

3) if something still looks bad, transform also the response (you should put your model formula inside car::powerTransform).

Then do the diagnostics (outliers and cook's distance, heteroschedasticity..);

if everything satisfies you, run the summary.

Answer 4

I'm not allowed to comment yet, but if I understood the variables correctly, then you are expecting an interaction term, y ~ x3 * x5, since benefits paid should be equal to avg benefit per week times number of weeks (or at least more equal to that than any linear combination).