Multivariate Data Analaysis using R
Setia Pramana, PhD
October 28, 2014
Matrix Operation in R
Row and column vectors are drawn the same way.
The vector in R printed row format but can really be regarded as a column
vector, cfr the convention above.
> a <- c(2,2)
> a
[1] 2 2
>
>
>
>
>
>
>
>
>
#draw the vector
plot(c(0,2),c(0,2), xlim=c(-4,4), ylim=c(-4,4))
arrows(0,0,2,2)
abline(h=0,v=0)
text(2.3,2.3,"a=(2,2)")
4
0
4
c(0, 2)
a=(2,2)
c(0, 2)
> a <- c(1:3)
> a
[1] 1 2 3
> # transpose
>
> t(a)
[1,]
[,1] [,2] [,3]
1
2
3
> # Multiplication by a number
>
> 8*a
[1]
8 16 24
> #Vector Addition
>
> a <- c(1,3,2)
> b <- c(6,8,1)
> a+b
2
[1]
7 11
> ## length of a Vector (Norm)
>
> sqrt(sum(a*a))
[1] 3.741657
>
>
>
>
>
>
x <- c(-1,5,2,-2)
y <- c(4,-3,0,1)
lx <- sqrt(sum(x^2))
ly <- sqrt(sum(y^2))
cosxy <- 1/lx * 1/ly * sum (x*y)
acos (cosxy)
[1] 2.355063
>
>
>
> A <- matrix(c(1,3,2,4,8,7),ncol=3)
> A
[1,]
[2,]
[,1] [,2] [,3]
1
2
8
3
4
7
> 8*A
[1,]
[2,]
[,1] [,2] [,3]
8
16
64
24
32
56
> t(A)
[1,]
[2,]
[3,]
[,1] [,2]
1
3
2
4
8
7
> ## Addition of matrices
>
> B <- matrix(c(5,8,3,4,2,7),ncol=3,byrow=T)
> A+B
[1,]
[2,]
[,1] [,2] [,3]
6
10
11
7
6
14
3
> ## Multiplication ##
>
> A%*%a
[1,]
[2,]
[,1]
23
29
> # Different with
> A*a
[1,]
[2,]
[,1] [,2] [,3]
1
4
24
9
4
14
> A <- matrix(c(1,3,2,2,8,9),ncol=2)
> B <- matrix(c(5,8,4,2), ncol=2)
> A
[1,]
[2,]
[3,]
[,1] [,2]
1
2
3
8
2
9
> B
[1,]
[2,]
[,1] [,2]
5
4
8
2
> A%*%B
[1,]
[2,]
[3,]
>
>
>
>
>
>
>
[,1] [,2]
21
8
79
28
82
26
# determinant ##
#det(A)
# Error in [Link](x, logarithm = TRUE, ...) :
# 'x' must be a square matrix
det(B)
[1] -22
> D <- matrix(c(5,8,4,2,5,6,7,8,9), ncol=3)
> det(D)
4
[1] 101
> ## DIagonal Matrix ##
> C <- diag(c(1,2,3,5))
> C
[1,]
[2,]
[3,]
[4,]
[,1] [,2] [,3] [,4]
1
0
0
0
0
2
0
0
0
0
3
0
0
0
0
5
> det(C)
[1] 30
> diag(1,3)
[1,]
[2,]
[3,]
[,1] [,2] [,3]
1
0
0
0
1
0
0
0
1
> diag(A)
[1] 1 8
> ## Inverse Matrix ###
>
> A <- matrix(c(1,3,2,4),ncol=2,byrow=T)
> A
[1,]
[2,]
[,1] [,2]
1
3
2
4
> #invers#
> B <- solve(A)
> B
[1,]
[2,]
[,1] [,2]
-2 1.5
1 -0.5
> A%*%B
[1,]
[2,]
[,1] [,2]
1
0
0
1
5
> A <- matrix(c(1,3,2,6),ncol=2)
> A
[1,]
[2,]
>
>
>
>
>
>
>
>
>
>
>
#invers#
#solve(A)
#Error in [Link](A) :
# Lapack routine dgesv: system is exactly singular: U[2,2] = 0
## Solving systems of linear equations ##
A <- matrix(c(1,2,3,4),ncol=2)
A
[1,]
[2,]
>
>
>
[,1] [,2]
1
2
3
6
[,1] [,2]
1
3
2
4
b <- c(7,10)
x <- solve(A)%*%b
x
[1,]
[2,]
[,1]
1
2
>
>
>
> A <- cbind(x1=c(42,52,48,58),x2=c(4,5,4,3))
> meanA <- colMeans(A)
> meanA
x1 x2
50 4
> var(A[,1])
[1] 45.33333
> var(A[,2])
[1] 0.6666667
6
> var(A)
x1
x2
x1 45.33333 -2.0000000
x2 -2.00000 0.6666667
> cor(A)
x1
x2
x1 1.0000000 -0.3638034
x2 -0.3638034 1.0000000
> t(A)-meanA
x1
x2
[,1] [,2] [,3] [,4]
-8
2
-2
8
0
1
0
-1
> a <- (t(A)-meanA)
> s1 <- sum(a[1,]^2)/4
> s1
[1] 34
> s2 <- sum(a[2,]^2)/4
> s2
[1] 0.5
> b <- t(a)
> b[,1]*b[,2]
[1]
0 -8
> s12 <- sum( b[,1]*b[,2] )/4
> s12
[1] -1.5
> r12 <- s12/(sqrt(s1)*sqrt(s2))
> r12
[1] -0.3638034
> cov(A)
x1
x2
x1 45.33333 -2.0000000
x2 -2.00000 0.6666667
7
> ## Visualization
>
> data(iris)
> head(iris)
1
2
3
4
5
6
[Link] [Link] [Link] [Link] Species
5.1
3.5
1.4
0.2 setosa
4.9
3.0
1.4
0.2 setosa
4.7
3.2
1.3
0.2 setosa
4.6
3.1
1.5
0.2 setosa
5.0
3.6
1.4
0.2 setosa
5.4
3.9
1.7
0.4 setosa
> # Plot #1: Basic scatterplot matrix of the four measurements
> pairs(~[Link]+[Link]+[Link]+[Link], data=iris)
>
>
3.0
4.0
0.5
1.5
2.5
6.5
7.5
2.0
4.0
4.5
5.5
[Link]
2.0
3.0
[Link]
1.5
2.5
[Link]
0.5
[Link]
4.5
5.5
6.5
7.5
> [Link] <- function(x, y, digits=2, prefix="", [Link], ...)
+ {
+
usr <- par("usr"); [Link](par(usr))
+
par(usr = c(0, 1, 0, 1))
+
r <- abs(cor(x, y))
8
+
txt <- format(c(r, 0.123456789), digits=digits)[1]
+
txt <- paste(prefix, txt, sep="")
+
if(missing([Link])) [Link] <- 0.8/strwidth(txt)
+
text(0.5, 0.5, txt, cex = [Link] * r)
+ }
> pairs(~[Link]+[Link]+[Link]+[Link], data=iris,
+
[Link]=[Link], [Link]=[Link],
+
pch=20, main=" Scatterplot Matrix")
>
>
Scatterplot Matrix
3.0
4.0
0.5
1.5
2.5
0.87
0.82
0.43
0.37
[Link]
0.96
5.5
0.12
4.0
4.5
[Link]
6.5
7.5
2.0
1.5
2.5
0.5
[Link]
4.5
>
>
>
>
5.5
6.5
7.5
#[Link]("corrplot")
library(corrplot)
M <- cor(mtcars)
corrplot(M, method = "circle")
1 2 3 4 5 6 7
2.0
3.0
[Link]
carb
gear
am
vs
qsec
wt
drat
hp
disp
cyl
mpg
mpg
0.8
cyl
0.6
disp
hp
0.4
drat
0.2
wt
qsec
0.2
vs
0.4
am
0.6
gear
0.8
carb
1
> corrplot(M, method = "ellipse")
10
carb
gear
am
vs
qsec
wt
drat
hp
disp
cyl
mpg
mpg
0.8
cyl
0.6
disp
hp
0.4
drat
0.2
wt
qsec
0.2
vs
0.4
am
0.6
gear
0.8
carb
1
> corrplot(M, type = "upper")
11
carb
gear
am
vs
qsec
wt
drat
hp
disp
cyl
mpg
mpg
0.8
cyl
0.6
disp
0.4
hp
0.2
drat
wt
qsec
0.2
vs
0.4
am
0.6
gear
0.8
carb
1
> corrplot(M, order = "hclust", addrect = 3)
12
gear
am
drat
mpg
vs
qsec
disp
cyl
hp
wt
carb
carb
0.8
wt
0.6
hp
cyl
0.4
disp
0.2
qsec
vs
0.2
mpg
0.4
drat
0.6
am
0.8
gear
1
>
>
>
>
>
>
>
>
>
### Larger Data
X1 <- rnorm(n=200,mean=0,sd=1)
X2 <- rnorm(n=200,mean=0,sd=1)
X3 <- rnorm(n=200,mean=0,sd=1)
hist(X1, prob=TRUE)
lines(density(X1))
dt <- cbind(X1,X2,X3)
colMeans(dt)
X1
0.035007409
X2
X3
0.001373853 -0.063052048
> #rowMeans(dt)
> cov(dt)
X1
X2
X3
X1 0.90863072 -0.1002754 -0.06591799
X2 -0.10027537 1.0821860 -0.10273225
X3 -0.06591799 -0.1027322 1.04013075
> cor(dt)
13
X1
X2
X3
X1 1.00000000 -0.10112292 -0.06780571
X2 -0.10112292 1.00000000 -0.09683043
X3 -0.06780571 -0.09683043 1.00000000
> dt2 <- cbind(X1,X2=-X1,X3=2*X1)
> cov(dt2)
X1
X2
X3
X1 0.9086307 -0.9086307 1.817261
X2 -0.9086307 0.9086307 -1.817261
X3 1.8172614 -1.8172614 3.634523
> cor(dt2)
X1 X2 X3
X1 1 -1 1
X2 -1 1 -1
X3 1 -1 1
>
>
>
0.2
0.1
0.0
Density
0.3
0.4
Histogram of X1
X1
14