1 Linear regression
Let’s say that we have \(n=10\) observations for a predictor (explanatory) variable \(x\) and a dependent variable \(y\):
set.seed(123)
(x <- round(rnorm(10, mean = 0, sd = 1), digits = 2))## [1] -0.56 -0.23 1.56 0.07 0.13 1.72 0.46 -1.27 -0.69 -0.45(y <- 2 + 3 * x + rnorm(10, mean = 0, sd = 1))## [1] 1.5440818 1.6698138 7.0807715 2.3206827 1.8341589 8.9469131
## [7] 3.8778505 -3.7766172 0.6313559 0.1772086ggplot(data.frame(x = x, y = y),
aes(x = x, y = y)) +
geom_point(size = 4) + theme_minimal()
Linear regression of \(y\) on \(x\) would imply using the model \[y_i = \beta_0 + \beta_1\cdot x_i + \varepsilon_i\] for \(i=1, ..., n\), and where \(\beta_0\) (intercept) and \(\beta_1\) (slope) are parameters that are unknown and will be estimated from the data. In other words, \[y_1 = \beta_0 + \beta_1\cdot x_1 + \varepsilon_1\\y_2 = \beta_0 + \beta_1\cdot x_2 + \varepsilon_2\\...\\y_n = \beta_0 + \beta_1\cdot x_n + \varepsilon_n\]
We can write this in matrix form as follows: \[\begin{pmatrix}y_1\\y_2\\...\\y_n\end{pmatrix} = \begin{pmatrix}1 & x_1\\1 & x_2\\... & ...\\1 & x_n\end{pmatrix}\begin{pmatrix}\beta_0\\\beta_1\end{pmatrix} + \begin{pmatrix}\varepsilon_1\\\varepsilon_2\\...\\\varepsilon_n\end{pmatrix}\]
The first matrix to the right of the equal sign is the design matrix.
In R, we would express this model in formula notation as y ~ x. This will by default include an intercept (unless we explicitly tell it not to), and since x is a continuous variable, it will include a single coefficient for it (what we have referred to as \(\beta_1\) above). The design matrix for a given model can be created using the model.matrix() function in R (note that this does not require knowledge about the response variable):
model.matrix(~ x)## (Intercept) x
## 1 1 -0.56
## 2 1 -0.23
## 3 1 1.56
## 4 1 0.07
## 5 1 0.13
## 6 1 1.72
## 7 1 0.46
## 8 1 -1.27
## 9 1 -0.69
## 10 1 -0.45
## attr(,"assign")
## [1] 0 1We can also fit a linear model in R, which will provide estimates for the coefficients \(\beta_0\) and \(\beta_1\):
summary(lm(y ~ x))##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.33130 -0.64381 -0.02429 0.49406 1.41357
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1622 0.2850 7.587 6.38e-05 ***
## x 3.6279 0.3131 11.587 2.80e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8982 on 8 degrees of freedom
## Multiple R-squared: 0.9438, Adjusted R-squared: 0.9367
## F-statistic: 134.3 on 1 and 8 DF, p-value: 2.799e-06Note that the coefficients are not displayed as \(\beta_0\) and \(\beta_1\), but rather as (Intercept) and x (to be interpreted as ‘the coefficient for the x column in the design matrix’).
For a given value of \(x\), the fitted/predicted value of \(y\) can now be obtained as \[\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1\cdot x_i\] or in matrix form \[\begin{pmatrix}\hat{y}_1\\\hat{y}_2\\...\\\hat{y}_n\end{pmatrix} = \begin{pmatrix}1 & x_1\\1 & x_2\\... & ...\\1 & x_n\end{pmatrix}\begin{pmatrix}\hat{\beta}_0\\\hat{\beta}_1\end{pmatrix}\] In other words, the fitted values are obtained by matrix multiplication of the design matrix and the estimated coefficients.
ggplot(data.frame(x = x, y = y),
aes(x = x, y = y)) +
geom_point(size = 4) + theme_minimal() +
geom_smooth(method = "lm", se = FALSE)## `geom_smooth()` using formula = 'y ~ x'
The ExploreModelMatrix Bioconductor package gives us a way to explore the design, either interactively or programmatically:
library(ExploreModelMatrix)
library(cowplot)
if (interactive()) {
ExploreModelMatrix(sampleData = data.frame(y = y, x = x),
designFormula = ~ x)
}
vd <- VisualizeDesign(sampleData = data.frame(y = y, x = x),
designFormula = ~ x)
cowplot::plot_grid(plotlist = vd$plotlist)
The displayed plot indicates that the fitted value for a given observation is obtained by taking the intercept + the value of \(x\) times the coefficient for x.
2 Categorical predictors
Let’s now assume that we have a categorical predictor treat:
set.seed(123)
(treat <- gl(n = 2, k = 5, labels = c("ctrl", "trt")))## [1] ctrl ctrl ctrl ctrl ctrl trt trt trt trt trt
## Levels: ctrl trty <- rnorm(n = 10, mean = rep(c(1, 3), each = 5))
ggplot(data.frame(treat = treat, y = y),
aes(x = treat, y = y)) +
geom_point(size = 4) + theme_minimal()
Of course we can’t multiply a coefficient with a character, so instead we recode the predictor using so called dummy variables. A dummy variable is 1 for all observations that belong to a certain category, and 0 otherwise. Let’s make a design matrix for this example:
model.matrix(~ treat)## (Intercept) treattrt
## 1 1 0
## 2 1 0
## 3 1 0
## 4 1 0
## 5 1 0
## 6 1 1
## 7 1 1
## 8 1 1
## 9 1 1
## 10 1 1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$treat
## [1] "contr.treatment"Note that we still have the intercept, and the treattrt column is 1 for all samples in treat group trt, and zero for all samples not in group trt. Also note that we don’t have a column which is 1 only for group ctrl - this would be overparametrizing the model, since that column would be identical to the (Intercept) column minus the treattrt column.
Fitting a linear model again allows us to estimate the coefficients:
summary(lm(y ~ treat))##
## Call:
## lm(formula = y ~ treat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2207 -0.5878 -0.2622 0.3629 1.7594
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1936 0.4485 2.661 0.0287 *
## treattrt 1.7621 0.6343 2.778 0.0240 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.003 on 8 degrees of freedom
## Multiple R-squared: 0.491, Adjusted R-squared: 0.4274
## F-statistic: 7.718 on 1 and 8 DF, p-value: 0.02399So how do we interpret this? Recall from above that the fitted/predicted value of \(y\) can be obtained by multiplying the design matrix by the estimated coefficients:
\[\begin{pmatrix}\hat{y}_1\\\hat{y}_2\\...\\\hat{y}_n\end{pmatrix} = \begin{pmatrix}1 & 0\\1 & 0\\... & ...\\1 & 1\end{pmatrix}\begin{pmatrix}\widehat{(Intercept)}\\\widehat{treattrt}\end{pmatrix}\]
So similarly to what we did for the linear regression, the fitted values for all samples in group ctrl (the first five), will be ‘1 times the intercept + 0 times the coefficient for treattrt’, i.e., just the intercept. The fitted value for all samples in group trt is ‘1 times the intercept + 1 times the coefficient for treattrt’, i.e., intercept + coefficient for treattrt. Thus, the coefficient for treattrt in this case directly represents the difference between the fitted values for the two groups, and the intercept represents the fitted value for the samples in group ctrl. We can again explore the design interactively using ExploreModelMatrix:
if (interactive()) {
ExploreModelMatrix(sampleData = data.frame(y = y, treat = treat),
designFormula = ~ treat)
}
vd <- VisualizeDesign(sampleData = data.frame(y = y, treat = treat),
designFormula = ~ treat)
cowplot::plot_grid(plotlist = vd$plotlist)
There are other ways of parametrizing this model - if we want the coefficients to directly give us the fitted values for each group, we can exclude the intercept:
model.matrix(~ 0 + treat)## treatctrl treattrt
## 1 1 0
## 2 1 0
## 3 1 0
## 4 1 0
## 5 1 0
## 6 0 1
## 7 0 1
## 8 0 1
## 9 0 1
## 10 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$treat
## [1] "contr.treatment"summary(lm(y ~ 0 + treat))##
## Call:
## lm(formula = y ~ 0 + treat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2207 -0.5878 -0.2622 0.3629 1.7594
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## treatctrl 1.1936 0.4485 2.661 0.028748 *
## treattrt 2.9557 0.4485 6.590 0.000171 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.003 on 8 degrees of freedom
## Multiple R-squared: 0.8633, Adjusted R-squared: 0.8291
## F-statistic: 25.26 on 2 and 8 DF, p-value: 0.0003494if (interactive()) {
ExploreModelMatrix(sampleData = data.frame(y = y, treat = treat),
designFormula = ~ 0 + treat)
}
vd <- VisualizeDesign(sampleData = data.frame(y = y, treat = treat),
designFormula = ~ 0 + treat)
cowplot::plot_grid(plotlist = vd$plotlist)
Both models are equally ‘valid’ - the main difference appears when we want to interpret the coefficients or perform statistical tests. For example, to perform a test comparing groups ctrl and trt with the first approach, we would directly test whether the coefficient treattrt is different from zero. In the second case, we would have to test whether the difference between the treattrt and treatctrl coefficients is different from zero (such linear combinations of coefficients are usually referred to as contrasts).
3 What if we have two explanatory variables?
Multiple explanatory variables, assumed to have additive effects on the response, can be acccommodated as well:
(gt <- rep(c("wt", "mut"), 5))## [1] "wt" "mut" "wt" "mut" "wt" "mut" "wt" "mut" "wt" "mut"model.matrix(~ gt + treat)## (Intercept) gtwt treattrt
## 1 1 1 0
## 2 1 0 0
## 3 1 1 0
## 4 1 0 0
## 5 1 1 0
## 6 1 0 1
## 7 1 1 1
## 8 1 0 1
## 9 1 1 1
## 10 1 0 1
## attr(,"assign")
## [1] 0 1 2
## attr(,"contrasts")
## attr(,"contrasts")$gt
## [1] "contr.treatment"
##
## attr(,"contrasts")$treat
## [1] "contr.treatment"Note how the intercept absorbs the reference level of each of the predictors.
if (interactive()) {
ExploreModelMatrix(sampleData = data.frame(y = y, gt = gt, treat = treat),
designFormula = ~ gt + treat)
}
vd <- VisualizeDesign(sampleData = data.frame(y = y, gt = gt, treat = treat),
designFormula = ~ gt + treat)
cowplot::plot_grid(plotlist = vd$plotlist)
Finally, if we don’t want to assume that the effects of the two predictors are additive, we can use an interaction term (effectively allowing the effect of one of the predictors to depend on the value of another):
model.matrix(~ gt * treat) ## equivalent to ~ gt + treat + gt:treat## (Intercept) gtwt treattrt gtwt:treattrt
## 1 1 1 0 0
## 2 1 0 0 0
## 3 1 1 0 0
## 4 1 0 0 0
## 5 1 1 0 0
## 6 1 0 1 0
## 7 1 1 1 1
## 8 1 0 1 0
## 9 1 1 1 1
## 10 1 0 1 0
## attr(,"assign")
## [1] 0 1 2 3
## attr(,"contrasts")
## attr(,"contrasts")$gt
## [1] "contr.treatment"
##
## attr(,"contrasts")$treat
## [1] "contr.treatment"if (interactive()) {
ExploreModelMatrix(sampleData = data.frame(y = y, gt = gt, treat = treat),
designFormula = ~ gt * treat)
}
vd <- VisualizeDesign(sampleData = data.frame(y = y, gt = gt, treat = treat),
designFormula = ~ gt * treat)
cowplot::plot_grid(plotlist = vd$plotlist)
4 Exercises
- Take a look at the ExploreModelMatrix vignette and explore some of the designs in there. The ExploreModelMatrix app also provides example designs, which you can load via the ‘Use example design’ dropdown menu. Some background information and a description of the app is found in the accompanying paper.
- Law et al (2020) provides an excellent guide to creating design matrices in R.