Correlation between variables and principal dimensions

Variables can be visualized as follow:

plot(res.mca, choix = "var")

Coordinates of variable categories

head(round(var$coord, 2))

         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
Abdo_n    1.32 -0.04 -0.01 -0.15 -0.07
Abdo_y   -0.64  0.02  0.00  0.07  0.03

Use the function fviz_mca_var() [in factoextra] to visualize only variable categories:

# Default plot
fviz_mca_var(res.mca)

It’s possible to change the color and the shape of the variable points using the arguments col.var and shape.var as follow:

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

# Hide individuals
plot(res.mca, invisible="ind") 

Contribution of variable categories to the dimensions

The contribution of the variable categories (in %) to the definition of the dimensions can be extracted as follow:

head(round(var$contrib,2))

         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
Abdo_n   15.42  0.03  0.00  0.73  0.18
Abdo_y    7.50  0.01  0.00  0.36  0.09

The variable categories with the larger value, contribute the most to the definition of the dimensions.

The different categories in the table are:

categories <- rownames(var$coord)
length(categories)

[1] 22

print(categories)

 [1] "Nausea_n"   "Nausea_y"   "Vomit_n"    "Vomit_y"    "Abdo_n"     "Abdo_y"     "Fever_n"   
 [8] "Fever_y"    "Diarrhea_n" "Diarrhea_y" "Potato_n"   "Potato_y"   "Fish_n"     "Fish_y"    
[15] "Mayo_n"     "Mayo_y"     "Courg_n"    "Courg_y"    "Cheese_n"   "Cheese_y"   "Icecream_n"
[22] "Icecream_y"

It’s possible to use the function corrplot to highlight the most contributing variables for each dimension:

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

The function fviz_contrib()[in factoextra] can be used to draw a bar plot of variable contributions:

# Contributions of variables on Dim.1
fviz_contrib(res.mca, choice = "var", axes = 1)

It can be seen that the categories Abdo_n, Diarrhea_n, Fever_n and Mayo_n are the most important in the definition of the first dimension.

# Contributions of rows on Dim.2
fviz_contrib(res.mca, choice = "var", axes = 2)

The row items Courg_n, Potato_n, Vomit_y and Icecream_n contribute the most to the dimension 2.

# Total contribution on Dim.1 and Dim.2
fviz_contrib(res.mca, choice = "var", axes = 1:2)

If your data contains many categories, the top contributing categories can be displayed as follow:

fviz_contrib(res.mca, choice = "var", axes = 1, top = 10)

Read more about fviz_contrib(): fviz_contrib

A second option is to draw a scatter plot of categories and to highlight categories according to the amount of their contributions. The function fviz_mca_var() is used.

Note that, using factoextra package, the color or the transparency of the variable categories can be automatically controlled by the value of their contributions, their cos2, their coordinates on x or y axis.

# Control category point colors using their contribution
# Possible values for the argument col.row are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_mca_var(res.mca, col.var = "contrib")

# Change the gradient color
fviz_mca_var(res.mca, col.var="contrib")+
scale_color_gradient2(low="white", mid="blue", 
                      high="red", midpoint=2)+theme_minimal()

It’s also possible to control automatically the transparency of variable categories by their contributions. The argument alpha.var is used:

# Control the transparency of categories using their contribution
# Possible values for the argument alpha.var are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_mca_var(res.mca, alpha.var="contrib")+
  theme_minimal()

It’s possible to select and display only the top contributing categories as illustrated in the R code below.

# Select the top 10 contributing categories
fviz_mca_var(res.mca, select.var=list(contrib=10))

Variable category/individual selections are discussed in details in the next sections

Read more about fviz_mca_var(): fviz_mca_var

Cos2 : The quality of representation of variable categories

The two dimensions 1 and 2 are sufficient to retain 46% of the total inertia contained in the data.

However, not all the points are equally well displayed in the two dimensions.

The quality of representation of the categories on the factor map is called the squared cosine (cos2) or the squared correlations.

The cos2 measures the degree of association between variable categories and a particular axis.

The cos2 of variable categories can be extracted as follow:

head(var$cos2)

             Dim 1        Dim 2        Dim 3       Dim 4       Dim 5
Nausea_n 0.2562007 0.0528025759 2.527485e-01 0.004084375 0.019466197
Nausea_y 0.2562007 0.0528025759 2.527485e-01 0.004084375 0.019466197
Vomit_n  0.3442016 0.2511603912 1.070855e-02 0.112294813 0.004126898
Vomit_y  0.3442016 0.2511603912 1.070855e-02 0.112294813 0.004126898
Abdo_n   0.8451157 0.0006215864 1.262496e-05 0.011479077 0.002374929
Abdo_y   0.8451157 0.0006215864 1.262496e-05 0.011479077 0.002374929

The values of the cos2 are comprised between 0 and 1.

The sum of the cos2 for rows on all the MCA dimensions is equal to one.

The quality of representation of a variable category or an individual in n dimensions is simply the sum of the squared cosine of that variable category or individual over the n dimensions.

If a variable category is well represented by two dimensions, the sum of the cos2 is closed to one.

For some of the categories, more than 2 dimensions are required to perfectly represent the data.

Visualize the cos2 of variable categories using corrplot:

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

The function fviz_cos2()[in factoextra] can be used to draw a bar plot of rows cos2:

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

Note that, variable categories Fish_n, Fish_y, Icecream_n and Icecream_y are not very well represented by the first two dimensions. This implies that the position of the corresponding points on the scatter plot should be interpreted with some caution. A higher dimensional solution is probably necessary.

Read more about fviz_cos2(): fviz_cos2

Coordinates of individuals

head(ind$coord)

       Dim 1       Dim 2       Dim 3       Dim 4       Dim 5
1 -0.4525811 -0.26415072  0.17151614  0.01369348 -0.11696806
2  0.8361700 -0.03193457 -0.07208249 -0.08550351  0.51978710
3 -0.4481892  0.13538726 -0.22484048 -0.14170168 -0.05004753
4  0.8803694 -0.08536230 -0.02052044 -0.07275873 -0.22935022
5 -0.4481892  0.13538726 -0.22484048 -0.14170168 -0.05004753
6 -0.3594324 -0.43604390 -1.20932223  1.72464616  0.04348157

Use the function fviz_mca_ind() [in factoextra] to visualize only column points:

fviz_mca_ind(res.mca)

Read more about fviz_mca_ind(): fviz_mca_ind

# Hide variable categories
plot(res.mca, invisible="var") 

Contribution of individuals to the dimensions

head(ind$contrib)

     Dim 1      Dim 2        Dim 3        Dim 4      Dim 5
1 1.110927 0.98238297  0.498254685  0.003555817 0.31554778
2 3.792117 0.01435818  0.088003703  0.138637089 6.23134138
3 1.089470 0.25806722  0.856229950  0.380768961 0.05776914
4 4.203611 0.10259105  0.007132055  0.100387990 1.21319013
5 1.089470 0.25806722  0.856229950  0.380768961 0.05776914
6 0.700692 2.67693398 24.769968729 56.404214518 0.04360547

Note that, you can use the previously mentioned corrplot() function to visualize the contribution of individuals.

Use the function fviz_contrib()[in factoextra] to visualize column contributions on dimensions 1+2:

fviz_contrib(res.mca, choice = "ind", axes = 1:2, top = 20)

Draw a scatter plot of individuals points and highlight individuals according to the amount of their contributions. The function fviz_mca_ind() [in factoextra] is used:

# Control individual colors using their contribution
# Possible values for the argument col.ind are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_mca_ind(res.mca, col.ind="contrib")+
scale_color_gradient2(low="white", mid="blue", 
                      high="red", midpoint=0.85)+theme_minimal()

# Control the transparency of individuals using their contribution
# Possible values for the argument alpha.ind are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_mca_ind(res.mca, alpha.ind="contrib")

Cos2 : The quality of representation of individuals

head(ind$cos2)

       Dim 1        Dim 2        Dim 3        Dim 4        Dim 5
1 0.34652591 0.1180447167 0.0497683175 0.0003172275 0.0231460846
2 0.55589562 0.0008108236 0.0041310808 0.0058126211 0.2148103098
3 0.54813888 0.0500176790 0.1379484860 0.0547920948 0.0068349171
4 0.74773962 0.0070299584 0.0004062504 0.0051072923 0.0507479873
5 0.54813888 0.0500176790 0.1379484860 0.0547920948 0.0068349171
6 0.02485357 0.0365775483 0.2813443706 0.5722083217 0.0003637178

Note that, the value of the cos2 is between 0 and 1. A cos2 closed to 1 corresponds to a variable categories/individuals that are well represented on the factor map.

The function fviz_cos2()[in factoextra] can be used to draw a bar plot of individuals cos2:

# Cos2 of individuals on Dim.1 and Dim.2
fviz_cos2(res.mca, choice = "ind", axes = 1:2, top = 20)

Change the color of individuals by groups

As mentioned above, our data contains supplementary qualitative variables: Columns 3 and 4 corresponding to the columns Sick and Sex, respectively. These factor variables will be used to color individuals by groups.

sick <- as.factor(poison$Sick)
head(sick)

[1] Sick_y Sick_n Sick_y Sick_n Sick_y Sick_y
Levels: Sick_n Sick_y

sex <- as.factor(poison$Sex)
head(sex)

[1] F F F F M M
Levels: F M

Individuals factor map :

# Default plot
fviz_mca_ind(res.mca, label ="none")

Change individual colors by groups using the levels of the variable sick. The argument habillage is used:

fviz_mca_ind(res.mca, label = "none", habillage=sick)

Add ellipses of point concentrations : the argument habillage is used to specify the factor variable for coloring the observations by groups.

fviz_mca_ind(res.mca, label="none", habillage = sick,
             addEllipses = TRUE, ellipse.level = 0.95)

Now, let’s :

• make a biplot of individuals and variable categories
• change the color of individuals by groups (sick levels)
• show only the labels for variables

fviz_mca_biplot(res.mca, 
  habillage = sick, addEllipses = TRUE,
  label = "var", shape.var = 15) +
  scale_color_brewer(palette="Dark2")+
  theme_minimal()

Note that, it’s possible to color the individuals using any of the qualitative variable in the initial data table (poison)

Let’s color the individuals by groups using the levels of the variable Vomiting:

fviz_mca_ind(res.mca, 
  habillage = poison$Vomiting, addEllipses = TRUE) +
  scale_color_brewer(palette="Dark2")+
  theme_minimal()

It’s also possible to use the index of the column as follow (habillage = 2):

fviz_mca_ind(res.mca, 
  habillage = 2, addEllipses = TRUE) +
  scale_color_brewer(palette="Dark2")+
  theme_minimal()

You can also use the function plotellipses() [in FactoMineR] to draw confidence ellipses around the categories. The simplified format is:

plotellipses(model, keepvar="all", axis =c(1,2))

• model: object of class MCA or PCA
• keppvar: a boolean or numeric vector of indexes of variables or a character vector of names of variables. If keepvar is “all”, “quali” or “quali.sup”, variables which are plotted are all the categorical variables, only those which are used to compute the dimensions (active variables) or only the supplementary categorical variables. If keepvar is a numeric vector of indexes or a character vector of names of variables, only relevant variables are plotted.

plotellipses(res.mca, keepvar=1)

plotellipses(res.mca, keepvar=1:4)

plotellipses(res.mca, keepvar="Vomiting")

plotellipses(res.mca, keepvar=c("Vomiting", "Fever"))

plotellipses(res.mca, keepvar="all")

Make a biplot of individuals and variable categories

FactomineR base graph:

plot(res.mca)

Use factoextra:

fviz_mca_biplot(res.mca) +
  theme_minimal()

Visualize supplementary variables

The graph below highlight the correlation between variables (active & supplementary) and dimensions:

plot(res.mca, choix ="var")

res.mca$quali.sup

$coord
             Dim 1         Dim 2       Dim 3        Dim 4       Dim 5
Sick_n  1.41809140  0.0020394048  0.13199139 -0.016036841 -0.08354663
Sick_y -0.63026284 -0.0009064021 -0.05866284  0.007127485  0.03713184
F      -0.03108147  0.1123143957  0.05033124 -0.055927173 -0.06832928
M       0.03356798 -0.1212995474 -0.05435774  0.060401347  0.07379562

$cos2
             Dim 1        Dim 2       Dim 3        Dim 4       Dim 5
Sick_n 0.893770319 1.848521e-06 0.007742990 0.0001143023 0.003102240
Sick_y 0.893770319 1.848521e-06 0.007742990 0.0001143023 0.003102240
F      0.001043342 1.362369e-02 0.002735892 0.0033780765 0.005042401
M      0.001043342 1.362369e-02 0.002735892 0.0033780765 0.005042401

$v.test
            Dim 1        Dim 2      Dim 3       Dim 4      Dim 5
Sick_n  6.7514655  0.009709509  0.6284047 -0.07635063 -0.3977615
Sick_y -6.7514655 -0.009709509 -0.6284047  0.07635063  0.3977615
F      -0.2306739  0.833551410  0.3735378 -0.41506855 -0.5071119
M       0.2306739 -0.833551410 -0.3735378  0.41506855  0.5071119

$eta2
           Dim 1        Dim 2       Dim 3        Dim 4       Dim 5
Sick 0.893770319 1.848521e-06 0.007742990 0.0001143023 0.003102240
Sex  0.001043342 1.362369e-02 0.002735892 0.0033780765 0.005042401

fviz_mca_var(res.mca) + theme_minimal()

# Hide active variables
fviz_mca_var(res.mca, invisible ="var") +
  theme_minimal()

# Hide supplementary qualitative variables
fviz_mca_var(res.mca, invisible ="quali.sup") +
  theme_minimal()

res.mca$quanti

$coord
            Dim 1       Dim 2       Dim 3       Dim 4       Dim 5
Age   0.003934896 -0.00741340 -0.26494536  0.20015501  0.02928483
Time -0.838158507 -0.08330586 -0.08718851 -0.08421599 -0.02316931

plot(res.mca, choix="quanti.sup")

Visualize supplementary individuals

The results for supplementary individuals can be extracted as follow :

res.mca$ind.sup

$coord
        Dim 1     Dim 2      Dim 3      Dim 4      Dim 5
53  1.0835684 0.5172478  0.5794063  0.5390903  0.4553650
54 -0.1249473 0.1417271 -0.1765234 -0.1526587 -0.2779565
55 -0.4315948 0.1270468 -0.2071580 -0.1186804 -0.1891760

$cos2
        Dim 1      Dim 2      Dim 3      Dim 4      Dim 5
53 0.36304957 0.08272764 0.10380536 0.08986204 0.06411692
54 0.03157652 0.04062716 0.06302535 0.04713607 0.15626590
55 0.50232519 0.04352713 0.11572730 0.03798314 0.09650827

Factor map for individuals:

fviz_mca_ind(res.mca) +
  theme_minimal()

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

Supplementary individuals are shown in darkblue.

Visual inspection

The above contingency table is not very large. Therefore, it’s easy to visually inspect and interpret row and column profiles:

• It’s evident that, the housetasks - Laundry, Main_Meal and Dinner - are more frequently done by the “Wife”.
• Repairs and driving are dominantly done by the husband
• Holidays are frequently associated with the column “jointly”

Visualize a contingency table using graphical matrix

It’s also possible to visualize a contingency table using the function balloonplot() [in gplots package]. This function draws a graphical matrix where each cell contains a dot whose size reflects the relative magnitude of the corresponding component.

To execute the R code below, you should install the package gplots: install.packages(“gplots”).

library("gplots")
# 1. convert the data as a table
dt <- as.table(as.matrix(housetasks))
# 2. Graph
balloonplot(t(dt), main ="housetasks", xlab ="", ylab="",
            label = FALSE, show.margins = FALSE)

Note that, row and column sums are printed by default in the bottom and right margins, respectively. These values can be hidden using the argument show.margins = FALSE.

Mosaic / association plots

The function mosaicplot() from the built-in R package garphics can be used also to visualize a contingency table.

library("graphics")
mosaicplot(dt, shade = TRUE, las=2,
           main = "housetasks")

• The argument shade is used to color the graph
• The argument las = 2 produces vertical labels

The surface of an element of the mosaic reflects the relative magnitude of its value.

• Blue color indicates that the observed value is higher than the expected value if the data were random
• Red color specifies that the observed value is lower than the expected value if the data were random

From this mosaic plot, it can be seen that the housetasks Laundry, Main_meal, Dinner and breakfeast (blue color) are mainly done by the wife in our example.

It’s also possible to use the package vcd to make a mosaic plot (function mosaic()) or an association plot (function assoc()).

# install.packages("vcd")
library("vcd")
# plot just a subset of the table
assoc(head(dt), shade = T, las=3)

Chi-square statistic

Another method to analyse a frequency table is to use the Chi-square test of independence. The Chi-square test evaluates whether there is a significant dependence between row and column categories.

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 = 1944.456, df = 36, p-value < 2.2e-16

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

Read more: correspondence analysis basics

Significance of the association between rows and columns

To interpret correspondence analysis, the first step is to evaluate whether there is a significant dependency between the rows and columns.

There are two methods to inspect the significance:

1. Using the trace
2. Using the Chi-square statistic

The trace is the the total inertia of the table (i.e, the sum of the eigenvalues). The square root of the trace is interpreted as the correlation coefficient between rows and columns.

The correlation coefficient is calculated as follow:

eig <- get_eigenvalue(res.ca)
trace <- sum(eig$eigenvalue) 
cor.coef <- sqrt(trace)
cor.coef

[1] 1.055907

Note that, as a rule of thumb 0.2 is the threshold above which the correlation can be considered as important (Bendixen 1995, 576; Healey 2013, 289-290).

In our example, the correlation coefficient is 1.0559074 indicating a strong association between row and column variables.

A more rigorous method is to use the chi-square statistic for examining the association. This appears at the top of the report generated by the function summary.CA(). A high chi-square statistic means strong link between row and column variables.

In our example, the association is highly significant (chi-square: 1944.456, p = 0).

# Chi-square statistics
chi2 <- trace*sum(as.matrix(housetasks))
chi2

[1] 1944.456

# Degree of freedom
df <- (nrow(housetasks) - 1) * (ncol(housetasks) - 1)
# P-value
pval <- pchisq(chi2, df = df, lower.tail = FALSE)
pval

[1] 0

Eigenvalues and scree plot

How many dimensions are sufficient for the data interpretation?

The number of dimensions to retain in the solution can be determined by examining the table of eigenvalues.

As mentioned above, trace is the total sum of eigenvalues. For a given axis, the ratio of the axis eigenvalue to the trace is called the percentage of variance (or total inertia or chi-square value) explained by that axis.

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

eigenvalues <- get_eigenvalue(res.ca)
head(round(eigenvalues, 2))

      eigenvalue variance.percent cumulative.variance.percent
Dim.1       0.54            48.69                       48.69
Dim.2       0.45            39.91                       88.60
Dim.3       0.13            11.40                      100.00
Dim.4       0.00             0.00                      100.00

There is no “rule of thumb” to choose the number of dimension to keep for the data interpretation. It depends on the research question and the researcher’s need. For example, if you are satisfied with 80% of the total inertia explained then use the number of dimensions necessary to achieve that.

Another method is to visually inspect the scree plot in which dimensions are ordered decreasingly according the amount of explained inertia.

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

fviz_screeplot(res.ca)

The point at which the scree plot shows a bend (so called “elbow”) can be considered as indicating an optimal dimensionality.

It’s also possible to calculate an average eigenvalue above which the axis should be kept in the solution.

Any axis with a contribution larger than the maximum of these two percentages should be considered as important and included in the solution for the interpretation of the data (see, Bendixen 1995, 577).

The R code below, draws the scree plot with a red dashed line specifying the average eigenvalue:

fviz_screeplot(res.ca) +
 geom_hline(yintercept=33.33, linetype=2, color="red")

According to the graph above, only dimensions 1 and 2 should be used in the solution. The dimension 3 explains only 11.4% of the total inertia which is below the average eigeinvalue (33.33%) and too little to be kept for further analysis.

Note that, you can use more than 2 dimensions. However, the supplementary dimensions are unlikely to contribute significantly to the interpretation of nature of the association between the rows and columns.

Dimensions 1 and 2 explain approximately 48.7% and 39.9% of the total inertia respectively. This corresponds to a cumulative total of 88.6% of total inertia retained by the 2 dimensions.

The higher the retention, the more subtlety in the original data is retained in the low-dimensional solution (Mike Bendixen, 2003).

Read more about eigenvalues and screeplot: Eigenvalues data visualization

CA scatter plot: Biplot of row and column variables

The function plot.CA()[in FactoMineR] can be used to plot the coordinates of rows and columns presented in the correspondence analysis output.

A simplified format is :

plot.CA(x, axes = c(1,2), col.row = "blue", col.col = "red")

FactoMineR base graph for CA:

plot(res.ca)

It’s also possible to use the function fviz_ca_biplot()[in factoextra package] to draw a nice looking plot:

fviz_ca_biplot(res.ca)

# Change the theme
fviz_ca_biplot(res.ca) +
  theme_minimal()

Read more about fviz_ca_biplot(): fviz_ca_biplot

The graph above is called symetric plot and shows a global pattern within the data. Rows are represented by blue points and columns 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.

The next step for the interpretation is to determine which row and column variables contribute the most in the definition of the different dimensions retained in the model.

Row variables

The function get_ca_row()[in factoextra] is used to extract the results for row variables. This function returns a list containing the coordinates, the cos2, the contribution and the inertia of row variables:

row <- get_ca_row(res.ca)
row

Correspondence Analysis - Results for rows
 ===================================================
  Name       Description                
1 "$coord"   "Coordinates for the rows" 
2 "$cos2"    "Cos2 for the rows"        
3 "$contrib" "contributions of the rows"
4 "$inertia" "Inertia of the rows"      

head(row$coord)

                Dim 1      Dim 2       Dim 3
Laundry    -0.9918368  0.4953220 -0.31672897
Main_meal  -0.8755855  0.4901092 -0.16406487
Dinner     -0.6925740  0.3081043 -0.20741377
Breakfeast -0.5086002  0.4528038  0.22040453
Tidying    -0.3938084 -0.4343444 -0.09421375
Dishes     -0.1889641 -0.4419662  0.26694926

# Default plot
fviz_ca_row(res.ca)

fviz_ca_row(res.ca, col.row="steelblue", shape.row = 15)

# Hide columns
plot(res.ca, invisible="col") 

head(row$contrib)

                Dim 1    Dim 2    Dim 3
Laundry    18.2867003 5.563891 7.968424
Main_meal  12.3888433 4.735523 1.858689
Dinner      5.4713982 1.321022 2.096926
Breakfeast  3.8249284 3.698613 3.069399
Tidying     1.9983518 2.965644 0.488734
Dishes      0.4261663 2.844117 3.634294

library("corrplot")
corrplot(row$contrib, is.corr=FALSE)

# Contributions of rows on Dim.1
fviz_contrib(res.ca, choice = "row", axes = 1)

# Contributions of rows on Dim.2
fviz_contrib(res.ca, choice = "row", axes = 2)

# Total contribution on Dim.1 and Dim.2
fviz_contrib(res.ca, choice = "row", axes = 1:2)

fviz_contrib(res.ca, choice = "row", axes = 1, top = 5)

# Control row point colors using their contribution
# Possible values for the argument col.row are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_ca_row(res.ca, col.row = "contrib")

# Change the gradient color
fviz_ca_row(res.ca, col.row="contrib")+
scale_color_gradient2(low="white", mid="blue", 
                      high="red", midpoint=10)+theme_minimal()

# Control the transparency of rows using their contribution
# Possible values for the argument alpha.var are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_ca_row(res.ca, alpha.row="contrib")+
  theme_minimal()

# Select the top 5 contributing rows
fviz_ca_row(res.ca, alpha.row="contrib", select.row=list(contrib=5))

head(row$cos2)

               Dim 1     Dim 2      Dim 3
Laundry    0.7399874 0.1845521 0.07546047
Main_meal  0.7416028 0.2323593 0.02603787
Dinner     0.7766401 0.1537032 0.06965666
Breakfeast 0.5049433 0.4002300 0.09482670
Tidying    0.4398124 0.5350151 0.02517249
Dishes     0.1181178 0.6461525 0.23572969

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

# Cos2 of rows on Dim.1 and Dim.2
fviz_cos2(res.ca, choice = "row", axes = 1:2)

Column varables

The function get_ca_col()[in factoextra] is used to extract the results for column variables. This function returns a list containing the coordinates, the cos2, the contribution and the inertia of columns variables:

col <- get_ca_col(res.ca)
col

Correspondence Analysis - Results for columns
 ===================================================
  Name       Description                   
1 "$coord"   "Coordinates for the columns" 
2 "$cos2"    "Cos2 for the columns"        
3 "$contrib" "contributions of the columns"
4 "$inertia" "Inertia of the columns"      

The result for columns gives the same information as described for rows. For this reason, I’ll just displayed the result for columns in this section without commenting.

head(col$coord)

                  Dim 1      Dim 2       Dim 3
Wife        -0.83762154  0.3652207 -0.19991139
Alternating -0.06218462  0.2915938  0.84858939
Husband      1.16091847  0.6019199 -0.18885924
Jointly      0.14942609 -1.0265791 -0.04644302

fviz_ca_col(res.ca)

# Hide rows
plot(res.ca, invisible="row") 

head(col$contrib)

                Dim 1     Dim 2      Dim 3
Wife        44.462018 10.312237 10.8220753
Alternating  0.103739  2.782794 82.5492464
Husband     54.233879 17.786612  6.1331792
Jointly      1.200364 69.118357  0.4954991

fviz_contrib(res.ca, choice = "col", axes = 1:2)

# Control column point colors using their contribution
# Possible values for the argument col.col are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_ca_col(res.ca, col.col="contrib")+
scale_color_gradient2(low="white", mid="blue", 
                      high="red", midpoint=24.5)+theme_minimal()

# Control the transparency of rows using their contribution
# Possible values for the argument alpha.col are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_ca_col(res.ca, alpha.col="contrib")

Cos2 : The quality of representation of columns

head(col$cos2)

                  Dim 1     Dim 2       Dim 3
Wife        0.801875947 0.1524482 0.045675847
Alternating 0.004779897 0.1051016 0.890118521
Husband     0.772026244 0.2075420 0.020431728
Jointly     0.020705858 0.9772939 0.002000236

Note that, the value of the cos2 is between 0 and 1. A cos2 closed to 1 corresponds to a column/row variables that are well represented on the factor map.

The function fviz_cos2() [in factoextra] can be used to draw a bar plot of columns cos2:

# Cos2 of columns on Dim.1 and Dim.2
fviz_cos2(res.ca, choice = "col", axes = 1:2)

Note that, only the column item Alternating is not very well displayed on the first two dimensions. The position of this item must be interpreted with caution in the space formed by dimensions 1 and 2.

Symmetric biplot

As mentioned above, the standard plot of correspondence analysis is a symmetric biplot in which both rows (blue points) and columns (red triangles) are represented in the same space using the principal coordinates. These coordinates represent the row and column profiles. In this case, only the distance between row points or the distance between column points can be really interpreted.

With symmetric plot, the inter-distance between rows and columns can’t be interpreted. Only a general statements can be made about the pattern.

fviz_ca_biplot(res.ca)+
  theme_minimal()

Remove the points from the graph, use texts only :

fviz_ca_biplot(res.ca, geom="text")

Note that, in order to interpret the distance between column points and row points, the simplest way is to make an asymmetric plot (Bendixen, 2003). This means that, the column profiles must be presented in row space or vice-versa.

Read more about fviz_ca_biplot(): fviz_ca_biplot

Asymmetric biplot for correspondence analysis

To make an asymetric plot, rows (or columns) points are plotted from the standard co-ordinates (S) and the profiles of the columns (or the rows) are plotted from the principale coordinates (P) (Bendixen 2003).

Depending on the situation, other types of display can be set using the argument map for the function fviz_ca_biplot()[in factoextra]. This is inspired from ca package (Michael Greenacre).

The allowed options for the argument map are:

1. “rowprincipal” or “colprincipal” - these are the so-called asymmetric biplots, with either rows in principal coordinates and columns in standard coordinates, or vice versa (also known as row-metric-preserving or column-metric-preserving respectively).

• “rowprincipal”: columns are represented in row space
• “colprincipal”: rows are represented in column space

2. “symbiplot” - both rows and columns are scaled to have variances equal to the singular values (square roots of eigenvalues), which gives a symmetric biplot but does not preserve row or column metrics.

3. “rowgab” or “colgab”: Asymetric maps proposed by Gabriel & Odoroff (1990):

• “rowgab”: rows in principal coordinates and columns in standard coordinates multiplied by the mass.
• “colgab”: columns in principal coordinates and rows in standard coordinates multiplied by the mass.

4. “rowgreen” or “colgreen”: The so-called contribution biplots showing visually the most contributing points (Greenacre 2006b).

• “rowgreen”: rows in principal coordinates and columns in standard coordinates multiplied by square root of the mass.
• “colgreen”: columns in principal coordinates and rows in standard coordinates multiplied by the square root of the mass.

The R code below draw a standard asymetric biplot:

fviz_ca_biplot(res.ca, map ="rowprincipal", arrow = c(TRUE, TRUE))

The argument arrows is a vector of two logicals specifying if the plot should contain points (FALSE, default) or arrows (TRUE). First value sets the rows and the second value sets the columns.

Contribution biplot

In correspondence analysis, biplot is a graphical display of rows and columns in 2 or 3 dimensions.

In the standard symmetric biplot (mentioned in the previous sections), it’s difficult to know the most contributing points to the solution of the CA.

Michael Greenacre proposed a new scaling displayed (called contribution biplot) which incorporates the contribution of points. In this display, points that contribute very little to the solution, are close to the center of the biplot and are relatively unimportant to the interpretation.

A contribution biplot can be drawn using the argument map = “rowgreen” or map = “colgreen”.

Firstly, you have to decide whether to analyse the contributions of rows or columns to the definition of the axes.

In our example we’ll interpret the contribution of rows to the axes. The argument map =“colgreen” is used. In this case, remember that columns are in principal coordinates and rows in standard coordinates multiplied by the square root of the mass. For a given row, the square of the new coordinate on an axis i is exactly the contribution of this row to the inertia of the axis i.

fviz_ca_biplot(res.ca, map ="colgreen",
               arrow = c(TRUE, FALSE))

In the graph above, the position of the column profile points is unchanged relative to that in the conventional biplot. However, the distances of the row points from the plot origin are related to their contributions to the two-dimensional factor map.

The closer an arrow is (in terms of angular distance) to an axis the greater is the contribution of the row category on that axis relative to the other axis. If the arrow is halfway between the two, its row category contributes to the two axes to the same extent.

Plot rows or columns only

It’s also possible to draw the rows or columns only using the function fviz_ca_biplot() (instead of using fviz_ca_row() and fviz_ca_col)

Plot rows only by hiding the columns (invisible =“col”):

fviz_ca_biplot(res.ca, invisible = "col")+
  theme_minimal()

Plot columns only by hiding the rows (invisible =“row”):

fviz_ca_biplot(res.ca, invisible = "row")+
  theme_minimal()

Data

We’ll use the data set children [in FactoMineR package]. It contains 18 rows and 8 columns:

data(children)
# head(children)

The data used here is a contingency table describing the answers given by different categories of people to the following question: What are the reasons that can make hesitate a woman or a couple to have children?

In CA terminology, our data contains :

CA with supplementary rows/columns

As mentioned above, supplementary rows and columns are not used for the definition of the principal dimensions. Their coordinates are predicted using only the informations provided by the performed CA on active rows/columns.

To specify supplementary rows/columns, the function CA()[in FactoMineR] can be used as follow :

CA(X,  ncp = 5, row.sup = NULL, col.sup = NULL,
   graph = TRUE)

Example of usage :

res.ca <- CA (children, row.sup = 15:18, col.sup = 6:8,
              graph = FALSE)

The summary of the CA is :

summary(res.ca, nb.dec = 2, ncp = 2)

Call:
rmarkdown::render("factominer-correspondance-analysis.Rmd", encoding = "UTF-8") 

The chi square of independence between the two variables is equal to 98.80159 (p-value =  9.748064e-05 ).

Eigenvalues
                      Dim.1  Dim.2  Dim.3  Dim.4  Dim.5
Variance               0.04   0.01   0.01   0.01   0.00
% of var.             57.04  21.13  11.76  10.06   0.00
Cumulative % of var.  57.04  78.17  89.94 100.00 100.00

Rows (the 10 first)
                      Dim.1   ctr  cos2   Dim.2   ctr  cos2  
money               | -0.12  4.55  0.43 |  0.02  0.37  0.01 |
future              |  0.18 17.57  0.72 | -0.10 14.59  0.22 |
unemployment        | -0.21 22.62  0.87 | -0.07  6.78  0.10 |
circumstances       |  0.40  6.27  0.58 |  0.33 11.54  0.40 |
hard                | -0.25  2.99  0.88 |  0.07  0.59  0.06 |
economic            |  0.35 12.00  0.48 |  0.32 26.60  0.40 |
egoism              |  0.06  0.68  0.07 | -0.03  0.34  0.01 |
employment          | -0.14  2.62  0.16 |  0.22 17.55  0.41 |
finances            | -0.24  2.79  0.28 | -0.21  5.69  0.21 |
war                 |  0.22  2.17  0.75 | -0.07  0.69  0.09 |

Columns
                      Dim.1   ctr  cos2   Dim.2   ctr  cos2  
unqualified         | -0.21 25.11  0.68 | -0.08 10.08  0.10 |
cep                 | -0.14 18.30  0.64 |  0.06  8.08  0.11 |
bepc                |  0.11  6.76  0.31 | -0.03  1.25  0.02 |
high_school_diploma |  0.27 37.98  0.76 | -0.12 20.10  0.15 |
university          |  0.23 11.86  0.31 |  0.32 60.49  0.59 |

Supplementary rows
                      Dim.1 cos2   Dim.2 cos2  
comfort             |  0.21 0.07 |  0.70 0.78 |
disagreement        |  0.15 0.13 |  0.12 0.09 |
world               |  0.52 0.88 |  0.14 0.07 |
to_live             |  0.31 0.14 |  0.50 0.37 |

Supplementary columns
                      Dim.1  cos2   Dim.2  cos2  
thirty              |  0.11  0.14 | -0.06  0.04 |
fifty               | -0.02  0.01 |  0.05  0.09 |
more_fifty          | -0.18  0.29 | -0.05  0.02 |

For the supplementary rows/columns, the coordinates and the quality of representation (cos2) on the factor maps are displayed. They don’t contribute to the dimensions.

Make a biplot of rows and columns

FactomineR base graph:

plot(res.ca)

Use factoextra:

fviz_ca_biplot(res.ca) +
  theme_minimal()

It’s also possible to hide supplementary rows and columns using the argument invisible:

fviz_ca_biplot(res.ca, invisible = c("row.sup", "col.sup") ) +
  theme_minimal()

The argument invisible is also available in FactoMineR base graph.

Visualize supplementary rows

All the results (coordinates and cos2) for the supplementary rows can be extracted as follow :

res.ca$row.sup

$coord
                 Dim 1     Dim 2      Dim 3      Dim 4
comfort      0.2096705 0.7031677 0.07111168  0.3071354
disagreement 0.1462777 0.1190106 0.17108916 -0.3132169
world        0.5233045 0.1429707 0.08399269 -0.1063597
to_live      0.3083067 0.5020193 0.52093397  0.2557357

$cos2
                  Dim 1      Dim 2       Dim 3      Dim 4
comfort      0.06892759 0.77524032 0.007928672 0.14790342
disagreement 0.13132177 0.08692632 0.179649183 0.60210272
world        0.87587685 0.06537746 0.022564054 0.03618163
to_live      0.13899699 0.36853645 0.396830367 0.09563620

Factor map for rows :

fviz_ca_row(res.ca) +
  theme_minimal()

Supplementary rows are shown in darkblue color.

Visualize supplementary columns

Factor map for columns:

fviz_ca_col(res.ca) +
  theme_minimal()

Supplementary columns are shown in darkred.

The results for supplementary columns can be extracted as follow :

res.ca$col.sup

$coord
                 Dim 1       Dim 2       Dim 3       Dim 4
thirty      0.10541339 -0.05969594 -0.10322613  0.06977996
fifty      -0.01706444  0.04907657 -0.01568923 -0.01306117
more_fifty -0.17706810 -0.04813788  0.10077299 -0.08517528

$cos2
               Dim 1      Dim 2       Dim 3       Dim 4
thirty     0.1375601 0.04411543 0.131910759 0.060278490
fifty      0.0108695 0.08990298 0.009188167 0.006367804
more_fifty 0.2860989 0.02114509 0.092666735 0.066200714

Extract the eigenvalues

Eigenvalues measure the amount of variation retained by a principal axis :

summary(res.ca)

Class: coa dudi
Call: dudi.coa(df = housetasks, scannf = FALSE, nf = 5)

Total inertia: 1.115

Eigenvalues:
    Ax1     Ax2     Ax3 
 0.5429  0.4450  0.1270 

Projected inertia (%):
    Ax1     Ax2     Ax3 
  48.69   39.91   11.40 

Cumulative projected inertia (%):
    Ax1   Ax1:2   Ax1:3 
  48.69   88.60  100.00 

You can also use the function get_eigenvalue() [in factoextra package] to extract the eigenvalues :

eig.val <- get_eigenvalue(res.ca)
head(eig.val)

      eigenvalue variance.percent cumulative.variance.percent
Dim.1  0.5428893         48.69222                    48.69222
Dim.2  0.4450028         39.91269                    88.60491
Dim.3  0.1270484         11.39509                   100.00000

Make a scree plot using ade4 base graphics

The function screeplot() can be used to draw the amount of inertia (variance) retained by the dimensions.

A simplified format is:

screeplot(x, ncps = length(x$eig), type = c("barplot", "lines"))

Example of usage :

screeplot(res.ca, main ="Screeplot - Eigenvalues")

~89% of the information contained in the data are retained by the first two dimensions.

Make the scree plot using factoextra

It’s also possible to use the function fviz_screeplot() [in factoextra] to make the scree plot. In the R code below, we’ll draw the percentage of variances retained by each component :

fviz_screeplot(res.ca, ncp=3)

Read more about eigenvalues and screeplot: Eigenvalues data visualization

Coordinates of rows

The coordinates of the rows on the factor map are :

head(res.ca$li)

               Axis1      Axis2       Axis3
Laundry    0.9918368 -0.4953220 -0.31672897
Main_meal  0.8755855 -0.4901092 -0.16406487
Dinner     0.6925740 -0.3081043 -0.20741377
Breakfeast 0.5086002 -0.4528038  0.22040453
Tidying    0.3938084  0.4343444 -0.09421375
Dishes     0.1889641  0.4419662  0.26694926

Use the function fviz_ca_row() [in factoextra package] to visualize only row points:

# Default plot
fviz_ca_row(res.ca)

s.label(res.ca$li, xax = 1, yax = 2)

Contribution of rows to the dimensions

The cos2 and the contributions of rows / columns are calculated using the function inertia.dudi() as follow :

inertia <- inertia.dudi(res.ca, row.inertia = TRUE,
                        col.inertia = TRUE)

Note that, the contributions and the cos2 are printed in 1/10 000. The sign is the sign of the coordinates.

The contributions can be printed in % as follow :

# absolute contribution of columns
contrib <- inertia$col.abs/100
head(contrib)

            Comp1 Comp2 Comp3
Wife        44.46 10.31 10.82
Alternating  0.10  2.78 82.55
Husband     54.23 17.79  6.13
Jointly      1.20 69.12  0.50

Recall that, as mentioned above, the simplest way is to use the function get_ca_row() [in factoextra package]. It provides a list of matrices containing all the results for the active rows(coordinates, squared cosine and contributions).

row <- get_ca_row(res.ca)
row

Correspondence Analysis - Results for rows
 ===================================================
  Name       Description                
1 "$coord"   "Coordinates for the rows" 
2 "$cos2"    "Cos2 for the rows"        
3 "$contrib" "contributions of the rows"
4 "$inertia" "Inertia of the rows"      

# Row contributions
row$contrib

           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

The row category with the largest value, contribute the most to the definition of the dimensions.

The function fviz_contrib()[in factoextra] can be used to visualize the most important row variables:

# Contributions of rows on Dim.1
fviz_contrib(res.ca, choice = "row", axes = 1)

The row items Repairs, Laundry, Main_meal and Driving contribute the most in the definition of the first axis.

# Contributions of rows on Dim.2
fviz_contrib(res.ca, choice = "row", axes = 2)

Read more about fviz_contrib(): fviz_contrib

Using factoextra package, the color of rows can be automatically controlled by the value of their contributions

fviz_ca_row(res.ca, col.row="contrib")+
scale_color_gradient2(low="white", mid="blue", 
                      high="red", midpoint=10)+theme_minimal()

The graph above highlight the most important rows in the correspondence analysis solution.

Read more about fviz_ca_row(): fviz_ca_row

Cos2 : quality of representation of rows on the factor map

• A high cos2 indicates a good representation of the rows on the factor map.
• A low cos2 indicates that the variable is not perfectly represented by the principal dimensions.

The cos2 of the rows are (factoextra code) :

head(row$cos2)

            Dim.1  Dim.2  Dim.3
Laundry    0.7400 0.1846 0.0755
Main_meal  0.7416 0.2324 0.0260
Dinner     0.7766 0.1537 0.0697
Breakfeast 0.5049 0.4002 0.0948
Tidying    0.4398 0.5350 0.0252
Dishes     0.1181 0.6462 0.2357

Note that, the ade4 code is:

# relative contributions of rows
cos2 <- abs(inertia$row.rel/10000)
head(cos2)

The values of the cos2 are comprised between 0 and 1.

The function fviz_cos2()[in factoextra] can be used to draw a bar plot of rows cos2:

# Cos2 of rows on Dim.1 and Dim.2
fviz_cos2(res.ca, choice = "row", axes = 1:2)

Note that, all row points except Official are well represented by the first two dimensions. The position of the point corresponding the item Official on the scatter plot should be interpreted with some caution.

Using factoextra package, the color of rows can be automatically controlled by the value of their cos2.

fviz_ca_row(res.ca, col.row="cos2")+
scale_color_gradient2(low="white", mid="blue", 
      high="red", midpoint=0.5) + theme_minimal()

Read more about fviz_cos2(): fviz_cos2

Coordinates of columns

The coordinates of the columns on the factor maps can be extracted as follow :

# ade4 code
head(res.ca$co)

                  Comp1      Comp2       Comp3
Wife         0.83762154 -0.3652207 -0.19991139
Alternating  0.06218462 -0.2915938  0.84858939
Husband     -1.16091847 -0.6019199 -0.18885924
Jointly     -0.14942609  1.0265791 -0.04644302

Use the function fviz_ca_col() [in factoextra] to visualize only column points:

fviz_ca_col(res.ca)

s.label(res.ca$co, xax = 1, yax = 2)

Contribution of columns

The contributions can be printed in % as follow :

# absolute contributions of columns
# ade4 code
contrib <- inertia$col.abs/100
head(contrib)

            Comp1 Comp2 Comp3
Wife        44.46 10.31 10.82
Alternating  0.10  2.78 82.55
Husband     54.23 17.79  6.13
Jointly      1.20 69.12  0.50

It’s simple to use the function get_ca_col() [from factoextra package]. factoextra provides, a list of matrices containing all the results for the active columns (coordinates, squared cosine and contributions)./span>

columns <- get_ca_col(res.ca)
columns

Correspondence Analysis - Results for columns
 ===================================================
  Name       Description                   
1 "$coord"   "Coordinates for the columns" 
2 "$cos2"    "Cos2 for the columns"        
3 "$contrib" "contributions of the columns"
4 "$inertia" "Inertia of the columns"      

# Contributions of columns
head(columns$contrib)

            Dim.1 Dim.2 Dim.3
Wife        44.46 10.31 10.82
Alternating  0.10  2.78 82.55
Husband     54.23 17.79  6.13
Jointly      1.20 69.12  0.50

Use the function fviz_contrib()[factoextra package] to visualize the most contributing columns :

# Contributions of columns on Dim.1
fviz_contrib(res.ca, choice = "col", axes = 1)

# Contributions of columns on Dim.2
fviz_contrib(res.ca, choice = "col", axes = 2)

Read more about fviz_contrib(): fviz_contrib

Draw a scatter plot of column points and highlight columns according to the amount of their contributions. The function fviz_ca_col() [in factoextra] is used:

# Control column point colors using their contribution
# Possible values for the argument col.col are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_ca_col(res.ca, col.col="contrib")+
scale_color_gradient2(low="white", mid="blue", 
                      high="red", midpoint=24.5)+theme_minimal()

Cos2 : The quality of representation of columns

# relative contributions of columns
cos2 <- abs(inertia$col.rel)/10000
head(cos2)

             Comp1  Comp2  Comp3 con.tra
Wife        0.8019 0.1524 0.0457  0.2700
Alternating 0.0048 0.1051 0.8901  0.1057
Husband     0.7720 0.2075 0.0204  0.3421
Jointly     0.0207 0.9773 0.0020  0.2823

The function fviz_cos2()[in factoextra] can be used to draw a bar plot of columns cos2:

# Cos2 of columns on Dim.1 and Dim.2
fviz_cos2(res.ca, choice = "col", axes = 1:2)

Note that, only the column item Alternating is not very well displayed on the first two dimensions. The position of this item must be interpreted with caution in the space formed by dimensions 1 and 2.

Read more about fviz_cos2(): fviz_cos2

Data

We’ll use the data set children available on STHDA website. It contains 18 rows and 8 columns:

ff <- "https://www.sthda.com/sthda/RDoc/data/ca-children.txt"
children <- read.table(file = ff, sep ="\t", 
                       header = TRUE, row.names = 1)

The data used here is a contingency table describing the answers given by different categories of people to the following question: What are the reasons that can make hesitate a woman or a couple to have children? (source of the data: FactoMineR package)

In CA terminology, our data contains :

R functions

The functions suprow() and supcol() [in ade4 package] are used to calculate the coordinates of supplementary rows and columns, respectively.

The simplified formats are :

# For supplementary rows
suprow(x, Xsup)

# For supplementary columns
supcol(x, Xsup)

Supplementary rows

# Data for the supplementary rows
row.sup <- children[15:18, 1:5, drop = FALSE]
head(row.sup)

             unqualified cep bepc high_school_diploma university
comfort                2   4    3                   1          4
disagreement           2   8    2                   5          2
world                  1   5    4                   6          3
to_live                3   3    1                   3          4

STEP 1/2 - CA using active rows/columns:

d.active <- children[1:14, 1:5]
res.ca <- dudi.coa(d.active, scannf = FALSE, nf =5)

STEP 2/2 - Predict the coordinates of the supplementary rows:

row.sup.ca <- suprow(res.ca, row.sup)
names(row.sup.ca)

[1] "tabsup" "lisup" 

# coordinates 
row.sup.coord <- row.sup.ca$lisup
head(row.sup.coord)

                 Axis1     Axis2      Axis3      Axis4
comfort      0.2096705 0.7031677 0.07111168  0.3071354
disagreement 0.1462777 0.1190106 0.17108916 -0.3132169
world        0.5233045 0.1429707 0.08399269 -0.1063597
to_live      0.3083067 0.5020193 0.52093397  0.2557357

How to visualize supplementary rows on the factor map?

The function fviz_add() is used :

# Plot of active rows
p <- fviz_ca_row(res.ca)
# Add supplementary rows
fviz_add(p, row.sup.coord, color ="darkgreen")

Supplementary columns

# Data for the supplementary quantitative variables
col.sup <- children[1:14, 6:8, drop = FALSE]
head(col.sup)

              thirty fifty more_fifty
money             59    66         70
future           115   117         86
unemployment      79    88        177
circumstances      9     8          5
hard               2    17         18
economic          18    19         17

Recall that, rows 15:18 are supplementary rows. We don’t want them in this current analysis. This is why, I extracted only rows 1:14.

Predict the coordinates of the supplementary columns :

col.sup.ca <- supcol(res.ca, col.sup)
names(col.sup.ca)

[1] "tabsup" "cosup" 

# coordinates 
col.sup.coord <- col.sup.ca$cosup
head(col.sup.coord)

                 Comp1       Comp2       Comp3       Comp4
thirty      0.10541339 -0.05969594 -0.10322613  0.06977996
fifty      -0.01706444  0.04907657 -0.01568923 -0.01306117
more_fifty -0.17706810 -0.04813788  0.10077299 -0.08517528

Visualize supplementary columns on the factor map using factoextra :

# Plot of active columns
p <- fviz_ca_col(res.ca)
# Add supplementary active variables
fviz_add(p, col.sup.coord , color ="darkgreen")

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:

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)

In the table above, the row TOTAL (in light blue) is called the average row profile (or marginal profile of columns or column margin)

The average row profile is computed as follow:

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.3440367   0.1456422   0.2184633   0.2918578 

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:

This distance is called Chi-square distance.

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.03684787

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.772028

Note that, the rows Laundry and Main_meal are very close (d2 ~ 0.036, similar profiles) compared to the rows Laundry and Driving (d2 ~ 3.77)

You can also compute the squared distance between each row profile and the average row profile in order to view rows that are the most similar or different to the average row.

Squared distance between each row profile and the average row 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

The rows Repairs, Holidays, Laundry and Driving have the most different profiles from the average profile.

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)

The size of the circle is proportional to the magnitude of the distance between row profiles.

When the data contains many categories, correspondence analysis is very useful to visualize the similarity between items.

Row mass and inertia

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

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.10091743 0.08772936 0.06192661 0.08027523 0.06995413 0.06479358 

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.inertia <- row.mass * d2.row
head(row.inertia)

   Laundry  Main_meal     Dinner Breakfeast    Tidying     Dishes 
0.13415976 0.09069235 0.03824633 0.04112368 0.02466697 0.01958732 

# Total inertia
sum(row.inertia)

[1] 1.11494

The total inertia corresponds to the amount of the information the data contains.

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:

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)

In the table above, the column TOTAL is called the average column profile (or marginale profile of rows)

The average column profile is calculated as follow:

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.10091743 0.08772936 0.06192661 0.08027523 0.06995413 0.06479358 

Distance (similarity) between column profiles

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

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.050311

You can also compute the squared distance between each column profile and the average column profile

Squared distance between each column profile and the average column 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.sum <- apply(housetasks, 2, sum)
grand.total <- sum(housetasks)
col.mass <- col.sum/grand.total
head(col.mass)

       Wife Alternating     Husband     Jointly 
  0.3440367   0.1456422   0.2184633   0.2918578 

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 * d2.col
head(col.inertia)

       Wife Alternating     Husband     Jointly 
  0.3010185   0.1178242   0.3813729   0.3147248 

# total inertia
sum(col.inertia)

[1] 1.11494

Recall that the total inertia corresponds to the amount of the information the data contains. Note that, the total inertia obtained using column profile is the same as the one obtained when analyzing row profile. That’s normal, because we are analyzing the same data with just a different angle of view.

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:

The Chi-square statistic is calculated as follow:

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.

Note that, Chi-square test should only be applied when the expected frequency of any cell is at least 5.

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 = 1944.456, df = 36, p-value < 2.2e-16

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

Note that, while Chi-square test can help to establish dependence between rows and the columns, the nature of the dependency is unknown.

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.55       25.63   38.45   51.37
Main_meal  52.64       22.28   33.42   44.65
Dinner     37.16       15.73   23.59   31.52
Breakfeast 48.17       20.39   30.58   40.86
Tidying    41.97       17.77   26.65   35.61
Dishes     38.88       16.46   24.69   32.98
Shopping   41.28       17.48   26.22   35.02
Official   33.03       13.98   20.97   28.02
Driving    47.82       20.24   30.37   40.57
Finances   38.88       16.46   24.69   32.98
Insurance  47.82       20.24   30.37   40.57
Repairs    56.77       24.03   36.05   48.16
Holidays   55.05       23.30   34.95   46.70

As mentioned above the Chi-square statistic is 1944.456196.

Which are the most contributing cells to the definition of the total Chi-square statistic?

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.266      -2.298  -5.878  -6.609
Main_meal   9.836      -0.484  -4.917  -6.084
Dinner      6.537      -1.192  -3.416  -3.299
Breakfeast  4.875       3.457  -2.818  -5.297
Tidying     1.702      -1.606  -4.969   3.585
Dishes     -1.103       1.859  -4.163   3.486
Shopping   -1.289       1.321  -3.362   3.376
Official   -3.659       8.563   0.443  -2.459
Driving    -5.469       6.836   8.100  -5.898
Finances   -4.150      -0.852  -0.742   5.750
Insurance  -5.758      -4.277   4.107   5.720
Repairs    -7.534      -4.290  20.646  -6.651
Holidays   -7.419      -4.620  -4.897  15.556

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.

Note that, correspondence analysis is just the singular value decomposition of the standardized residuals. This will be explained in the next section.

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

• 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

Chi-square statistic and the 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:

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

[1] 1.11494

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)

Likelihood ratio test in R

The functions likelihood.test()[in Deducer package] or G.test()[in 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() work as chisq.test():

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

    G-test

data:  as.matrix(housetasks)
G = 1907.658, df = 36, p-value < 2.2e-16

Interpret the association between rows and columns using likelihood ratio

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

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.5428893 48.69222    48.69222
2 0.4450028 39.91269    88.60491
3 0.1270484 11.39509   100.00000

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 :

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.83762154 -0.3652207 -0.19991139
Alternating -0.06218462 -0.2915938  0.84858939
Husband      1.16091847 -0.6019199 -0.18885924
Jointly      0.14942609  1.0265791 -0.04644302

# 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:

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 :

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.46 10.31 10.82
Alternating  0.10  2.78 82.55
Husband     54.23 17.79  6.13
Jointly      1.20 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.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:

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.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.1889641 0.4419662 0.2669493

## 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 = "/")

Coordinates of variables on the principal components

The correlation between variables and principal components is used as coordinates. It can be calculated as follow :

Variable correlations with PCs = loadings * the component standard deviations.

# Helper function : 
# Correlation between variables and principal components
var_cor_func <- function(var.loadings, comp.sdev){
  var.loadings*comp.sdev
  }

# Variable correlation/coordinates
loadings <- res.pca$rotation
sdev <- res.pca$sdev

var.coord <- var.cor <- t(apply(loadings, 1, var_cor_func, sdev))
head(var.coord[, 1:4])

                    PC1         PC2        PC3         PC4
X100m        -0.8506257  0.17939806 -0.3015564  0.03357320
Long.jump     0.7941806 -0.28085695  0.1905465 -0.11538956
Shot.put      0.7339127 -0.08540412 -0.5175978  0.12846837
High.jump     0.6100840  0.46521415 -0.3300852  0.14455012
X400m        -0.7016034 -0.29017826 -0.2835329  0.43082552
X110m.hurdle -0.7641252  0.02474081 -0.4488873 -0.01689589

Graph of variables using R base graph

# Plot the correlation circle
a <- seq(0, 2*pi, length = 100)
plot( cos(a), sin(a), type = &#39;l&#39;, col="gray",
      xlab = "PC1",  ylab = "PC2")

abline(h = 0, v = 0, lty = 2)

# Add active variables
arrows(0, 0, var.coord[, 1], var.coord[, 2], 
      length = 0.1, angle = 15, code = 2)

# Add labels
text(var.coord, labels=rownames(var.coord), cex = 1, adj=1)

Graph of variables using factoextra

fviz_pca_var(res.pca)

Read more about the function fviz_pca_var() : Graph of variables - Principal Component Analysis

How to interpret the correlation plot?

The graph of variables shows the relationships between all variables :

• Positively correlated variables are grouped together.
• Negatively correlated variables are positioned on opposite sides of the plot origin (opposed quadrants).
• The distance between variables and the origine measures the quality of the variables on the factor map. Variables that are away from the origin are well represented on the factor map.

Cos2 : quality of representation for variables on the factor map

The cos2 of variables are calculated as the squared coordinates : var.cos2 = var.coord * var.coord

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

                   PC1          PC2        PC3          PC4
X100m        0.7235641 0.0321836641 0.09093628 0.0011271597
Long.jump    0.6307229 0.0788806285 0.03630798 0.0133147506
Shot.put     0.5386279 0.0072938636 0.26790749 0.0165041211
High.jump    0.3722025 0.2164242070 0.10895622 0.0208947375
X400m        0.4922473 0.0842034209 0.08039091 0.1856106269
X110m.hurdle 0.5838873 0.0006121077 0.20149984 0.0002854712

Using factoextra package, the color of variables can be automatically controlled by the value of their cos2.

fviz_pca_var(res.pca, col.var="contrib")+
scale_color_gradient2(low="white", mid="blue", 
      high="red", midpoint=55) + theme_minimal()

Contributions of the variables to the principal components

The contribution of a variable to a given principal component is (in percentage) : (var.cos2 * 100) / (total cos2 of the component)

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.544293  1.7505098  7.338659  0.13755240
Long.jump    15.293168  4.2904162  2.930094  1.62485936
Shot.put     13.060137  0.3967224 21.620432  2.01407269
High.jump     9.024811 11.7715838  8.792888  2.54987951
X400m        11.935544  4.5799296  6.487636 22.65090599
X110m.hurdle 14.157544  0.0332933 16.261261  0.03483735

Highlight the most important (i.e, contributing) variables :

fviz_pca_var(res.pca, col.var="contrib") +
scale_color_gradient2(low="white", mid="blue", 
      high="red", midpoint=50) + theme_minimal()

You can also use the function fviz_contrib() described here : Principal Component Analysis: How to reveal the most important variables in your data?

Coordinates of individuals on the principal components

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

                 PC1        PC2        PC3         PC4
SEBRLE     0.1912074 -1.5541282 -0.6283688  0.08205241
CLAY       0.7901217 -2.4204156  1.3568870  1.26984296
BERNARD   -1.3292592 -1.6118687 -0.1961500 -1.92092203
YURKOV    -0.8694134  0.4328779 -2.4739822  0.69723814
ZSIVOCZKY -0.1057450  2.0233632  1.3049312 -0.09929630
McMULLEN   0.1185550  0.9916237  0.8435582  1.31215266

Cos2 : quality of representation for individuals on the principal components

To calculate the cos2 of individuals, 2 simple steps are required :

1. 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 + …+..

2. Calculate the cos2 = ind.coord^2/d2

# Compute the square of the distance between an individual and the
# 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)

# Compute the cos2
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.007530179 0.49747323 0.081325232 0.001386688
CLAY      0.048701249 0.45701660 0.143628117 0.125791741
BERNARD   0.197199804 0.28996555 0.004294015 0.411819183
YURKOV    0.096109800 0.02382571 0.778230322 0.061812637
ZSIVOCZKY 0.001574385 0.57641944 0.239754152 0.001388216
McMULLEN  0.002175437 0.15219499 0.110137872 0.266486530

The sum of each row is 1, if we consider the 10 components

Contribution of individuals to the princial components

The contribution of individuals (in percentage) to the principal components can be computed as follow :

100 * (1 / number_of_individuals)*(ind.coord^2 / comp_sdev^2)

# 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.03854254  5.7118249  1.3854184  0.03572215
CLAY      0.65814114 13.8541889  6.4600973  8.55568792
BERNARD   1.86273218  6.1441319  0.1349983 19.57827284
YURKOV    0.79686310  0.4431309 21.4755770  2.57939100
ZSIVOCZKY 0.01178829  9.6816398  5.9748485  0.05231437
McMULLEN  0.01481737  2.3253860  2.4967890  9.13531719

Note that the sum of all the contributions per column is 100

Graph of individuals : base graph

plot(ind.coord[,1], ind.coord[,2], pch = 19,  
     xlab="PC1 - 41.2%",ylab="PC2 - 18.4%")
abline(h=0, v=0, lty = 2)
text(ind.coord[,1], ind.coord[,2], labels=rownames(ind.coord),
        cex=0.7, pos = 3)

Biplot of individuals and variables :

biplot(res.pca, cex = 0.8, col = c("black", "red") )

Graph of individuals : factoextra

ind <- get_pca_ind(res.pca)
ind

Principal Component 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 for individuals
head(ind$coord[, 1:4])

               Dim.1      Dim.2      Dim.3       Dim.4
SEBRLE     0.1912074 -1.5541282 -0.6283688  0.08205241
CLAY       0.7901217 -2.4204156  1.3568870  1.26984296
BERNARD   -1.3292592 -1.6118687 -0.1961500 -1.92092203
YURKOV    -0.8694134  0.4328779 -2.4739822  0.69723814
ZSIVOCZKY -0.1057450  2.0233632  1.3049312 -0.09929630
McMULLEN   0.1185550  0.9916237  0.8435582  1.31215266

fviz_pca_ind(res.pca)

fviz_pca_ind(res.pca, col.ind="cos2") +
scale_color_gradient2(low="white", mid="blue", 
    high="red", midpoint=0.50) + theme_minimal()

fviz_pca_biplot(res.pca,  geom = "text") +
  theme_minimal()

Supplementary quantitative variables

As described above, the data sets decathlon2 contain some supplementary continuous variables at columns 11 and 12 corresponding respectively to the rank and the points of athletes.

# Data for the supplementary quantitative variables
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

Recall that, rows 24:27 are supplementary individuals. We don’t want them in this current analysis. This is why, I extracted only rows 1:23.

In this section we’ll see how to calculate the predicted coordinates of these two variables using the information provided by the previously performed principal component analysis.

2 simples steps are required :

1. Calculate the correlation between each supplementary quantitative variables and the principal components
2. Make a factor map of all variables (active and supplementary ones) to visualize the position of the supplementary variables

The R code below can be used :

# Calculate the correlations between supplementary variables
# and the principal components
ind.coord <- res.pca$x
quanti.coord <- cor(quanti.sup, ind.coord)
head(quanti.coord[, 1:4])

              PC1         PC2        PC3         PC4
Rank   -0.7014777  0.24519443  0.1834294  0.05575186
Points  0.9637075 -0.07768262 -0.1580225 -0.16623092

# Variable factor maps
#++++++++++++++++++
# Plot the correlation circle
a <- seq(0, 2*pi, length = 100)
plot( cos(a), sin(a), type = &#39;l&#39;, col="gray",
      xlab = "PC1",  ylab = "PC2")
abline(h = 0, v = 0, lty = 2)
# Add active variables
var.coord <- get_pca_var(res.pca)$coord
arrows(0 ,0, x1=var.coord[,1], y1 = var.coord[,2], 
       col="black", length = 0.09)
text(var.coord[,1], var.coord[,2],
     labels=rownames(var.coord), cex=0.8)
# Add supplementary quantitative variables
arrows(0 ,0, x1= quanti.coord[,1], y1 = quanti.coord[,2], 
       col="blue", lty =2, length = 0.09)
text(quanti.coord[,1], quanti.coord[,2],
     labels=rownames(quanti.coord), cex=0.8, col ="blue")

It’s also possible to make the graph of variables using factoextra:

# Plot of active variables
p <- fviz_pca_var(res.pca)
# Add supplementary active variables
fviz_add(p, quanti.coord, color ="blue", geom="arrow")

# get the cos2 of the supplementary quantitative variables
(quanti.coord^2)[, 1:4]

             PC1         PC2        PC3        PC4
Rank   0.4920710 0.060120310 0.03364635 0.00310827
Points 0.9287322 0.006034589 0.02497110 0.02763272

Supplementary qualitative variables

The data sets decathlon2 contain a supplementary qualitative variable at columns 13 corresponding to the type of competitions.

Qualitative variable can be helpful for interpreting the data and for coloring individuals by groups :

# Data for the supplementary qualitative variables
quali.sup <- as.factor(decathlon2[1:23, 13])
head(quali.sup)

[1] Decastar Decastar Decastar Decastar Decastar Decastar
Levels: Decastar OlympicG

Color individuals by groups :

fviz_pca_ind(res.pca, 
  habillage = quali.sup, addEllipses = TRUE, ellipse.level = 0.68) +
  theme_minimal()

Note that, the argument habillage is used to specify the variable containing the groups of individuals

It’s very easy to get the coordinates for the levels of a supplementary qualitative variables. The helper function below can be used :

# Return the coordinates of a group levels
# x : coordinate of individuals on x axis
# y : coordinate of indiviuals on y axis
get_coord_quali<-function(x, y, groups){
  data.frame(
    x= tapply(x, groups, mean),
    y = tapply(y, groups, mean)
  )
}

Calculate the coordinates on components 1 and 2 :

coord.quali <- get_coord_quali(ind.coord[,1], ind.coord[,2],
                               groups = quali.sup)
coord.quali

                 x          y
Decastar -1.313921 -0.1191322
OlympicG  1.204428  0.1092046

Supplementary individuals

The data sets decathlon2 contain some supplementary individuals from row 24 to 27.

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

        X100m Long.jump Shot.put High.jump X400m X110m.hurdle
KARPOV  11.02      7.30    14.77      2.04 48.37        14.09
WARNERS 11.11      7.60    14.31      1.98 48.68        14.23
Nool    10.80      7.53    14.26      1.88 48.81        14.80
Drews   10.87      7.38    13.07      1.88 48.51        14.01

Remember that, columns 11:13 are supplementary variables. We don’t want them in this current analysis. This is why, I extracted only columns 1:10. I used also the argument drop = FALSE to preserve the type of the data (which is a data.frame).

In this section we’ll see how to predict the coordinates of the supplementary individuals using only the information provided by the previously performed principal component analysis.

A simple function to predict the coordinates of new individuals data

One simple approach is to use the function predict() from the built-in R stats package :

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

               PC1         PC2       PC3        PC4
KARPOV   0.7772521 -0.76237804 1.5971253  1.6863286
WARNERS -0.3779697  0.11891968 1.7005146 -0.6908084
Nool    -0.5468405 -1.93402211 0.4724184 -2.2283706
Drews   -1.0848227 -0.01703198 2.9818031 -1.5006207

Calculate the predicted coordinates by hand

2 simples steps are required :

1. Center and scale the values for the supplementary individuals 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
scale_func <- function(ind_row, center, scale){
  (ind_row-center)/scale
}

ind.scaled <- t(apply(ind.sup, 1, scale_func, res.pca$center, res.pca$scale))

# Coordinates of the individividuals
pca.loadings <- res.pca$rotation
coord_func <- function(ind, loadings){
  r <- loadings*ind
  r <- apply(r, 2, sum)
  r
}

ind.sup.coord <- t(apply(ind.scaled, 1, coord_func, pca.loadings ))
ind.sup.coord[, 1:4]

               PC1         PC2       PC3        PC4
KARPOV   0.7772521 -0.76237804 1.5971253  1.6863286
WARNERS -0.3779697  0.11891968 1.7005146 -0.6908084
Nool    -0.5468405 -1.93402211 0.4724184 -2.2283706
Drews   -1.0848227 -0.01703198 2.9818031 -1.5006207

Make a factor map including the supplementary individuals using factoextra

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

Extract the eigenvalues

Eigenvalues measure the amount of variation retained by a principal component :

summary(res.pca)

Class: pca dudi
Call: dudi.pca(df = decathlon2.active, scannf = FALSE, nf = 5)

Total inertia: 10

Eigenvalues:
    Ax1     Ax2     Ax3     Ax4     Ax5 
 4.1242  1.8385  1.2391  0.8194  0.7016 

Projected inertia (%):
    Ax1     Ax2     Ax3     Ax4     Ax5 
 41.242  18.385  12.391   8.194   7.016 

Cumulative projected inertia (%):
    Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
  41.24   59.63   72.02   80.21   87.23 

(Only 5 dimensions (out of 10) are shown)

You can also use the package factoextra to extract the eigenvalues :

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

      eigenvalue variance.percent cumulative.variance.percent
Dim 1  4.1242133        41.242133                    41.24213
Dim 2  1.8385309        18.385309                    59.62744
Dim 3  1.2391403        12.391403                    72.01885
Dim 4  0.8194402         8.194402                    80.21325
Dim 5  0.7015528         7.015528                    87.22878
Dim 6  0.4228828         4.228828                    91.45760

Make a scree plot using ade4 base graphics

The function scree plot() can be used to represent the amount of inertia (variance) associated with each principal component (PC).

A simplified format is :

screeplot(x, ncps = length(x$eig), type = c("barplot", "lines"))

Example of usage :

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

You can also customize the plot using the standard barplot() function. In the R code below, we’ll draw the percentage of variances retained by each component :

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

~60% of the information (variances) contained in the data are retained by the first two principal components.

Make the scree plot using the package factoextra

fviz_screeplot(res.pca, ncp=10)

Coordinates of variables on the principal components

The coordinates of the variables on the factor map are :

# Column coordinates
head(res.pca$co)

                  Comp1       Comp2      Comp3       Comp4      Comp5
X100m         0.8506257 -0.17939806 -0.3015564  0.03357320  0.1944440
Long.jump    -0.7941806  0.28085695  0.1905465 -0.11538956 -0.2331567
Shot.put     -0.7339127  0.08540412 -0.5175978  0.12846837  0.2488129
High.jump    -0.6100840 -0.46521415 -0.3300852  0.14455012 -0.4027002
X400m         0.7016034  0.29017826 -0.2835329  0.43082552 -0.1039085
X110m.hurdle  0.7641252 -0.02474081 -0.4488873 -0.01689589 -0.2242200

Graph of variables using ade4 base graph

The function s.corcircle() can be used to plot the correlation circle. A simplified format is :

s.corcircle(dfxy, label = row.names(dfxy), grid = TRUE,
            box = FALSE)

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

Graph of variables using factoextra

The function fviz_pca_var() is used to visualize variables :

# Default plot
fviz_pca_var(res.pca)

# Change color and theme
fviz_pca_var(res.pca, col.var="steelblue")+
  theme_minimal()

Read more about the function fviz_pca_var() : Graph of variables - Principal Component Analysis

How to calculate the cos2 and the contribution of variables?

The cos2 and the contributions of variables (columns) / individuals (rows) are calculated using the function inertia.dudi() as follow :

inertia <- inertia.dudi(res.pca, row.inertia = TRUE,
                        col.inertia = TRUE)

Note that, the contributions and the cos2 are printed in 1/10 000. The sign is the sign of the coordinates.

Cos2 : quality of the representation for variables on the factor map

The squared coordinates of variables are called cos2.

• A high cos2 indicates a good representation of the variable on the principal component. In this case the variable is positioned close to the circumference of the correlation circle.

• A low cos2 indicates that the variable is not perfectly represented by the PCs. In this case the variable is close to the center of the circle.

The cos2 of the variables are :

# relative contributions of columns
var.cos2 <- abs(inertia$col.rel/10000)
head(var.cos2)

              Comp1  Comp2  Comp3  Comp4  Comp5 con.tra
X100m        0.7236 0.0322 0.0909 0.0011 0.0378     0.1
Long.jump    0.6307 0.0789 0.0363 0.0133 0.0544     0.1
Shot.put     0.5386 0.0073 0.2679 0.0165 0.0619     0.1
High.jump    0.3722 0.2164 0.1090 0.0209 0.1622     0.1
X400m        0.4922 0.0842 0.0804 0.1856 0.0108     0.1
X110m.hurdle 0.5839 0.0006 0.2015 0.0003 0.0503     0.1

It can also be calculated as follow :

# squared coordinates
head(res.pca$co^2)

                 Comp1        Comp2      Comp3        Comp4      Comp5
X100m        0.7235641 0.0321836641 0.09093628 0.0011271597 0.03780845
Long.jump    0.6307229 0.0788806285 0.03630798 0.0133147506 0.05436203
Shot.put     0.5386279 0.0072938636 0.26790749 0.0165041211 0.06190783
High.jump    0.3722025 0.2164242070 0.10895622 0.0208947375 0.16216747
X400m        0.4922473 0.0842034209 0.08039091 0.1856106269 0.01079698
X110m.hurdle 0.5838873 0.0006121077 0.20149984 0.0002854712 0.05027463

Using factoextra package, the color of variables can be automatically controlled by the value of their cos2.

fviz_pca_var(res.pca, col.var="contrib")+
scale_color_gradient2(low="white", mid="blue", 
      high="red", midpoint=55) + theme_minimal()

Contributions of the variables to the principal components

The contributions can be printed in % as follow :

# absolute contribution of columns
var.contrib <- inertia$col.abs/100
head(var.contrib)

             Comp1 Comp2 Comp3 Comp4 Comp5
X100m        17.54  1.75  7.34  0.14  5.39
Long.jump    15.29  4.29  2.93  1.62  7.75
Shot.put     13.06  0.40 21.62  2.01  8.82
High.jump     9.02 11.77  8.79  2.55 23.12
X400m        11.94  4.58  6.49 22.65  1.54
X110m.hurdle 14.16  0.03 16.26  0.03  7.17

Note that, You can also use the function get_pca_var() [from factoextra package]. It provides a list of matrices containing all the results for the active variables (coordinates, correlation between variables and axes, squared cosine and contributions).

var <- get_pca_var(res.pca)
names(var)

[1] "coord"   "cor"     "cos2"    "contrib"

# Contributions of variables
head(var$contrib)

                 Dim.1      Dim.2     Dim.3       Dim.4     Dim.5
X100m        17.544293  1.7505098  7.338659  0.13755240  5.389252
Long.jump    15.293168  4.2904162  2.930094  1.62485936  7.748815
Shot.put     13.060137  0.3967224 21.620432  2.01407269  8.824401
High.jump     9.024811 11.7715838  8.792888  2.54987951 23.115504
X400m        11.935544  4.5799296  6.487636 22.65090599  1.539012
X110m.hurdle 14.157544  0.0332933 16.261261  0.03483735  7.166193

Using factoextra package, the color of variables can be automatically controlled by the value of their contributions

fviz_pca_var(res.pca, col.var="contrib") +
scale_color_gradient2(low="white", mid="blue", 
      high="red", midpoint=50) + theme_minimal()

This is helpful to highlight the most important variables for the principal components.

The most important variables for a given PC can be visualized using the function fviz_pca_contrib()[factoextra package] :

(factoextra >= 1.0.1 is required)

# Contributions of variables on PC1
fviz_pca_contrib(res.pca, choice = "var", axes = 1)

# Contributions of variables on PC2
fviz_pca_contrib(res.pca, choice = "var", axes = 2)

Read more about fviz_pca_contrib() : Principal Component Analysis: How to reveal the most important variables in your data?

Coordinates of individuals on the principal components

The coordinates of the individuals on the factor maps can be extracted as follow :

# The row coordinates
head(res.pca$li)

                Axis1      Axis2      Axis3       Axis4       Axis5
SEBRLE     -0.1955047  1.5890567 -0.6424912  0.08389652 -1.16829387
CLAY       -0.8078795  2.4748137  1.3873827  1.29838232  0.82498206
BERNARD     1.3591340  1.6480950 -0.2005584 -1.96409420 -0.08419345
YURKOV      0.8889532 -0.4426067 -2.5295843  0.71290837 -0.40782264
ZSIVOCZKY   0.1081216 -2.0688377  1.3342591 -0.10152796  0.20145217
McMULLEN   -0.1212195 -1.0139102  0.8625170  1.34164291 -1.62151286

Cos2 : quality of the representation for individuals on the principal components

# relative contributions of rows
ind.cos2 <- abs(inertia$row.rel)/10000
head(ind.cos2)

            Axis1  Axis2  Axis3  Axis4  Axis5 con.tra
SEBRLE     0.0075 0.4975 0.0813 0.0014 0.2689  0.0221
CLAY       0.0487 0.4570 0.1436 0.1258 0.0508  0.0583
BERNARD    0.1972 0.2900 0.0043 0.4118 0.0008  0.0407
YURKOV     0.0961 0.0238 0.7782 0.0618 0.0202  0.0357
ZSIVOCZKY  0.0016 0.5764 0.2398 0.0014 0.0055  0.0323
McMULLEN   0.0022 0.1522 0.1101 0.2665 0.3893  0.0294

Contribution of the individuals to the princial components

The contributions can be printed in % as follow :

# absolute contributions of rows
ind.contrib <- inertia$row.abs/100
head(ind.contrib)

           Axis1 Axis2 Axis3 Axis4 Axis5
SEBRLE      0.04  5.97  1.45  0.04  8.46
CLAY        0.69 14.48  6.75  8.94  4.22
BERNARD     1.95  6.42  0.14 20.47  0.04
YURKOV      0.83  0.46 22.45  2.70  1.03
ZSIVOCZKY   0.01 10.12  6.25  0.05  0.25
McMULLEN    0.02  2.43  2.61  9.55 16.29

It’s also possible to use the function get_pca_ind() [from factoextra package]. factoextra provides, a list of matrices containing all the results for the active individuals (coordinates, squared cosine and contributions)./span>

ind <- get_pca_ind(res.pca)
names(ind)

[1] "coord"   "cos2"    "contrib"

# Contributions of individuals
head(ind$contrib)

                Dim.1      Dim.2      Dim.3       Dim.4       Dim.5
SEBRLE     0.04029447  5.9714533  1.4483919  0.03734589  8.45894063
CLAY       0.68805664 14.4839248  6.7537381  8.94458283  4.21794385
BERNARD    1.94740183  6.4234107  0.1411345 20.46819433  0.04393073
YURKOV     0.83308415  0.4632733 22.4517396  2.69663605  1.03075263
ZSIVOCZKY  0.01232413 10.1217143  6.2464325  0.05469230  0.25151025
McMULLEN   0.01549089  2.4310854  2.6102794  9.55055888 16.29493304

Use the function fviz_pca_contrib()[factoextra package] to visualize the most contributing individuals :

(factoextra >= 1.0.1 is required)

# Contributions of variables on PC1
fviz_pca_contrib(res.pca, choice = "ind", axes = 1)

# Contributions of variables on PC2
fviz_pca_contrib(res.pca, choice = "ind", axes = 2)

Read more about fviz_pca_contrib() : Principal Component Analysis: How to reveal the most important variables in your data?

Graph of individuals using ade4 base graph

The function s.label() can be used. A simplified format is :

s.label(dfxy, xax = 1, yax = 2)

Factor map of individuals :

s.label(res.pca$li, xax = 1, yax = 2)

Biplot of individuals and variables using ade4

Biplot can be drawn using the combination of the two functions below :

• s.label() to plot individuals
• s.arrow() to add variables

# Plot of individuals
s.label(res.pca$li, xax = 1, yax = 2)
# Add variables
s.arrow(7*res.pca$c1, add.plot = TRUE)

It’s also possible to use the function scatter() or biplot() :

scatter(res.pca)

# Remove the scree plot (posieig ="none")
# Remove row labels (clab.row = 0)
scatter(res.pca,  posieig = "none", clab.row = 0)

NULL

Note that, to remove variable labels the argument clab.col = 0 can be used.

Graph of individuals using factoextra

The function fviz_pca_ind() is used to visualize individuals :

fviz_pca_ind(res.pca)

Control automatically the color of individuals using the cos2 values (the quality of the individuals on the factor map) :

fviz_pca_ind(res.pca, col.ind="cos2") +
scale_color_gradient2(low="white", mid="blue", 
                  high="red", midpoint=0.50)

Change the theme :

fviz_pca_ind(res.pca, col.ind="cos2") +
scale_color_gradient2(low="white", mid="blue", 
    high="red", midpoint=0.50) + theme_minimal()

Read more about fviz_pca_ind() : Graph of individuals - principal component analysis

Make a biplot of individuals and variables :

fviz_pca_biplot(res.pca, geom = "text") +
  theme_minimal()

Read more about fviz_pca_biplot() : Biplot of individuals and variables - principal component analysis

Change the color of individuals by groups

The data sets decathlon2 contain a supplementary qualitative variable at columns 13 corresponding to the type of competitions.

Qualitative variable can be helpful for interpreting the data and for coloring individuals by groups :

# Data for the supplementary qualitative variables
quali.sup <- as.factor(decathlon2[1:23, 13])
head(quali.sup)

[1] Decastar Decastar Decastar Decastar Decastar Decastar
Levels: Decastar OlympicG

The function s.class() can be used to visualize the classes (groups) of points :

s.class(dfxy, fac, xax = 1, yax = 2, col)

Color individuals by groups :

s.class(res.pca$li, fac = quali.sup, xax = 1, yax = 2)

# Change the colors
s.class(res.pca$li, fac = quali.sup, col = c("blue", "red"))

# Make a biplot
# clab.row : hide the label for rows (individuals)
res <- scatter(res.pca, clab.row = 0, posieig = "none")
s.class(res.pca$li, fac = quali.sup, col = c("blue", "red"),
        add.plot = TRUE)

# Customize the biplot
# - remove row labels (clab.row = 0)
# - hide the scree plot (posieig = 0)
# - remove stars (cstar = 0)
# - remove ellipse (cellipse = 0)
res <- scatter(res.pca, clab.row = 0, posieig = "none")
s.class(res.pca$li, fac = quali.sup, col = c("blue", "red"),
        add.plot = TRUE, cstar = 0, cellipse = 0)

# remove labels for classes (clabel = 0)
res <- scatter(res.pca, clab.row = 0, posieig = "none")
s.class(res.pca$li, fac = quali.sup, col = c("blue", "red"),
        add.plot = TRUE, cstar = 0, cellipse = 0, clabel = 0)

It’s also possible to use factoextra :

fviz_pca_ind(res.pca, habillage = quali.sup,
     addEllipses =TRUE, ellipse.level = 0.68) +
  theme_minimal()

Elegant biplot using factoextra and iris data :

data(iris)

head(iris)

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

# The variable Species (index = 5) is removed
# before PCA analysis
iris.pca <- dudi.pca(iris[,-5], scannf = FALSE, nf = 2)

Now, let’s :

• make a biplot of individuals and variables
• change the color of individuals by groups
• change the transparency of variable colors by their contribution values
• show only the labels for variables

fviz_pca_biplot(iris.pca, 
  habillage = iris$Species, addEllipses = TRUE,
  col.var = "red", alpha.var ="cos2",
  label = "var") +
  scale_color_brewer(palette="Dark2")+
  theme_minimal()

Supplementary individuals

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

         X100m Long.jump Shot.put High.jump X400m X110m.hurdle
KARPOV   11.02      7.30    14.77      2.04 48.37        14.09
WARNERS  11.11      7.60    14.31      1.98 48.68        14.23
Nool     10.80      7.53    14.26      1.88 48.81        14.80
Drews    10.87      7.38    13.07      1.88 48.51        14.01

Predict the coordinates of the supplementary individuals :

ind.sup.pca <- suprow(res.pca, ind.sup)
names(ind.sup.pca)

[1] "tabsup" "lisup" 

# coordinates 
ind.sup.coord <- ind.sup.pca$lisup
head(ind.sup.coord)

              Axis1       Axis2     Axis3      Axis4      Axis5
KARPOV   -0.7947206  0.77951227 1.6330203  1.7242283 0.75070396
WARNERS   0.3864645 -0.12159237 1.7387332 -0.7063341 0.03230011
Nool      0.5591306  1.97748871 0.4830358 -2.2784526 0.25461493
Drews     1.1092038  0.01741477 3.0488182 -1.5343468 0.32642192

How to visualize supplementary individuals on the factor map?

The function fviz_add() is used :

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

How to calculate the cos2 (quality of the representation) for supplementary individuals?

cos2.func <-function(x){x^2/sum(x^2)}
ind.sup.cos2 <- t(apply(ind.sup.coord, 1, cos2.func))
head(ind.sup.cos2)

              Axis1        Axis2      Axis3     Axis4        Axis5
KARPOV   0.08486144 8.164458e-02 0.35831467 0.3994579 0.0757214366
WARNERS  0.04050537 4.009646e-03 0.81989704 0.1353050 0.0002829447
Nool     0.03218782 4.026179e-01 0.02402281 0.5344967 0.0066747159
Drews    0.09473792 2.335268e-05 0.71575477 0.1812793 0.0082046453

Supplementary quantitative variables

# Data for the supplementary quantitative variables
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

Remember that, rows 24:27 are supplementary individuals. We don’t want them in this current analysis. This is why, I extracted only rows 1:23.

Predict the coordinates of the supplementary variables :

(You have to scale the supplementary variables before the analysis as the PCA has been performed on scaled data.)

quanti.pca <- supcol(res.pca, scale(quanti.sup))
names(quanti.pca)

[1] "tabsup" "cosup" 

# coordinates 
quanti.coord <- quanti.pca$cosup
head(quanti.coord)

            Comp1      Comp2      Comp3      Comp4      Comp5
Rank    0.6860587 -0.2398049  0.1793975  0.0545264 0.07220371
Points -0.9425246  0.0759751 -0.1545490 -0.1625770 0.03046248

Visualize supplementary variables on the factor map using factoextra :

# Plot of active variables
p <- fviz_pca_var(res.pca)
# Add supplementary active variables
fviz_add(p, quanti.coord, geom="arrow", color ="blue")

# Get the cos2 of the supplementary quantitative variables
(quanti.coord^2)[, 1:4]

           Comp1       Comp2      Comp3       Comp4
Rank   0.4706766 0.057506383 0.03218347 0.002973128
Points 0.8883526 0.005772216 0.02388540 0.026431296

Descriptive statistics

decathlon2.active_stats <- data.frame(
  Min = apply(decathlon2.active, 2, min), # minimum
  Q1 = apply(decathlon2.active, 2, quantile, 1/4), # First quartile
  Med = apply(decathlon2.active, 2, median), # median
  Mean = apply(decathlon2.active, 2, mean), # mean
  Q3 = apply(decathlon2.active, 2, quantile, 3/4), # Third quartile
  Max = apply(decathlon2.active, 2, max) # Maximum
  )
decathlon2.active_stats <- round(decathlon2.active_stats, 1)
head(decathlon2.active_stats)

              Min   Q1  Med Mean   Q3  Max
X100m        10.4 10.8 11.0 11.0 11.2 11.6
Long.jump     6.8  7.2  7.3  7.3  7.5  8.0
Shot.put     12.7 14.2 14.7 14.6 15.1 16.4
High.jump     1.9  1.9  2.0  2.0  2.1  2.1
X400m        46.8 49.0 49.4 49.4 50.0 51.2
X110m.hurdle 14.0 14.2 14.4 14.5 14.9 15.7

Note that, you can also use the built-in R function summary() for the descriptive statistics but I don’t like the format of the output on data frame.

Correlation matrix

The correlation between variables can be calculated as follow :

cor.mat <- round(cor(decathlon2.active),2)
head(cor.mat[, 1:6])

             X100m Long.jump Shot.put High.jump X400m X110m.hurdle
X100m         1.00     -0.76    -0.45     -0.40  0.59         0.73
Long.jump    -0.76      1.00     0.44      0.34 -0.51        -0.59
Shot.put     -0.45      0.44     1.00      0.53 -0.31        -0.38
High.jump    -0.40      0.34     0.53      1.00 -0.37        -0.25
X400m         0.59     -0.51    -0.31     -0.37  1.00         0.58
X110m.hurdle  0.73     -0.59    -0.38     -0.25  0.58         1.00

Visualize the correlation matrix using a correlogram : the package corrplot is required.

# install.packages("corrplot")
library("corrplot")
corrplot(cor.mat, type="upper", order="hclust", 
         tl.col="black", tl.srt=45)

Read more about visualizing correlation matrix : Correlation matrix visualization

Make a scatter plot matrix showing the correlation coefficients between variables and the significance levels : the package PerformanceAnalytics is required.

# install.packages("PerformanceAnalytics")
library("PerformanceAnalytics")
chart.Correlation(decathlon2.active[, 1:6], histogram=TRUE, pch=19)

You can read more about this plot here : Correlation matrix visualization

Variances of the principal components

The proportion of variances retained by the principal components can be extracted as follow :

eigenvalues <- res.pca$eig
head(eigenvalues[, 1:2])

       eigenvalue percentage of variance
comp 1  4.1242133              41.242133
comp 2  1.8385309              18.385309
comp 3  1.2391403              12.391403
comp 4  0.8194402               8.194402
comp 5  0.7015528               7.015528
comp 6  0.4228828               4.228828

Make a scree plot using base graphics : A scree plot is a graph of the eigenvalues/variances associated with components.

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

~60% of the informations (variances) contained in the data are retained by the first two principal components.

Make the scree plot using the package factoextra :

fviz_screeplot(res.pca, ncp=10)

Coordinates of variables on the principal components

head(res.pca$var$coord)

                  Dim.1       Dim.2      Dim.3       Dim.4      Dim.5
X100m        -0.8506257 -0.17939806  0.3015564  0.03357320 -0.1944440
Long.jump     0.7941806  0.28085695 -0.1905465 -0.11538956  0.2331567
Shot.put      0.7339127  0.08540412  0.5175978  0.12846837 -0.2488129
High.jump     0.6100840 -0.46521415  0.3300852  0.14455012  0.4027002
X400m        -0.7016034  0.29017826  0.2835329  0.43082552  0.1039085
X110m.hurdle -0.7641252 -0.02474081  0.4488873 -0.01689589  0.2242200

Cos2 : quality of variables on the factor map

The quality of representation of the variables of the principal components are called the cos2.

head(res.pca$var$cos2)

                 Dim.1        Dim.2      Dim.3        Dim.4      Dim.5
X100m        0.7235641 0.0321836641 0.09093628 0.0011271597 0.03780845
Long.jump    0.6307229 0.0788806285 0.03630798 0.0133147506 0.05436203
Shot.put     0.5386279 0.0072938636 0.26790749 0.0165041211 0.06190783
High.jump    0.3722025 0.2164242070 0.10895622 0.0208947375 0.16216747
X400m        0.4922473 0.0842034209 0.08039091 0.1856106269 0.01079698
X110m.hurdle 0.5838873 0.0006121077 0.20149984 0.0002854712 0.05027463

Contributions of the variables to the principal components

Variable contributions in the determination of a given principal component are (in percentage) : (var.cos2 * 100) / (total cos2 of the component)

head(res.pca$var$contrib)

                 Dim.1      Dim.2     Dim.3       Dim.4     Dim.5
X100m        17.544293  1.7505098  7.338659  0.13755240  5.389252
Long.jump    15.293168  4.2904162  2.930094  1.62485936  7.748815
Shot.put     13.060137  0.3967224 21.620432  2.01407269  8.824401
High.jump     9.024811 11.7715838  8.792888  2.54987951 23.115504
X400m        11.935544  4.5799296  6.487636 22.65090599  1.539012
X110m.hurdle 14.157544  0.0332933 16.261261  0.03483735  7.166193

Graph of variables using FactoMineR base graph

plot(res.pca, choix = "var")

Graph of variables using factoextra

The function fviz_pca_var() is used to visualize variables :

# Default plot
fviz_pca_var(res.pca)

# Change color and theme
fviz_pca_var(res.pca, col.var="steelblue")+
  theme_minimal()

Note that, using factoextra package, the color or the transparency of variables can be automatically controlled by the value of their contributions, their cos2, their coordinates on x or y axis.

# Control variable colors using their contribution
# Possible values for the argument col.var are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_pca_var(res.pca, col.var="contrib")

# Change the gradient color
fviz_pca_var(res.pca, col.var="contrib")+
scale_color_gradient2(low="white", mid="blue", 
                      high="red", midpoint=55)+theme_bw()

This is helpful to highlight the most important variables in the determination of the principal components.

It’s also possible to control automatically the transparency of variables by their contributions :

# Control the transparency of variables using their contribution
# Possible values for the argument alpha.var are :
  # "cos2", "contrib", "coord", "x", "y"
fviz_pca_var(res.pca, alpha.var="contrib")+
  theme_minimal()

Read more about ggplot2 and colors here : ggplot2 colors - How to change colors automatically and manually?

Coordinates of individuals on the principal components

head(res.pca$ind$coord)

                Dim.1      Dim.2      Dim.3       Dim.4       Dim.5
SEBRLE      0.1955047  1.5890567  0.6424912  0.08389652  1.16829387
CLAY        0.8078795  2.4748137 -1.3873827  1.29838232 -0.82498206
BERNARD    -1.3591340  1.6480950  0.2005584 -1.96409420  0.08419345
YURKOV     -0.8889532 -0.4426067  2.5295843  0.71290837  0.40782264
ZSIVOCZKY  -0.1081216 -2.0688377 -1.3342591 -0.10152796 -0.20145217
McMULLEN    0.1212195 -1.0139102 -0.8625170  1.34164291  1.62151286

Cos2 : quality of representation of individuals on the principal components

head(res.pca$ind$cos2)

                 Dim.1      Dim.2       Dim.3       Dim.4        Dim.5
SEBRLE     0.007530179 0.49747323 0.081325232 0.001386688 0.2689026575
CLAY       0.048701249 0.45701660 0.143628117 0.125791741 0.0507850580
BERNARD    0.197199804 0.28996555 0.004294015 0.411819183 0.0007567259
YURKOV     0.096109800 0.02382571 0.778230322 0.061812637 0.0202279796
ZSIVOCZKY  0.001574385 0.57641944 0.239754152 0.001388216 0.0054654972
McMULLEN   0.002175437 0.15219499 0.110137872 0.266486530 0.3892621478

Contribition of individuals to the princial components

head(res.pca$ind$contrib)

                Dim.1      Dim.2      Dim.3       Dim.4       Dim.5
SEBRLE     0.04029447  5.9714533  1.4483919  0.03734589  8.45894063
CLAY       0.68805664 14.4839248  6.7537381  8.94458283  4.21794385
BERNARD    1.94740183  6.4234107  0.1411345 20.46819433  0.04393073
YURKOV     0.83308415  0.4632733 22.4517396  2.69663605  1.03075263
ZSIVOCZKY  0.01232413 10.1217143  6.2464325  0.05469230  0.25151025
McMULLEN   0.01549089  2.4310854  2.6102794  9.55055888 16.29493304

Graph of individuals using FactoMineR base graph

plot(res.pca, choix = "ind")

Graph of individuals using factoextra

The function fviz_pca_ind() is used to visualize individuals :

fviz_pca_ind(res.pca)

Remove the points from the graph, use texts only :

fviz_pca_ind(res.pca, geom="text")

Control automatically the color of individuals using the cos2 values (the quality of the individuals on the factor map) :

fviz_pca_ind(res.pca, col.ind="cos2") +
scale_color_gradient2(low="white", mid="blue", 
                      high="red", midpoint=0.50)

Read more about ggplot2 and colors here : ggplot2 colors - How to change colors automatically and manually?

Change the theme :

fviz_pca_ind(res.pca,  col.ind="cos2") +
scale_color_gradient2(low="white", mid="blue", 
                      high="red", midpoint=0.50)+
  theme_minimal()

Read more about ggplot2 themes here : ggplot2 themes and background colors

Make a biplot of individuals and variables :

fviz_pca_biplot(res.pca,  geom = "text")

Change the color of individuals by groups

We will use iris data sets in this section :

data(iris)

head(iris)

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa

# The variable Species (index = 5) is removed
# before PCA analysis
iris.pca <- PCA(iris[,-5], graph = FALSE)

Individuals factor map :

# Default plot
fviz_pca_ind(iris.pca, label="none")

Change individual colors by groups :

fviz_pca_ind(iris.pca,  label="none", habillage=iris$Species)

Add ellipses of point concentrations : the argument habillage is used to specify the factor variable for coloring the observations by groups.

fviz_pca_ind(iris.pca, label="none", habillage=iris$Species,
             addEllipses=TRUE, ellipse.level=0.95)

Now, let’s :

• make a biplot of individuals and variables
• change the color of individuals by groups
• change the transparency of variable colors by their contribution values
• show only the labels for variables

fviz_pca_biplot(iris.pca, 
  habillage = iris$Species, addEllipses = TRUE,
  col.var = "red", alpha.var ="cos2",
  label = "var") +
  scale_color_brewer(palette="Dark2")+
  theme_minimal()

Visualize supplementary quantitative variables

All the results (coordinates, correlation and cos2) for the supplementary quantitative variables can be extracted as follow :

res.pca$quanti.sup

$coord
            Dim.1       Dim.2      Dim.3       Dim.4       Dim.5
Rank   -0.7014777 -0.24519443 -0.1834294  0.05575186 -0.07382647
Points  0.9637075  0.07768262  0.1580225 -0.16623092 -0.03114711

$cor
            Dim.1       Dim.2      Dim.3       Dim.4       Dim.5
Rank   -0.7014777 -0.24519443 -0.1834294  0.05575186 -0.07382647
Points  0.9637075  0.07768262  0.1580225 -0.16623092 -0.03114711

$cos2
           Dim.1       Dim.2      Dim.3      Dim.4        Dim.5
Rank   0.4920710 0.060120310 0.03364635 0.00310827 0.0054503477
Points 0.9287322 0.006034589 0.02497110 0.02763272 0.0009701427

Variables factor map using FactoMineR base graph :

plot(res.pca, choix = "var")

Supplementary quantitative variables are shown in blue color and dashed lines.

It’s also possible to make the variables factor map using factoextra :

fviz_pca_var(res.pca)

Visualize supplementary individuals

The data sets decathlon2 contain some supplementary individuals from row 24 to 27.

# 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.02      7.30    14.77      2.04 48.37        14.09
WARNERS  11.11      7.60    14.31      1.98 48.68        14.23
Nool     10.80      7.53    14.26      1.88 48.81        14.80
Drews    10.87      7.38    13.07      1.88 48.51        14.01

Individuals factor map using FactoMineR base graph :

plot(res.pca, choix="ind")

Supplementary individuals are shown in blue. The levels of the supplementary qualitative variable are shown in magnenta color.

The results for supplementary individuals can be extracted as follow :

res.pca$ind.sup

$coord
              Dim.1       Dim.2      Dim.3      Dim.4       Dim.5
KARPOV    0.7947206  0.77951227 -1.6330203  1.7242283 -0.75070396
WARNERS  -0.3864645 -0.12159237 -1.7387332 -0.7063341 -0.03230011
Nool     -0.5591306  1.97748871 -0.4830358 -2.2784526 -0.25461493
Drews    -1.1092038  0.01741477 -3.0488182 -1.5343468 -0.32642192

$cos2
              Dim.1        Dim.2      Dim.3      Dim.4        Dim.5
KARPOV   0.05104677 4.911173e-02 0.21553730 0.24028620 0.0455487744
WARNERS  0.02422707 2.398250e-03 0.49039677 0.08092862 0.0001692349
Nool     0.02897149 3.623868e-01 0.02162236 0.48108780 0.0060077529
Drews    0.09207094 2.269527e-05 0.69560547 0.17617609 0.0079736753

$dist
 KARPOV  WARNERS     Nool    Drews  
3.517470 2.482899 3.284943 3.655527 

Supplementary qualitative variables

The data sets decathlon2 contain a supplementary qualitative variable at columns 13 corresponding to the type of competitions.

Qualitative variable can be helpful for interpreting the data and for coloring individuals by groups.

The argument habillage is used to specify the index of the supplementary qualitative variable :

plot(res.pca, choix = "ind", habillage = 13)

It’s also possible to use factoextra :

fviz_pca_ind(res.pca, habillage = 13,
  addEllipses =TRUE, ellipse.level = 0.68) +
  scale_color_brewer(palette="Dark2") +
  theme_minimal()

Supplementary individuals are shown in blue color

The results concerning the supplementary qualitative variable are :

res.pca$quali

$coord
             Dim.1      Dim.2       Dim.3      Dim.4      Dim.5
Decastar -1.343451  0.1218097 -0.03789524  0.1808357  0.1343364
OlympicG  1.231497 -0.1116589  0.03473730 -0.1657661 -0.1231417

$cos2
             Dim.1       Dim.2        Dim.3      Dim.4       Dim.5
Decastar 0.9051233 0.007440939 0.0007201669 0.01639956 0.009050062
OlympicG 0.9051233 0.007440939 0.0007201669 0.01639956 0.009050062

$v.test
             Dim.1      Dim.2      Dim.3      Dim.4      Dim.5
Decastar -2.970766  0.4034256 -0.1528767  0.8971036  0.7202457
OlympicG  2.970766 -0.4034256  0.1528767 -0.8971036 -0.7202457

$dist
Decastar OlympicG 
1.412108 1.294433 

$eta2
                Dim 1      Dim 2       Dim 3      Dim 4      Dim 5
Competition 0.4011568 0.00739783 0.001062332 0.03658159 0.02357972