What's the justification for comparing two separate models built on subsets of data versus using one model that uses the whole dataset?

I've noticed that there are some data analysis being done in some scientific field where the authors would split out an entire dataset into subsets based on a particular property. One classic example would be splitting an entire dataset into subsets containing control versus disease subjects. What got me a little confused is that the authors would then proceed to build two separate models, one for each subset of the data, and then proceed to compare the models head to head.

Some authors would simply look at the coefficients and say there are "distinct" signatures that separate the regression, others could use a Z test (credit to the post: Testing equality of coefficients from two different regressions) where the formula is as follows:

$$Z = \frac{( \beta_1 - \beta_2)}{\sqrt((SE\beta_1)^2 + (SE\beta_2)^2)}$$

Where $$\beta_1$$ and $$\beta_2$$ are the coefficients of each respective model, and $$SE\beta_1$$ and $$SE\beta_2$$ are the standard errors of their respective models.

What confuses me is what's the distinction between using one linear model that interrogates the entire dataset by having the groups (in this example disease versus control) as a fixed effect versus building two separate models and then comparing their coefficients instead?

I believe there is a difference in the math, especially since if it is a linear regression, the OLS method would result in different coefficients if the entire dataset was used (versus building two separate models based on their respective subsets). This is even more alarming to me since if the residuals are different, then the interpretationof one large model versus comapring two separate models are not the same, and then so what is the difference in the interpretation?

Running two models lets you compare the model for the two subsets -- e.g. in your example, control and disease subsets. Also, when you have specifically "control" and "disease" that is often the dependent variable, but that's only for that particular case. Another common splitting variable is sex: Male and female or maybe male, female, other, or some other list.

Running one model does not let you compare, unless you add interaction terms, but they can be hard to interpret and, if you have a lot of independent variables, can substantially reduce the degrees of freedom. (As evidence of the difficulty, just look at how many questions there are on this site about interpreting them!)

I'm not sure why this alarms you. As long as you are clear on what you are doing, there shouldn't be any confusion or alarm.

EDIT: An example in R

set.seed(1234)

x <- rnorm(1000)
sex <- c(rep("M", 500), rep("F", 500))
y <- 3*x + as.numeric(sex == "M") + as.numeric(sex == "M")*x + rnorm(1000, 0, 3)

mod1 <- lm(y~x + sex + x*sex)
summary(mod1) # parameter estimate x 3.10, se 0.14, R^2 = 0.62

modfemale <- lm(y~x, subset = sex == "F")
summary(modfemale) # parameter estimate x 3.10, se 0.14 (same), R^2 = 0.48

modmale <- lm(y~x, subset = sex == "M")
summary(modmale) # parameter estimate x 4.21, se 0.1 (different) R^2 = 0.71

• Would it be fair to say then that I should get the same result if I was to run one model with interaction versus running two separate models? My confusion stems from the fact that the residuals of the model would be different (two models vs one), which could affect the coefficients, which then alter the interpretation of the data. If they are the same, then why do these two methods exist? If not, then what's the difference? Feb 19 at 21:47
• If you add the interaction term, the parameter estimate for the other independent variable will be the same as for the subset that codes the class variable to 0. But not for the other subset. I will add an example. But R^2 will be different for all 3, as will df and other things. Feb 19 at 22:45
• Thank you for the example! When I look at it, my interpretation is that the male model is a better fit, and the female model somehow has the same coefficient for the regression line as the combined model, just with a lower goodness of fit metric. Shouldn't there be three different parameter estiamtes for mod1? When I ran the example in Rstudio, I had three parameter estimates (x, sexM, interaction between x and sexM) and what's surprising to me is that the estimate for sexM is nothing like the estimate for modmale. Could you help explain why? I'm still confused about that. Feb 20 at 2:29
• Yes, there are three. I just put one in my answer. There is no reason why the estimate for sexM ought to be the like anything in modmale. The first is looking at the relationship of sex to the DV, the second is looking at the relationship of x to y among males. Feb 20 at 12:07
• Since the grouped model is looking at the effect of sex whilst the second is looking at the relationship of x to y among males, then why is the model looking at the relationship of x to y among females have the same coefficient for x as the grouped model? I'm interpretating it as the relationship of x on y for females only (the modfemale) is the same as the relationship of x on y as the grouped model, just with a lower degree of "goodness of fit". Feb 20 at 20:25