How can I deal with a covariate being defined only for a subset of my sample?

Question

My model looks like: ROI_size ~ diagnosis + medication_dose + sex + age. Specifically, I want to find the effect of disease (1 or 0), adjusted for current medication dose (measured in mg) on brain region size, while controlling for age and sex. Both of the latter variables are available for the entire dataset. However, the former (medication dose) is defined only for the patients (diagnosis = 1).

I have thought of two solutions. One would be to assign medication_dose = 0 to all healthy controls. However, this creates 1. collinearity with diagnosis 2. zero-inflated data and is per se not correct (I think), since the difference between 0 and 1 mg is not the same as the difference between 1 mg and 2 mg.

The other solution would be to adjust the ROI_size for medication_dose within the patient cohort, i.e., to run the model ROI_size ~ medication_dose and then by calculating its beta and then adding beta*medication_dose to the ROI_size for all patients.

My question in this case is if I should have age and sex in the second model, since including them will yield a more accurate beta for medication_dose. In that case, do I need to correct for age and sex separately for the healthy controls as well and then to find the effect of disease to just do ROI_size ~ Diagnosis?

Histogram for current medication use in patients

Histogram for current medication use in all

Answer 1

By (arbitrary) convention, people say you have 'a problem' with collinearity when the variance inflation factor ($VIF$) is greater than or equal to $10$. When considering a pair of variables, that means the Pearson correlation is $\gtrapprox .95$. You report that the Spearman correlation is $\approx .6$. If the Pearson correlation were the same (and the other variables were uncorrelated), that would make the $VIF = 1.56$. In other words, the variance of the estimated sampling distributions would be $\approx 1.6\times$ larger than what you could have ideally achieved, or the standard errors would be $1.25\times$ as wide. That's hardly any difference at all.

Assuming there is sufficient variability in doses, and that the relationship is sufficiently linear from a dose of $0$ to your maximum dose, the coefficient on disease gives you an estimate of the size of the ROI for a patient with the disease but who isn't taking any medication. The coefficient on the dose gives you the estimated association between changes in size and a one unit increase in dose. That is potentially useful information. I would start with that model (i.e., assigning ${\rm dose}=0$ for patients not on the medication) and assess to see if the assumptions appear tenable (cf., EdM's answer). Also bear in mind that these data are observational and the dose is endogenous (physicians prescribed higher doses based on their assessments of the patients), so you cannot assume the coefficient is an unbiased estimate of the causal effect of the medication. But that's typical of biomedical research, and the results still have some value (properly understood).

Answer 2

... there is no healthy control that has medication > 0 and no patient that has medication = 0.

That suggests a different two-step modeling approach from your second solution. Start without the medication.dose term in the model. That evaluates the critical issue of diagnosis directly, controlling for Sex and Age. Then add the medication.dose term and see whether it improves the fit; e.g., anova() comparison of those two nested models.

If the fit is not improved, then you have documented that adjusting for medication.dose doesn't matter (much) in terms of ROI_size when diagnosis, Sex and Age are taken into account.

If the fit is improved, then you have an estimate of the association between ROI_size and medication.dose. The model is then equivalent to what @whuber recommends in a different context for a predictor whose value is necessarily 0 for some cases. Here, the diagnosis predictor serves the role of the "loan indicator" in that question, with a corresponding interpretation of intercepts and regression coefficients.

There's an important caution about how to handle medication.dose, however. If you treat it as an untransformed continuous variable then you are imposing a linear association between medication.dose and ROI_size. If there are only a few different values of medication.dose then you might be best off treating it as a multi-level categorical predictor, perhaps ordinal. If there's a wide range of values, you might need to transform appropriately.

You would have a complete collinearity problem if there was only one value of medication.dose for all those with a diagnosis; in that case you would not be able to fit a model with both diagnosis and medication.dose uniquely. The partial collinearity that you fear with more than one medication.dose might tend to make the error estimates of the regression coefficients larger. If you are interested specifically in the association of medication.dose with ROI_size for those with diagnosis you might be better off restricting analysis to those with diagnosis. The question (which I don't think can be answered except by looking at the data) is whether having the correction for Age and Sex from extra cases without diagnosis makes up for that.