Correspondence Analysis in R , with Two-and Three-dimensional Graphics: The ca Package

We describe an implementation of simple, multiple and joint correspondence analysis in R . The resulting package comprises two parts, one for simple correspondence analysis and one for multiple and joint correspondence analysis. Within each part, functions for computation, summaries and visualization in two and three dimensions are provided, including options to display supplementary points and perform subset analyses. Special emphasis has been put on the visualization functions that oﬀer features such as diﬀerent scaling options for biplots and three-dimensional maps using the rgl package. Graphical options include shading and sizing plot symbols for the points according to their contributions to the map and masses respectively.


Introduction
The geometric interpretation of correspondence analysis originated in the research and teaching of Jean-Paul Benzécri in France -the classic reference is Benzécri (1973).Interest in correspondence analysis increased in the late 1980s and 1990s, and simple and multiple correspondence analysis were introduced into most of the mainstream statistical software packages.In R (R Development Core Team 2007) the functions corresp() and mca() (from the MASS package, Venables and Ripley 2002) provide a facility for the computation of CA and MCA.However, the implementation in these functions is kept to a minimum.For example, the simple CA function does not offer the inclusion of supplementary rows or columns.In the case of MCA, the implementation offers neither supplementary points, nor more recent developments such as adjustment of eigenvalues for improved fit and corresponding adjustments of contributions, joint correspondence analysis (JCA) and subset analysis.Apart from the XL-STAT software (Addinsoft 2007), none of the major programs offers these recent developments either.
This paper illustrates the implementation of CA and MCA in the R package ca (as of version 0.1-9, October 2006).The package comprises functions for simple, multiple and joint CA with support for subset analyses and the inclusion of supplementary variables.Furthermore, it offers functions for the graphical display of the results in two and three dimensions.
The package comprises the following components: • Simple CA: -Computation: ca() -Printing and Summaries: print() and summary() methods for ca objects -Plotting: plot.ca() and plot3d.ca() • MCA and JCA: -Computation: mjca() -Printing and Summaries: print() and summary() methods for ca objects -Plotting: plot.mjca() and plot3d.mjca() • Datasets: smoke, author and wg93 The package contains further functions, such as iterate.mjca()for the updating of the Burt matrix in JCA.The remaining sections describe the functions for the various forms of CA along with selected examples.Since the visualization functions are very similar for both cases, they are jointly covered in one section.

Simple correspondence analysis
As in principal component analysis, the idea in CA is to reduce the dimensionality of a data matrix and visualize it in a subspace of low-dimensionality, commonly two-or threedimensional.The data of interest in simple CA are usually a two-way contingency table or any other table of nonnegative ratio-scale data for which relative values are of primary interest.The CA solution was shown by (Greenacre 1984, Chapter 2 and Appendix) to be neatly encapsulated in the singular-value decomposition (SVD) of a suitably transformed matrix.To summarize the theory, first divide the I × J data matrix, denoted by N, by its grand total n to obtain the so-called correspondence matrix P = N/n.Let the row and column marginal totals of P be the vectors r and c respectively, that is the vectors of row and column masses, and D r and D c be the diagonal matrices of these matrices.The computational algorithm to obtain coordinates of the row and column profiles with respect to principal axes, using the SVD, is as follows: 1. Calculate the matrix of standarized residuals: The rows of the coordinate matrices in (3)-( 6) above refer to the rows or columns, as the case may be, of the original table, while the columns of these matrices refer to the principal axes, or dimensions, of the solution.Notice that the row and column principal coordinates are scaled in such a way that FD r F = GD c G = D 2 α , i.e. the weighted sum-of-squares of the coordinates on the k-th dimension (i.e., their inertia in the direction of this dimension) is equal to the principal inertia (or eigenvalue) α 2 k , the square of the k-th singular value, whereas the standard coordinates have weighted sum-of-squares equal to 1: XD r X = YD c Y = I.
The implementation of the algorithm follows Blasius and Greenacre (1994).The function ca() computes a simple CA, for example R> data("smoke") R> ca(smoke) performs a simple CA on the provided smoke dataset (Greenacre 1984), this dataset has become a test example for CA in various software packages, and is also discussed in Greenacre (1993Greenacre ( , 2007)).This dataset contains frequencies of smoking habits (none, light, medium and heavy) for staff groups (senior managers, junior managers, senior employees, junior employees and secretaries) in a fictional company.
The output of ca() is controlled by the printing method for CA, i.e. the following output is given by default: The output contains the eigenvalues and percentages of explained inertia for all possible dimensions.Additionally, values for the rows and columns (masses, chi-squared distances of points to their average, inertias and standard coordinates are given -again, see Greenacre (2007) for more details about these concepts).However, these values are restricted to two dimensions where applicable, e.g. for the standard coordinates.A list of all available entries that are returned by ca() is obtained with names():

R> names(ca(smoke))
The output of ca is structured as a list-object, for example, the row standard coordinates are obtained with

R> ca(smoke)$rowcoord
Optional arguments for the ca() function include an option for setting the dimensionality of the solution (nd), options for marking selected rows and/or columns as supplementary ones (suprow and supcol, respectively) and options for setting subset rows and/or columns (subsetrow and subsetcol, respectively) for subset CA (Greenacre and Pardo 2006a,b).The subset option restricts the analysis to the selected subset(s) while maintaining the original margins of the table.As an extension to the printing method, a summary method is also provided.This gives a more detailed output in the classic style of Tabet (1973).

R> summary(ca(smoke))
returns the summary of the CA: Again, eigenvalues and relative percentages of explained inertia are given for all available dimensions.Additionally, cumulated percentages and a scree plot are shown.The items given in Rows and Columns include the principal coordinates for the first two dimensions (k = 1 and k = 2).Squared correlations (cor) and contributions (ctr) for the points are displayed next to the coordinates.Notice that the quantities in these tables are multiplied by 1000 (e.g., the coordinates and masses) which for the cor and ctr quantities means they are expressed in permills.The total quality (qlt) is given with respect to the dimensionality of the solution, i.e. in this case it is the sum of the squared correlations over the two included dimensions.In the case of supplementary variables, an asterisk is appended to the supplementary variable names in the output.For example, the summary of a CA on the smoke data where the "none" category (i.e. the first column) is treated as a supplementary one is given by: R> summary(ca(smoke, supcol = 1)) In the corresponding section of the output the following is given: Supplementary variables have no impact on the computation.They are projected onto the solution space afterwards.Thus, contributions are not applicable for this case.Squared correlations (cor) as a measure of how well a point is represented by the axes are still meaningful for the case of supplementary variables and thus are included in the output.

Multiple and joint correspondence analysis
Multiple and Joint Correspondence Analysis (MCA and JCA, respectively) are extensions of simple CA of a single cross-tabulation to more than two categorical variables.More details about these methods can be found in Greenacre (2007), Chapters 18 and 19 respectively, while the computation of MCA and JCA is described in detail in the appendix of Nenadić and Greenacre (2006).Essentially, four approaches for the computation are considered.
The classic approach to MCA is to perform a simple CA on the indicator matrix, i.e. by performing a SVD on the matrix of standardized residuals, as shown previously, calculated on the indicator matrix.The indicator matrix Z is the cases×variables matrix with columns being dummy variables (with values only 0 or 1), a dummy variable for each category of the set of categorical variables.This approach yields principal inertias and principal coordinates equal to the eigenvalues and scale values in homogeneity analysis in the Gifi system (Michailidis and de Leeuw 1998).
An almost equivalent and more preferable approach from a CA point of view is given by performing an eigenvalue-eigenvector decomposition based on the Burt matrix, which is equal to the cross-product of the indicator matrix, Z Z: i.e., the matrix which concatenates all twoway cross-tabulations between pairs of variables.Due to the structure of the Burt matrix, with submatrices on the main diagonal that are cross-tabulations of each variable with itself, the solution overestimates the total inertia.
In order to overcome this problem, two alternative approaches are considered, namely the adjustment of inertias (Greenacre 1993) and joint correspondence analysis (Greenacre 1988).The adjustment approach improves the MCA solution by rescaling the coordinates of the solution to best fit the pairwise cross-tabulations off the main diagonal of the Burt matrix.JCA is a different iterative algorithm which finds the optimal weighted least-squares fit to these off-diagonal tables.
MCA and JCA are performed with the function mjca().The structure of the function is kept similar to its counterpart from simple CA.The two most striking differences are the format of the input data and the restriction to columns for the analyses.The function mjca() takes a response pattern matrix as input.In R terms this is a data frame with the columns containing factors.Within the function, the response pattern matrix is converted to an indicator matrix and a Burt matrix, depending on the type of analysis.The restriction to columns means that by default, only values for the columns are given in the output.Also, the specification of supplementary variables or of a subset is limited to columns.The "approach" to MCA is specified by the lambda option in mjca(): • lambda = "indicator": Analysis based on a simple CA of the indicator matrix • lambda = "Burt": Analysis based on an eigenvalue-decomposition of the Burt matrix • lambda = "adjusted": Analysis based on the Burt matrix with an adjustment of inertias • lambda = "JCA": Joint correspondence analysis By default, mjca() performs an adjusted analysis, i.e. lambda = "adjusted", which we believe to be the best default option, since the optimal scaling properties of MCA are conserved while raising the percentages of inertia and squared correlations to be usually very close to those one would get with a JCA.In the case of a full-blooded JCA, which involves an updating of the Burt matrix by iteratively weighted least squares, the auxiliary function iterate.mjca() is internally used for the updating of the Burt matrix.The updating function has two convergence criteria, namely epsilon and maxit.The option epsilon sets a convergence criterion by means of maximum absolute difference of the Burt matrix in an iteration step compared to the Burt matrix of the previous step.The maximum number of iterations is given by the option maxit.This way, the program iterates until any one of the two conditions is satisfied.Setting one option to NA results in that criterion being ignored -for example, exactly 50 iterations without considering convergence are performed with maxit=50 and epsilon=NA.
As with simple CA, the solution is restricted by the nd option to two dimensions.However, eigenvalues are given for all possible dimensions, which is equal to (J − Q) for the "indicator" and "Burt" case.In the case of an adjusted analysis or a JCA, the eigenvalues are given only for those dimensions k, where the singular values from the Burt matrix λ k (i.e., the principal inertias of the indicator matrix) satisfy the condition λ k > 1/Q.For example, a MCA (based on an adjusted analysis) is performed with R> data("wg93") R> mjca(wg93[,1:4]) In this case a MCA of the first four columns of the provided dataset wg93 (taken from the International Social Survey Programme 1993, see http://www.issp.org/) is performed.These columns contain attitudes of 871 individuals towards science and the environment (see Greenacre 2006a, Chapter 2, for more details on these data).Each category contains five possible answers (strongly agree, somewhat agree, neither agree nor disagree, somewhat disagree, strongly disagree, coded as 1 to 5).Thus, the output labels are given by appending the level names to the category names: Eigenvalues Notice that the percentages of inertia do not add up to 100% in the adjusted analysis.By entering R> summary(mjca(wg93[,1:4], lambda = "Burt")) a summary of a MCA based on the Burt matrix is given in the same style as for simple CA:  The graphical representation of results from CA and MCA is commonly done with so-called symmetric maps.In that case, the row and column coordinates on each axis are scaled to have inertias equal to the principal inertia along that axis: these are the principal row and column coordinates.Depending on the situation, other types of display are appropriate.This can be set with the scaling option map in the plotting functions for CA and MCA.Table 1 gives a brief overview over the available options and their meanings.option description "symmetric" Rows and columns in principal coordinates (default) "rowprincipal" Rows in principal and columns in standard coordinates "colprincipal" Rows in standard and columns in principal coordinates "symbiplot" Row and column coordinates are scaled to have variances equal to the singular values "rowgab" Rows in principal coordinates and columns in standard coordinates times mass "colgab" Columns in principal coordinates and rows in standard coordinates times mass (according to a proposal by Gabriel and Odoroff 1990) "rowgreen" Rows in principal coordinates and columns in standard coordinates times the square root of the mass "colgreen" Columns in principal coordinates and rows in standard coordinates times the square root of the mass (according to a proposal by Greenacre 2006b) Table 1: Scaling options in plot.ca and plot.mjca By default, a symmetric map is plotted.For example, a symmetric map of a CA of the smoke data (with the "none" smoking category set as a supplementary one) is created with the plot method for CA: R> plot(ca(smoke, supcol = 1)) The result is shown in Figure 4.The symmetric map is not a true biplot but all other options are.
By default, supplementary variables are added to the plot with a different symbol.The symbols can be defined with the pch option in plot.ca().This option takes four values in the following order: plotting point character or symbol for i.) active rows, ii.) supplementary rows, iii.) active columns and iv.) supplementary columns.As a general rule, options that contain entries for rows and for columns contain the entries for the rows first and than the entries for the columns.(For example, the colour of the symbols is specified with the col option, by default it is col = c("#000000", "#FF0000"), i.e. black for the rows and red for the columns.)The option what controls the content of the plot.It can be set to "all", "active", "passive" or "none" for the rows and for the columns.For example, a plot of only the active (i.e.not supplementary) points is created by using what=c("active", "active").
Apart from the scaling option with map, the plotting methods offer options for the inclusion Similarly, relative or absolute contributions can be indicated by the colour intensity in the plot by using the contrib option.Figure 4 shows the resulting plot of R> plot(ca(smoke), mass = TRUE, contrib = "absolute", + map = "rowgreen", arrows = c(FALSE, TRUE)) Greenacre (2006b) justifies the biplot options implemented here as "rowgreen" and "colgreen" calling them "standard biplot" because they give displays which function well for low and high inertia examples.The option dim controls which dimensions to plot.The default value is dim = c(1, 2), i.e. the first two dimensions are plotted.A plot of e.g. the second and third dimensions is obtained by setting dim = c(2, 3).Another possibility for adding the third dimension to the plot is given with the functions plot3d.ca()and plot3d.mjca().These two functions rely on the rgl package (Adler, Nenadić, and Zucchini 2003;Adler and Murdoch 2006) for a three-dimensional display in R. Their structure is kept similar to their counterparts for two dimensions, for example: R> plot3d.ca(ca(smoke,nd=3)) creates a three-dimensional display of the CA.The resulting display is shown in Figure 4: This type of display offers the advantage that one can zoom and navigate using the mouse.This way, maps of two dimensions are revealed by navigating to the appropriate viewpoint, e.g. to (0 o , 0 o ) in azimuthal coordinates for the first two dimensions or (−90 o , 0 o ) for the second and third dimension.

Summary
We have presented the R package ca for simple, multiple and joint correspondence analysis.This package contains all the features of present commercially available software packages as well as various new features that are not available elsewhere.Amongst these new features are the inclusion of adjustments in MCA, the corresponding adjustments in the percentages of inertia and squared correlations, the facility for subset analysis, joint correspondence analysis and provision of fully integrated three-dimensional graphics.