# Recreate lm Categorical Regression

Consider the code, which contains regression using lm of two categorical and one continuous variables without interaction using data from the correct model:

set.seed(12345)

n <- 100
X1 <- as.factor(sample(c("a", "b", "c"), n, replace = TRUE))
X2 <- as.factor(sample(c("d", "e"), n, replace = TRUE))
X3 <- rnorm(n, mean = 0, sd = 1)

effeta <- 1
effetb <- 2
effetc <- 3
effetd <- 4
effete <- 5

Y <- 10 + 3*X3 + effeta*(X1 == "a") + effetb*(X1 == "b") + effetc*(X1 == "c") + effetd*(X2 == "d") + effete*(X2 == "e") + rnorm(n, mean = 0, sd = 1)

model <- lm(Y ~ X1 + X2 + X3)
print(summary(model))


The output is :

Call:
lm(formula = Y ~ X1 + X2 + X3)

Residuals:
Min       1Q   Median       3Q      Max
-2.74679 -0.60571  0.06521  0.67873  1.85143

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.06741    0.20151  74.773  < 2e-16 ***
X1b          0.95913    0.23994   3.997 0.000127 ***
X1c          1.96010    0.23637   8.292 7.26e-13 ***
X2e          1.08521    0.18887   5.746 1.10e-07 ***
X3           3.03324    0.09537  31.804  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9327 on 95 degrees of freedom
Multiple R-squared:  0.9194,    Adjusted R-squared:  0.916
F-statistic: 270.8 on 4 and 95 DF,  p-value: < 2.2e-16


The objective is to understand MATHEMATICALLY what exactly is going on and then to reproduce it in R. More specifically :

1. One can say that X3 coefficient is correctly estimated. Sources online indicate that lm will take one level as reference, set its effect to zero. In a generic situation, how to predict in advance which variable will be taken as reference? In alphabetical order perhaps? Once the effects are estimated, how do we see that they are correct ?
2. Although there is no reference, I would believe that R writes the categorical observations using dummy variables to combine everything in a matrix like mentionned in Searle and then do multiple linear regression. How do I write down such a matrix in this case?

Let us write $$x_3 \in \mathbb{R}^n$$ in place of x3. Then :

$$\begin{equation*} x_3 = \begin{bmatrix} x_{1,3} \\ x_{2,3} \\ \vdots \\ x_{n,3} \end{bmatrix} \end{equation*}$$

To do it blindly for the categorical variables, one writes :

$$\begin{equation*} x_1 = \begin{bmatrix} x_{1,1a} && x_{1,1b} && x_{1,1c} \\ x_{2,1a} && x_{2,1b} && x_{2,1c} \\ \vdots && \vdots && \vdots\\ x_{n,1a} && x_{n,1b} && x_{n,1c} \\ \end{bmatrix}; x_2 = \begin{bmatrix} x_{1,1d} && x_{1,1e} \\ x_{2,1d} && x_{2,1e}\\ \vdots && \vdots \\ x_{n,1d} && x_{n,1e} \\ \end{bmatrix} \end{equation*}$$

All entries in the two matrices are $$0$$ or $$1$$. And then put :

$$\begin{equation*} X = \begin{bmatrix} 1 &&x_{1,1a} && x_{1,1b} && x_{1,1c} && x_{1,1d} && x_{1,1e} && x_{1,3} \\ 1 && x_{2,1a} && x_{2,1b} && x_{2,1c} && x_{2,1d} && x_{2,1e} &&x_{2,3} \\\\ \vdots && \vdots && \vdots && \vdots && \vdots && \vdots && \vdots \\ 1 && x_{n,1a} && x_{n,1b} && x_{n,1c} && x_{n,1d} && x_{n,1e} && x_{n,3} \end{bmatrix} \end{equation*}$$

So that we can do regression lm(Y~X). But this is obviously not what happens. How do I take away some columns to get what actually happens?

• If you drop the second and fifth column in $X$, you get the design matrix that R passes to the internal QR solver. Check with model.matrix(model). By default, factors are dummy encoded, dropping the "first" level (as being linearly dependent). Commented Jul 18 at 20:06

There actually is a reference given, Statistical Models in S by Chambers and Hastie.

R does not necessarily use dummy variables. The contrasts option, or contrasts set with the function C describe how factors are coded. The factory default is contr.treatment, which does give dummy variables for all levels except the first. Here's the coding for a four-level variable. Note that it has three columns.

> contr.treatment(4)
2 3 4
1 0 0 0
2 1 0 0
3 0 1 0
4 0 0 1


There are other options

> contr.sum(4)
[,1] [,2] [,3]
1    1    0    0
2    0    1    0
3    0    0    1
4   -1   -1   -1
> contr.helmert(4)
[,1] [,2] [,3]
1   -1   -1   -1
2    1   -1   -1
3    0    2   -1
4    0    0    3
> contr.SAS(4)
1 2 3
1 1 0 0
2 0 1 0
3 0 0 1
4 0 0 0
> contr.poly(4)
.L   .Q         .C
[1,] -0.6708204  0.5 -0.2236068
[2,] -0.2236068 -0.5  0.6708204
[3,]  0.2236068 -0.5 -0.6708204
[4,]  0.6708204  0.5  0.2236068


This gives blocks of columns of the design matrix, which just get pasted together sideways. If you have interactions things get a bit more complicated: you multiply each of the columns for one variable by each of the columns for the other variable (easy) and then decide which columns need to be dropped (annoying).

In your example, if you do model.matrix(model) you get the model matrix. Here's the top of it

> head(model.matrix(model))
(Intercept) X1b X1c X2e          X3
1           1   1   0   0  1.39160685
2           1   0   1   0 -0.30383538
3           1   1   0   1  1.59350961
4           1   1   0   1  0.78706366
5           1   0   0   1  0.45555145
6           1   0   1   1  0.08665505


and, importantly, here is the bottom of it

99            1   1   0   1 -0.49347541
100           1   0   1   1  0.16951075
attr(,"assign")
[1] 0 1 1 2 3
attr(,"contrasts")
attr(,"contrasts")$X1 [1] "contr.treatment" attr(,"contrasts")$X2
[1] "contr.treatment"


The assign attribute tells you which formula terms each column came from. The first column is the intercept, the next two are X1, then one for X2 and one for X3. It also tells you that the two factor variables did use treatment 'contrasts'.

R gets this information by processing the model formula with terms

> terms(model)
Y ~ X1 + X2 + X3
attr(,"variables")
list(Y, X1, X2, X3)
attr(,"factors")
X1 X2 X3
Y   0  0  0
X1  1  0  0
X2  0  1  0
X3  0  0  1
attr(,"term.labels")
[1] "X1" "X2" "X3"
attr(,"order")
[1] 1 1 1
attr(,"intercept")
[1] 1
attr(,"response")
[1] 1
attr(,".Environment")
<environment: R_GlobalEnv>
attr(,"predvars")
list(Y, X1, X2, X3)
attr(,"dataClasses")
Y        X1        X2        X3
"numeric"  "factor"  "factor" "numeric"


The "term.labels" attribute specifies the terms in the model. model.matrix goes through and makes a set of columns for each term using the specified contrasts and accounting appropriately for interactions.

In this case there aren't any interactions, so you get 1 for the intercept, you get coding by

> contr.treatment(3)
2 3
1 0 0
2 1 0
3 0 1


for X1, you get coding by

> contr.treatment(2)
2
1 0
2 1


for X2, and X3 just gets put in as a column on its own.

• Thank you very much! This is extremely helpful. Just one more thing though. Do you mind answering the first question in the post as well? Namely, given lm approximation of effects, how do you determine mathematically what the effects are in the original format ?
– 温泽海
Commented Jul 19 at 14:30
• The contrasts matrix specifies which contrasts are estimated -- eg, if you're using treatment contrasts, they are differences between the each other level of the factor and the first level. The ordering is the ordering of the factor levels. If you didn't tell R what ordering to use it will be alphabetic in your current locale Commented Jul 21 at 21:40