Best way to model a reponse when a covariate only applies to one group I was using a GLM but am happy with other models.
Fake example: Lets say I have a response such as income. I want to include GED as a factor (0 or 1), but I also have 'parental income' that I only want to apply to people that have a GED. I am mainly modelling against 'time being in employment', with other factors to consider.
What is the best way to do this? Would it be stratify into two subgroups and then model them? I feel like I would lose the relationship between the two groups. Or if there a way to alter model for this? I am an R user if that helps...
 A: You can include a main effect of GED and an interaction between GED and parental income without including parental income. This is equivalent to encoding the restriction that for those without a GED, there is no relationship between parental income (PI) and income, which is a strong assumption to make. Here's the model you would fit:
$$
E[income|GED, PI] = \beta_0 + \beta_1 GED + \beta_2 GED \times PI
$$
You would use the mean-centered version of PI. The interpretation of the coefficients would be the following:

*

*$\beta_0$ - the expected income for those without a GED

*$\beta_1$ - the difference between the expected income for those with a GED and those without a GED for those with an average level of PI

*$\beta_2$ - the slope of PI on income for those with a GED

Typically, a model would also have a main effect term for $PI$, which represents the effect of $PI$ on income for those without a GED, and which changes the interpretation of $\beta_2$ to be the difference between the effect of $PI$ on income for those with a GED and the effect of $PI$ on income for those without a GED.
In R, you would fit this model using the following syntax:
fit <- lm(income ~ GED + GED:PI)

You must use a robust standard error for the tests to be valid because you are by definition explaining more variance in income for those with a GED than for those without, so the residual variance differs between groups. To estimate the robust standard errors, you could use the following syntax:
lmtest::coeftest(fit, vcov. = sandwich::vcovHC)

I don't recommend you fit a model like this unless you were unable to measure PI for those without a GED. The assumption you make is a strong one and almost certainly false, so the rest of the model may not be valid.
