# 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$$?

• I assume by "feature importance" in the linear case, you assume the size of the coefficients. However, even in this case this is not the case. perhaps try to identify how to define feature importance, then consider taking non linear combinations of $x_i$ in linear model... Commented Jan 20, 2022 at 0:11
• I've tried to clarify my question and removed the vague "feature importance" :)
– pir
Commented Jan 23, 2022 at 0:30
• Instead of the controlled directed effect, you may also consider the natural direct effect of $x_1$,...,$x_n$ on $y$, taking into account that $x_c$ is a mediating variable. There are several methods for investigating this direct effect (as opposed to the total effect, which is computed by ignoring the influence of $x_c$) that also work in the nonlinear setting, e.g., see "Pearl’s more general approach" in bmcmedresmethodol.biomedcentral.com/articles/10.1186/… Commented Jan 28, 2022 at 7:42
• This question is really too broad and vague to be suitable here: calling something "nonlinear" is like describing a shape as "not square:" it covers an awful lot of territory! Could you instead give us some details of your specific problem?
– whuber
Commented Jan 28, 2022 at 17:38
• This is a problem that I’ve faced with multiple modeling tasks so I was hoping to find a generalized solution. Sounds like that might not be possible, which is reasonable.
– pir
Commented Jan 28, 2022 at 17:40

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.

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).

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!

• I added my own answer as I think this approach works, but I'm not 100% sure and would love input. However, I'd also really appreciate other completely unrelated answers on how to tackle this modeling problem.
– pir
Commented Jan 19, 2022 at 22:18
• I think it would be better to post this in your question as a proposed solution, rather than giving it an air of authority by posting an answer.
– Dave
Commented Jan 19, 2022 at 22:18
• I was initially planning to do that, but was fearing that the answers might then spend a bunch of time explaining why this approach works / doesn't work. I'd worry that it'd disencourage others from answering with completely unrelated approaches to this. Perhaps it'd be fine to do this if I make the warning even more clear?
– pir
Commented Jan 19, 2022 at 22:21