10.3 More than two categories
Residence has two categories. What about a variable with seven?
The data also record religion, which takes seven values.
#>
#> Buddhist Christian Hindu Jain Muslim Sikh Tribal
#> 174 546 16588 22 2143 411 116The temptation is to code these 1 through 7 and enter that as a variable. It would run, and it would be meaningless.
Such a coding tells the regression that the categories are equally spaced and ordered. It would treat “Christian minus Buddhist” as the same distance as “Muslim minus Jain”, and would place Jain exactly midway between Hindu and Muslim. None of those statements means anything. The numbers are labels, not quantities.
The correct approach is one dummy per category — but not quite all of them.
The dummy variable trap
Suppose we build a dummy for urban and another for rural, and include both.
#> (Intercept) urban rural
#> 9.818 1.076 NANA again, and for the reason established in Section 8.4:
\(\text{urban} + \text{rural} = 1\) for every worker, which is exactly the
intercept’s column. One column is an exact linear combination of the others, so
CLRM7 fails.
The dummy variable trap. With \(g\) categories, include \(g-1\) dummies, not \(g\).
The omitted category is the base group. Every coefficient is a comparison with it, and the intercept describes it.
This is not a statistical subtlety but the same arithmetic as Section 8.4: a full set of category dummies always sums to one, and the regression already has a column of ones.
An alternative is to drop the intercept instead and keep all \(g\) dummies, in which case each coefficient is a group mean rather than a difference. Both are defensible; what is not defensible is keeping both the intercept and all \(g\) dummies.
Letting R do it
R handles this automatically if the variable is a factor.
religion_fit <- lm(lw ~ factor(religion), data = wages)
round(coef(summary(religion_fit))[, 1:2], 4)#> Estimate Std. Error
#> (Intercept) 10.1923 0.0951
#> factor(religion)Christian 0.4595 0.1092
#> factor(religion)Hindu -0.1222 0.0956
#> factor(religion)Jain 0.3452 0.2838
#> factor(religion)Muslim -0.0616 0.0989
#> factor(religion)Sikh 0.4263 0.1134
#> factor(religion)Tribal 0.0981 0.1504Six coefficients for seven categories. The missing one is Buddhist, which R chose because it comes first alphabetically — that is the only reason.
factor() is not doing anything clever.
It builds the dummy variables for you, drops one, and hands the rest to the same
lm() that has been used all along. We can see the columns it made.
#> (Intercept) factor(religion)Christian factor(religion)Hindu
#> 1 1 0 1
#> 2 1 0 0
#> 3 1 0 1
#> factor(religion)Jain
#> 1 0
#> 2 0
#> 3 0#> (Intercept) factor(religion)Christian factor(religion)Hindu
#> 20000 546 16588
#> factor(religion)Jain factor(religion)Muslim factor(religion)Sikh
#> 22 2143 411
#> factor(religion)Tribal
#> 116Those are ones and zeros, exactly as if they had been typed by hand. The column of ones is the intercept, and the six others are the dummies — note that their totals are the category counts from the table above, and that Buddhist has no column at all.
Doing it manually gives the same answer:
wages$christian <- as.integer(wages$religion == "Christian")
wages$hindu <- as.integer(wages$religion == "Hindu")
round(coef(lm(lw ~ christian + hindu, data = wages)), 4)#> (Intercept) christian hindu
#> 10.2140 0.4378 -0.1439Here the base group is everyone else, since only two dummies were built. That is
the point: factor() saves typing and prevents the trap, but the model it fits
is one you could have assembled yourself.
The base group is arbitrary, and the default is rarely the best choice.
Here it is a category with 174 observations out of 20,000. Every coefficient in the table is a comparison with a small group, and every standard error is inflated by that group’s imprecision. Nothing is wrong — the fitted values, \(R^2\) and \(F\) statistic are identical whatever the base — but the table is harder to read than it needs to be.
A large or otherwise natural category makes a better reference.
wages$religion <- relevel(factor(wages$religion), ref = "Hindu")
round(coef(summary(lm(lw ~ religion, data = wages)))[, 1:2], 4)#> Estimate Std. Error
#> (Intercept) 10.0701 0.0097
#> religionBuddhist 0.1222 0.0956
#> religionChristian 0.5817 0.0546
#> religionJain 0.4674 0.2676
#> religionMuslim 0.0606 0.0288
#> religionSikh 0.5484 0.0626
#> religionTribal 0.2203 0.1169With the largest group as the base, each coefficient is a comparison with it and the standard errors are smaller.
Notice that changing the base group changed every coefficient without changing the model at all. The fit is identical; only the comparison being reported has moved.
A coefficient on a dummy is never meaningful on its own. It is meaningful relative to a base group, and a reader who does not know which group that is cannot interpret the table.