I have a model $Y = \alpha + \beta_1X_1 + \beta_2X_2$.
$Y$ has a bimodal normal(ish) distribution, so I'm looking to see if the relationship between the predictors and the response is different for low values of $Y$ compared with high values of $Y$.
Thus far I've recoded $Y$ as a dummy variable $D$ (0 if Y < median(Y), else 1) and added it as a predictor yielding $Y = \alpha + \beta_1X_1 + \beta_2X_2 + \beta_3D$ (plus interactions).
$D$ is significant as both a main effect (obviously) and in interaction with $X_1$. This suggests to me that, indeed, the relationship between $X_1$ and $Y$ differs for high values of $Y$ compared with low values of $Y$.
But somehow this method feels off.
Adding the dummy increased $R^2$ from 20% to 80% (caused by the fact that the dummy is derived from $Y$). I'm getting a large amount of multicollinearity between my dummy and the dummy in interaction with $X_1$ and $X_2$, which, again makes sense, but signals to me this method is not what I'm looking for.
I can't find any useful information on recoding a response variable as a dummy and adding it as a predictor (again, for me, that's a cause for concern), so I'm hoping someone here will point out the flaw and explain why this method is not cricket.
EDIT: For clarification, I'm trying to avoid splitting my data and running the analysis separately for each mode.