Installing FactoMineR

The FactoMineR package can be installed and loaded as follow:

# Install
install.packages("FactoMineR")

# Load
library("FactoMineR")

Installing factoextra

• factoextra can be installed from CRAN as follow:

install.packages("factoextra")

• Or, install the latest developmental version from Github

if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/factoextra")

• Load factoextra as follow :

library("factoextra")

Main functions in FactoMineR

Functions for computing principal component methods and clustering:

Main functions in factoextra

factoextra functions covered in this book are listed in the table below. See the online documentation (https://www.sthda.com/english/rpkgs/factoextra) for a complete list.

• Visualizing principal component method outputs

• Extracting data from principal component method outputs. The following functions extract all the results (coordinates, squared cosine, contributions) for the active individuals/variables from the analysis outputs.

• Clustering analysis and visualization

R functions

• cmdscale() [stats package]: Compute classical (metric) multidimensional scaling.
• isoMDS() [MASS package]: Compute Kruskal’s non-metric multidimensional scaling (one form of non-metric MDS).
• sammon() [MASS package]: Compute sammon’s non-linear mapping (one form of non-metric MDS).

All these functions take a distance object as the main argument and k is the desired number of dimensions in the scaled output. By default, they return two dimension solutions, but we can change that through the parameter k which defaults to 2.

Demo data

swiss data that contains fertility and socio-economic data on 47 French speaking provinces in Switzerland.

data("swiss")
head(swiss)

##              Fertility Agriculture Examination Education Catholic
## Courtelary        80.2        17.0          15        12     9.96
## Delemont          83.1        45.1           6         9    84.84
## Franches-Mnt      92.5        39.7           5         5    93.40
## Moutier           85.8        36.5          12         7    33.77
## Neuveville        76.9        43.5          17        15     5.16
## Porrentruy        76.1        35.3           9         7    90.57
##              Infant.Mortality
## Courtelary               22.2
## Delemont                 22.2
## Franches-Mnt             20.2
## Moutier                  20.3
## Neuveville               20.6
## Porrentruy               26.6

Classical MDS

# Load required packages
library(magrittr)
library(dplyr)
library(ggpubr)

# Cmpute MDS
mds <- swiss %>%
  dist() %>%          
  cmdscale() %>%
  as_tibble()
colnames(mds) <- c("Dim.1", "Dim.2")

# Plot MDS
ggscatter(mds, x = "Dim.1", y = "Dim.2", 
          label = rownames(swiss),
          size = 1,
          repel = TRUE)

Create 3 groups using k-means clustering. Color points by groups

# K-means clustering
clust <- kmeans(mds, 3)$cluster %>%
  as.factor()

mds <- mds %>%
  mutate(groups = clust)

# Plot and color by groups
ggscatter(mds, x = "Dim.1", y = "Dim.2", 
          label = rownames(swiss),
          color = "groups",
          palette = "jco",
          size = 1, 
          ellipse = TRUE,
          ellipse.type = "convex",
          repel = TRUE)

Non-metric MDS

• Load general packages:

library(magrittr)
library(dplyr)
library(ggpubr)

• Kruskal’s non-metric multidimensional scaling

# Cmpute MDS
library(MASS)
mds <- swiss %>%
  dist() %>%          
  isoMDS() %>%
  .$points %>%
  as_tibble()
colnames(mds) <- c("Dim.1", "Dim.2")

# Plot MDS
ggscatter(mds, x = "Dim.1", y = "Dim.2", 
          label = rownames(swiss),
          size = 1,
          repel = TRUE)

• Sammon’s non-linear mapping:

# Cmpute MDS
library(MASS)
mds <- swiss %>%
  dist() %>%          
  sammon() %>%
  .$points %>%
  as_tibble()
colnames(mds) <- c("Dim.1", "Dim.2")

# Plot MDS
ggscatter(mds, x = "Dim.1", y = "Dim.2", 
          label = rownames(swiss),
          size = 1,
          repel = TRUE)

Row profiles

The profile of a given row is calculated by taking each row point and dividing by its margin (i.e, the sum of all row points). The formula is:

\[ row.profile = \frac{row}{row.sum} \]

For example the profile of the row point Laundry/wife is P = 156/176 = 88.6%.

The R code below can be used to compute row profiles:

row.profile <- housetasks/row.sum
# head(row.profile)

The average row profile is computed as follow:

\[ average.rp = \frac{column.sum}{grand.total} \]

For example, the average row profile is : (600/1744, 254/1744, 381/1744, 509/1744). It can be computed in R as follow:

# Column sums
col.sum <- apply(housetasks, 2, sum)
# average row profile = Column sums / grand total
average.rp <- col.sum/n 
average.rp

##        Wife Alternating     Husband     Jointly 
##       0.344       0.146       0.218       0.292

Distance (or similarity) between row profiles

If we want to compare 2 rows (row1 and row2), we need to compute the squared distance between their profiles as follow:

\[ d^2(row_1, row_2) = \sum{\frac{(row.profile_1 - row.profile_2)^2}{average.profile}} \]

For example the distance between the rows Laundry and Main_meal are:

\[ d^2(Laundry, Main\_meal) = \frac{(0.886-0.810)^2}{0.344} + \frac{(0.0795-0.131)^2}{0.146} + ... = 0.036 \]

The distance between Laundry and Main_meal can be calculated as follow in R:

# Laundry and Main_meal profiles
laundry.p <- row.profile["Laundry",]
main_meal.p <- row.profile["Main_meal",]
# Distance between Laundry and Main_meal
d2 <- sum(((laundry.p - main_meal.p)^2) / average.rp)
d2

## [1] 0.0368

The distance between Laundry and Driving is:

# Driving profile
driving.p <- row.profile["Driving",]
# Distance between Laundry and Driving
d2 <- sum(((laundry.p - driving.p)^2) / average.rp)
d2

## [1] 3.77

Squared distance between each row profile and the average row profile

\[ d^2(row_i, average.profile) = \sum{\frac{(row.profile_i - average.profile)^2}{average.profile}} \]

The R code below computes the distance from the average profile for all the row variables:

d2.row <- apply(row.profile, 1, 
        function(row.p, av.p){sum(((row.p - av.p)^2)/av.p)}, 
        average.rp)
as.matrix(round(d2.row,3))

##             [,1]
## Laundry    1.329
## Main_meal  1.034
## Dinner     0.618
## Breakfeast 0.512
## Tidying    0.353
## Dishes     0.302
## Shopping   0.218
## Official   0.968
## Driving    1.274
## Finances   0.456
## Insurance  0.727
## Repairs    3.307
## Holidays   2.140

Distance matrix

In this section the squared distance is computed between each row profile and the other rows in the contingency table.

The result is a distance matrix (a kind of correlation or dissimilarity matrix).

The custom R function below is used to compute the distance matrix:

## data: a data frame or matrix; 
## average.profile: average profile
dist.matrix <- function(data, average.profile){
   mat <- as.matrix(t(data))
    n <- ncol(mat)
    dist.mat<- matrix(NA, n, n)
    diag(dist.mat) <- 0
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            d2 <- sum(((mat[, i] - mat[, j])^2) / average.profile)
            dist.mat[i, j] <- dist.mat[j, i] <- d2
        }
    }
  colnames(dist.mat) <- rownames(dist.mat) <- colnames(mat)
  dist.mat
}

Compute and visualize the distance between row profiles. The package corrplot is required for the visualization. It can be installed as follow: install.packages("corrplot").

# Distance matrix
dist.mat <- dist.matrix(row.profile, average.rp)
dist.mat <-round(dist.mat, 2)
# Visualize the matrix
library("corrplot")
corrplot(dist.mat, type="upper",  is.corr = FALSE)

Row mass and inertia

The Row mass (or row weight) is the total frequency of a given row. It’s calculated as follow:

\[ row.mass = \frac{row.sum}{grand.total} \]

row.sum <- apply(housetasks, 1, sum)
grand.total <- sum(housetasks)
row.mass <- row.sum/grand.total
head(row.mass)

##    Laundry  Main_meal     Dinner Breakfeast    Tidying     Dishes 
##     0.1009     0.0877     0.0619     0.0803     0.0700     0.0648

The row inertia is calculated as the row mass multiplied by the squared distance between the row and the average row profile:

\[ row.inertia = row.mass * d^2(row) \]

# Row inertia
row.inertia <- row.mass * d2.row
head(row.inertia)

##    Laundry  Main_meal     Dinner Breakfeast    Tidying     Dishes 
##     0.1342     0.0907     0.0382     0.0411     0.0247     0.0196

# Total inertia
sum(row.inertia)

## [1] 1.11

Row summary

The result for rows can be summarized as follow:

row <- cbind.data.frame(d2 = d2.row, mass = row.mass, inertia = row.inertia)
round(row,3)

##               d2  mass inertia
## Laundry    1.329 0.101   0.134
## Main_meal  1.034 0.088   0.091
## Dinner     0.618 0.062   0.038
## Breakfeast 0.512 0.080   0.041
## Tidying    0.353 0.070   0.025
## Dishes     0.302 0.065   0.020
## Shopping   0.218 0.069   0.015
## Official   0.968 0.055   0.053
## Driving    1.274 0.080   0.102
## Finances   0.456 0.065   0.030
## Insurance  0.727 0.080   0.058
## Repairs    3.307 0.095   0.313
## Holidays   2.140 0.092   0.196

Column profiles

These are calculated in the same way as the row profiles table.

The profile of a given column is computed as follow:

\[ col.profile = \frac{col}{col.sum} \]

The R code below can be used to compute column profile:

col.profile <- t(housetasks)/col.sum
col.profile <- as.data.frame(t(col.profile))
# head(col.profile)

The average column profile is calculated as follow:

\[ average.cp = row.sum/grand.total \]

For example, the average column profile is : (176/1744, 153/1744, 108/1744, 140/1744, …). It can be computed in R as follow:

# Row sums
row.sum <- apply(housetasks, 1, sum)
# average column profile= row sums/grand total
average.cp <- row.sum/n 
head(average.cp)

##    Laundry  Main_meal     Dinner Breakfeast    Tidying     Dishes 
##     0.1009     0.0877     0.0619     0.0803     0.0700     0.0648

Distance (similarity) between column profiles

If we want to compare columns, we need to compute the squared distance between their profiles as follow:

\[ d^2(col_1, col_2) = \sum{\frac{(col.profile_1 - col.profile_2)^2}{average.profile}} \]

For example the distance between the columns Wife and Husband are:

\[ d^2(Wife, Husband) = \frac{(0.26-0.005)^2}{0.10} + \frac{(0.21-0.013)^2}{0.09} + ... + ... = 4.05 \]

The distance between Wife and Husband can be calculated as follow in R:

# Wife and Husband profiles
wife.p <- col.profile[, "Wife"]
husband.p <- col.profile[, "Husband"]
# Distance between Wife and Husband
d2 <- sum(((wife.p - husband.p)^2) / average.cp)
d2

## [1] 4.05

Squared distance between each column profile and the average column profile

\[ d^2(col_i, average.profile) = \sum{\frac{(col.profile_i - average.profile)^2}{average.profile}} \]

The R code below computes the distance from the average profile for all the column variables

d2.col <- apply(col.profile, 2, 
        function(col.p, av.p){sum(((col.p - av.p)^2)/av.p)}, 
        average.cp)
round(d2.col,3)

##        Wife Alternating     Husband     Jointly 
##       0.875       0.809       1.746       1.078

Distance matrix

# Distance matrix
dist.mat <- dist.matrix(t(col.profile), average.cp)
dist.mat <-round(dist.mat, 2)
dist.mat

##             Wife Alternating Husband Jointly
## Wife        0.00        1.71    4.05    2.93
## Alternating 1.71        0.00    2.67    2.58
## Husband     4.05        2.67    0.00    3.70
## Jointly     2.93        2.58    3.70    0.00

# Visualize the matrix
library("corrplot")
corrplot(dist.mat, type="upper", order="hclust", is.corr = FALSE)

column mass and inertia

The column mass(or column weight) is the total frequency of each column. It’s calculated as follow:

\[ col.mass = \frac{col.sum}{grand.total} \]

col.sum <- apply(housetasks, 2, sum)
grand.total <- sum(housetasks)
col.mass <- col.sum/grand.total
head(col.mass)

##        Wife Alternating     Husband     Jointly 
##       0.344       0.146       0.218       0.292

The column inertia is calculated as the column mass multiplied by the squared distance between the column and the average column profile:

\[ col.inertia = col.mass * d^2(col) \]

col.inertia <- col.mass * d2.col
head(col.inertia)

##        Wife Alternating     Husband     Jointly 
##       0.301       0.118       0.381       0.315

# total inertia
sum(col.inertia)

## [1] 1.11

Column summary

The result for rows can be summarized as follow:

col <- cbind.data.frame(d2 = d2.col, mass = col.mass, 
                        inertia = col.inertia)
round(col,3)

##                d2  mass inertia
## Wife        0.875 0.344   0.301
## Alternating 0.809 0.146   0.118
## Husband     1.746 0.218   0.381
## Jointly     1.078 0.292   0.315

Chi-square test

Chi-square test issued to examine whether rows and columns of a contingency table are statistically significantly associated.

• Null hypothesis (H0): the row and the column variables of the contingency table are independent.
• Alternative hypothesis (H1): row and column variables are dependent

For each cell of the table, we have to calculate the expected value under null hypothesis.

For a given cell, the expected value is calculated as follow:

\[ e = \frac{row.sum * col.sum}{grand.total} \]

The Chi-square statistic is calculated as follow:

\[ \chi^2 = \sum{\frac{(o - e)^2}{e}} \]

• o is the observed value
• e is the expected value

This calculated Chi-square statistic is compared to the critical value (obtained from statistical tables) with \(df = (r - 1)(c - 1)\) degrees of freedom and p = 0.05.

• r is the number of rows in the contingency table
• c is the number of column in the contingency table

If the calculated Chi-square statistic is greater than the critical value, then we must conclude that the row and the column variables are not independent of each other. This implies that they are significantly associated.

Chi-square statistic can be easily computed using the function chisq.test() as follow:

chisq <- chisq.test(housetasks)
chisq

## 
##  Pearson's Chi-squared test
## 
## data:  housetasks
## X-squared = 2000, df = 40, p-value <2e-16

In our example, the row and the column variables are statistically significantly associated (p-value = 0)

The observed and the expected counts can be extracted from the result of the test as follow:

# Observed counts
chisq$observed

##            Wife Alternating Husband Jointly
## Laundry     156          14       2       4
## Main_meal   124          20       5       4
## Dinner       77          11       7      13
## Breakfeast   82          36      15       7
## Tidying      53          11       1      57
## Dishes       32          24       4      53
## Shopping     33          23       9      55
## Official     12          46      23      15
## Driving      10          51      75       3
## Finances     13          13      21      66
## Insurance     8           1      53      77
## Repairs       0           3     160       2
## Holidays      0           1       6     153

# Expected counts
round(chisq$expected,2)

##            Wife Alternating Husband Jointly
## Laundry    60.5        25.6    38.5    51.4
## Main_meal  52.6        22.3    33.4    44.6
## Dinner     37.2        15.7    23.6    31.5
## Breakfeast 48.2        20.4    30.6    40.9
## Tidying    42.0        17.8    26.6    35.6
## Dishes     38.9        16.5    24.7    33.0
## Shopping   41.3        17.5    26.2    35.0
## Official   33.0        14.0    21.0    28.0
## Driving    47.8        20.2    30.4    40.6
## Finances   38.9        16.5    24.7    33.0
## Insurance  47.8        20.2    30.4    40.6
## Repairs    56.8        24.0    36.0    48.2
## Holidays   55.0        23.3    35.0    46.7

As mentioned above the Chi-square statistic is 1944.456.

If you want to know the most contributing cells to the total Chi-square score, you just have to calculate the Chi-square statistic for each cell:

\[ r = \frac{o - e}{\sqrt{e}} \]

The above formula returns the so-called Pearson residuals (r) for each cell (or standardized residuals). Cells with the highest absolute standardized residuals contribute the most to the total Chi-square score.

Pearson residuals can be easily extracted from the output of the function chisq.test():

round(chisq$residuals, 3)

##             Wife Alternating Husband Jointly
## Laundry    12.27      -2.298  -5.878   -6.61
## Main_meal   9.84      -0.484  -4.917   -6.08
## Dinner      6.54      -1.192  -3.416   -3.30
## Breakfeast  4.88       3.457  -2.818   -5.30
## Tidying     1.70      -1.606  -4.969    3.58
## Dishes     -1.10       1.859  -4.163    3.49
## Shopping   -1.29       1.321  -3.362    3.38
## Official   -3.66       8.563   0.443   -2.46
## Driving    -5.47       6.836   8.100   -5.90
## Finances   -4.15      -0.852  -0.742    5.75
## Insurance  -5.76      -4.277   4.107    5.72
## Repairs    -7.53      -4.290  20.646   -6.65
## Holidays   -7.42      -4.620  -4.897   15.56

Let’s visualize Pearson residuals using the package corrplot:

library(corrplot)
corrplot(chisq$residuals, is.cor = FALSE)

For a given cell, the size of the circle is proportional to the amount of the cell contribution.

The sign of the standardized residuals is also very important to interpret the association between rows and columns as explained in the block below:

1. Positive residuals are in blue. Positive values in cells specify an attraction (positive association) between the corresponding row and column variables.
   • In the image above, it’s evident that there are an association between the column Wife and, the rows Laundry and Main_meal.
   • There is a strong positive association between the column Husband and the row Repair
2. Negative residuals are in red. This implies a repulsion (negative association) between the corresponding row and column variables. For example the column Wife are negatively associated (~ “not associated”) with the row Repairs. There is a repulsion between the column Husband and, the rows Laundry and Main_meal

The contribution (in %) of a given cell to the total Chi-square score is calculated as follow:

\[ contrib = \frac{r^2}{\chi^2} \]

• r is the residual of the cell

# Contibution in percentage (%)
contrib <- 100*chisq$residuals^2/chisq$statistic
round(contrib, 3)

##             Wife Alternating Husband Jointly
## Laundry    7.738       0.272   1.777   2.246
## Main_meal  4.976       0.012   1.243   1.903
## Dinner     2.197       0.073   0.600   0.560
## Breakfeast 1.222       0.615   0.408   1.443
## Tidying    0.149       0.133   1.270   0.661
## Dishes     0.063       0.178   0.891   0.625
## Shopping   0.085       0.090   0.581   0.586
## Official   0.688       3.771   0.010   0.311
## Driving    1.538       2.403   3.374   1.789
## Finances   0.886       0.037   0.028   1.700
## Insurance  1.705       0.941   0.868   1.683
## Repairs    2.919       0.947  21.921   2.275
## Holidays   2.831       1.098   1.233  12.445

# Visualize the contribution
corrplot(contrib, is.cor = FALSE)

The relative contribution of each cell to the total Chi-square score give some indication of the nature of the dependency between rows and columns of the contingency table.

It can be seen that:

1. The column “Wife” is strongly associated with Laundry, Main_meal, Dinner
2. The column “Husband” is strongly associated with the row Repairs
3. The column jointly is frequently associated with the row Holidays

Total inertia

As mentioned above, the total inertia is the amount of the information contained in the data table.

It’s called \(\phi^2\) (squared phi) and is calculated as follow:

\[ \phi^2 = \frac{\chi^2}{grand.total} \]

phi2 <- as.numeric(chisq$statistic/sum(housetasks))
phi2

## [1] 1.11

The square root of \(\phi^2\) are called trace and may be interpreted as a correlation coefficient (Bendixen, 2003). Any value of the trace > 0.2 indicates a significant dependency between rows and columns (Bendixen M., 2003)

Mosaic plot

Mosaic plot is used to visualize a contingency table in order to examine the association between categorical variables.

The function mosaicplot() [graphics package] can be used.

library("graphics")
# Mosaic plot of observed values
mosaicplot(housetasks,  las=2, col="steelblue",
           main = "housetasks - observed counts")

# Mosaic plot of expected values
mosaicplot(chisq$expected,  las=2, col = "gray",
           main = "housetasks - expected counts")

In these plots, column variables are firstly spitted (vertical split) and then row variables are spited(horizontal split). For each cell, the height of bars is proportional to the observed relative frequency it contains:

\[ \frac{cell.value}{column.sum} \]

The blue plot, is the mosaic plot of the observed values. The gray one is the mosaic plot of the expected values under null hypothesis.

It’s also possible to color the mosaic plot according to the value of the standardized residuals:

mosaicplot(housetasks, shade = TRUE,   # Color the graph
           las = 2,                    # produces vertical labels
           main = "housetasks")

• This plot clearly show you that Laundry, Main_meal, Dinner and Breakfeast are more often done by the “Wife”.
• Repairs are done by the Husband

G-test: Likelihood ratio test

The G–test of independence is an alternative to the chi-square test of independence, and they will give approximately the same conclusion.

The test is based on the likelihood ratio defined as follow:

\[ ratio = \frac{o}{e} \]

• o is the observed value
• e is the expected value under null hypothesis

This likelihood ratio, or its logarithm, can be used to compute a p-value. When the logarithm of the likelihood ratio is used, the statistic is known as a log-likelihood ratio statistic.

This test is called G-test (or likelihood ratio test or maximum likelihood statistical significance test) and can be used in situations where Chi-square tests were previously recommended.

The G-test is generally defined as follow:

\[ G = 2 * \sum{o * log(\frac{o}{e})} \]

• o is the observed frequency in a cell
• e is the expected frequency under the null hypothesis
• log is the natural logarithm
• The sum is taken over all non-empty cells.

The distribution of G is approximately a chi-squared distribution, with the same number of degrees of freedom as in the corresponding chi-squared test:

\[df = (r - 1)(c - 1)\]

• r is the number of rows in the contingency table
• c is the number of column in the contingency table

The commonly used Pearson Chi-square test is, in fact, just an approximation of the log-likelihood ratio on which the G-tests are based.

Remember that, the Chi-square formula is:

\[ \chi^2 = \sum{\frac{(o - e)^2}{e}} \]

The functions likelihood.test() [Deducer package] or G.test() [RVAideMemoire] can be used to perform a G-test on a contingency table.

We’ll use the package RVAideMemoire which can be installed as follow : install.packages("RVAideMemoire").

The function G.test() works as chisq.test():

library("RVAideMemoire")
gtest <- G.test(as.matrix(housetasks))
gtest

## 
##  G-test
## 
## data:  as.matrix(housetasks)
## G = 2000, df = 40, p-value <2e-16

Interpret the association between rows and columns

To interpret the association between the rows and the columns of the contingency table, the likelihood ratio can be used as an index (i):

\[ ratio = \frac{o}{e} \]

For a given cell,

• If ratio > 1, there is an “attraction” (association) between the corresponding column and row
• If ratio < 1, there is a “repulsion” between the corresponding column and row

The ratio can be calculated as follow:

ratio <- chisq$observed/chisq$expected
round(ratio,3)

##             Wife Alternating Husband Jointly
## Laundry    2.576       0.546   0.052   0.078
## Main_meal  2.356       0.898   0.150   0.090
## Dinner     2.072       0.699   0.297   0.412
## Breakfeast 1.702       1.766   0.490   0.171
## Tidying    1.263       0.619   0.038   1.601
## Dishes     0.823       1.458   0.162   1.607
## Shopping   0.799       1.316   0.343   1.570
## Official   0.363       3.290   1.097   0.535
## Driving    0.209       2.519   2.470   0.074
## Finances   0.334       0.790   0.851   2.001
## Insurance  0.167       0.049   1.745   1.898
## Repairs    0.000       0.125   4.439   0.042
## Holidays   0.000       0.043   0.172   3.276

Note that, you can also use the R code : gtest\$observed/gtest\$expected

The package corrplot can be used to make a graph of the likelihood ratio:

corrplot(ratio, is.cor = FALSE)

The image above confirms our previous observations:

• The rows Laundry, Main_meal and Dinner are associated with the column Wife
• Repairs are done more often by the Husband
• Holidays are taken Jointly

Let’s take the log(ratio) to see the attraction and the repulsion in different colors:

• If ratio < 1 => log(ratio) < 0 (negative values) => red color
• If ratio > 1 = > log(ratio) > 0 (positive values) => blue color

We’ll also add a small value (0.5) to all cells to avoid log(0):

corrplot(log2(ratio + 0.5), is.cor = FALSE)

Eigenvalues and screeplot

# Eigenvalues
eig <- sv^2
# Variances in percentage
variance <- eig*100/sum(eig)
# Cumulative variances
cumvar <- cumsum(variance)

eig<- data.frame(eig = eig, variance = variance,
                     cumvariance = cumvar)
head(eig)

##     eig variance cumvariance
## 1 0.543     48.7        48.7
## 2 0.445     39.9        88.6
## 3 0.127     11.4       100.0

barplot(eig[, 2], names.arg=1:nrow(eig), 
       main = "Variances",
       xlab = "Dimensions",
       ylab = "Percentage of variances",
       col ="steelblue")
# Add connected line segments to the plot
lines(x = 1:nrow(eig), eig[, 2], 
      type="b", pch=19, col = "red")

How many dimensions to retain?:

1. The maximum number of axes in the CA is :

\[ nb.axes = min( r-1, c-1) \]

r and c are respectively the number of rows and columns in the table.

2. Use elbow method

Row coordinates

We can use the function apply to perform arbitrary operations on the rows and columns of a matrix.

A simplified format is:

apply(X, MARGIN, FUN, ...)

• x: a matrix
• MARGIN: allowed values can be 1 or 2. 1 specifies that we want to operate on the rows of the matrix. 2 specifies that we want to operate on the column.
• FUN: the function to be applied
• ...: optional arguments to FUN

# row sum
row.sum <- apply(housetasks, 1, sum)
# row mass
row.mass <- row.sum/n

# row coord = sv * u /sqrt(row.mass)
cc <- t(apply(u, 1, &#39;*&#39;, sv)) # each row X sv
row.coord <- apply(cc, 2, &#39;/&#39;, sqrt(row.mass))
rownames(row.coord) <- rownames(housetasks)
colnames(row.coord) <- paste0("Dim.", 1:nb.axes)
round(row.coord,3)

##             Dim.1  Dim.2  Dim.3
## Laundry    -0.992 -0.495 -0.317
## Main_meal  -0.876 -0.490 -0.164
## Dinner     -0.693 -0.308 -0.207
## Breakfeast -0.509 -0.453  0.220
## Tidying    -0.394  0.434 -0.094
## Dishes     -0.189  0.442  0.267
## Shopping   -0.118  0.403  0.203
## Official    0.227 -0.254  0.923
## Driving     0.742 -0.653  0.544
## Finances    0.271  0.618  0.035
## Insurance   0.647  0.474 -0.289
## Repairs     1.529 -0.864 -0.472
## Holidays    0.252  1.435 -0.130

# plot
plot(row.coord, pch=19, col = "blue")
text(row.coord, labels =rownames(row.coord), pos = 3, col ="blue")
abline(v=0, h=0, lty = 2)

Column coordinates

# Coordinates of columns
col.sum <- apply(housetasks, 2, sum)
col.mass <- col.sum/n
# coordinates sv * v /sqrt(col.mass)
cc <- t(apply(v, 1, &#39;*&#39;, sv))
col.coord <- apply(cc, 2, &#39;/&#39;, sqrt(col.mass))
rownames(col.coord) <- colnames(housetasks)
colnames(col.coord) <- paste0("Dim", 1:nb.axes)
head(col.coord)

##                Dim1   Dim2    Dim3
## Wife        -0.8376 -0.365 -0.1999
## Alternating -0.0622 -0.292  0.8486
## Husband      1.1609 -0.602 -0.1889
## Jointly      0.1494  1.027 -0.0464

# plot
plot(col.coord, pch=17, col = "red")
text(col.coord, labels =rownames(col.coord), pos = 3, col ="red")
abline(v=0, h=0, lty = 2)

Biplot of rows and columns to view the association

xlim <- range(c(row.coord[,1], col.coord[,1]))*1.1
ylim <- range(c(row.coord[,2], col.coord[,2]))*1.1
# Plot of rows
plot(row.coord, pch=19, col = "blue", xlim = xlim, ylim = ylim)
text(row.coord, labels =rownames(row.coord), pos = 3, col ="blue")
# plot off columns
points(col.coord, pch=17, col = "red")
text(col.coord, labels =rownames(col.coord), pos = 3, col ="red")
abline(v=0, h=0, lty = 2)

You can interpret the distance between rows points or between column points but the distance between column points and row points are not meaningful.

Diagnostic

Recall that, the total inertia contained in the data is:

\[ \phi^2 = \frac{\chi^2}{n} = 1.11494 \]

Our two-dimensional plot captures about 88% of the total inertia of the table.

Contribution of rows and columns

The contributions of a rows/columns to the definition of a principal axis are :

\[ row.contrib = \frac{row.mass * row.coord^2}{eigenvalue} \]

\[ col.contrib = \frac{col.mass * col.coord^2}{eigenvalue} \]

Contribution of rows in %

# contrib <- row.mass * row.coord^2/eigenvalue
cc <- apply(row.coord^2, 2, "*", row.mass)
row.contrib <- t(apply(cc, 1, "/", eig[1:nb.axes,1])) *100
round(row.contrib, 2)

##            Dim.1 Dim.2 Dim.3
## Laundry    18.29  5.56  7.97
## Main_meal  12.39  4.74  1.86
## Dinner      5.47  1.32  2.10
## Breakfeast  3.82  3.70  3.07
## Tidying     2.00  2.97  0.49
## Dishes      0.43  2.84  3.63
## Shopping    0.18  2.52  2.22
## Official    0.52  0.80 36.94
## Driving     8.08  7.65 18.60
## Finances    0.88  5.56  0.06
## Insurance   6.15  4.02  5.25
## Repairs    40.73 15.88 16.60
## Holidays    1.08 42.45  1.21

corrplot(row.contrib, is.cor = FALSE)

Contribution of columns in %

# contrib <- col.mass * col.coord^2/eigenvalue
cc <- apply(col.coord^2, 2, "*", col.mass)
col.contrib <- t(apply(cc, 1, "/", eig[1:nb.axes,1])) *100
round(col.contrib, 2)

##             Dim1  Dim2  Dim3
## Wife        44.5 10.31 10.82
## Alternating  0.1  2.78 82.55
## Husband     54.2 17.79  6.13
## Jointly      1.2 69.12  0.50

corrplot(col.contrib, is.cor = FALSE)

Quality of the representation

The quality of the representation is called COS2.

The quality of the representation of a row on an axis is:

\[ row.cos2 = \frac{row.coord^2}{d^2} \]

• row.coord is the coordinate of the row on the axis
• \(d^2\) is the squared distance from the average profile

Recall that the distance between each row profile and the average row profile is:

\[ d^2(row_i, average.profile) = \sum{\frac{(row.profile_i - average.profile)^2}{average.profile}} \]

row.profile <- housetasks/row.sum
head(round(row.profile, 3))

##             Wife Alternating Husband Jointly
## Laundry    0.886       0.080   0.011   0.023
## Main_meal  0.810       0.131   0.033   0.026
## Dinner     0.713       0.102   0.065   0.120
## Breakfeast 0.586       0.257   0.107   0.050
## Tidying    0.434       0.090   0.008   0.467
## Dishes     0.283       0.212   0.035   0.469

average.profile <- col.sum/n
head(round(average.profile, 3))

##        Wife Alternating     Husband     Jointly 
##       0.344       0.146       0.218       0.292

The R code below computes the distance from the average profile for all the row variables

d2.row <- apply(row.profile, 1, 
                function(row.p, av.p){sum(((row.p - av.p)^2)/av.p)}, 
                average.rp)
head(round(d2.row,3))

##    Laundry  Main_meal     Dinner Breakfeast    Tidying     Dishes 
##      1.329      1.034      0.618      0.512      0.353      0.302

The cos2 of rows on the factor map are:

row.cos2 <- apply(row.coord^2, 2, "/", d2.row)
round(row.cos2, 3)

##            Dim.1 Dim.2 Dim.3
## Laundry    0.740 0.185 0.075
## Main_meal  0.742 0.232 0.026
## Dinner     0.777 0.154 0.070
## Breakfeast 0.505 0.400 0.095
## Tidying    0.440 0.535 0.025
## Dishes     0.118 0.646 0.236
## Shopping   0.064 0.748 0.189
## Official   0.053 0.066 0.881
## Driving    0.432 0.335 0.233
## Finances   0.161 0.837 0.003
## Insurance  0.576 0.309 0.115
## Repairs    0.707 0.226 0.067
## Holidays   0.030 0.962 0.008

visualize the cos2:

corrplot(row.cos2, is.cor = FALSE)

Cos2 of columns

\[ col.cos2 = \frac{col.coord^2}{d^2} \]

col.profile <- t(housetasks)/col.sum
col.profile <- t(col.profile)
#head(round(col.profile, 3))

average.profile <- row.sum/n
#head(round(average.profile, 3))

The R code below computes the distance from the average profile for all the column variables

d2.col <- apply(col.profile, 2, 
        function(col.p, av.p){sum(((col.p - av.p)^2)/av.p)}, 
        average.profile)
#round(d2.col,3)

The cos2 of columns on the factor map are:

col.cos2 <- apply(col.coord^2, 2, "/", d2.col)
round(col.cos2, 3)

##              Dim1  Dim2  Dim3
## Wife        0.802 0.152 0.046
## Alternating 0.005 0.105 0.890
## Husband     0.772 0.208 0.020
## Jointly     0.021 0.977 0.002

visualize the cos2:

corrplot(col.cos2, is.cor = FALSE)

Supplementary rows/columns

# Supplementary row
sup.row <- as.data.frame(housetasks["Dishes",, drop = FALSE])
# Supplementary row profile
sup.row.sum <- apply(sup.row, 1, sum)
sup.row.profile <- sweep(sup.row, 1, sup.row.sum, "/")
# V/sqrt(col.mass)
vv <- sweep(v, 1, sqrt(col.mass), FUN = "/")
# Supplementary row coord
sup.row.coord <- as.matrix(sup.row.profile) %*% vv
sup.row.coord

##          [,1]  [,2]  [,3]
## Dishes -0.189 0.442 0.267

## COS2 = coor^2/Distance from average profile
d2.row <- apply(sup.row.profile, 1, 
        function(row.p, av.p){sum(((row.p - av.p)^2)/av.p)}, 
        average.rp)
sup.row.cos2 <- sweep(sup.row.coord^2, 1, d2.row, FUN = "/")

Visualize using factoextra

The factoextra R package creates ggplot2-based visualization.

1. Visualize eigenvalues (scree plot). Show the percentage of variances explained by each principal component.

fviz_eig(res.pca)

2. Graph of individuals. Individuals with a similar profile are grouped together.

fviz_pca_ind(res.pca,
             col.ind = "cos2", # Color by the quality of representation
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             )

3. Graph of variables. Positive correlated variables point to the same side of the plot. Negative correlated variables point to opposite sides of the graph.

fviz_pca_var(res.pca,
             col.var = "contrib", # Color by contributions to the PC
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             )

4. Biplot of individuals and variables

fviz_pca_biplot(res.pca, repel = TRUE,
                col.var = "#2E9FDF", # Variables color
                col.ind = "#696969"  # Individuals color
                )

Visualize using ade4

The ade4 package creates R base plots.

# Scree plot
screeplot(res.pca, main = "Screeplot - Eigenvalues")

# Correlation circle of variables
s.corcircle(res.pca$co)

# Graph of individuals
s.label(res.pca$li, 
        xax = 1,     # Dimension 1
        yax = 2)     # Dimension 2

# Biplot of individuals and variables
scatter(res.pca,
        posieig = "none", # Hide the scree plot
        clab.row = 0      # Hide row labels
        )

## NULL

Supplementary individuals

1. Data: rows 24 to 27 and columns 1 to to 10 [in decathlon2 data sets]. The new data must contain columns (variables) with the same names and in the same order as the active data used to compute PCA.

# Data for the supplementary individuals
ind.sup <- decathlon2[24:27, 1:10]
ind.sup[, 1:6]

##         X100m Long.jump Shot.put High.jump X400m X110m.hurdle
## KARPOV   11.0      7.30     14.8      2.04  48.4         14.1
## WARNERS  11.1      7.60     14.3      1.98  48.7         14.2
## Nool     10.8      7.53     14.3      1.88  48.8         14.8
## Drews    10.9      7.38     13.1      1.88  48.5         14.0

2. Predict the coordinates of new individuals data.

ind.sup.coord <- suprow(res.pca, ind.sup) %>%
  .$lisup
ind.sup.coord[, 1:4]

##          Axis1   Axis2 Axis3  Axis4
## KARPOV  -0.795  0.7795 1.633  1.724
## WARNERS  0.386 -0.1216 1.739 -0.706
## Nool     0.559  1.9775 0.483 -2.278
## Drews    1.109  0.0174 3.049 -1.534

3. Graph of individuals including the supplementary individuals:

# Plot of active individuals
p <- fviz_pca_ind(res.pca, repel = TRUE)
# Add supplementary individuals
fviz_add(p, ind.sup.coord, color ="blue")

Supplementary variables

groups <- as.factor(decathlon2$Competition[1:23])
fviz_pca_ind(res.pca,
             col.ind = groups, # color by groups
             palette = c("#00AFBB",  "#FC4E07"),
             addEllipses = TRUE, # Concentration ellipses
             ellipse.type = "confidence",
             legend.title = "Groups",
             repel = TRUE
             )

groups <- as.factor(decathlon2$Competition[1:23])
s.class(res.pca$li,
        fac = groups,  # color by groups
        col = c("#00AFBB",  "#FC4E07")
        )

# Biplot
res <- scatter(res.pca, clab.row = 0, posieig = "none")
s.class(res.pca$li, 
        fac = groups,
        col = c("#00AFBB",  "#FC4E07"),
        add.plot = TRUE,         # Add onto the scatter plot
        cstar = 0,               # Remove stars
        cellipse = 0             # Remove ellipses
        )

quanti.sup <- decathlon2[1:23, 11:12, drop = FALSE]
head(quanti.sup)

##           Rank Points
## SEBRLE       1   8217
## CLAY         2   8122
## BERNARD      4   8067
## YURKOV       5   8036
## ZSIVOCZKY    7   8004
## McMULLEN     8   7995

# Predict coordinates and compute cos2
quanti.coord <- supcol(res.pca, scale(quanti.sup)) %>%
  .$cosup
quanti.cos2 <- quanti.coord^2
# Graph of variables including supplementary variables
p <- fviz_pca_var(res.pca)
fviz_add(p, quanti.coord, color ="blue", geom="arrow")

Supplementary individuals

1. Data: rows 24 to 27 and columns 1 to to 10 [in decathlon2 data sets]. The new data must contain columns (variables) with the same names and in the same order as the active data used to compute PCA.

# Data for the supplementary individuals
ind.sup <- decathlon2[24:27, 1:10]
ind.sup[, 1:6]

##         X100m Long.jump Shot.put High.jump X400m X110m.hurdle
## KARPOV   11.0      7.30     14.8      2.04  48.4         14.1
## WARNERS  11.1      7.60     14.3      1.98  48.7         14.2
## Nool     10.8      7.53     14.3      1.88  48.8         14.8
## Drews    10.9      7.38     13.1      1.88  48.5         14.0

2. Predict the coordinates of new individuals data. Use the R base function predict():

ind.sup.coord <- predict(res.pca, newdata = ind.sup)
ind.sup.coord[, 1:4]

##            PC1    PC2   PC3    PC4
## KARPOV   0.777 -0.762 1.597  1.686
## WARNERS -0.378  0.119 1.701 -0.691
## Nool    -0.547 -1.934 0.472 -2.228
## Drews   -1.085 -0.017 2.982 -1.501

3. Graph of individuals including the supplementary individuals:

# Plot of active individuals
p <- fviz_pca_ind(res.pca, repel = TRUE)
# Add supplementary individuals
fviz_add(p, ind.sup.coord, color ="blue")

The predicted coordinates of individuals can be manually calculated as follow:

1. Center and scale the new individuals data using the center and the scale of the PCA
2. Calculate the predicted coordinates by multiplying the scaled values with the eigenvectors (loadings) of the principal components.

The R code below can be used :

# Centering and scaling the supplementary individuals
ind.scaled <- scale(ind.sup, 
                    center = res.pca$center,
                    scale = res.pca$scale)
# Coordinates of the individividuals
coord_func <- function(ind, loadings){
  r <- loadings*ind
  apply(r, 2, sum)
}
pca.loadings <- res.pca$rotation
ind.sup.coord <- t(apply(ind.scaled, 1, coord_func, pca.loadings ))
ind.sup.coord[, 1:4]

##            PC1    PC2   PC3    PC4
## KARPOV   0.777 -0.762 1.597  1.686
## WARNERS -0.378  0.119 1.701 -0.691
## Nool    -0.547 -1.934 0.472 -2.228
## Drews   -1.085 -0.017 2.982 -1.501

Supplementary variables

groups <- as.factor(decathlon2$Competition[1:23])
fviz_pca_ind(res.pca,
             col.ind = groups, # color by groups
             palette = c("#00AFBB",  "#FC4E07"),
             addEllipses = TRUE, # Concentration ellipses
             ellipse.type = "confidence",
             legend.title = "Groups",
             repel = TRUE
             )

library(magrittr) # for pipe %>%
library(dplyr)   # everything else
# 1. Individual coordinates
res.ind <- get_pca_ind(res.pca)
# 2. Coordinate of groups
coord.groups <- res.ind$coord %>%
  as_data_frame() %>%
  select(Dim.1, Dim.2) %>%
  mutate(competition = groups) %>%
  group_by(competition) %>%
  summarise(
    Dim.1 = mean(Dim.1),
    Dim.2 = mean(Dim.2)
    )
coord.groups

## # A tibble: 2 x 3
##   competition Dim.1  Dim.2
##           
## 1    Decastar -1.31 -0.119
## 2    OlympicG  1.20  0.109

quanti.sup <- decathlon2[1:23, 11:12, drop = FALSE]
head(quanti.sup)

##           Rank Points
## SEBRLE       1   8217
## CLAY         2   8122
## BERNARD      4   8067
## YURKOV       5   8036
## ZSIVOCZKY    7   8004
## McMULLEN     8   7995

# Predict coordinates and compute cos2
quanti.coord <- cor(quanti.sup, res.pca$x)
quanti.cos2 <- quanti.coord^2
# Graph of variables including supplementary variables
p <- fviz_pca_var(res.pca)
fviz_add(p, quanti.coord, color ="blue", geom="arrow")

PCA results for variables

Here we’ll show how to calculate the PCA results for variables: coordinates, cos2 and contributions:

• var.coord = loadings * the component standard deviations
• var.cos2 = var.coord^2
• var.contrib. The contribution of a variable to a given principal component is (in percentage) : (var.cos2 * 100) / (total cos2 of the component)

# Helper function 
#::::::::::::::::::::::::::::::::::::::::
var_coord_func <- function(loadings, comp.sdev){
  loadings*comp.sdev
}
# Compute Coordinates
#::::::::::::::::::::::::::::::::::::::::
loadings <- res.pca$rotation
sdev <- res.pca$sdev
var.coord <- t(apply(loadings, 1, var_coord_func, sdev)) 
head(var.coord[, 1:4])

##                 PC1     PC2    PC3     PC4
## X100m        -0.851  0.1794 -0.302  0.0336
## Long.jump     0.794 -0.2809  0.191 -0.1154
## Shot.put      0.734 -0.0854 -0.518  0.1285
## High.jump     0.610  0.4652 -0.330  0.1446
## X400m        -0.702 -0.2902 -0.284  0.4308
## X110m.hurdle -0.764  0.0247 -0.449 -0.0169

# Compute Cos2
#::::::::::::::::::::::::::::::::::::::::
var.cos2 <- var.coord^2
head(var.cos2[, 1:4])

##                PC1      PC2    PC3      PC4
## X100m        0.724 0.032184 0.0909 0.001127
## Long.jump    0.631 0.078881 0.0363 0.013315
## Shot.put     0.539 0.007294 0.2679 0.016504
## High.jump    0.372 0.216424 0.1090 0.020895
## X400m        0.492 0.084203 0.0804 0.185611
## X110m.hurdle 0.584 0.000612 0.2015 0.000285

# Compute contributions
#::::::::::::::::::::::::::::::::::::::::
comp.cos2 <- apply(var.cos2, 2, sum)
contrib <- function(var.cos2, comp.cos2){var.cos2*100/comp.cos2}
var.contrib <- t(apply(var.cos2,1, contrib, comp.cos2))
head(var.contrib[, 1:4])

##                PC1     PC2   PC3     PC4
## X100m        17.54  1.7505  7.34  0.1376
## Long.jump    15.29  4.2904  2.93  1.6249
## Shot.put     13.06  0.3967 21.62  2.0141
## High.jump     9.02 11.7716  8.79  2.5499
## X400m        11.94  4.5799  6.49 22.6509
## X110m.hurdle 14.16  0.0333 16.26  0.0348

PCA results for individuals

• ind.coord = res.pca$x
• Cos2 of individuals. Two steps:
  • Calculate the square distance between each individual and the PCA center of gravity: d2 = [(var1_ind_i - mean_var1)/sd_var1]^2 + …+ [(var10_ind_i - mean_var10)/sd_var10]^2 + …+..
  • Calculate the cos2 as ind.coord^2/d2
• Contributions of individuals to the principal components: 100 * (1 / number_of_individuals)*(ind.coord^2 / comp_sdev^2). Note that the sum of all the contributions per column is 100

# Coordinates of individuals
#::::::::::::::::::::::::::::::::::
ind.coord <- res.pca$x
head(ind.coord[, 1:4])

##              PC1    PC2    PC3     PC4
## SEBRLE     0.191 -1.554 -0.628  0.0821
## CLAY       0.790 -2.420  1.357  1.2698
## BERNARD   -1.329 -1.612 -0.196 -1.9209
## YURKOV    -0.869  0.433 -2.474  0.6972
## ZSIVOCZKY -0.106  2.023  1.305 -0.0993
## McMULLEN   0.119  0.992  0.844  1.3122

# Cos2 of individuals
#:::::::::::::::::::::::::::::::::
# 1. square of the distance between an individual and the
# PCA center of gravity
center <- res.pca$center
scale<- res.pca$scale
getdistance <- function(ind_row, center, scale){
  return(sum(((ind_row-center)/scale)^2))
  }
d2 <- apply(decathlon2.active,1,getdistance, center, scale)
# 2. Compute the cos2. The sum of each row is 1
cos2 <- function(ind.coord, d2){return(ind.coord^2/d2)}
ind.cos2 <- apply(ind.coord, 2, cos2, d2)
head(ind.cos2[, 1:4])

##               PC1    PC2     PC3     PC4
## SEBRLE    0.00753 0.4975 0.08133 0.00139
## CLAY      0.04870 0.4570 0.14363 0.12579
## BERNARD   0.19720 0.2900 0.00429 0.41182
## YURKOV    0.09611 0.0238 0.77823 0.06181
## ZSIVOCZKY 0.00157 0.5764 0.23975 0.00139
## McMULLEN  0.00218 0.1522 0.11014 0.26649

# Contributions of individuals
#:::::::::::::::::::::::::::::::
contrib <- function(ind.coord, comp.sdev, n.ind){
  100*(1/n.ind)*ind.coord^2/comp.sdev^2
}
ind.contrib <- t(apply(ind.coord, 1, contrib, 
                       res.pca$sdev, nrow(ind.coord)))
head(ind.contrib[, 1:4])

##              PC1    PC2    PC3     PC4
## SEBRLE    0.0385  5.712  1.385  0.0357
## CLAY      0.6581 13.854  6.460  8.5557
## BERNARD   1.8627  6.144  0.135 19.5783
## YURKOV    0.7969  0.443 21.476  2.5794
## ZSIVOCZKY 0.0118  9.682  5.975  0.0523
## McMULLEN  0.0148  2.325  2.497  9.1353

Case 1: Continuous variables

In the situation where you have a multidimensional data set containing multiple continuous variables, the principal component analysis (PCA) can be used to reduce the dimension of the data into few continuous variables containing the most important information in the data. Next, you can perform cluster analysis on the PCA results.

The PCA step can be considered as a denoising step which can lead to a more stable clustering. This might be very useful if you have a large data set with multiple variables, such as in gene expression data.

Case 2: Clustering on categorical data

In order to perform clustering analysis on categorical data, the correspondence analysis (CA, for analyzing contingency table) and the multiple correspondence analysis (MCA, for analyzing multidimensional categorical variables) can be used to transform categorical variables into a set of few continuous variables (the principal components). The cluster analysis can be then applied on the (M)CA results.

In this case, the (M)CA method can be considered as pre-processing steps which allow to compute clustering on categorical data.

Case 3: Clustering on mixed data

When you have a mixed data of continuous and categorical variables, you can first perform FAMD (factor analysis of mixed data) or MFA (multiple factor analysis). Next, you can apply cluster analysis on the FAMD/MFA outputs.

R packages

We’ll use two R packages: i) FactoMineR for computing HCPC and ii) factoextra for visualizing the results.

To install the packages, type this:

install.packages(c("FactoMineR", "factoextra"))

After the installation, load the packages as follow:

library(factoextra)
library(FactoMineR)

R function

The function HCPC() [in FactoMineR package] can be used to compute hierarchical clustering on principal components.

A simplified format is:

HCPC(res, nb.clust = 0, min = 3, max = NULL, graph = TRUE)

• res: Either the result of a factor analysis or a data frame.
• nb.clust: an integer specifying the number of clusters. Possible values are:
  • 0: the tree is cut at the level the user clicks on
  • -1: the tree is automatically cut at the suggested level
  • Any positive integer: the tree is cut with nb.clusters clusters
• min, max: the minimum and the maximum number of clusters to be generated, respectively
• graph: if TRUE, graphics are displayed

Case of continuous variables

We start by computing again the principal component analysis (PCA). The argument ncp = 3 is used in the function PCA() to keep only the first three principal components. Next, the HCPC is applied on the result of the PCA.

library(FactoMineR)
# Compute PCA with ncp = 3
res.pca <- PCA(USArrests, ncp = 3, graph = FALSE)
# Compute hierarchical clustering on principal components
res.hcpc <- HCPC(res.pca, graph = FALSE)

To visualize the dendrogram generated by the hierarchical clustering, we’ll use the function fviz_dend() [factoextra package]:

fviz_dend(res.hcpc, 
          cex = 0.7,                     # Label size
          palette = "jco",               # Color palette see ?ggpubr::ggpar
          rect = TRUE, rect_fill = TRUE, # Add rectangle around groups
          rect_border = "jco",           # Rectangle color
          labels_track_height = 0.8      # Augment the room for labels
          )

It’s possible to visualize individuals on the principal component map and to color individuals according to the cluster they belong to. The function fviz_cluster() [in factoextra] can be used to visualize individuals clusters.

fviz_cluster(res.hcpc,
             repel = TRUE,            # Avoid label overlapping
             show.clust.cent = TRUE, # Show cluster centers
             palette = "jco",         # Color palette see ?ggpubr::ggpar
             ggtheme = theme_minimal(),
             main = "Factor map"
             )

You can also draw a three dimensional plot combining the hierarchical clustering and the factorial map using the R base function plot():

# Principal components + tree
plot(res.hcpc, choice = "3D.map")

The function HCPC() returns a list containing:

• data.clust: The original data with a supplementary column called class containing the partition.
• desc.var: The variables describing clusters
• desc.ind: The more typical individuals of each cluster
• desc.axes: The axes describing clusters

To display the original data with cluster assignments, type this:

head(res.hcpc$data.clust, 10)

##             Murder Assault UrbanPop Rape clust
## Alabama       13.2     236       58 21.2     3
## Alaska        10.0     263       48 44.5     4
## Arizona        8.1     294       80 31.0     4
## Arkansas       8.8     190       50 19.5     3
## California     9.0     276       91 40.6     4
## Colorado       7.9     204       78 38.7     4
## Connecticut    3.3     110       77 11.1     2
## Delaware       5.9     238       72 15.8     2
## Florida       15.4     335       80 31.9     4
## Georgia       17.4     211       60 25.8     3

In the table above, the last column contains the cluster assignments.

To display quantitative variables that describe the most each cluster, type this:

res.hcpc$desc.var$quanti

Here, we show only some columns of interest: “Mean in category”, “Overall Mean”, “p.value”

## $`1`
##          Mean in category Overall mean  p.value
## UrbanPop             52.1        65.54 9.68e-05
## Murder                3.6         7.79 5.57e-05
## Rape                 12.2        21.23 5.08e-05
## Assault              78.5       170.76 3.52e-06
## 
## $`2`
##          Mean in category Overall mean p.value
## UrbanPop            73.88        65.54 0.00522
## Murder               5.66         7.79 0.01759
## 
## $`3`
##          Mean in category Overall mean  p.value
## Murder               13.9         7.79 1.32e-05
## Assault             243.6       170.76 6.97e-03
## UrbanPop             53.8        65.54 1.19e-02
## 
## $`4`
##          Mean in category Overall mean  p.value
## Rape                 33.2        21.23 8.69e-08
## Assault             257.4       170.76 1.32e-05
## UrbanPop             76.0        65.54 2.45e-03
## Murder               10.8         7.79 3.58e-03

Similarly, to show principal dimensions that are the most associated with clusters, type this:

res.hcpc$desc.axes$quanti

## $`1`
##       Mean in category Overall mean  p.value
## Dim.1            -1.96    -5.64e-16 2.27e-07
## 
## $`2`
##       Mean in category Overall mean  p.value
## Dim.2            0.743    -5.37e-16 0.000336
## 
## $`3`
##       Mean in category Overall mean  p.value
## Dim.1            1.061    -5.64e-16 3.96e-02
## Dim.3            0.397     3.54e-17 4.25e-02
## Dim.2           -1.477    -5.37e-16 5.72e-06
## 
## $`4`
##       Mean in category Overall mean  p.value
## Dim.1             1.89    -5.64e-16 6.15e-07

Finally, representative individuals of each cluster can be extracted as follow:

res.hcpc$desc.ind$para

## Cluster: 1
##         Idaho  South Dakota         Maine          Iowa New Hampshire 
##         0.367         0.499         0.501         0.553         0.589 
## -------------------------------------------------------- 
## Cluster: 2
##         Ohio     Oklahoma Pennsylvania       Kansas      Indiana 
##        0.280        0.505        0.509        0.604        0.710 
## -------------------------------------------------------- 
## Cluster: 3
##        Alabama South Carolina        Georgia      Tennessee      Louisiana 
##          0.355          0.534          0.614          0.852          0.878 
## -------------------------------------------------------- 
## Cluster: 4
##   Michigan    Arizona New Mexico   Maryland      Texas 
##      0.325      0.453      0.518      0.901      0.924

Case of categorical variables

For categorical variables, compute CA or MCA and then apply the function HCPC() on the results as described above.

Here, we’ll use the tea data [in FactoMineR] as demo data set: Rows represent the individuals and columns represent categorical variables.

We start, by performing an MCA on the individuals. We keep the first 20 axes of the MCA which retain 87% of the information.

# Loading data
library(FactoMineR)
data(tea)

# Performing MCA
res.mca <- MCA(tea, 
               ncp = 20,            # Number of components kept
               quanti.sup = 19,     # Quantitative supplementary variables
               quali.sup = c(20:36), # Qualitative supplementary variables
               graph=FALSE)

Next, we apply hierarchical clustering on the results of the MCA:

res.hcpc <- HCPC (res.mca, graph = FALSE, max = 3)

The results can be visualized as follow:

# Dendrogram
fviz_dend(res.hcpc, show_labels = FALSE)

# Individuals facor map
fviz_cluster(res.hcpc, geom = "point", main = "Factor map")

As mentioned above, clusters can be described by i) variables and/or categories, ii) principal axes and iii) individuals. In the example below, we display only a subset of the results.

• Description by variables and categories

# Description by variables
res.hcpc$desc.var$test.chi2

##          p.value df
## where   8.47e-79  4
## how     3.14e-47  4
## price   1.86e-28 10
## tearoom 9.62e-19  2

# Description by variable categories
res.hcpc$desc.var$category

## $`1`
##                     Cla/Mod Mod/Cla Global  p.value
## where=chain store      85.9    93.8   64.0 2.09e-40
## how=tea bag            84.1    81.2   56.7 1.48e-25
## tearoom=Not.tearoom    70.7    97.2   80.7 1.08e-18
## price=p_branded        83.2    44.9   31.7 1.63e-09
## 
## $`2`
##                 Cla/Mod Mod/Cla Global  p.value
## where=tea shop     90.0    84.4   10.0 3.70e-30
## how=unpackaged     66.7    75.0   12.0 5.35e-20
## price=p_upscale    49.1    81.2   17.7 2.39e-17
## Tea=green          27.3    28.1   11.0 4.44e-03
## 
## $`3`
##                            Cla/Mod Mod/Cla Global  p.value
## where=chain store+tea shop    85.9    72.8   26.0 5.73e-34
## how=tea bag+unpackaged        67.0    68.5   31.3 1.38e-19
## tearoom=tearoom               77.6    48.9   19.3 1.25e-16
## pub=pub                       63.5    43.5   21.0 1.13e-09

The variables that characterize the most the clusters are the variables “where” and “how”. Each cluster is characterized by a category of the variables “where” and “how”. For example, individuals who belong to the first cluster buy tea as tea bag in chain stores.

• Description by principal components

res.hcpc$desc.axes

• Description by individuals

res.hcpc$desc.ind$para

R packages

Install FactoMineR and factoextra as follow:

install.packages(c("FactoMineR", "factoextra"))

Load the packages:

library("FactoMineR")
library("factoextra")

Data format

We’ll use the demo data sets wine available in FactoMineR package. This data set is about a sensory evaluation of wines by different judges.

library("FactoMineR")
data(wine)
colnames(wine)

##  [1] "Label"                         "Soil"                         
##  [3] "Odor.Intensity.before.shaking" "Aroma.quality.before.shaking" 
##  [5] "Fruity.before.shaking"         "Flower.before.shaking"        
##  [7] "Spice.before.shaking"          "Visual.intensity"             
##  [9] "Nuance"                        "Surface.feeling"              
## [11] "Odor.Intensity"                "Quality.of.odour"             
## [13] "Fruity"                        "Flower"                       
## [15] "Spice"                         "Plante"                       
## [17] "Phenolic"                      "Aroma.intensity"              
## [19] "Aroma.persistency"             "Aroma.quality"                
## [21] "Attack.intensity"              "Acidity"                      
## [23] "Astringency"                   "Alcohol"                      
## [25] "Balance"                       "Smooth"                       
## [27] "Bitterness"                    "Intensity"                    
## [29] "Harmony"                       "Overall.quality"              
## [31] "Typical"

An image of the data is shown below:

(Image source, FactoMineR, http://factominer.free.fr)

The data contains 21 rows (wines, individuals) and 31 columns (variables):

• The first two columns are categorical variables: label (Saumur, Bourgueil or Chinon) and soil (Reference, Env1, Env2 or Env4).
• The 29 next columns are continuous sensory variables. For each wine, the value is the mean score for all the judges.

The variables are organized in groups as follow:

1. First group - A group of categorical variables specifying the origin of the wines, including the variables label and soil corresponding to the first 2 columns in the data table. In FactoMineR terminology, the arguments group = 2 is used to define the first 2 columns as a group.

2. Second group - A group of continuous variables, describing the odor of the wines before shaking, including the variables: Odor.Intensity.before.shaking, Aroma.quality.before.shaking, Fruity.before.shaking, Flower.before.shaking and Spice.before.shaking. These variables corresponds to the next 5 columns after the first group. FactoMineR terminology: group = 5.

3. Third group - A group of continuous variables quantifying the visual inspection of the wines, including the variables: Visual.intensity, Nuance and Surface.feeling. These variables corresponds to the next 3 columns after the second group. FactoMineR terminology: group = 3.

4. Fourth group - A group of continuous variables concerning the odor of the wines after shaking, including the variables: Odor.Intensity, Quality.of.odour, Fruity, Flower, Spice, Plante, Phenolic, Aroma.intensity, Aroma.persistency and Aroma.quality. These variables corresponds to the next 10 columns after the third group. FactoMineR terminology: group = 10.

5. Fith group - A group of continuous variables evaluating the taste of the wines, including the variables Attack.intensity, Acidity, Astringency, Alcohol, Balance, Smooth, Bitterness, Intensity and Harmony. These variables corresponds to the next 9 columns after the fourth group. FactoMineR terminology: group = 9.

6. Sixth group - A group of continuous variables concerning the overall judgement of the wines, including the variables Overall.quality and Typical. These variables corresponds to the next 2 columns after the fith group. FactoMineR terminology: group = 2.

R code

The function MFA()[FactoMiner package] can be used. A simplified format is :

MFA (base, group, type = rep("s",length(group)), ind.sup = NULL, 
    name.group = NULL, num.group.sup = NULL, graph = TRUE)

• base : a data frame with n rows (individuals) and p columns (variables)
• group: a vector with the number of variables in each group.
• type: the type of variables in each group. By default, all variables are quantitative and scaled to unit variance. Allowed values include:
  • “c” or “s” for quantitative variables. If “s”, the variables are scaled to unit variance.
  • “n” for categorical variables.
  • “f” for frequencies (from a contingency tables).
• ind.sup: a vector indicating the indexes of the supplementary individuals.
• name.group: a vector containing the name of the groups (by default, NULL and the group are named group.1, group.2 and so on).
• num.group.sup: the indexes of the illustrative groups (by default, NULL and no group are illustrative).
• graph : a logical value. If TRUE a graph is displayed.

The R code below performs the MFA on the wines data using the groups: odor, visual, odor after shaking and taste. These groups are named active groups. The remaining group of variables - origin (the first group) and overall judgement (the sixth group) - are named supplementary groups; num.group.sup = c(1, 6):

library(FactoMineR)
data(wine)
res.mfa <- MFA(wine, 
               group = c(2, 5, 3, 10, 9, 2), 
               type = c("n", "s", "s", "s", "s", "s"),
               name.group = c("origin","odor","visual",
                              "odor.after.shaking", "taste","overall"),
               num.group.sup = c(1, 6),
               graph = FALSE)

The output of the MFA() function is a list including :

print(res.mfa)

## **Results of the Multiple Factor Analysis (MFA)**
## The analysis was performed on 21 individuals, described by 31 variables
## *Results are available in the following objects :
## 
##    name                
## 1  "$eig"              
## 2  "$separate.analyses"
## 3  "$group"            
## 4  "$partial.axes"     
## 5  "$inertia.ratio"    
## 6  "$ind"              
## 7  "$quanti.var"       
## 8  "$quanti.var.sup"   
## 9  "$quali.var.sup"    
## 10 "$summary.quanti"   
## 11 "$summary.quali"    
## 12 "$global.pca"       
##    description                                           
## 1  "eigenvalues"                                         
## 2  "separate analyses for each group of variables"       
## 3  "results for all the groups"                          
## 4  "results for the partial axes"                        
## 5  "inertia ratio"                                       
## 6  "results for the individuals"                         
## 7  "results for the quantitative variables"              
## 8  "results for the quantitative supplementary variables"
## 9  "results for the categorical supplementary variables" 
## 10 "summary for the quantitative variables"              
## 11 "summary for the categorical variables"               
## 12 "results for the global PCA"

Eigenvalues / Variances

The proportion of variances retained by the different dimensions (axes) can be extracted using the function get_eigenvalue() [factoextra package] as follow:

library("factoextra")
eig.val <- get_eigenvalue(res.mfa)
head(eig.val)

##       eigenvalue variance.percent cumulative.variance.percent
## Dim.1      3.462            49.38                        49.4
## Dim.2      1.367            19.49                        68.9
## Dim.3      0.615             8.78                        77.7
## Dim.4      0.372             5.31                        83.0
## Dim.5      0.270             3.86                        86.8
## Dim.6      0.202             2.89                        89.7

The function fviz_eig() or fviz_screeplot() [factoextra package] can be used to draw the scree plot:

fviz_screeplot(res.mfa)

Graph of variables

group <- get_mfa_var(res.mfa, "group")
group

## Multiple Factor Analysis results for variable groups 
##  ===================================================
##   Name           Description                                          
## 1 "$coord"       "Coordinates"                                        
## 2 "$cos2"        "Cos2, quality of representation"                    
## 3 "$contrib"     "Contributions"                                      
## 4 "$correlation" "Correlation between groups and principal dimensions"

# Coordinates of groups
head(group$coord)

# Cos2: quality of representation on the factore map
head(group$cos2)

# Contributions to the  dimensions
head(group$contrib)

fviz_mfa_var(res.mfa, "group")

# Contribution to the first dimension
fviz_contrib(res.mfa, "group", axes = 1)

# Contribution to the second dimension
fviz_contrib(res.mfa, "group", axes = 2)

quanti.var <- get_mfa_var(res.mfa, "quanti.var")
quanti.var 

## Multiple Factor Analysis results for quantitative variables 
##  ===================================================
##   Name       Description                      
## 1 "$coord"   "Coordinates"                    
## 2 "$cos2"    "Cos2, quality of representation"
## 3 "$contrib" "Contributions"

# Coordinates
head(quanti.var$coord)

# Cos2: quality on the factore map
head(quanti.var$cos2)

# Contributions to the dimensions
head(quanti.var$contrib)

fviz_mfa_var(res.mfa, "quanti.var", palette = "jco", 
             col.var.sup = "violet", repel = TRUE)

fviz_mfa_var(res.mfa, "quanti.var", palette = "jco", 
             col.var.sup = "violet", repel = TRUE,
             geom = c("point", "text"), legend = "bottom")

# Contributions to dimension 1
fviz_contrib(res.mfa, choice = "quanti.var", axes = 1, top = 20,
             palette = "jco")

# Contributions to dimension 2
fviz_contrib(res.mfa, choice = "quanti.var", axes = 2, top = 20,
             palette = "jco")

fviz_mfa_var(res.mfa, "quanti.var", col.var = "contrib", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             col.var.sup = "violet", repel = TRUE,
             geom = c("point", "text"))

# Color by cos2 values: quality on the factor map
fviz_mfa_var(res.mfa, col.var = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             col.var.sup = "violet", repel = TRUE)

fviz_cos2(res.mfa, choice = "quanti.var", axes = 1)

Graph of individuals

To get the results for individuals, type this:

ind <- get_mfa_ind(res.mfa)
ind

## Multiple Factor Analysis results for individuals 
##  ===================================================
##   Name                      Description                      
## 1 "$coord"                  "Coordinates"                    
## 2 "$cos2"                   "Cos2, quality of representation"
## 3 "$contrib"                "Contributions"                  
## 4 "$coord.partiel"          "Partial coordinates"            
## 5 "$within.inertia"         "Within inertia"                 
## 6 "$within.partial.inertia" "Within partial inertia"

To plot individuals, use the function fviz_mfa_ind() [in factoextra]. By default, individuals are colored in blue. However, like variables, it’s also possible to color individuals by their cos2 values:

fviz_mfa_ind(res.mfa, col.ind = "cos2", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE)

Individuals with similar profiles are close to each other on the factor map. The first axis, mainly opposes the wine 1DAM and, the wines 1VAU and 2ING. As described in the previous section, the first dimension represents the harmony and the intensity of wines. Thus, the wine 1DAM (positive coordinates) was evaluated as the most “intense” and “harmonious” contrary to wines 1VAU and 2ING (negative coordinates) which are the least “intense” and “harmonious”. The second axis is essentially associated with the two wines T1 and T2 characterized by a strong value of the variables Spice.before.shaking and Odor.intensity.before.shaking.

Most of the supplementary qualitative variable categories are close to the origin of the map. This result indicates that the concerned categories are not related to the first axis (wine “intensity” & “harmony”) or the second axis (wine T1 and T2).

The category Env4 has high coordinates on the second axis related to T1 and T2.

The category “Reference” is known to be related to an excellent wine-producing soil. As expected, our analysis demonstrates that the category “Reference” has high coordinates on the first axis, which is positively correlated with wines “intensity” and “harmony”.

Note that, it’s possible to color the individuals using any of the qualitative variables in the initial data table. To do this, the argument habillage is used in the fviz_mfa_ind() function. For example, if you want to color the wines according to the supplementary qualitative variable “Label”, type this:

fviz_mfa_ind(res.mfa, 
             habillage = "Label", # color by groups 
             palette = c("#00AFBB", "#E7B800", "#FC4E07"),
             addEllipses = TRUE, ellipse.type = "confidence", 
             repel = TRUE # Avoid text overlapping
             ) 

If you want to color individuals using multiple categorical variables at the same time, use the function fviz_ellipses() [in factoextra] as follow:

fviz_ellipses(res.mfa, c("Label", "Soil"), repel = TRUE)

Alternatively, you can specify categorical variable indices:

fviz_ellipses(res.mca, 1:2, geom = "point")

Graph of partial individuals

The results for individuals obtained from the analysis performed with a single group are named partial individuals. In other words, an individual considered from the point of view of a single group is called partial individual.

In the default fviz_mfa_ind() plot, for a given individual, the point corresponds to the mean individual or the center of gravity of the partial points of the individual. That is, the individual viewed by all groups of variables.

For a given individual, there are as many partial points as groups of variables.

The graph of partial individuals represents each wine viewed by each group and its barycenter. To plot the partial points of all individuals, type this:

fviz_mfa_ind(res.mfa, partial = "all") 

If you want to visualize partial points for wines of interest, let say c(“1DAM”, “1VAU”, “2ING”), use this:

fviz_mfa_ind(res.mfa, partial = c("1DAM", "1VAU", "2ING")) 

Red color represents the wines seen by only the odor variables; violet color represents the wines seen by only the visual variables, and so on.

The wine 1DAM has been described in the previous section as particularly “intense” and “harmonious”, particularly by the odor group: It has a high coordinate on the first axis from the point of view of the odor variables group compared to the point of view of the other groups.

From the odor group’s point of view, 2ING was more “intense” and “harmonious” than 1VAU but from the taste group’s point of view, 1VAU was more “intense” and “harmonious” than 2ING.

Graph of partial axes

The graph of partial axes shows the relationship between the principal axes of the MFA and the ones obtained from analyzing each group using either a PCA (for groups of continuous variables) or a MCA (for qualitative variables).

fviz_mfa_axes(res.mfa)

R packages

Install required packages as follow:

install.packages(c("FactoMineR", "factoextra"))

Load the packages:

library("FactoMineR")
library("factoextra")

Data format

We’ll use a subset of the wine data set available in FactoMineR package:

library("FactoMineR")
data(wine)
df <- wine[,c(1,2, 16, 22, 29, 28, 30,31)]
head(df[, 1:7], 4)

##           Label Soil Plante Acidity Harmony Intensity Overall.quality
## 2EL      Saumur Env1   2.00    2.11    3.14      2.86            3.39
## 1CHA     Saumur Env1   2.00    2.11    2.96      2.89            3.21
## 1FON Bourgueuil Env1   1.75    2.18    3.14      3.07            3.54
## 1VAU     Chinon Env2   2.30    3.18    2.04      2.46            2.46

To see the structure of the data, type this:

str(df)

The data contains 21 rows (wines, individuals) and 8 columns (variables):

• The first two columns are factors (categorical variables): label (Saumur, Bourgueil or Chinon) and soil (Reference, Env1, Env2 or Env4).
• The remaining columns are numeric (continuous variables).

R code

The function FAMD() [FactoMiner package] can be used to compute FAMD. A simplified format is :

FAMD (base, ncp = 5, sup.var = NULL, ind.sup = NULL, graph = TRUE)

• base : a data frame with n rows (individuals) and p columns (variables).
• ncp: the number of dimensions kept in the results (by default 5)
• sup.var: a vector indicating the indexes of the supplementary variables.
• ind.sup: a vector indicating the indexes of the supplementary individuals.
• graph : a logical value. If TRUE a graph is displayed.

To compute FAMD, type this:

library(FactoMineR)
res.famd <- FAMD(df, graph = FALSE)

The output of the FAMD() function is a list including :

print(res.famd)

## *The results are available in the following objects:
## 
##   name          description                             
## 1 "$eig"        "eigenvalues and inertia"               
## 2 "$var"        "Results for the variables"             
## 3 "$ind"        "results for the individuals"           
## 4 "$quali.var"  "Results for the qualitative variables" 
## 5 "$quanti.var" "Results for the quantitative variables"

Eigenvalues / Variances

The proportion of variances retained by the different dimensions (axes) can be extracted using the function get_eigenvalue() [factoextra package] as follow:

library("factoextra")
eig.val <- get_eigenvalue(res.famd)
head(eig.val)

##       eigenvalue variance.percent cumulative.variance.percent
## Dim.1      4.832            43.92                        43.9
## Dim.2      1.857            16.88                        60.8
## Dim.3      1.582            14.39                        75.2
## Dim.4      1.149            10.45                        85.6
## Dim.5      0.652             5.93                        91.6

The function fviz_eig() or fviz_screeplot() [factoextra package] can be used to draw the scree plot (the percentages of inertia explained by each FAMD dimensions):

fviz_screeplot(res.famd)

Graph of variables

var <- get_famd_var(res.famd)
var

## FAMD results for variables 
##  ===================================================
##   Name       Description                      
## 1 "$coord"   "Coordinates"                    
## 2 "$cos2"    "Cos2, quality of representation"
## 3 "$contrib" "Contributions"

# Coordinates of variables
head(var$coord)

# Cos2: quality of representation on the factore map
head(var$cos2)

# Contributions to the  dimensions
head(var$contrib)

# Plot of variables
fviz_famd_var(res.famd, repel = TRUE)

# Contribution to the first dimension
fviz_contrib(res.famd, "var", axes = 1)

# Contribution to the second dimension
fviz_contrib(res.famd, "var", axes = 2)

quanti.var <- get_famd_var(res.famd, "quanti.var")
quanti.var 

## FAMD results for quantitative variables 
##  ===================================================
##   Name       Description                      
## 1 "$coord"   "Coordinates"                    
## 2 "$cos2"    "Cos2, quality of representation"
## 3 "$contrib" "Contributions"

fviz_famd_var(res.famd, "quanti.var", repel = TRUE,
              col.var = "black")

fviz_famd_var(res.famd, "quanti.var", col.var = "contrib", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE)

# Color by cos2 values: quality on the factor map
fviz_famd_var(res.famd, "quanti.var", col.var = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE)

quali.var <- get_famd_var(res.famd, "quali.var")
quali.var 

## FAMD results for qualitative variable categories 
##  ===================================================
##   Name       Description                      
## 1 "$coord"   "Coordinates"                    
## 2 "$cos2"    "Cos2, quality of representation"
## 3 "$contrib" "Contributions"

fviz_famd_var(res.famd, "quali.var", col.var = "contrib", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
             )

Graph of individuals

To get the results for individuals, type this:

ind <- get_famd_ind(res.famd)
ind

## FAMD results for individuals 
##  ===================================================
##   Name       Description                      
## 1 "$coord"   "Coordinates"                    
## 2 "$cos2"    "Cos2, quality of representation"
## 3 "$contrib" "Contributions"

To plot individuals, use the function fviz_mfa_ind() [in factoextra]. By default, individuals are colored in blue. However, like variables, it’s also possible to color individuals by their cos2 and contribution values:

fviz_famd_ind(res.famd, col.ind = "cos2", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE)

In the plot above, the qualitative variable categories are shown in black. Env1, Env2, Env3 are the categories of the soil. Saumur, Bourgueuil and Chinon are the categories of the wine Label. If you don’t want to show them on the plot, use the argument invisible = "quali.var".

Individuals with similar profiles are close to each other on the factor map. For the interpretation, read more at Chapter @ref(multiple-correspondence-analysis) (MCA) and Chapter @ref(multiple-factor-analysis) (MFA).

Note that, it’s possible to color the individuals using any of the qualitative variables in the initial data table. To do this, the argument habillage is used in the fviz_famd_ind() function. For example, if you want to color the wines according to the supplementary qualitative variable “Label”, type this:

fviz_mfa_ind(res.famd, 
             habillage = "Label", # color by groups 
             palette = c("#00AFBB", "#E7B800", "#FC4E07"),
             addEllipses = TRUE, ellipse.type = "confidence", 
             repel = TRUE # Avoid text overlapping
             ) 

If you want to color individuals using multiple categorical variables at the same time, use the function fviz_ellipses() [in factoextra] as follow:

fviz_ellipses(res.famd, c("Label", "Soil"), repel = TRUE)

Alternatively, you can specify categorical variable indices:

fviz_ellipses(res.famd, 1:2, geom = "point")

R packages

Several functions from different packages are available in the R software for computing multiple correspondence analysis. These functions/packages include:

• MCA() function [FactoMineR package]
• dudi.mca() function [ade4 package]
• and epMCA() [ExPosition package]

No matter what function you decide to use, you can easily extract and visualize the MCA results using R functions provided in the factoextra R package.

Here, we’ll use FactoMineR (for the analysis) and factoextra (for ggplot2-based elegant visualization). To install the two packages, type this:

install.packages(c("FactoMineR", "factoextra"))

Load the packages:

library("FactoMineR")
library("factoextra")

Data format

We’ll use the demo data sets poison available in FactoMineR package:

data(poison)
head(poison[, 1:7], 3)

##   Age Time   Sick Sex   Nausea Vomiting Abdominals
## 1   9   22 Sick_y   F Nausea_y  Vomit_n     Abdo_y
## 2   5    0 Sick_n   F Nausea_n  Vomit_n     Abdo_n
## 3   6   16 Sick_y   F Nausea_n  Vomit_y     Abdo_y

This data is a result from a survey carried out on children of primary school who suffered from food poisoning. They were asked about their symptoms and about what they ate.

The data contains 55 rows (individuals) and 15 columns (variables). We’ll use only some of these individuals (children) and variables to perform the multiple correspondence analysis. The coordinates of the remaining individuals and variables on the factor map will be predicted from the previous MCA results.

In MCA terminology, our data contains :

Subset only active individuals and variables for multiple correspondence analysis:

poison.active <- poison[1:55, 5:15]
head(poison.active[, 1:6], 3)

##     Nausea Vomiting Abdominals   Fever   Diarrhae   Potato
## 1 Nausea_y  Vomit_n     Abdo_y Fever_y Diarrhea_y Potato_y
## 2 Nausea_n  Vomit_n     Abdo_n Fever_n Diarrhea_n Potato_y
## 3 Nausea_n  Vomit_y     Abdo_y Fever_y Diarrhea_y Potato_y

Data summary

The R base function summary() can be used to compute the frequency of variable categories. As the data table contains a large number of variables, we’ll display only the results for the first 4 variables.

Statistical summaries:

# Summary of the 4 first variables
summary(poison.active)[, 1:4]

##       Nausea      Vomiting   Abdominals     Fever   
##  Nausea_n:43   Vomit_n:33   Abdo_n:18   Fever_n:20  
##  Nausea_y:12   Vomit_y:22   Abdo_y:37   Fever_y:35

The summary() functions return the size of each variable category.

It’s also possible to plot the frequency of variable categories. The R code below, plots the first 4 columns:

for (i in 1:4) {
  plot(poison.active[,i], main=colnames(poison.active)[i],
       ylab = "Count", col="steelblue", las = 2)
  }

R code

The function MCA()[FactoMiner package] can be used. A simplified format is :

MCA(X, ncp = 5, graph = TRUE)

• X: a data frame with n rows (individuals) and p columns (categorical variables)
• ncp: number of dimensions kept in the final results.
• graph: a logical value. If TRUE a graph is displayed.

In the R code below, the MCA is performed only on the active individuals/variables :

res.mca <- MCA(poison.active, graph = FALSE)

The output of the MCA() function is a list including :

print(res.mca)

## **Results of the Multiple Correspondence Analysis (MCA)**
## The analysis was performed on 55 individuals, described by 11 variables
## *The results are available in the following objects:
## 
##    name              description                       
## 1  "$eig"            "eigenvalues"                     
## 2  "$var"            "results for the variables"       
## 3  "$var$coord"      "coord. of the categories"        
## 4  "$var$cos2"       "cos2 for the categories"         
## 5  "$var$contrib"    "contributions of the categories" 
## 6  "$var$v.test"     "v-test for the categories"       
## 7  "$ind"            "results for the individuals"     
## 8  "$ind$coord"      "coord. for the individuals"      
## 9  "$ind$cos2"       "cos2 for the individuals"        
## 10 "$ind$contrib"    "contributions of the individuals"
## 11 "$call"           "intermediate results"            
## 12 "$call$marge.col" "weights of columns"              
## 13 "$call$marge.li"  "weights of rows"

Eigenvalues / Variances

The proportion of variances retained by the different dimensions (axes) can be extracted using the function get_eigenvalue() [factoextra package] as follow:

library("factoextra")
eig.val <- get_eigenvalue(res.mca)
 # head(eig.val)

To visualize the percentages of inertia explained by each MCA dimensions, use the function fviz_eig() or fviz_screeplot() [factoextra package]:

fviz_screeplot(res.mca, addlabels = TRUE, ylim = c(0, 45))

Biplot

The function fviz_mca_biplot() [factoextra package] is used to draw the biplot of individuals and variable categories:

fviz_mca_biplot(res.mca, 
               repel = TRUE, # Avoid text overlapping (slow if many point)
               ggtheme = theme_minimal())

The plot above shows a global pattern within the data. Rows (individuals) are represented by blue points and columns (variable categories) by red triangles.

The distance between any row points or column points gives a measure of their similarity (or dissimilarity). Row points with similar profile are closed on the factor map. The same holds true for column points.

Graph of variables

var <- get_mca_var(res.mca)
var

## Multiple Correspondence Analysis Results for variables
##  ===================================================
##   Name       Description                  
## 1 "$coord"   "Coordinates for categories" 
## 2 "$cos2"    "Cos2 for categories"        
## 3 "$contrib" "contributions of categories"

# Coordinates
head(var$coord)
# Cos2: quality on the factore map
head(var$cos2)
# Contributions to the principal components
head(var$contrib)

fviz_mca_var(res.mca, choice = "mca.cor", 
            repel = TRUE, # Avoid text overlapping (slow)
            ggtheme = theme_minimal())

head(round(var$coord, 2), 4)

##          Dim 1 Dim 2 Dim 3 Dim 4 Dim 5
## Nausea_n  0.27  0.12 -0.27  0.03  0.07
## Nausea_y -0.96 -0.43  0.95 -0.12 -0.26
## Vomit_n   0.48 -0.41  0.08  0.27  0.05
## Vomit_y  -0.72  0.61 -0.13 -0.41 -0.08

fviz_mca_var(res.mca, 
             repel = TRUE, # Avoid text overlapping (slow)
             ggtheme = theme_minimal())

fviz_mca_var(res.mca, col.var="black", shape.var = 15,
             repel = TRUE)

head(var$cos2, 4)

##          Dim 1  Dim 2  Dim 3   Dim 4   Dim 5
## Nausea_n 0.256 0.0528 0.2527 0.00408 0.01947
## Nausea_y 0.256 0.0528 0.2527 0.00408 0.01947
## Vomit_n  0.344 0.2512 0.0107 0.11229 0.00413
## Vomit_y  0.344 0.2512 0.0107 0.11229 0.00413

# Color by cos2 values: quality on the factor map
fviz_mca_var(res.mca, col.var = "cos2",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE, # Avoid text overlapping
             ggtheme = theme_minimal())

# Change the transparency by cos2 values
fviz_mca_var(res.mca, alpha.var="cos2",
             repel = TRUE,
             ggtheme = theme_minimal())

library("corrplot")
corrplot(var$cos2, is.corr=FALSE)

# Cos2 of variable categories on Dim.1 and Dim.2
fviz_cos2(res.mca, choice = "var", axes = 1:2)

head(round(var$contrib,2), 4)

##          Dim 1 Dim 2 Dim 3 Dim 4 Dim 5
## Nausea_n  1.52  0.81  4.67  0.08  0.49
## Nausea_y  5.43  2.91 16.73  0.30  1.76
## Vomit_n   3.73  7.07  0.36  4.26  0.19
## Vomit_y   5.60 10.61  0.54  6.39  0.29

# Contributions of rows to dimension 1
fviz_contrib(res.mca, choice = "var", axes = 1, top = 15)
# Contributions of rows to dimension 2
fviz_contrib(res.mca, choice = "var", axes = 2, top = 15)

# Total contribution to dimension 1 and 2
fviz_contrib(res.mca, choice = "var", axes = 1:2, top = 15)

fviz_mca_var(res.mca, col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE, # avoid text overlapping (slow)
             ggtheme = theme_minimal()
             )

# Change the transparency by contrib values
fviz_mca_var(res.mca, alpha.var="contrib",
             repel = TRUE,
             ggtheme = theme_minimal())

Graph of individuals

ind <- get_mca_ind(res.mca)
ind

## Multiple Correspondence Analysis Results for individuals
##  ===================================================
##   Name       Description                       
## 1 "$coord"   "Coordinates for the individuals" 
## 2 "$cos2"    "Cos2 for the individuals"        
## 3 "$contrib" "contributions of the individuals"

# Coordinates of column points
head(ind$coord)
# Quality of representation
head(ind$cos2)
# Contributions
head(ind$contrib)

fviz_mca_ind(res.mca, col.ind = "cos2", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE, # Avoid text overlapping (slow if many points)
             ggtheme = theme_minimal())

# Cos2 of individuals
fviz_cos2(res.mca, choice = "ind", axes = 1:2, top = 20)
# Contribution of individuals to the dimensions
fviz_contrib(res.mca, choice = "ind", axes = 1:2, top = 20)

Color individuals by groups

The R code below colors the individuals by groups using the levels of the variable Vomiting. The argument habillage is used to specify the factor variable for coloring the individuals by groups. A concentration ellipse can be also added around each group using the argument addEllipses = TRUE. If you want a confidence ellipse around the mean point of categories, use ellipse.type = "confidence" The argument palette is used to change group colors.

fviz_mca_ind(res.mca, 
             label = "none", # hide individual labels
             habillage = "Vomiting", # color by groups 
             palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, ellipse.type = "confidence",
             ggtheme = theme_minimal()) 

Note that, to specify the value of the argument habillage, it’s also possible to use the index of the column as follow (habillage = 2). Additionally, you can provide an external grouping variable as follow: habillage = poison$Vomiting. For example:

# habillage = index of the column to be used as grouping variable
fviz_mca_ind(res.mca, habillage = 2, addEllipses = TRUE)
# habillage = external grouping variable
fviz_mca_ind(res.mca, habillage = poison$Vomiting, addEllipses = TRUE)

If you want to color individuals using multiple categorical variables at the same time, use the function fviz_ellipses() [in factoextra] as follow:

fviz_ellipses(res.mca, c("Vomiting", "Fever"),
              geom = "point")

Alternatively, you can specify categorical variable indices:

fviz_ellipses(res.mca, 1:4, geom = "point")

Dimension description

The function dimdesc() [in FactoMineR] can be used to identify the most correlated variables with a given dimension:

res.desc <- dimdesc(res.mca, axes = c(1,2))
# Description of dimension 1
res.desc[[1]]
# Description of dimension 2
res.desc[[2]]

Definition and types

As described above (section @ref(mca-data-format)), the data set poison contains:

• supplementary continuous variables (quanti.sup = 1:2, columns 1 and 2 corresponding to the columns age and time, respectively)
• supplementary qualitative variables (quali.sup = 3:4, corresponding to the columns Sick and Sex, respectively). This factor variables are used to color individuals by groups

The data doesn’t contain supplementary individuals. However, for demonstration, we’ll use the individuals 53:55 as supplementary individuals.

Specification in MCA

To specify supplementary individuals and variables, the function MCA() can be used as follow :

MCA(X,  ind.sup = NULL, quanti.sup = NULL, quali.sup=NULL,
    graph = TRUE, axes = c(1,2))

• X: a data frame. Rows are individuals and columns are variables.
• ind.sup: a numeric vector specifying the indexes of the supplementary individuals.
• quanti.sup, quali.sup: a numeric vector specifying, respectively, the indexes of the quantitative and qualitative variables.
• graph: a logical value. If TRUE a graph is displayed.
• axes: a vector of length 2 specifying the components to be plotted.

For example, type this:

res.mca <- MCA(poison, ind.sup = 53:55, 
               quanti.sup = 1:2, quali.sup = 3:4,  graph=FALSE)

Results

The predicted results for supplementary individuals/variables can be extracted as follow:

# Supplementary qualitative variable categories
res.mca$quali.sup
# Supplementary quantitative variables
res.mca$quanti
# Supplementary individuals
res.mca$ind.sup

Plots

To make a biplot of individuals and variable categories, type this:

# Biplot of individuals and variable categories
fviz_mca_biplot(res.mca, repel = TRUE,
                ggtheme = theme_minimal())

If you want to highlight the correlation between variables (active & supplementary) and dimensions, use the function fviz_mca_var() with the argument choice = “mca.cor”:

fviz_mca_var(res.mca, choice = "mca.cor",
             repel = TRUE)

The R code below plots qualitative variable categories (active & supplementary variables):

fviz_mca_var(res.mca, repel = TRUE,
             ggtheme= theme_minimal())

For supplementary quantitative variables, type this:

fviz_mca_var(res.mca, choice = "quanti.sup",
             ggtheme = theme_minimal())

To visualize supplementary individuals, type this:

fviz_mca_ind(res.mca, 
             label = "ind.sup", #Show the label of ind.sup only
             ggtheme = theme_minimal())

Export plots to PDF/PNG files

Two steps:

1. Create the plot of interest as an R object:

# Scree plot
scree.plot <- fviz_eig(res.mca)
# Biplot of row and column variables
biplot.mca <- fviz_mca_biplot(res.mca)

2. Export the plots into a single pdf file as follow (one plot per page):

library(ggpubr)
ggexport(plotlist = list(scree.plot, biplot.mca), 
         filename = "MCA.pdf")

More options at: Chapter @ref(principal-component-analysis) (section: Exporting results).

Export results to txt/csv files

Easy to use R function: write.infile() [in FactoMineR] package.

# Export into a TXT file
write.infile(res.mca, "mca.txt", sep = "\t")
# Export into a CSV file
write.infile(res.mca, "mca.csv", sep = ";")