regressing out confounders from negative binomial data I have a very large dataset of count data. I know that these count data are discrete and follow a negative binomial distribution.
I want to control the effects of two confounders. To do this, I plan to use a regression model adjusting for these confounders, then perform other analyses with the residuals.
To "regress out" these two confounders, would I use a linear model or GLM (poisson/negative binomial)? Are there other choices? And in what ways would I identify a single modeling approach to be better than the others?
 A: This is not the correct approach.
A "confounding" relationship is a three way relationship. There is an exposure ($X$), what is also called a "predictor", "regressor", or an "independent variable". Then there is an outcome ($Y$), equivalently "response", or "dependent variable". The confounder ($W$) is the third element which is causal of both the exposure and the outcome. 
While you haven't mentioned an exposure, we have to assume it's there. If not, "confounder" is the wrong term.
Your suggestion is first to fit a model:
$$ E[Y|W] = \gamma_0 + \gamma_1 W$$ 
Then use the residuals $r_w = Y - \hat{\gamma}_0 - \hat{\gamma}_1 W$ to model the association 
$$E[Y | X, W] = \beta_0 + \beta_1 X + \text{offset}(\hat{\gamma}_0 - \hat{\gamma}_1 W)$$.
It turns out this approach will be biased because the marginal effect of $W$ is different from the conditional one. The first regression model $Y$ on $W$ models the total effect of $W$ on $Y$ which is its direct effect, but also the indirect effect of $W$ on $X$ then $X$ on $Y$. The resulting inference for the association between $X$ and $Y$ will be biased toward the null if a) $W$ and $X$ have the same direction of effect on $Y$ and they are positively correlated, b) away from the null if either direction is different or if $W$ and $X$ are negatively correlated, or c) biased again toward the null if they have opposite direction of effect and are negatively correlated. In all cases, inference is underpowered.
The correct approach is to use the full conditional model in one fell swoop:
$$E[Y |X,W] = \eta_0 + \eta_1 X + \eta_2 W$$
