The goal of this script is to help understand matrix algebra and notation in linear and generalized linear models.
Imagine there are \(n=10\) data points and you want to fit a linear regression model and predict \(\hat{y}\) from the estimated coeffients (\(\boldsymbol{\beta}\)) and covariates with the design matrix (\(\textbf{X}\)).
We will ignore the estimation of coeffients and just pretend that we know what they are.
Lets setup the coeffients as one intercept and two slopes (\(\boldsymbol{\beta}\))
betas = matrix(c(2,1,5),nrow=3)
betas## [,1]
## [1,] 2
## [2,] 1
## [3,] 5Next, create a design matrix that has one intercept and two continuous covariates
set.seed(34524254)
X = matrix(cbind(rep(1,10),rnorm(10,2,1),rpois(10,10)),nrow=10)
head(X)## [,1] [,2] [,3]
## [1,] 1 2.91595970 9
## [2,] 1 -0.05838379 7
## [3,] 1 0.10943754 9
## [4,] 1 2.60919391 8
## [5,] 1 3.75040827 12
## [6,] 1 2.39877244 11If I want to multiple the coeffients by the covariates and intercept to predict \(\hat{y}\), how do I relate \(\textbf{X}\) and (\(\boldsymbol{\beta}\))?
First, what are the dimensions of the matrices…
dim(X)## [1] 10 3 dim(betas)## [1] 3 1Can we multiply them together and what does that even mean?
X*betas## Error in `X * betas`:
## ! non-conformable arrays betas*X## Error in `betas * X`:
## ! non-conformable arraysNo, we can’t do that.
Rather, we need to do matrix multiplication (called a dot product (\(\cdot\))), as
X%*%betas## [,1]
## [1,] 49.91596
## [2,] 36.94162
## [3,] 47.10944
## [4,] 44.60919
## [5,] 65.75041
## [6,] 59.39877
## [7,] 68.79571
## [8,] 64.08380
## [9,] 45.07889
## [10,] 55.24953Notice that the dimensions of the output is the number of samples, n = 10
But what is the dot product exactly doing?
Lets take one row at a time of (\(\boldsymbol{\beta}\)). Start with the first row
X[1,]## [1] 1.00000 2.91596 9.00000We want to multiply betas[1] by 1, and betas[2] by 2.91596 and betas[3] by 9
X[1,]*betas## [,1]
## [1,] 2.00000
## [2,] 2.91596
## [3,] 45.00000And then we want to sum them…
sum(X[1,]*betas)## [1] 49.91596This is the same as the first element of
( X%*%betas)[1]## [1] 49.91596Remember what the \(\boldsymbol{\beta}\) are. They are the partial effects of a known variable or covariate. The coeffients in columns 2 and 3 are slopes of a line or magnitude of effect that all together sum (linear approximation) to relate to the data, \(y\).
Go ahead and calculate manually the matrix multiplication of row 7 of \(\textbf{X}\) and the \(\boldsymbol{\beta}\)
Is your answer the same as (X%*%betas)[7]?
Summary: So X%*%betas is multiplying all the elements of coeffients by each row of \(\textbf{X}\) and then summing it.
Transpose is an important matrix operation. It takes a matrix and turns it on its side.
X## [,1] [,2] [,3]
## [1,] 1 2.91595970 9
## [2,] 1 -0.05838379 7
## [3,] 1 0.10943754 9
## [4,] 1 2.60919391 8
## [5,] 1 3.75040827 12
## [6,] 1 2.39877244 11
## [7,] 1 1.79571183 13
## [8,] 1 2.08380215 12
## [9,] 1 3.07889390 8
## [10,] 1 3.24953278 10 t(X)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1.00000 1.00000000 1.0000000 1.000000 1.000000 1.000000 1.000000
## [2,] 2.91596 -0.05838379 0.1094375 2.609194 3.750408 2.398772 1.795712
## [3,] 9.00000 7.00000000 9.0000000 8.000000 12.000000 11.000000 13.000000
## [,8] [,9] [,10]
## [1,] 1.000000 1.000000 1.000000
## [2,] 2.083802 3.078894 3.249533
## [3,] 12.000000 8.000000 10.000000Note the dimension change
dim(X)## [1] 10 3 dim(t(X))## [1] 3 10What about betas?
dim(betas)## [1] 3 1 dim(t(betas))## [1] 1 3If you multiply two matrices together, you should know the dimensions of the output matrix before the multiplication.
X%*%t(X)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 90.50282 63.82976 82.31912 80.60830 119.93604 106.99472 123.23622
## [2,] 63.82976 50.00341 63.99361 56.84767 84.78104 77.85995 91.89516
## [3,] 82.31912 63.99361 82.01198 73.28554 109.41044 100.26252 118.19652
## [4,] 80.60830 56.84767 73.28554 71.80789 106.78554 95.25886 109.68536
## [5,] 119.93604 84.78104 109.41044 106.78554 159.06556 141.99638 163.73465
## [6,] 106.99472 77.85995 100.26252 95.25886 141.99638 127.75411 148.30750
## [7,] 123.23622 91.89516 118.19652 109.68536 163.73465 148.30750 173.22458
## [8,] 115.07628 84.87834 109.22805 102.43704 152.81511 137.99857 160.74191
## [9,] 81.97793 56.82024 73.33695 73.03343 108.54711 96.38557 110.52881
## [10,] 100.47551 70.81028 91.35562 89.47866 133.18707 118.79489 136.83522
## [,8] [,9] [,10]
## [1,] 115.07628 81.97793 100.47551
## [2,] 84.87834 56.82024 70.81028
## [3,] 109.22805 73.33695 91.35562
## [4,] 102.43704 73.03343 89.47866
## [5,] 152.81511 108.54711 133.18707
## [6,] 137.99857 96.38557 118.79489
## [7,] 160.74191 110.52881 136.83522
## [8,] 149.34223 103.41581 127.77138
## [9,] 103.41581 74.47959 91.00497
## [10,] 127.77138 91.00497 111.55946\(\textbf{X}\) by \(\textbf{X}^{t}\) is matrix multiplying a 10x3 matrix by a 3x10 matrix.
We know we can multiply them because the columns of the first matrix (3) is the same number as the rows of the second matrix (3). The matrix that is produced has the dimensions of rows of the first matrix (10) and the number of columns of the second matrix (10). Thus, the outcome is a 10x10 matrix
dim(X)## [1] 10 3 dim(t(X))## [1] 3 10 dim(X%*%t(X))## [1] 10 10We can not multiple matrices that don’t align in terms of the dimensions
X%*%X## Error in `X %*% X`:
## ! non-conformable argumentsWhy is X%*%X non-conformable?
Can we always multiply a square matrix by itself? Try it.
The Identity matrix is a special matrix, in which the diagonals are all one, and all other elements (off diagonals) are zeros.
diag(10)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 0 0 0 0 0 0 0 0 0
## [2,] 0 1 0 0 0 0 0 0 0 0
## [3,] 0 0 1 0 0 0 0 0 0 0
## [4,] 0 0 0 1 0 0 0 0 0 0
## [5,] 0 0 0 0 1 0 0 0 0 0
## [6,] 0 0 0 0 0 1 0 0 0 0
## [7,] 0 0 0 0 0 0 1 0 0 0
## [8,] 0 0 0 0 0 0 0 1 0 0
## [9,] 0 0 0 0 0 0 0 0 1 0
## [10,] 0 0 0 0 0 0 0 0 0 1 diag(3)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1What is the outcome of multiplying \(\boldsymbol{\beta}\) by an identity matrix
diag(3)%*%(betas)## [,1]
## [1,] 2
## [2,] 1
## [3,] 5Nothing changed!
In GLM’s an identity link is the identity matrix- basically, no transformation. In contrast to a log-transformation for the Poisson or the logit for the Binomial.
If I multiply a matrix that is 8x10 times a matrix that is 8x10, what is the dimension of the produced matrix? See if you can figure it out.
If I multiply a matrix that is 9x3 times a matrix that is 3x1000, what is the dimension of the produced matrix? See if you can figure it out.
Can we multiply a matrix by a scalar?
scalar=3
diag(3)*scalar## [,1] [,2] [,3]
## [1,] 3 0 0
## [2,] 0 3 0
## [3,] 0 0 3# OR, how about
X*100## [,1] [,2] [,3]
## [1,] 100 291.595970 900
## [2,] 100 -5.838379 700
## [3,] 100 10.943754 900
## [4,] 100 260.919391 800
## [5,] 100 375.040827 1200
## [6,] 100 239.877244 1100
## [7,] 100 179.571183 1300
## [8,] 100 208.380215 1200
## [9,] 100 307.889390 800
## [10,] 100 324.953278 1000Yes, we can do that, no problem.
Can we divide, add, and substract by a scalar?
diag(3)/scalar## [,1] [,2] [,3]
## [1,] 0.3333333 0.0000000 0.0000000
## [2,] 0.0000000 0.3333333 0.0000000
## [3,] 0.0000000 0.0000000 0.3333333#OR
diag(3)+1000## [,1] [,2] [,3]
## [1,] 1001 1000 1000
## [2,] 1000 1001 1000
## [3,] 1000 1000 1001 diag(3)-1000## [,1] [,2] [,3]
## [1,] -999 -1000 -1000
## [2,] -1000 -999 -1000
## [3,] -1000 -1000 -999Yes, not a problem
Can we do these in a different order?
3/diag(3)## [,1] [,2] [,3]
## [1,] 3 Inf Inf
## [2,] Inf 3 Inf
## [3,] Inf Inf 3 1000+diag(3)## [,1] [,2] [,3]
## [1,] 1001 1000 1000
## [2,] 1000 1001 1000
## [3,] 1000 1000 1001 1000-diag(3)## [,1] [,2] [,3]
## [1,] 999 1000 1000
## [2,] 1000 999 1000
## [3,] 1000 1000 999Yes, but the order of operations matters, so the output is different. We also shouldn’t divide by zero!
Does it matter what order of operations we do our matrix multiplication?
A=matrix(c(10,2,1,5),nrow=2,ncol=2)
B=matrix(c(1,5,2,10),nrow=2,ncol=2)
A## [,1] [,2]
## [1,] 10 1
## [2,] 2 5 B## [,1] [,2]
## [1,] 1 2
## [2,] 5 10 A%*%B## [,1] [,2]
## [1,] 15 30
## [2,] 27 54 B%*%A## [,1] [,2]
## [1,] 14 11
## [2,] 70 55Are the output matrices the same? No!
We can’t assume \(A\cdot B\) is equal to \(B\cdot A\) when dealing with matrices, which we can assume with regular multiplication: \(3\times2 = 2\times3\)
What if we don’t have continiuos variables, but categorical?
Imagine a categorical variable (habtat) with three levels: grassland, forest, shrub
We want to relate elk density (\(y\); our data) to habitat with these three levels. We thus need a design matrix (\(\textbf{X}\)) that has 3 columns, to translate which row of data correspond to each level.
Lets consider 12 elk densities in which 4 of each correspond to each of the habitat levels
Here is the variable
habitat=factor(c(rep("Forest",4),rep("Grass",4),rep("Shrub",4)))
habitat## [1] Forest Forest Forest Forest Grass Grass Grass Grass Shrub Shrub
## [11] Shrub Shrub
## Levels: Forest Grass ShrubLet’s turn this into a design matrix
X=model.matrix(~habitat)
X## (Intercept) habitatGrass habitatShrub
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 1 0 0
## 5 1 1 0
## 6 1 1 0
## 7 1 1 0
## 8 1 1 0
## 9 1 0 1
## 10 1 0 1
## 11 1 0 1
## 12 1 0 1
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$habitat
## [1] "contr.treatment"This transformation is a really really important part of understanding all linear models.
Notice that the first column is the intercept, the second and third columns are called partial intercepts or dummy codes to indicate the different levels.
For each row, there is a unique combination of 1’s and 0’s that indicate which is forest, grassland, and shrub. The first row is 1 0 0, and we know looking at ‘habitat’ that the first element is ‘forest’. The 5th row is 1 1 0, which from ‘habitat’, we know is grassland. Lastly, row 9 is 1 0 1, which we can see is Shrub.
How do we interpret the intercept?
The intercept alone ( 1 0 0 ), indicates Forest because no matter what the other coefficient values are we will always only be adding 1beta1 + 0beta2 + 0beta3, which is equal to just 1beta1. Thus, the intercept and beta1 indicate the density of elk in forest.
What does the 2nd column indicate? The beta2 corresponds to the difference in elk density between forest and grassland, so if grassland has less elk than forests beta2 will be negative, and if grassland has more elk than forests it will be positive. If there is no difference, we would expect beta2 to be zero.
The same is true for column three. beta3 is the difference from the intercept (forest) because anytime we multiply the betas by the rows of X, we only ever add the intercept and one of the other habitat types.
habitat2=relevel(habitat, ref="Shrub")
X2=model.matrix(~habitat2)
X2## (Intercept) habitat2Forest habitat2Grass
## 1 1 1 0
## 2 1 1 0
## 3 1 1 0
## 4 1 1 0
## 5 1 0 1
## 6 1 0 1
## 7 1 0 1
## 8 1 0 1
## 9 1 0 0
## 10 1 0 0
## 11 1 0 0
## 12 1 0 0
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$habitat2
## [1] "contr.treatment"Notice how the design matrix has changed. Now, we have the 1 0 0 be the last part of the data, thus the intercept corresponds to “Shrub”,
Let use a set of betas and get the prediction of y (elk density) using the new design matrix
betas2=c(0.5,1,2)
X2%*%betas2## [,1]
## 1 1.5
## 2 1.5
## 3 1.5
## 4 1.5
## 5 2.5
## 6 2.5
## 7 2.5
## 8 2.5
## 9 0.5
## 10 0.5
## 11 0.5
## 12 0.5By changing the design matrix ( (\(\textbf{X}\))) you can manipulate what betas are estimate and what they mean, but not necessarily the predictions of y-hat. Below is an example.
#This is our data…
y=c(3, 3, 3, 3, 1.5, 1.5, 1.5, 1.5, 1, 1, 1, 1)Let’s relate our data to two different betas and sets of design matrices
y.hat.1=t(X%*%c(3,-1.5,-2))
y.hat.2=t(X2%*%c(1,2,0.5))
y==y.hat.1## 1 2 3 4 5 6 7 8 9 10 11 12
## [1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE y==y.hat.2## 1 2 3 4 5 6 7 8 9 10 11 12
## [1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE y.hat.1==y.hat.2## 1 2 3 4 5 6 7 8 9 10 11 12
## [1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEIs it an issue that our betas are not unique solutions, but can change depending on our design matrix?
No, not at all. We just need to be certain that we understand the design matrix so that we can understand the meaning of out betas. But, irregardless of their meaning, we still predict the same y-hat’s.
When using a design matrix (\(\textbf{X}\)) what is the definition of the intercept (beta1)? This indicates the elk desnity in Forests
When using design matrix (\(\textbf{X2}\)) what is the definition of the intercept (beta1)? This indicates the elk desnity in Shrubs
When using design matrix (\(\textbf{X}\)) what is the definition of slope 1 (beta2)? It is the difference in elk density between forests and …..
When using design matrix (\(\textbf{X2}\)) what is the definition of slope 1 (beta2)? It is the difference in elk density between shrubs and …..
When using design matrix (\(\textbf{X}\)) what is the definition of slope 2 (beta3)?
When using design matrix (\(\textbf{X2}\)) what is the definition of slope 2 (beta3)?
How does this all relate to linear modeling and generalized linear modeling?
Does using variable habitat or habitat2 change the inference on the coefficients? Does the variable habitat or habitat2 change the predictions from the model?
Lets fit a linear regression model with the variable habitat
fit1=lm(y~habitat)
fit1##
## Call:
## lm(formula = y ~ habitat)
##
## Coefficients:
## (Intercept) habitatGrass habitatShrub
## 3.0 -1.5 -2.0 fit1$model## y habitat
## 1 3.0 Forest
## 2 3.0 Forest
## 3 3.0 Forest
## 4 3.0 Forest
## 5 1.5 Grass
## 6 1.5 Grass
## 7 1.5 Grass
## 8 1.5 Grass
## 9 1.0 Shrub
## 10 1.0 Shrub
## 11 1.0 Shrub
## 12 1.0 ShrubDo the predictions match y?
predict(fit1)## 1 2 3 4 5 6 7 8 9 10 11 12
## 3.0 3.0 3.0 3.0 1.5 1.5 1.5 1.5 1.0 1.0 1.0 1.0 y## [1] 3.0 3.0 3.0 3.0 1.5 1.5 1.5 1.5 1.0 1.0 1.0 1.0#Now use habitat2 instead of habitat. Did the coefficients change? What about the predictions?
fit2=lm(y~habitat2)
fit2##
## Call:
## lm(formula = y ~ habitat2)
##
## Coefficients:
## (Intercept) habitat2Forest habitat2Grass
## 1.0 2.0 0.5 predict(fit2)## 1 2 3 4 5 6 7 8 9 10 11 12
## 3.0 3.0 3.0 3.0 1.5 1.5 1.5 1.5 1.0 1.0 1.0 1.0 y## [1] 3.0 3.0 3.0 3.0 1.5 1.5 1.5 1.5 1.0 1.0 1.0 1.0Instead of doing y~habitat, lets specify the design matrix we want to use, X1 and X2.
If we fit the y data with the design matrix, do they match with above? Meaning do we get the same results by using habitat or by using X? #Do we also get the same results by using habitat2 or by using X2?
lm(y~X+0)##
## Call:
## lm(formula = y ~ X + 0)
##
## Coefficients:
## X(Intercept) XhabitatGrass XhabitatShrub
## 3.0 -1.5 -2.0 lm(y~X2+0)##
## Call:
## lm(formula = y ~ X2 + 0)
##
## Coefficients:
## X2(Intercept) X2habitat2Forest X2habitat2Grass
## 1.0 2.0 0.5Do the same as above but for glm..
fit3=glm(y~habitat,family=gaussian)
fit3##
## Call: glm(formula = y ~ habitat, family = gaussian)
##
## Coefficients:
## (Intercept) habitatGrass habitatShrub
## 3.0 -1.5 -2.0
##
## Degrees of Freedom: 11 Total (i.e. Null); 9 Residual
## Null Deviance: 8.667
## Residual Deviance: 1.42e-29 AIC: -784.9 predict(fit3)## 1 2 3 4 5 6 7 8 9 10 11 12
## 3.0 3.0 3.0 3.0 1.5 1.5 1.5 1.5 1.0 1.0 1.0 1.0 y## [1] 3.0 3.0 3.0 3.0 1.5 1.5 1.5 1.5 1.0 1.0 1.0 1.0Same results as y~habitat?
glm(y~X+0,family=gaussian)##
## Call: glm(formula = y ~ X + 0, family = gaussian)
##
## Coefficients:
## X(Intercept) XhabitatGrass XhabitatShrub
## 3.0 -1.5 -2.0
##
## Degrees of Freedom: 12 Total (i.e. Null); 9 Residual
## Null Deviance: 49
## Residual Deviance: 1.42e-29 AIC: -784.9 fit4=glm(y~habitat2,family=gaussian)
fit4##
## Call: glm(formula = y ~ habitat2, family = gaussian)
##
## Coefficients:
## (Intercept) habitat2Forest habitat2Grass
## 1.0 2.0 0.5
##
## Degrees of Freedom: 11 Total (i.e. Null); 9 Residual
## Null Deviance: 8.667
## Residual Deviance: 8.677e-30 AIC: -790.8 predict(fit4)## 1 2 3 4 5 6 7 8 9 10 11 12
## 3.0 3.0 3.0 3.0 1.5 1.5 1.5 1.5 1.0 1.0 1.0 1.0 y## [1] 3.0 3.0 3.0 3.0 1.5 1.5 1.5 1.5 1.0 1.0 1.0 1.0Same results as y~habitat2?
glm(y~X2+0,family=gaussian)##
## Call: glm(formula = y ~ X2 + 0, family = gaussian)
##
## Coefficients:
## X2(Intercept) X2habitat2Forest X2habitat2Grass
## 1.0 2.0 0.5
##
## Degrees of Freedom: 12 Total (i.e. Null); 9 Residual
## Null Deviance: 49
## Residual Deviance: 8.677e-30 AIC: -790.8