I would like to regress total energy expenditure on age, gender, height and weight.But how do I check whether this relationship is linear?
Even if that relationship is not linear, a regression model provides an estimate of the "first order trend" or averaged linear response. So we can still extract useful information from a linear model set to curvilinear data. This is a common approach in inference and inferring associations. As another point, I am actually very surprised exactly how effectively linear trends usually model data.
For the purposes of prediction, modeling the trend as closely as possible is important. When doing model selection and we are concerned about curvilinear trends, we may fit higher order effects and structurally test models for parsimony. Higher order effects include quadratic terms. This usually can provide a reasonable approximation for slightly exponential or logarithmic trends. Cubic terms can be better at capturing S shaped trends.
A more general form of doing this is fitting splines with cross validation.
Lastly, a cautionary note, bivariate scatter plots of regressors against the outcome of interest are not suited to demonstrate whether multivariate linearity holds. Thinking in 3 or 4 dimensions, you can visualize linear planes in which data are distributed in such a way that their projection appears nonlinear. A conditioning plot is a very nice way to take care of that and display multidimensional data using bivariate scatter plots.
You could simply look at some sort of multiple nonparametric smoothing regression (e.g. lowess, running line, splines, GAM, etc.) to directly assess whether linearity is a good assumption (as opposed to trying to do it through residuals). An added depth to the interpretability of this approach is that any non-linear realtionships you view are often suggestive of a specific functional form (e.g. a saturation or threshold hinge-function, quadratic relationship, etc.) which you might then specify in a subsequent parametric regression (e.g. nonlinear least squares).
One way to assess the suitability of the linearity assumptions for each variable would be to look at added variable (partial regression) plots for each variable; this removes the effect of all the other variables (as long as they're correctly specified).
This will identify any one variable being nonlinearly related to $y$ (though generally you'll still see problems even if several relationships are curved).
Here's an example - (done using
Alternatively (and perhaps better) there's partial residual or component+residual plots:
If there's a chance multiple variables may be non-linear, you might be better fitting an additive model ($y=\alpha+s_1(x_1)+s_2(x_2)+...+s_p(x_p)+\varepsilon$) and checking if each one is reasonably well approximated by a line.