# What would be the simplest way to do a linear regression separated by factors? [closed]

I have this data in a dataframe, df:

>head(df)
sx   rk yr        dg yd    sl
1   male full 25 doctorate 35 36350
2   male full 13 doctorate 22 35350
3   male full 10 doctorate 23 28200
4 female full  7 doctorate 27 26775
5   male full 19   masters 30 33696
6   male full 16 doctorate 21 28516
...


Where sx is "sex", yd is "years after highest degree received" and sl is "salary." I want to produce a regression for sl~yd, but I want to split across sx. In other words, I want to see how males and female salaries change over their career.

Would there be a way to do this using a single lm() call? Thanks.

• I, personally, would be a little more interested in a loess vs. age. It might rapidly accelerate, or cross-over, or get stabilized difference and have equal slope with constant separation. Nice question. – EngrStudent Nov 15 '17 at 23:08
• I'm not sure whether this question is on-topic on CV or would be better on StackOverflow - the way it is phrased ("using a single lm() call") makes it sound like a code question, but there are also statistical issues involved, e.g. if this can be achieved by using a different model. Perhaps it can be rephrased somehow to make it a clearer fit for the sit (i.e. more clearly "statistical" and less "codey"). – Silverfish Nov 15 '17 at 23:39
• – Scortchi - Reinstate Monica Nov 16 '17 at 0:22

> fit <- lm(sl ~ 0 + sx + sx:rd, data=df)
> summary(fit)


This fits one regression model with separate intercepts and slopes for the two sexes (four parameters in total).

The intercepts and slopes from the above fit will be the same as you would get from fitting simple linear regressions to the males and females separately. However the residual variance and standard errors will be different. Statistical inference from the combined fit will be more powerful than doing entirely separate fits for males and females because the error variance is pooled over both sexes, giving more residual df and hence more statistical power than the separate regressions would provide.

Another advantage of fitting it all with one model is that the combined approach allows you to conduct a statistical tests of whether the slopes or intercepts are different between the sexes. Surely you would want to do that.

• Ding ding ding! We have a winner! I have designated you as the accepted answer. The only downside is that p-values seem to change--lm() must be doing some multiple hypothesis testing behind the scenes. – Jeff Nov 15 '17 at 23:10
• lm does not adjust for multiple comparisons behind the scenes. If you fit models on different subsets or a different design matrix you would expect different p-values. – alexpghayes Nov 15 '17 at 23:17
• What model do you intend to fit with that code? – Scortchi - Reinstate Monica Nov 16 '17 at 0:31

Why not use these two lines rather than a single lm call?

lm_male <- lm(sl ~ yd, data=subset(df, sx=='male'))
lm_female <- lm(sl ~ yd, data=subset(df, sx=='female'))


EDIT: I've edited the answer to make it clear I'm not asking a question about the code.

• I think it adds to readability and reproducibility. – EngrStudent Nov 15 '17 at 23:09
• Nothing is wrong with two lines. I just want to see if I can be clever and only use one call to lm and end up with only one object to deal with. If I were to call lm on a model with many factors, I'd prefer to only have to make one call to get my model. – Jeff Nov 15 '17 at 23:12
• If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. - From Review – Ferdi Nov 15 '17 at 23:25
• @Ferdi It isn't asking for debugging, it's a suggested alternative. I think it would benefit from a little rewrite: "Why not use these two lines rather than a single lm call?" or similar, so that it doesn't resemble a question so much! – Silverfish Nov 15 '17 at 23:41
• Nothing, except that you might want to use a common estimate of the error variance. – Scortchi - Reinstate Monica Nov 16 '17 at 0:25