I believe effects like these are frequently caused by collinearity (see [this question][1]). I think the book on multilevel modeling by Gelman and Hill talks about it. The problem is that `IV1` is correlated with one or more of the other predictors, and when they are all included in the model, their estimation becomes erratic. If the coefficient flipping is due to collinearity, then it's not really interesting to report, because it's not due to the relationship between your predictors to the outcome, but really due to the relationship between predictors. What I've seen suggested to resolve this problem is residualization. First, you fit a model for `IV2 ~ IV1`, then take the residuals of that model as `rIV2`. If all of your variables are correlated, you should really residualize all of them. You may choose do to so like this rIV2 <- resid(IV2 ~ IV1) rIV3 <- resid(IV3 ~ IV1 + rIV2) rIV4 <- resid(IV4 ~ IV1 + rIV2 + rIV3) Now, fit the final model with DV ~ IV1 + rIV2 + rIV3 + rIV4 Now, the coefficient for `rIV2` represents the independent effect of `IV2` given its correlation with `IV1`. I've heard you won't get the same result if you residualized in a different order, and that choosing the residualization order is really a judgment call within your research. [1]: https://stats.stackexchange.com/questions/1149/is-there-an-intuitive-explanation-why-multicollinearity-is-a-problem-in-linear-re