# How can statistical inference be done when using multicollinear data?

I want to model a variable $$y$$, which is known to be affected by a set of 5 variables: {$$x_1$$, $$x_2$$, $$x_3$$, $$x_4$$, $$x_5$$}. These variables are known to be correlated with one another and $$y$$ to some extent. From the literature it also seems probable that some polynomial (e.g. $$x_1^2$$) and interaction effects (e.g. $$x_1x_2$$) may also be present, but I have no way of telling beforehand which (if any) of these 5 variables exhibits this behaviour.

What I want to do is to infer the effect of $$x_1$$ on $$y$$; specifically I want to know the coefficient and its standard error. However, with the potential polynomial and interaction terms from all 5 variables I don't think that a simple OLS fit will work here.

From what I understand the ideal model of $$y$$ should be both parsimonious and have a high prediction score in order for any information about $$x_1$$ can be derived. I've seen in the past that LASSO, Ridge, and Elastic Net regularisations have been proposed to perform feature selection on multicollinear data, but I can't calculate standard errors or p-values using these techniques because these are biased estimators. A lot of material online tackles multicollinearity from the perspective of prediction-only models, but I've seen little discussion on how inference can be performed using such data.

What should my approach here be? Can regularisation methods still be used for inference somehow? What are the best practices for statistical inference using multicollinear data?

If interactions are important, there is no single "effect of $$x_1$$ on $$y$$." Rather, that "effect" (by the definition of an "interaction") will depend on the values of all of the other predictors that interact with $$x_1$$. Although you can't get a single coefficient and standard error for that "effect," you have a couple of ways to get some estimate of the overall importance of $$x_1$$ in the model.
One is to do the OLS multiple regression, including all the interaction terms that your knowledge of the subject matter suggests are important (if you can do so without overfitting). You will get coefficients for $$x_1$$ itself (representing the situation with other predictors at reference values or at 0) and coefficients for all interaction terms involving $$x_1$$. You then evaluate the significance of the overall contribution of $$x_1$$ to the model. The anova() function in the R rms package provides one generally applicable test, a Wald test on all of the coefficients involving $$x_1$$ that takes into account the correlations among the coefficient estimates.
If you have absolutely no idea about what interactions to include, you could consider gradient-boosted regression with enough tree depth to allow for a reasonable number of interactions. Such models can report the importance of each predictor, based on how much each predictor improves the overall model quality. There still won't be a single coefficient for $$x_1$$, however, and there is no easily interpretable way to evaluate which interactions are important in those models.
For reference, there is work on inference with LASSO and related models; see Statistical Learning with Sparsity for a detailed explanation and an introduction to later literature. That might not, however, give you what you are looking for here. There will still be issues about predictor selection, what happens if LASSO selects an interaction involving $$x_1$$ but not the main $$x_1$$ coefficient, etc. See many pages on this site, for example this page on predictor selection and this page on LASSO and interaction terms.
• Thank you very much for your answer. I have some more questions, if you don't mind: 1) You mention interactions only when discussing anova and gradient-boosting. Would the polynomial terms require similar consideration as well? 2) Would it be theoretically correct to take the first-order derivative of my fitted model in order to estimate the overall effect of $x_1$ (i.e. $dy/dx_1$)? For instance, if an interaction term $x_1x_2$ had a fitted coefficient $a$, then the derivative would include the term: $ax_2$ Dec 1 '20 at 15:02
• @ElectronicAnt (1) polynomial terms are interactions of a continuous predictor with itself, so treat them accordingly; (2) as the values of those derivatives include the values of the interacting predictors, there still is no single "overall effect" of $x_1$ independent of the values of the interacting predictors. You can, however, do things like generate examples of "overall effects" (and associated standard errors) conditional on particular values of the interacting predictors. Be careful, however, not to lead yourself or your audience astray into thinking that it is the overall effect.