## Check the data mtcars with head and save a new data as `mtcars.subset`

after dropping two non-numeric (binary) variables for PCA analysis

```
data <- mtcars
head(data)
```

```
## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
```

`str(data)`

```
## 'data.frame': 32 obs. of 11 variables:
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
## $ disp: num 160 160 108 258 360 ...
## $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
## $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec: num 16.5 17 18.6 19.4 17 ...
## $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
## $ am : num 1 1 1 0 0 0 0 0 0 0 ...
## $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
## $ carb: num 4 4 1 1 2 1 4 2 2 4 ...
```

In our data vs and am are binary variable so I drop them here.

```
library(dplyr)
mtcars.subset <- data[,-c(8,9)]
```

## Fit PCA in the data by `mtcars.pca`

`matcars.subset`

with `cor = TRUE`

and `scores = TRUE`

`mtcars.pca<-prcomp(mtcars.subset, cor = TRUE, scores = TRUE)`

```
## Warning: In prcomp.default(mtcars.subset, cor = TRUE, scores = TRUE) :
## extra arguments 'cor', 'scores' will be disregarded
```

## Get summary of `mtcars.pca`

`summary(mtcars.pca)`

```
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 136.532 38.14735 3.06642 1.27492 0.90474 0.64734 0.3054
## Proportion of Variance 0.927 0.07237 0.00047 0.00008 0.00004 0.00002 0.0000
## Cumulative Proportion 0.927 0.99938 0.99985 0.99993 0.99997 0.99999 1.0000
## PC8 PC9
## Standard deviation 0.2859 0.2159
## Proportion of Variance 0.0000 0.0000
## Cumulative Proportion 1.0000 1.0000
```

From above summary we see when standard deviation is greater propertion

of variance is high similarly when standard deviation is low propertion

of variation is also low.

## Get eigenvalue of the components using standard deviation of `mtcars.pca`

and chose the number of components based on Kaiser’s criteria

`mtcars.pca$sdev ^2`

```
## [1] 1.864106e+04 1.455220e+03 9.402948e+00 1.625431e+00 8.185525e-01
## [6] 4.190430e-01 9.327903e-02 8.175127e-02 4.660443e-02
```

## Get scree plot and chose the number of components best on “first bend” of this plot

```
#Calculating total variance explained by each principal component
var_explained<-mtcars.pca$sdev^2/sum(mtcars.pca$sdev^2)
##Creating scree plot
library(ggplot2)
qplot(c(1:9), var_explained) + geom_line() + xlab("Principal Component") + ylab("Variance explained") + ggtitle("Scree Plot")
```

## How many components must be retained based on Kaiser’s rule and/or scree plot?

Kisher’s rule suggaest us to use 2 components and Scree plot

suggest us to retain 4 components for the problem.

## Fit the final PCA model based on the retained components

`library(psych)`

`## Warning: package 'psych' was built under R version 4.1.2`

```
mtcars.pca<- psych::principal(mtcars.subset, nfactors = 4, rotate = "none")
mtcars.pca
```

```
## Principal Components Analysis
## Call: psych::principal(r = mtcars.subset, nfactors = 4, rotate = "none")
## Standardized loadings (pattern matrix) based upon correlation matrix
## PC1 PC2 PC3 PC4 h2 u2 com
## mpg -0.93 0.04 -0.16 0.00 0.90 0.0995 1.1
## cyl 0.96 0.02 -0.18 0.02 0.95 0.0504 1.1
## disp 0.94 -0.13 -0.06 0.17 0.94 0.0569 1.1
## hp 0.87 0.39 -0.01 0.04 0.91 0.0854 1.4
## drat -0.74 0.49 0.11 0.44 0.99 0.0062 2.5
## wt 0.89 -0.25 0.32 0.10 0.96 0.0360 1.5
## qsec -0.53 -0.70 0.45 -0.02 0.97 0.0283 2.6
## gear -0.50 0.79 0.15 -0.15 0.92 0.0775 1.8
## carb 0.58 0.70 0.33 -0.11 0.95 0.0525 2.5
##
## PC1 PC2 PC3 PC4
## SS loadings 5.66 2.08 0.50 0.27
## Proportion Var 0.63 0.23 0.06 0.03
## Cumulative Var 0.63 0.86 0.92 0.95
## Proportion Explained 0.66 0.24 0.06 0.03
## Cumulative Proportion 0.66 0.91 0.97 1.00
##
## Mean item complexity = 1.7
## Test of the hypothesis that 4 components are sufficient.
##
## The root mean square of the residuals (RMSR) is 0.02
## with the empirical chi square 0.96 with prob < 0.99
##
## Fit based upon off diagonal values = 1
```

## Get the head of the saved loadings of `mtcars.pca`

`head(mtcars.pca)`

```
## $values
## [1] 5.65593947 2.08210029 0.50421482 0.26502753 0.18315864 0.12379319 0.10506192
## [8] 0.05851375 0.02219038
##
## $rotation
## [1] "none"
##
## $n.obs
## [1] 32
##
## $communality
## mpg cyl disp hp drat wt qsec gear
## 0.9004692 0.9495949 0.9430951 0.9146434 0.9938301 0.9639755 0.9716949 0.9224616
## carb
## 0.9475174
##
## $loadings
##
## Loadings:
## PC1 PC2 PC3 PC4
## mpg -0.935 -0.157
## cyl 0.957 -0.179
## disp 0.945 -0.128 0.175
## hp 0.873 0.389
## drat -0.742 0.493 0.106 0.435
## wt 0.888 -0.248 0.322
## qsec -0.534 -0.698 0.446
## gear -0.498 0.795 0.147 -0.145
## carb 0.582 0.699 0.330 -0.110
##
## PC1 PC2 PC3 PC4
## SS loadings 5.656 2.082 0.504 0.265
## Proportion Var 0.628 0.231 0.056 0.029
## Cumulative Var 0.628 0.860 0.916 0.945
##
## $fit
## [1] 0.9982615
```

## Retain two components, get their loadings

```
mtcars.pca_2 <- principal(mtcars.subset, nfactors = 2, rotate="none")
mtcars.pca_2$loadings
```

```
##
## Loadings:
## PC1 PC2
## mpg -0.935
## cyl 0.957
## disp 0.945 -0.128
## hp 0.873 0.389
## drat -0.742 0.493
## wt 0.888 -0.248
## qsec -0.534 -0.698
## gear -0.498 0.795
## carb 0.582 0.699
##
## PC1 PC2
## SS loadings 5.656 2.082
## Proportion Var 0.628 0.231
## Cumulative Var 0.628 0.860
```

Principal components (PCs) are constructed by the linear combination of

the original variables, where PCA loading are the coefficients. Here,

cyl has the weights of 0.957 on PC1 computation but not in PC2. Positive

loading in above data indicates a variable and a component are

positively correlated. Negative loading indicate a negative correlation

between the variable and component. Similarly disp has positive loading

with PC1 and negative loading with PC2. Large (either positive or

negative) loading indicate that a variable has a strong effect on that

principal component. The larger value of cyl indicates the strong effect

on PC1.

## Get biplot of these two component loadings

`biplot(mtcars.pca, col = c("blue", "black"), cex = c(0.5, 1.3))`

## Get the head of the saved scores of `mtcars.pca`

`head(mtcars.pca$scores)`

```
## PC1 PC2 PC3 PC4
## Mazda RX4 -0.27929417 0.8132290 -0.2877380 -0.2447853
## Mazda RX4 Wag -0.26793045 0.6770476 0.1560073 -0.1664252
## Datsun 710 -0.96699807 -0.2263347 -0.2959516 -0.2110014
## Hornet 4 Drive -0.09052843 -1.3699797 -0.4639869 -0.5984019
## Hornet Sportabout 0.66729435 -0.5743299 -1.4547534 0.2862895
## Valiant 0.02085807 -1.6956141 0.1574155 -1.6929588
```

The original ddataset is projected into four principal components.

## Get the head of the scores of first two components of `mtcars.pca`

`head(mtcars.pca_2$scores)`

```
## PC1 PC2
## Mazda RX4 -0.27929417 0.8132290
## Mazda RX4 Wag -0.26793045 0.6770476
## Datsun 710 -0.96699807 -0.2263347
## Hornet 4 Drive -0.09052843 -1.3699797
## Hornet Sportabout 0.66729435 -0.5743299
## Valiant 0.02085807 -1.6956141
```

## Get biplot of these two component scores

`biplot(mtcars.pca_2, col = c("blue", "red"))`

Here we can observe that hp, cyl, disp and wt contribute to PC1 with

higher values. And mpg which has negative loadings is in opposite

direction to PC1 with higher values. Gear and carb has higher

contribution to PC2 with positive values and qsec has negative value.

## Get dissimilar distance of all the variables of `mtcars`

data as `mtcars.dist`

```
#Distance calculation
mtcars.dist<- dist(mtcars.subset)
```

## Fit classical multi-dimensional scaling model with the `mtcars.dist`

in 2-dimensional state as `cars.mds.2d`

```
# Fitting classical two- dimensional scaling model
cars.mds.2d<-cmdscale(mtcars.dist)
summary(cars.mds.2d)
```

```
## V1 V2
## Min. :-181.07 Min. :-139.047
## 1st Qu.:-116.69 1st Qu.: -10.373
## Median : -43.99 Median : 2.144
## Mean : 0.00 Mean : 0.000
## 3rd Qu.: 132.85 3rd Qu.: 29.375
## Max. : 242.81 Max. : 52.503
```

## Plot the `cars.mds.2d`

and compare it with the biplot of `mtcars.pca`

```
plot(cars.mds.2d, pch = 19)
abline(h = 0, v = 0, lty =2)
mtcars.subset<-mtcars[, 1:2] %>% scale
text(cars.mds.2d, pos = 4, labels = rownames(mtcars.subset), col = "tomato")
```

Hornet 4 Drive, Pontiac Fire bird etc lies on the positive orthant which

means they have positive contribution to the first and second

components. However, Lotus Europa and Ferari has opposite but highest

weight component 1 and component 2.

## Fit classical multi-dimensional scaling model with the mtcars.dist in 3-dimensional state as `cars.mds.3d`

```
#Fiting multi-dimensional scaling model with mtcars.dist in 3 - dimensional state
cars.mds.3d<-cmdscale(mtcars.dist, k = 3)
summary(cars.mds.3d)
```

```
## V1 V2 V3
## Min. :-181.07 Min. :-139.047 Min. :-6.8611
## 1st Qu.:-116.69 1st Qu.: -10.373 1st Qu.:-1.8374
## Median : -43.99 Median : 2.144 Median : 0.8492
## Mean : 0.00 Mean : 0.000 Mean : 0.0000
## 3rd Qu.: 132.85 3rd Qu.: 29.375 3rd Qu.: 2.2806
## Max. : 242.81 Max. : 52.503 Max. : 5.0029
```

## Create a 3-d scatterplot of `cars.mds.3d`

with `type = “h”`

, `pch=20`

and `lty.hplot=2`

and interpret it carefully

```
library(scatterplot3d)
cars.mds.3d <- data.frame(cmdscale(mtcars.dist, k = 3))
scatterplot3d(cars.mds.3d, type = "h", pch = 19, lty.hplot = 2)
```

## Create a 3-d scatterplot of `cars.mds.3d`

with `type = “h”, pch=20, lty.hplot=2`

and `color=mtcars$cyl`

```
library(scatterplot3d)
cars.mds.3d <- data.frame(cmdscale(mtcars.dist, k = 3))
scatterplot3d(cars.mds.3d, type = "h", pch = 19, lty.hplot = 2, color = mtcars$cyl)
```

We plotted the principal components in 3- dimensional scatter plot which is

distinguished by color cyl. We can use higher dimensions by changing the

k argument in the cmdscale() function to a higher value for eg. k = 3

for 3 dimension.

## Write a summary comparing PCA and MDS fits done above for `mtcars`

data

The input to PCA is the original vectors in n-dimensional

space. Similarly, input to MDS is the pairwise distances between points.

PCA behaves as an algorithm but MDS is a visualization technique for any

factor analysis. MDS applies PCA for the dimensionality reduction. For

the mtcars, it shows that two or more but less than or equal to 5 latent

features can be generated from the given dataset.

You must log in to post a comment.