lm versus lmFit

We have seen two ways to use linear models

1. To model a single variable (CT of one gene, number of cells, methylation level, etc.) we use just R

fitted.model <- lm(formula, data)

2. To model many genes in an array or RNAseq, we use the limma package

fitted.model <- lmFit(design.matrix, data)

The math is the same, the code is different

lm versus lmFit

• Data for lm is a data frame with one colum for each variable
  • for example: diet/stress, tissue, age
• Data for lmFit is a matrix with genes/miRNA in the rows, experiments in the columns, and the expression level in each cell
  • We need an extra table describing the conditions of each experiment

From formula to design.matrix

The main difference between lm and lmFit is how to describe the model

• formula
• matrix

If we have formula, we can get the design matrix using

design.matrix <- model.matrix(formula)

It is not hard to make your own design matrix if you know the math

Math to matrix

\[y_i = β_0 + β_1 x_i + e_i\]

\[\begin{pmatrix} y_1\\ ⋮ \\ y_i\\ ⋮ \\ y_n \end{pmatrix} = \begin{pmatrix} 1 & x_1\\ ⋮ & ⋮ \\ 1 & x_i\\ ⋮ & ⋮ \\ 1 & x_n\\ \end{pmatrix} \begin{pmatrix}β_0\\β_1\end{pmatrix} + \begin{pmatrix} e_1\\ ⋮ \\ e_i\\ ⋮ \\ e_n \end{pmatrix}\]

The fitted model has the smallest value for \(∑_i e_i^2\)

Two independent variables

\[y_i = β_0 + β_1 x_{i,1} + β_2 x_{i,2} + e_i\]

\[\begin{pmatrix} y_1\\ ⋮ \\ y_i\\ ⋮ \\ y_n \end{pmatrix} = \begin{pmatrix} 1 & x_{1,1} & x_{1,2}\\ ⋮ & ⋮ & ⋮ \\ 1 & x_{i,1} & x_{i,2}\\ ⋮ & ⋮ & ⋮ \\ 1 & x_{n,1} & x_{n,2}\\ \end{pmatrix} \begin{pmatrix}β_0\\β_1\\β_2\end{pmatrix} + \begin{pmatrix} e_1\\ ⋮ \\ e_i\\ ⋮ \\ e_n \end{pmatrix}\] \(β_0\) (intercept) is the value of \(y_i\) when all \(x_{i,j}=0\)

General case

\[y_i = β_0 + ∑_{j=1}^k β_j x_{i,j} + e_i\]

\[\begin{pmatrix} y_1\\ ⋮ \\ y_i\\ ⋮ \\ y_n \end{pmatrix} = \begin{pmatrix} 1 & x_{1,1} &\cdots& x_{1,j} &\cdots& x_{1,k}\\ ⋮ & ⋮ & ⋱ & ⋮ & ⋱ & ⋮ \\ 1 & x_{i,1} &\cdots& x_{i,j} &\cdots& x_{i,k}\\ ⋮ & ⋮ & ⋱ & ⋮ & ⋱ & ⋮ \\ 1 & x_{n,1} &\cdots& x_{n,j} &\cdots& x_{n,k}\\ \end{pmatrix} \begin{pmatrix}β_0\\β_1\\⋮\\β_j\\⋮\\β_k\end{pmatrix} + \begin{pmatrix} e_1\\ ⋮ \\ e_i\\ ⋮ \\ e_n \end{pmatrix}\] \(β_1,…,β_k\) are the effects of each independent variables

Model without intercept

\[y_i = β_0 + β_1 x_{i,1} + β_2 x_{i,2} + e_i\]

\[\begin{pmatrix} y_1\\ ⋮ \\ y_i\\ ⋮ \\ y_n \end{pmatrix} = \begin{pmatrix} x_{1,1} & x_{1,2}\\ ⋮ & ⋮ \\ x_{i,1} & x_{i,2}\\ ⋮ & ⋮ \\ x_{n,1} & x_{n,2}\\ \end{pmatrix} \begin{pmatrix}β_1\\β_2\end{pmatrix} + \begin{pmatrix} e_1\\ ⋮ \\ e_i\\ ⋮ \\ e_n \end{pmatrix}\] This is useful when \(x_{i,j}\) represent a factor, but may not be realistic when \(x_{i,j}\) is a number

General case without intercept

\[y_i = ∑_{j=1}^k β_j x_{i,j} + e_i\]

\[\begin{pmatrix} y_1\\ ⋮ \\ y_i\\ ⋮ \\ y_n \end{pmatrix} = \begin{pmatrix} x_{1,1} & x_{1,j} & x_{1,k}\\ ⋮ & ⋮ & ⋮ \\ x_{i,1} & x_{i,j} & x_{i,k}\\ ⋮ & ⋮ & ⋮ \\ x_{n,1} & x_{n,j} & x_{n,k}\\ \end{pmatrix} \begin{pmatrix}β_1\\⋮\\β_j\\⋮\\β_k\end{pmatrix} + \begin{pmatrix} e_1\\ ⋮ \\ e_i\\ ⋮ \\ e_n \end{pmatrix}\]

Factors

Let’s say that \(f\) is a factor with \(k\) levels \((f_i ∈\{l_1,…,l_k\})\)

We represent it with \(k\) columns

\[x_{i,j} = \begin{cases}1\text{ if }f_i\text{ is equal to }l_j\\ 0\text{ if }f_i\text{ is not equal to }l_j\end{cases}\]

Example: sex

\[y_i = β_1 x_{i,1} + β_2 x_{i,2} + e_i\]

\[\begin{pmatrix} y_1\\ ⋮ \\ y_n \end{pmatrix} = \begin{pmatrix} x_{1,1} & x_{1,2}\\ ⋮ & ⋮ \\ x_{n,1} & x_{n,2}\\ \end{pmatrix} \begin{pmatrix}β_1\\β_2\end{pmatrix} + \begin{pmatrix} e_1\\ ⋮ \\ e_n \end{pmatrix}\] If, for each female we have \(x_{i,1}=1\) and \(x_{i,2}=0,\)
then \(β_1\) will be the mean of \(y_i\) for all female

Formula to math

Let’s say that x is a numeric variable

• y ~ x is \(y_i=β_0 + β_1 x_i + e_i\)

• y ~ x+0 is \(y_i=β_1 x_i + e_i\)

• y ~ x1 + x2 + x3 is \(y_i=β_0 + β_1 x_{i,1} + β_2 x_{i,2} + β_3 x_{i,3} + e_i\)

• y ~ x1 + x2 + x3 + 0 is \(y_i=β_1 x_{i,1} + β_2 x_{i,2} + β_3 x_{i,3} + e_i\)

Formulas with factors

Let’s say that f is a factor with \(k\) levels \((f_i ∈\{l_1,…,l_k\})\)

y ~ f+0 is \(y_i=\sum_{j=1}^k β_j x_{i,j} +e_i,\) where \[x_{i,j} = \begin{cases}1\text{ if }f_i\text{ is equal to }l_j\\ 0\text{ if }f_i\text{ is not equal to }l_j\end{cases}\]

f <- factor(c("A","B","C"))
model.matrix(~f+0) |> as.data.frame()

  fA fB fC
1  1  0  0
2  0  1  0
3  0  0  1

Formulas with factors and intercept

Now the first level is not encoded

y ~ x is \(y_i=β_0 + \sum_{j=2}^k β_j x_{i,j} +e_i\)

f <- factor(c("A","B","C"))
model.matrix(~f) |> as.data.frame()

  (Intercept) fB fC
1           1  0  0
2           1  1  0
3           1  0  1

It is critical to choose well the first level
It is the baseline

Interpretation

\[y_i=β_0 + \sum_{j=2}^k β_j x_{i,j} +e_i\]

• \(β_0\) is the mean of \(y_i\) for all cases \(f_i\) is \(l_1\) \[β_0=\text{mean}(y_i \vert f_i=l_1)\]
• \(β_j\) is the difference of means between \(l_j\) and \(l_1\) \[β_j=\text{mean}(y_i | f_i=l_j)-\text{mean}(y_i | f_i=l_1)\]