When a predictor variable of one model needs to be adjusted for other measures before it is interpretable, is it better to fit a model for that predictor and then feed the residuals of that model into the outcome model? Or is it better to add the predictor-related measures into the outcome model?

As an example, say I am trying to investigate the relationship between whole blood WGS-derived mitochondrial DNA copy number and overall mortality.

Overall risk of mortality may need to be adjusted for some covariates, so I might write the model in R as something like this:

Model 1: coxph(Surv(Follow_up_time, Death) ~ Age + Sex + Smoking_Status + BMI + Education + Ethnicity + mtDNA_CN, DF)

However, because in this case mitochondrial DNA copy number is a tissue-specific relative value, it also needs to be adjusted for a number of technical covariates before it is interpretable, such as sex, age, white blood cell count, platelet count, blood cell subtype parameters, autosomal coverage, etc.

Model 2: lm(mtDNA_CN ~ Age + Sex + WBC_n + Platelet_n + Blood_Parameters + Autosomal_Coverage, DF)

Should I take the residuals of Model 2 and feed them into Model 1 instead of mtDNA_CN? e.g.

Model 3: coxph(Surv(Follow_up_time, Death) ~ Age + Sex + Smoking_Status + BMI + Education + Ethnicity + Model_2_Residuals, DF)

Or would it be better to simply add the technical covariates for mtDNA_CN into the survival model even though they may not themselves be related to mortality? e.g.

Model 4: coxph(Surv(Follow_up_time, Death) ~ Age + Sex + WBC_n + Platelet_n + Blood_Parameters + Autosomal_Coverage + Smoking_Status + BMI + Education + Ethnicity + mtDNA_CN, DF)

I have asked a few researchers I know and I've received conflicting answers. I have also seen similar questions asked in other threads, but was not able find one that answered this question specifically.


1 Answer 1


Model 4 provides the advantage of allowing for associations of WBC_n, Platelet_n, Blood_Parameters and Autosomal_Coverage with outcome beyond their association with mtDNA_CN (in Model 1) or its deviation from what would be predicted by Model 2 (in Model 3). In effect, Model 1 isn't substantively different from Model 3, except in interpretation. The disadvantage in moving to Model 4 is that you need to estimate more Cox coefficients and thus need a larger study size to avoid overfitting.

To illustrate, consider a set of simpler models only including Age, mtDNA_CN and WBC_n.

Then we have:

mod1a <- coxph(Surv(time, status) ~ Age + mtDNA_CN)
mod2a <- lm(mtDNA_CN ~ Age + WBC_n)
mod3a <- coxph(Surv(time, status) ~ Age + mod2resid)
mod4a <- coxph(Surv(time, status) ~ Age + mtDNA_CN + WBC_n)

Say that mod2a has an intercept b0 and coefficients bA for Age and bW for WBC_n. Then:

mod2resid = mtDNA_CN - b0 - bA*Age - bW*WBC_n

and you can re-write mod1a and mod3a with flipped predictors associated with mtDNA_CN:

mod1a' <- coxph(Surv(time, status) ~ Age + mod2resid + b0 + bA*Age + bW*WBC_n)
mod3a' <- coxph(Surv(time, status) ~ Age +  mtDNA_CN - b0 - bA*Age - bW*WBC_n)

The point is that (with these simple linear associations) mod1a and mod3a aren't fundamentally different. You will get different estimates for the Cox coefficient for Age depending on your choice, as Age affects outcome both directly and via its association with mtDNA_CN. You don't see a coefficient for WBC_n in mod1a or mod3a; it provides a hidden fixed-coefficient offset due to its involvement in modeling the relation between mtDNA_CN and mod2resid.

What you lose in mod1a and mod3a, versus mod4a, is the ability to capture any association of WBC_n with outcome except insofar as it affects mtDNA_CN or mod2resid.

The above simplicity can disappear if there are nonlinear terms or interactions in the models, but the principles are the same. If you omit anything that both predicts mtDNA_CN and might be associated with outcome, you won't pick that up and you are risk of omitted-variable bias unless you include it in something like your Model 4.

For example, mtDNA_CN can come either from white cells or platelets (or even cell-free DNA, depending on how samples are prepared). The association of mtDNA_CN might well differ depending on how much is from white cells versus platelets. I also suspect that WBC_n and Platelet_n have additional associations with outcome beyond their role in determining mtDNA_CN. For example, Platelet_n could be related to cardiovascular deaths and WBC_n to deaths from immune dysfunction or infection. It would seem wise to include them, at least, explicitly in your model.

If there are purely technical factors affecting the values of mtDNA_CN (perhaps Autosomal_Coverage?) you should consider correcting the mtDNA_CN directly for them first.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.