I believe effects like these are frequently caused by collinearity (see this question). 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.