Trouble with finding equation for predicting y based on data provided

Question

I am being asked to produce an equation for some characteristic $y$ based on measurements provided for explanatory variables $x_1,x_2,x_3,x_4$. The first thing I did was (naively) fit the linear regression, $y=\beta_0+\beta_1x_1 +\beta_2x_2 +\beta_3x_3 +\beta_4x_4 + \varepsilon$. I got the following results:

$$ \begin{array}{|l|c|} \hline & y \\ \hline x_1 & 0.227^{**} \\ & (0.100) \\ x_2 & 0.554 \\ & (0.370) \\ x_3 & -0.150^{***} \\ & (0.029) \\ x_4 & 0.155^{***} \\ & (0.006) \\ \text{const.} & -6.821 \\ & (10.123) \\ \hline \text{Observations} & 32 \\ \text{R}^{2} & 0.962 \\ \text{Adjusted R}^{2} & 0.957 \\ \text{Residual Std. Error} & 2.234 (\mathrm{df} = 27) \\ \text{F Statistic} & 171.713^{***} (\mathrm{df} = 4; 27) \\ \hline \hline \end{array} $$

Overall decent results. But when I check the fitted values against residuals I get the quadratic-ish relation.

Since the residuals are not randomly distributed and the trend reminds some quadratic function, it means that the way I assumed linearity at the very beginning was not actually correct.

In the previous example, changing $y$ to $log(y)$ solved the problem. However, here when I tried a similar approach I ended up with a switched direction of the curvature.

My question is what else can be done in order to find the "correct" non-linear model? Should I try squaring and logging each term to see how it affects the model? Or is there a more concise way of doing so? What should I do if I cannot find a combination that would solve this problem?

---

EDIT.

I have added the matrix of scatterplots so to show the relationship between each of the variables. It is hard for me to observe anything specific, to be honest. The only thing that is somehow transparent to me is the positive relation between $y$ and $x_4$. What I also find quite interesting is the relation between $y$ and $x_1$ but not entirely sure what that tells me.

Answer 1

It is important not to over-interpret these plots. The first plot of residuals vs fitted values is a little misleading in my opinion if you only focus on the red line, partly due to the fairly small sample size.

Yes, the red line has a curved shape, but looking at the data points, it is not clear at all that there is nonlinearity. This type of pattern can easily occur simply through random variation. We would like the line to be perfectly horizontal, but in practice this will never happen.

Here is a plot from a very simple simulation of a linear model:

set.seed(2)
X <- rnorm(30)
Y <- 10 + X + rnorm(30)

plot(lm(Y~X))

As you can see, even when sampling from a normal distribution with a linear relationship, the plot of residuals vs fitted values can still appear nonlinear if we only focus on the red line.

If you have reason to believe that all the $x$ variables are related to $y$ and there is no interdependence among the $x$s (as it appears from the corelation plots), then I would stick with the fist model.