I am working on a problem where I have a "target" variable $Z$ that I know for sure is influenced by a "predictor" $Y$. I also have a second predictor $X$ that is correlated with $Y$ (about -.3), and I would like to know whether is has an influence on $Z$
If I fit a linear model for $Z$ on $X$, I get a very significatively non-zero coefficient for the slope, but of course, this tells nothing, as any variable that is correlated with $Y$ would do so.
Next thing I tried was to fit a linear model for $Z$ on both $X$ and $Y$. The coefficients for both variables is significant, but, since the variables are correlated, I am not sure about the power of the model to "distinguish" between the effects of each of the predictors.
So I finally went for a third approach: I will build the model in two steps: First, I will fit a model for $Z$ on $Y$ to eliminate its effects. Then, I'll take its residuals $\epsilon_i$ and fit a second linear model for $\epsilon$ on $X$
Now my question is: Is this a valid method to achieve what I actually want to achieve? Also, should I add more complex terms in the "first-step" model, to eliminate other non-linear effects of $Y$ on $Z$?
NOTE: When I say "significant", I don't mean one of those $0.04$ $p$-values that fill mediocre research papers. I mean something more in the roum of $10^{-10}$, as sample and effect sizes are both quite big