Controlling for non-linear variable in non-linear modeling of response I need to model a continuous response variable $y$ based on continuous features $x_1, ..., x_n$ while controlling for another continuous feature $x_c$. The intent is to understand how much an increase of 1 in one of the features (e.g. $x_1$) would on average impact $y$ while controlling for $x_c$.
If the relationships were linear, then including $x_c$ in a normal linear regression would be enough to control for $x_c$, but that's not the case. The data has the following patterns that makes it much more difficult to model:

*

*The relationship between $y$ and $x_c$ is highly non-linear

*The relationship between $y$ and $x_1, ..., x_n$ is highly non-linear

*Some of the features in $x_1, ..., x_n$ are highly correlated with $x_c$
How would I best model $y$ using features $x_1, ..., x_n$ while controlling for $x_c$?
 A: All you've really told us is that the relationships are nonlinear, and you are asking us to tell you how to get the best model.  It's not really possible to do that, since "nonlinear" is a huge class of relationships.  Nevertheless, here are some things to think about:

*

*Ideally, there would be some aspect of the variables that would suggest a class of nonlinear functions that might be plausible for the regression relationships.  Sometimes the nature of the variables and the context from which they are obtained give us some hints on the types of nonlinear relationships they might have.


*If this is not the case, you might need to proceed on a purely empirical basis.  For non-periodic functions you can usually approximate them locally by polynomials (based on the theory of Taylor series), and for periodic functions you can usually approximate them locally by sums of sinusoidal functions (based on the theory of Fourier series).  These approximations work fairly well on a wide class of nonlinear functions, so they are commonly used as a starting point in nonlinear regression in cases where there is no prior information about the nature of the nonlinear relationship.


*Once you have formulated a reasonable starting point for your model, fit the model and produce its diagnostic plots, paying particular attention to the added variable plots.  These latter plots will allow you to scrutinise the chosen regression relationshops in your model based on the data.  This can alert you to erroneous assumptions in your model, and allow you to choose a more appropriate regression function to fit your data.  (But there is a limit to this; beware of the phenomenon of over-fitting.)


*Finally, the issue of multicollinearity of the explanatory variables is a common issue in prediction.  Unless you are conducting a controlled trial and you have control over your explanatory variables, there is nothing much you can do about this, except to be aware of the effect that multicollinearity has on the accuracy of your regression estimates, and to make sure that you are clear on what type of predictive inferences you want to make (i.e., exactly what do you want to condition on).
A: 
The intent is to understand how much an increase of 1 in one of the features (e.g. $_1$) would on average impact $y$ while controlling for $x_c$. (Emphasis added.)

First, given nonlinear associations of predictors with outcome, there isn't a unique answer. You have to specify a particular value of $x_1$ from which to evaluate the change in $y$ or a range of $x_1$ values over which you would average. If the nonlinearities involve interactions with other predictors, you would need to specify the levels of the interacting predictors too. Keep that in mind.
Second, nonlinear associations of predictors with outcome can often be analyzed empirically with a linear regression model if you have no theoretical model in mind. A particular form of polynomial approximation, restricted cubic regression splines, is a common choice. The regression is then still linear in the coefficients, so once the general form of the spline is specified (via methods in standard statistical software) linear regression fitting is all that is required.
Chapter 2 of Frank Harrell's course notes outlines that approach to modeling nonlinear relationships among variables (Section 2.4), including how to evaluate model fit and handle interactions among such predictors (Section 2.7). There are related approaches with penalized splines and generalized additive models, discussed in this thread.
Finally, as the comments indicate a potential interest in "feature importance," see Section 5.4 of Harrell's notes. The anova() function in his rms package can provide a measure of predictor importance that combines all nonlinear and interaction terms for a predictor, the difference between the partial $\chi^2$ for the predictor and the number of degrees of freedom. He uses analysis of multiple bootstrap samples to illustrate how unreliable such a measure can be.
A: You could train a non-linear model using only $x_c$ to predict $y$, i.e. $f(x_c) = \hat{y}$. Then compute the residual $\epsilon = \hat{y} - y$. The residual $\epsilon$ will have been controlled for with regards to $x_c$ as $f$ has captured the relationship between $x_c$ and $y$ as well as it possibly could.
You can now use the residual $\epsilon$ as input for all further analysis such as any modeling or feature importance algorithms. You’ll want to include $x_c$ as a feature in your further analysis in case of any interaction effects between $x_c$ and $x_1, …, x_n$.
In terms of the choice of the non-linear model for $f$, random forest is a good and easy-to-use method, but many other non-linear models (e.g. neural networks) could be used as well.
Warning: I think this approach works, but would appreciate references to confirm it!
