MULTIVARIATE ANALYSIS USING R
B N Mandal
I.A.S.R.I., Library Avenue, New Delhi – 110 012
bnmandal @[Link]
1. Introduction
This article gives an exposition of how to use the R statistical software for multivariate
analyses, with a focus on principal components analysis (PCA) and cluster analysis. We
will use the Iris data set of Fisher for all purposes unless otherwise specified.
2. Reading Data into R
There are a number of ways of reading data into R. For our purpose, we will use
[Link]() function to read the data into R.
>iris=[Link](“d:/[Link]”, header=T)
Here, the data is loaded in R as iris. After reading the data in R, attach the data so that
variables can be directly accessed.
>attach(data)
3. Plotting Multivariate Data
Once you have read a multivariate data set into R, the next step is usually to make a plot
of the data.
A matrix scatter plot
One common way of plotting multivariate data is to make a “matrix scatterplot”, showing
each pair of variables plotted against each other. We can use the “scatterplotMatrix()”
function from the “car” R package to do this. To use this function, first install the “car” R
package. Once installed the “car” R package, load the “car” R package by using:
> library("car")
Now, you can then use the “scatterplotMatrix()” function to plot the multivariate data.
Alternative, one can use pairs() function.
>pairs(iris[,1:4])
The output is like
Multivariate Analysis using R
One can see the correlations of four variables here.
4. Calculating Summary Statistics for Multivariate Data
One may be interested to calculate summary statistics such as the mean and standard
deviation for each of the variables in the multivariate data set.
This can be done using the “mean()”, “sd()”, “cv()”, “var()”functions in R. For example,
to calculate the mean and standard deviations of each of four variables in the iris data set
> mean(Sepal_length)
5.843333
>var(Sepal_length)
0.6856935
Similarly use for other variables.
This tells that the mean of variable Sepal_length is 5.843333 and the variance is
0.6856935.
Multivariate Analysis using R
Summary() gives some commonly used descriptive statistics of the variables in the data
set. For example,
> summary(iris)
gives the following output.
Sepal_length Sepal_width Petal_length Petal_width Species
Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100 I. setosa :50
1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300 I. versicolor:50
Median :5.800 Median :3.000 Median :4.350 Median :1.300 I. virginica :50
Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
5. Principal Component Analysis
The purpose of principal component analysis is to find the best low-dimensional
representation of the variation in a multivariate data set. To carry out a principal
component analysis (PCA) on a multivariate data set, first standardize the variables under
study using the “scale()” function. This is necessary if the input variables have very
different variances.
Once you have standardized the variables, carry out a principal component analysis using
the “prcomp()” function in R.
For example, to carry out a principal components analysis on the standardized variables,
use:
>pc1=prcomp(iris[,1:4],scale=T)
Here, scale option in the function prcomp scales all the variable in the data set such that
variables have zero mean and variance 1.
Standard deviations:
[1] 1.7083611 0.9560494 0.3830886 0.1439265
Rotation:
PC1 PC2 PC3 PC4
Sepal_length 0.5210659 -0.37741762 0.7195664 0.2612863
Sepal_width -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal_length 0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal_width 0.5648565 -0.06694199 -0.6342727 0.5235971
Multivariate Analysis using R
A summary of the principal component analysis results can be obtained using the
“summary()” function on the output of “prcomp()”.
> summary(pc1)
Importance of components:
PC1 PC2 PC3 PC4
Standard deviation 1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion 0.7296 0.9581 0.99482 1.00000
This gives us the standard deviation of each component, and the proportion of variance
explained by each component. The standard deviation of the components is stored in a
named element called “sdev” of the output variable made by “prcomp”:
> pc1$sdev
1.7083611 0.9560494 0.3830886 0.1439265
In order to decide how many principal components should be retained, it is common to
summarize the results of a principal components analysis by making a scree plot, which
can be done in R using the “screeplot()” function:
> screeplot(pc1, type="lines")
Multivariate Analysis using R
From above, it may be seen that either first 2 or 3 components should be retained.
Another way to decide how many principal components to retain is to decide to keep the
number of components required to explain at least some minimum amount of the total
variance. For example, if it is important to explain at least 90% of the variance, we would
retain the first two principal components, which can be seen from the output of
“summary(pc1)” that the first two principal components explain 95.81% of the variance.
Loadings of the principal components
The loadings for the principal components are stored in a named element “rotation” of the
variable returned by “prcomp()”. This contains a matrix with the loadings of each
principal component, where the first column in the matrix contains the loadings for the
first principal component, the second column contains the loadings for the second
principal component, and so on.
Therefore, to obtain the loadings for the first principal component, use:
> pc1$rotation
PC1 PC2 PC3 PC4
Sepal_length 0.5210659 -0.37741762 0.7195664 0.2612863
Sepal_width -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal_length 0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal_width 0.5648565 -0.06694199 -0.6342727 0.5235971
This means that the first principal component is a linear combination of the variables:
0.52Z1 –0.27Z2 + 0.58Z3+0.56Z4 where Z1, Z2, Z3...Z4 are the standardized versions of
the variables in the data set (that each have mean of 0 and variance of 1).
Note that the square of the loadings sum to 1, as this is a constraint used in calculating the
loadings:
> sum((pc1$rotation[,1])^2)
[1] 1
6. Cluster Analysis
R has a number of functions for cluster analysis. In this section, k-means and hierarchical
agglomerative clustering using R are discussed.
Data Preparation
Prior to clustering data, you may want to remove or estimate missing data and rescale
variables for comparability.
# Prepare Data
>iris = [Link](iris) # listwise deletion of missing
>iris1=iris[,1:4]
>iris1= scale(iris1) # standardize variables
Multivariate Analysis using R
k-means clustering
K-means clustering is the most popular partitioning method. It requires the analyst to
specify the number of clusters to extract. Use following functions to perform a K-means
clustering.
>fit = kmeans(iris1, 5) # 5 cluster solution
To get the cluster means, use
>aggregate(iris1,by=list(fit$cluster),FUN=mean)
Group.1 Sepal_length Sepal_width Petal_length Petal_width
1 1 -1.1924784 0.4015443 -1.3090939 -1.26089024
2 2 -0.3516137 -1.3285553 0.1026061 0.01228268
3 3 -0.6259564 1.8042613 -1.2826445 -1.22905673
4 4 1.3926646 0.2323817 1.1567451 1.21327591
5 5 0.3804044 -0.3896455 0.6067908 0.56390985
# append cluster assignment
>iris1 <- [Link](iris1, fit$cluster)
Hierarchical Agglomerative
There are a wide range of hierarchical clustering approaches. To perform hierarchical
clustering the function hclust () function can be used with a number of options. For
example, Ward's method of clustering using R can be implemented as
# Ward Hierarchical Clustering
>d = dist(iris1, method = "euclidean") # distance matrix
>fit = hclust(d, method="ward")
>plot(fit) # display dendogram
Multivariate Analysis using R
>groups =cutree(fit, k=5) # cut tree into 5 clusters
# draw dendogram with red borders around the 5 clusters
>[Link](fit, k=5, border="blue")
From the above dendogram, five clusters separated by blue boundary lines can be
identified.