The Matrix
The matrix function: R wants the data to be entered by columns starting with column one
1st arg: c(2,3,-2,1,2,2) the values of the elements filling the columns
2nd arg: 3 the number of rows
3rd arg: 2 the number of columns
Define matrix A
A <- matrix(c(2, 3, -2, 1, 2, 2), 3, 2) A
## [,1] [,2] ## [1,] 2 1 ## [2,] 3 2 ## [3,] -2 2
Is something a matrix
is.matrix(A)
## [1] TRUE
Is something a vector
is.vector(A)
## [1] FALSE
Multiplication by a scalar
Multiply each element of the matrix by the scalar.
c <- 3 c * A
## [,1] [,2] ## [1,] 6 3 ## [2,] 9 6 ## [3,] -6 6
Matrix addition & subtraction
Add (or subtract) element by element.
Matrix addition and subtract require the matrices to have the same dimensions.
A
## [,1] [,2] ## [1,] 2 1 ## [2,] 3 2 ## [3,] -2 2
B <- matrix(c(1, 4, -2, 1, 2, 1), 3, 2) B
## [,1] [,2] ## [1,] 1 1 ## [2,] 4 2 ## [3,] -2 1
C <- A + B C
## [,1] [,2] ## [1,] 3 2 ## [2,] 7 4 ## [3,] -4 3
D <- A - B D
## [,1] [,2] ## [1,] 1 0 ## [2,] -1 0 ## [3,] 0 1
Matrix multiplication
Matrix multiplication proceeds row by column. Matrices must be conformable, i.e., the number of colums of the first matrix must equal the number of rows of the second matrix. The dimensions of the product will have the same number of rows as the first matrix and the same number of columns as the second matrix.
A(8,4)*B(4,3) = C(8,3) Mathematical notation: AD or A*D. R notation: A%*%
F <- matrix(c(3, 9, 7, 1), 4, 1) F
## [,1] ## [1,] 3 ## [2,] 9 ## [3,] 7 ## [4,] 1
G <- matrix(c(6, 1, 3, 6), 1, 4) G
## [,1] [,2] [,3] [,4] ## [1,] 6 1 3 6
G %*% F
## [,1] ## [1,] 54
F %*% G
## [,1] [,2] [,3] [,4] ## [1,] 18 3 9 18 ## [2,] 54 9 27 54 ## [3,] 42 7 21 42 ## [4,] 6 1 3 6
D <- matrix(c(2, -2, 1, 2, 3, 1), 2, 3) D
## [,1] [,2] [,3] ## [1,] 2 1 3 ## [2,] -2 2 1
C <- D %*% A C
## [,1] [,2] ## [1,] 1 10 ## [2,] 0 4
C <- A %*% D C
## [,1] [,2] [,3] ## [1,] 2 4 7 ## [2,] 2 7 11 ## [3,] -8 2 -4
D <- matrix(c(2, 1, 3), 1, 3) D
## [,1] [,2] [,3] ## [1,] 2 1 3
C <- D %*% A C
## [,1] [,2] ## [1,] 1 10
# non-conformable error C <- A %*% D
## Error: non-conformable arguments
Transpose of a matrix
Mathematical notation: A´ or A';. R notation: t(A) function.
Note: The trasnpose of a transpose yields the original matrix.
AT <- t(A) AT
## [,1] [,2] [,3] ## [1,] 2 3 -2 ## [2,] 1 2 2
ATT <- t(AT) ATT
## [,1] [,2] ## [1,] 2 1 ## [2,] 3 2 ## [3,] -2 2
Common vectors
Unit vector
U <- matrix(1, 3, 1) U
## [,1] ## [1,] 1 ## [2,] 1 ## [3,] 1
Zero vector
Z <- matrix(0, 3, 1) Z
## [,1] ## [1,] 0 ## [2,] 0 ## [3,] 0
Common matrices
Unit matrix
U <- matrix(1, 3, 2) U
## [,1] [,2] ## [1,] 1 1 ## [2,] 1 1 ## [3,] 1 1
Zero matrix
Z <- matrix(0, 3, 2) Z
## [,1] [,2] ## [1,] 0 0 ## [2,] 0 0 ## [3,] 0 0
Diagonal matrix
S <- matrix(c(2, 3, -2, 1, 2, 2, 4, 2, 3), 3, 3) S
## [,1] [,2] [,3] ## [1,] 2 1 4 ## [2,] 3 2 2 ## [3,] -2 2 3
D <- diag(S) D
## [1] 2 2 3
D <- diag(diag(S)) D
## [,1] [,2] [,3] ## [1,] 2 0 0 ## [2,] 0 2 0 ## [3,] 0 0 3
Identity matrix
A diagonal matrix of ones.
Note: A*I = A; I*B = B;
I <- diag(3) I
## [,1] [,2] [,3] ## [1,] 1 0 0 ## [2,] 0 1 0 ## [3,] 0 0 1
Symmetric matrix
Note: If matrix A is symmetric then A' = a.;
C <- matrix(c(2, 1, 5, 1, 3, 4, 5, 4, -2), 3, 3) C
## [,1] [,2] [,3] ## [1,] 2 1 5 ## [2,] 1 3 4 ## [3,] 5 4 -2
CT <- t(C) CT
## [,1] [,2] [,3] ## [1,] 2 1 5 ## [2,] 1 3 4 ## [3,] 5 4 -2
Inverse of a matrix
The matrix analog of the scalar reciprocal.
Note: A-1*A = A*A-1 = I
A <- matrix(c(4, 4, -2, 2, 6, 2, 2, 8, 4), 3, 3) A
## [,1] [,2] [,3] ## [1,] 4 2 2 ## [2,] 4 6 8 ## [3,] -2 2 4
AI <- solve(A) AI
## [,1] [,2] [,3] ## [1,] 1.0 -0.5 0.5 ## [2,] -4.0 2.5 -3.0 ## [3,] 2.5 -1.5 2.0
A %*% AI
## [,1] [,2] [,3] ## [1,] 1 0 0 ## [2,] 0 1 0 ## [3,] 0 0 1
AI %*% A
## [,1] [,2] [,3] ## [1,] 1 0 0 ## [2,] 0 1 0 ## [3,] 0 0 1
Determinant of a matrix
C <- matrix(c(2, 1, 6, 1, 3, 4, 6, 4, -2), 3, 3) C
## [,1] [,2] [,3] ## [1,] 2 1 6 ## [2,] 1 3 4 ## [3,] 6 4 -2
d <- det(C) d
## [1] -102
Rank of a matrix
The column rank of a matrix is the number of linearly independent columns.
A <- matrix(c(2, 3, -2, 1, 2, 2, 4, 7, 0), 3, 3) A
## [,1] [,2] [,3] ## [1,] 2 1 4 ## [2,] 3 2 7 ## [3,] -2 2 0
matA <- qr(A) matA$rank
## [1] 3
A <- matrix(c(2, 3, -2, 1, 2, 2, 4, 6, -4), 3, 3) A
## [,1] [,2] [,3] ## [1,] 2 1 4 ## [2,] 3 2 6 ## [3,] -2 2 -4
matA <- qr(A) matA$rank
## [1] 2
# note column 3 is 2 times column 1
Number of rows & columns
X <- matrix(c(3, 2, 4, 3, 2, -2, 6, 1), 4, 2) X
## [,1] [,2] ## [1,] 3 2 ## [2,] 2 -2 ## [3,] 4 6 ## [4,] 3 1
dim(X)
## [1] 4 2
r <- nrow(X) r
## [1] 4
c <- ncol(X) c
## [1] 2
Computing column & row sums
# note the uppercase S A <- matrix(c(2, 3, -2, 1, 2, 2), 3, 2) A
## [,1] [,2] ## [1,] 2 1 ## [2,] 3 2 ## [3,] -2 2
c <- colSums(A) c
## [1] 3 5
r <- rowSums(A) r
## [1] 3 5 0
# sum of all tthe elements a <- sum(A) a
## [1] 8
Computing column & row means
# note the uppercase M cm <- colMeans(A) cm
## [1] 1.000 1.667
rm <- rowMeans(A) rm
## [1] 1.5 2.5 0.0
# mean of all the elments m <- mean(A) m
## [1] 1.333
Horizontal concatenation
The two matrices to be concatenated must have the same number of rows.
A
## [,1] [,2] ## [1,] 2 1 ## [2,] 3 2 ## [3,] -2 2
B <- matrix(c(1, 3, 2, 1, 4, 2), 3, 2) B
## [,1] [,2] ## [1,] 1 1 ## [2,] 3 4 ## [3,] 2 2
C <- cbind(A, B) C
## [,1] [,2] [,3] [,4] ## [1,] 2 1 1 1 ## [2,] 3 2 3 4 ## [3,] -2 2 2 2
Vertical concatenation (appending)
The two matrices to be appended must have the same number of columns.
C <- rbind(A, B) C
## [,1] [,2] ## [1,] 2 1 ## [2,] 3 2 ## [3,] -2 2 ## [4,] 1 1 ## [5,] 3 4 ## [6,] 2 2
Matrix examples
Compute column means manually, aka compute mean value vector
B
## [,1] [,2] ## [1,] 1 1 ## [2,] 3 4 ## [3,] 2 2
o <- matrix(1, 1, nrow(B)) o
## [,1] [,2] [,3] ## [1,] 1 1 1
# column totals o %*% B
## [,1] [,2] ## [1,] 6 7
# number of rows o %*% t(o)
## [,1] ## [1,] 3
# column means solve(o %*% t(o)) %*% (o %*% B)
## [,1] [,2] ## [1,] 2 2.333
Seven important matrices for data analysis
1) Raw score matrix
# define raw score matrix X X <- matrix(c(25, 20, 19, 26, 23, 10, 9, 10, 11, 7, 55, 53, 50, 60, 57), 5, 3) X
## [,1] [,2] [,3] ## [1,] 25 10 55 ## [2,] 20 9 53 ## [3,] 19 10 50 ## [4,] 26 11 60 ## [5,] 23 7 57
2) Raw score sums of squares and crossproducts
t(X) %*% X
## [,1] [,2] [,3] ## [1,] 2591 1067 6256 ## [2,] 1067 451 2586 ## [3,] 6256 2586 15183
3) Deviation score matrix
# Vector of column means cm <- colMeans(X) cm
## [1] 22.6 9.4 55.0
# Create mean matrix MM <- matrix(cm, 5, 3, byrow = TRUE) MM
## [,1] [,2] [,3] ## [1,] 22.6 9.4 55 ## [2,] 22.6 9.4 55 ## [3,] 22.6 9.4 55 ## [4,] 22.6 9.4 55 ## [5,] 22.6 9.4 55
# Deviation score matrix D <- X - MM D
## [,1] [,2] [,3] ## [1,] 2.4 0.6 0 ## [2,] -2.6 -0.4 -2 ## [3,] -3.6 0.6 -5 ## [4,] 3.4 1.6 5 ## [5,] 0.4 -2.4 2
4) Deviation score sums of squares and cross products
dsscp <- t(D) %*% D dsscp
## [,1] [,2] [,3] ## [1,] 37.2 4.8 41 ## [2,] 4.8 9.2 1 ## [3,] 41.0 1.0 58
5) Covariance matrix
nr <- nrow(D) cov <- 1/(nr - 1) * dsscp cov
## [,1] [,2] [,3] ## [1,] 9.30 1.20 10.25 ## [2,] 1.20 2.30 0.25 ## [3,] 10.25 0.25 14.50
6) Standard score matrix
# Diagonal matrix of variances var <- diag(diag(cov)) var
## [,1] [,2] [,3] ## [1,] 9.3 0.0 0.0 ## [2,] 0.0 2.3 0.0 ## [3,] 0.0 0.0 14.5
# Diagonal matrix of standard deviations sd <- sqrt(var) sd
## [,1] [,2] [,3] ## [1,] 3.05 0.000 0.000 ## [2,] 0.00 1.517 0.000 ## [3,] 0.00 0.000 3.808
# Standard score matrix Z <- D %*% solve(sd) Z
## [,1] [,2] [,3] ## [1,] 0.7870 0.3956 0.0000 ## [2,] -0.8526 -0.2638 -0.5252 ## [3,] -1.1805 0.3956 -1.3131 ## [4,] 1.1149 1.0550 1.3131 ## [5,] 0.1312 -1.5825 0.5252
7) Standard score sums of squares and cross products, aka correlation matrix
1/(nr - 1) * t(Z) %*% Z
## [,1] [,2] [,3] ## [1,] 1.0000 0.25946 0.88267 ## [2,] 0.2595 1.00000 0.04329 ## [3,] 0.8827 0.04329 1.00000
Compute regression coefficients: Regression of Y on X
Mathematical notation: b = (X'X)-1X'Y
# Define Y Y <- matrix(c(30, 25, 20, 25, 28), 5, 1) Y
## [,1] ## [1,] 30 ## [2,] 25 ## [3,] 20 ## [4,] 25 ## [5,] 28
# Create vector of ones o <- matrix(1, 5, 1) o
## [,1] ## [1,] 1 ## [2,] 1 ## [3,] 1 ## [4,] 1 ## [5,] 1
# Concatenate ones and X X <- cbind(o, X) X
## [,1] [,2] [,3] [,4] ## [1,] 1 25 10 55 ## [2,] 1 20 9 53 ## [3,] 1 19 10 50 ## [4,] 1 26 11 60 ## [5,] 1 23 7 57
b <- solve(t(X) %*% X) %*% t(X) %*% Y b
## [,1] ## [1,] 43.2528 ## [2,] 1.9919 ## [3,] -1.7301 ## [4,] -0.8437
Principal Components Analysis and Factor Analysis
# read 5x5 correlation matrix C <- matrix(c(1, 0.5968, 0.6623, 0.6302, 0.6215, 0.5968, 1, 0.6174, 0.5704, 0.6048, 0.6623, 0.6174, 1, 0.6307, 0.5445, 0.6302, 0.5704, 0.6307, 1, 0.4651, 0.6215, 0.6048, 0.5445, 0.4651, 1), 5, 5) # compute eigenvalues and eigenvectors # eigenvalues aka characteristic values or latent values # eigenvectors aka characteristic vectors or latent vectors r <- eigen(C) eval <- r$values evec <- r$vectors # principal components analysis eval
## [1] 3.3808 0.5574 0.4068 0.3562 0.2988
evec
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.4664 -0.02744 -0.530988 -0.02020 0.7066 ## [2,] 0.4484 0.20776 0.806508 0.05575 0.3197 ## [3,] 0.4588 -0.26083 -0.000789 -0.78020 -0.3358 ## [4,] 0.4356 -0.61095 -0.006783 0.58930 -0.2995 ## [5,] 0.4257 0.71749 -0.259906 0.20126 -0.4427
# principal component factor laodings # scale the eigenvectors by the square root of the eigenvalues d <- sqrt(diag(eval)) d
## [,1] [,2] [,3] [,4] [,5] ## [1,] 1.839 0.0000 0.0000 0.0000 0.0000 ## [2,] 0.000 0.7466 0.0000 0.0000 0.0000 ## [3,] 0.000 0.0000 0.6378 0.0000 0.0000 ## [4,] 0.000 0.0000 0.0000 0.5968 0.0000 ## [5,] 0.000 0.0000 0.0000 0.0000 0.5467
pcf_ld <- evec %*% d pcf_ld
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.8576 -0.02048 -0.3386581 -0.01206 0.3863 ## [2,] 0.8244 0.15511 0.5143817 0.03327 0.1748 ## [3,] 0.8435 -0.19473 -0.0005032 -0.46563 -0.1836 ## [4,] 0.8009 -0.45612 -0.0043259 0.35170 -0.1637 ## [5,] 0.7827 0.53566 -0.1657653 0.12011 -0.2420
# principal axes factor analysis # use squared multiple correlations (smc) instead of ones in the diagonal of # the correlation matrix # compute squared smcs: smc = 1 - 1/r^jj invC <- solve(C) smc <- diag(5) - solve(diag(diag(invC))) smc
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.5886 0.0000 0.0000 0.0000 0.0000 ## [2,] 0.0000 0.5208 0.0000 0.0000 0.0000 ## [3,] 0.0000 0.0000 0.5596 0.0000 0.0000 ## [4,] 0.0000 0.0000 0.0000 0.5003 0.0000 ## [5,] 0.0000 0.0000 0.0000 0.0000 0.4772
# replace ones in correlation matrix with smcs # R is the reduced correlation matrix R <- C - diag(5) + smc R
## [,1] [,2] [,3] [,4] [,5] ## [1,] 0.5886 0.5968 0.6623 0.6302 0.6215 ## [2,] 0.5968 0.5208 0.6174 0.5704 0.6048 ## [3,] 0.6623 0.6174 0.5596 0.6307 0.5445 ## [4,] 0.6302 0.5704 0.6307 0.5003 0.4651 ## [5,] 0.6215 0.6048 0.5445 0.4651 0.4772
# compute eigenvalues and eigenvectors of reduced correlation matrix r <- eigen(R) eval <- r$values eval
## [1] 2.91328 0.05254 -0.04770 -0.10451 -0.16715
# only first two eigenvalues are positive # subset eval: first two eigenvalues d <- sqrt(diag(eval[1:2])) d
## [,1] [,2] ## [1,] 1.707 0.0000 ## [2,] 0.000 0.2292
# subset evec: first two eigenvectors evec <- r$vectors[, 1:2] evec
## [,1] [,2] ## [1,] 0.4751 0.02411 ## [2,] 0.4464 -0.28499 ## [3,] 0.4630 0.32391 ## [4,] 0.4312 0.58545 ## [5,] 0.4179 -0.68595
# compute principal axes factor loadings # loadings are correlations of variables with latent factors paf_ld <- evec %*% d paf_ld
## [,1] [,2] ## [1,] 0.8110 0.005525 ## [2,] 0.7620 -0.065325 ## [3,] 0.7903 0.074245 ## [4,] 0.7360 0.134196 ## [5,] 0.7132 -0.157230