Linear Model with Categorical Predictor Variable: What's the Math behind it?

Question

I do understand how linear regression models work with numerical data (finding the line that has the least squares). Having read Grolemund's and Wickham's R for Data Science, I was wondering about one of their examples which goes like this (edit: This is only a fraction of the data they use but I wanted to create a small and reproducable example without tidyverse, as Roland pointed out below):

(daily <- data.frame(n = c(842,943,914,915,720,832,933,899,902,932,930,690,828,928,894,901,927,924,674,786), wday = c("Tue","Wed", "Thu", "Fri", "Sat", "Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun")))

mod <- lm(n ~ wday, data = daily)

Obviously, the cut could be regarded as an ordinal variable so that it can easily converted to the number but as the plot shows, this will not result in a linear line. I did consult R's lm() manual but did not find any information why the use of lm() here is ok.

Answer 1

First of all, days are not ordinal, but nominal. (The elegant way to include them in a regression is to convert them into a factor, you haven't done that, but in this particular case it doesn't matter, because right now they have character type, which will be automatically converted to factor.)

So the question is, what goes on under the hood with nominal variables? (R will call them "unordered factor".) This: $$Y_i = \beta_0 + \beta_{Tue} D_{i,Tue} + \beta_{Wed} D_{i,Wed} + \ldots + \beta_{Sun} D_{i,Sun}+\varepsilon.$$ Here $D_{i,Tue}$ is a variable which takes the value of 1, if the $i$th day is Tuesday, zero otherwise. Likewise for the remaining days. These are usually called dummy, or indicator variables. In R, they're generated totally automatically, you don't have to do anything, except for making sure that R handles the variable as factor. Note, that there is no indicator for Monday.

So, the important point is: you can now see that it is indeed a usual linear model!

How it works? Just think over the interpretation of the coefficients. If the day is Monday, then every dummy takes the value of zero, every product will be zero, so we have $Y_i = \beta_0 +\varepsilon$, i.e., $\beta_0$ will be the estimated value for Monday. (Simply the average in this case.) If we are on Tuesday, we have $Y_i = \beta_0 + \beta_{Tue} +\varepsilon$, so $\beta_{Tue}$ is the difference between the estimated values of Tuesday and Monday. (As $Y_i-\beta_0 = \beta_{Tue} +\varepsilon$, but we already know that $\beta_0$ is the estimated value for Monday.) And so on!

So you can see why we included one dummy less: the omitted level will be the so-called reference, what you see in the intercept, and every other one will have an interpretation of "compared to" the reference. One important remark: what will be the reference level is something that you can set yourself, by explicitly calling factor or using relevel on a variable that is already a factor, or it will be decided by R automatically, but if the original variable is a character, then it'll be the first in alphabetical order (which is not necessarily meaningful!).

If we were to include all dummys, we would run into perfect collinearity as the sum of the dummys would be a constant 1 - just as the intercept. Of course we could write $Y_i = \beta_{Mon} D_{i,Mon} + \beta_{Tue} D_{i,Tue} + \beta_{Wed} D_{i,Wed} + \ldots + \beta_{Sun} D_{i,Sun}+\varepsilon$, in this case every coefficient will take the value of the estimated response on that day - use the logic introduced above to track this! - but we don't really like this approach (just think over what would happen if we had to include two factors!).

Note, that these are just the most basic ways to encode categorical variables. ?contr.treatment or this page gives you more information on the other possible options.

Answer 2

summary(mod)

gives you a nice overview over what's going on under the hood. Basicallly what happens is that the categorical variable is automatically turned into a number of binary variables, i.e. is it a monday (0 or 1), is it a tuesday, and so on. One thing to note is that there are seven days, but only six binary variables, so that linear dependency is avoided.

Answer 3

How to treat "day of week" is an interesting question and I am not sure we can just automatically say that treating it as nominal is best in all cases. It will depend on context (and none was offered here). One commonly useful approach would be weekday and weekend - and, indeed, these data show huge differences between weekdays and weekends. But I could easily see cases where some other categorization would be useful, e.g.

"How much time do you spend in church?" Sunday vs. other seems sensible but

"How much time do you spend in worship?" might need special treatment depending on the religion of the person.

Or, getting to my old area of research - "How much alcohol do you drink on a typical day?" is often a bad question as it varies hugely by day.

Do you eat dinner out? is likely to be different on Fridays and Saturdays from other days (at least in the parts of the world where Sat and Sunday are typically weekend).

In addition, R (by default) chose Friday as the reference (because it is first, alphabetically). This may or may not be ideal. Yes, if you look carefully, the results mean essentially the same thing regardless of choice of reference, but people who get hung up on statistical sig. may find that changing the reference changes what is significant.