In Chapter 7, we saw how to summarize rectangular matrices whose columns were continuous variables. The maps we made used unsupervised dimensionality reduction techniques such as principal component analysis aimed at isolating the most important signal component in a matrix \(X\) when all the columns have meaningful variances.
Here we extend these ideas to more complex heterogeneous data where continuous and categorical variables are combined. Indeed, sometimes our observations cannot be easily described by sets of individual variables or coordinates – but it is possible to determine distances or (dis)similarities between them, or to describe relationships between them using a graph or a tree. Examples include species in a species tree or biological sequences. Outside of biology, examples include text documents or movie files, where we may have a reasonable method to determine (dis)similarity between them, but no obvious variables or coordinates.
This chapter contains more advanced techniques, for which we often omit technical details. Having come this far, we hope that by giving you some hands-on experience with examples, and extensive references, to enable you to understand and use some of the more `cutting edge’ techniques in nonlinear multivariate analysis.
In this chapter, we will:
Extend linear dimension reduction methods to cases when the distances between observations are available, known as multidimensional scaling (MDS) or principal coordinates analysis.
Find modifications of MDS that are nonlinear and robust to outliers.
Encode combinations of categorical data and continuous data as well as so-called ‘supplementary’ information. We will see that this enables us to deal with batch effects.
Use chi-square distances and correspondence analysis (CA) to see where categorical data (contingency tables) contain notable dependencies.
Generalize clustering methods that can uncover latent variables that are not categorical. This will allow us to detect gradients, “pseudotime” and hidden nonlinear effects in our data.
Generalize the notion of variance and covariance to the study of tables of data from multiple different data domains.
Sometimes, data are not represented as points in a feature space. This can occur when we are provided with (dis)similarity matrices between objects such as drugs, images, trees or other complex objects, which have no obvious coordinates in \({\mathbb R}^n\).
In Chapter 5 we saw how to produce clusters from distances. Here our goal is to visualize the data in maps in low dimensional spaces (e.g., planes) reminiscent of the ones we make from the first few principal axes in PCA.
We start with an intuitive example using geography data. In Figure 9.1, a heatmap and clustering of the distances between cities and places in Ukraine1 are shown.
1 The provenance of these data are described in the script ukraine-dists.R in the data folder.
library("pheatmap")
data("ukraine_dists", package = "MSMB")
as.matrix(ukraine_dists)[1:4, 1:4] Kyiv Odesa Sevastopol Chernihiv
Kyiv 0.0000 441.2548 687.7551 128.1287
Odesa 441.2548 0.0000 301.7482 558.6483
Sevastopol 687.7551 301.7482 0.0000 783.6561
Chernihiv 128.1287 558.6483 783.6561 0.0000pheatmap(as.matrix(ukraine_dists),
color = colorRampPalette(c("#0057b7", "#ffd700"))(50),
breaks = seq(0, max(ukraine_dists)^(1/2), length.out = 51)^2,
treeheight_row = 10, treeheight_col = 10)
ukraine_dists distance matrix. Distances are measured in kilometres. The function has re-arranged the order of the cities, and grouped the closest ones.
Besides ukraine_dists, which contains the pairwise distances, the RData file that we loaded above also contains the dataframe ukraine_coords with the longitudes and latitudes; we will use this later as a ground truth. Given the distances, multidimensional scaling (MDS) provides a “map” of their relative locations. It will not be possible to arrange the cities such that their Euclidean distances on a 2D plane exactly reproduce the given distance matrix: the cities lie on the curved surface of the Earth rather than in a plane. Nevertheless, we can expect to find a two dimensional embedding that represents the data well. With biological data, our 2D embeddings are likely to be much less clearcut. We call the function with:
ukraine_mds = cmdscale(ukraine_dists, eig = TRUE)We make a function that we will reuse several times in this chapter to make a screeplot from the result of a call to the cmdscale function:
library("dplyr")
library("ggplot2")
plotscree = function(x, m = length(x$eig)) {
ggplot(tibble(eig = x$eig[seq_len(m)], k = seq_along(eig)),
aes(x = k, y = eig)) + theme_minimal() +
scale_x_discrete("k", limits = as.factor(seq_len(m))) +
geom_bar(stat = "identity", width = 0.5, fill = "#ffd700", col = "#0057b7")
}plotscree(ukraine_mds, m = 4)
If you execute:
ukraine_mds$eig |> signif(3) [1] 3.91e+06 1.08e+06 3.42e+02 4.84e-01 2.13e-01 3.83e-05 5.90e-06
[8] 5.83e-07 8.79e-08 4.95e-08 7.70e-10 2.41e-10 1.82e-10 1.64e-10
[15] 9.17e-11 5.76e-11 -3.66e-11 -4.63e-11 -6.30e-11 -6.86e-11 -1.36e-10
[22] -1.44e-10 -1.78e-10 -3.75e-10 -4.24e-10 -2.35e-09 -1.52e-08 -9.60e-05
[29] -2.51e-04 -1.41e-02 -1.19e-01 -3.58e+02 -8.85e+02plotscree(ukraine_mds)
you will note that unlike in PCA, there are some negative eigenvalues. These are due to the way cmdscale works.
The main output from the cmdscale function are the coordinates of the two-dimensional embedding, which we show in Figure 9.4 (we will discuss how the algorithm works in the next section).
ukraine_mds_df = tibble(
PCo1 = ukraine_mds$points[, 1],
PCo2 = ukraine_mds$points[, 2],
labs = rownames(ukraine_mds$points)
)
library("ggrepel")
g = ggplot(ukraine_mds_df, aes(x = PCo1, y = PCo2, label = labs)) +
geom_point() + geom_text_repel(col = "#0057b7") + coord_fixed()
g
Note that while relative positions are correct, the orientation of the map is unconventional: Crimea is at the top. This is a common phenomenon with methods that reconstruct planar embeddings from distances. Since the distances between the points are invariant under rotations and reflections (axis flips), any solution is as good as any other solution that relates to it via rotation or reflection. Functions like cmdscale will pick one of the equally optimal solutions, and the particular choice can depend on minute details of the data or the computing platform being used. Here, we can transform our result into a more conventional orientation by reversing the sign of the \(y\)-axis. We redraw the map in Figure 9.5 and compare this to the true longitudes and latitudes from the ukraine_coords dataframe (Figure 9.6).
g + mutate(ukraine_mds_df, PCo1 = PCo1, PCo2 = -PCo2)
data("ukraine_coords", package = "MSMB")
print.data.frame(ukraine_coords[1:4, c("city", "lat", "lon")]) city lat lon
1 Kyiv 50.45003 30.52414
2 Odesa 46.48430 30.73229
3 Sevastopol 44.60544 33.52208
4 Chernihiv 51.49410 31.29433ggplot(ukraine_coords, aes(x = lon, y = lat, label = city)) +
geom_point() + geom_text_repel(col = "#0057b7")
ukraine_coords dataframe.
Note: MDS creates similar output as PCA, however there is only one ‘dimension’ to the data (the sample points). There is no ‘dual’ dimension, there are no biplots and no loading vectors. This is a drawback when coming to interpreting the maps. Interpretation can be facilitated by examining carefully the extreme points and their differences.
Let’s take a look at what would happen if we really started with points whose coordinates were known2. We put these coordinates into the two columns of a matrix with as many rows as there are points. Now we compute the distances between points based on these coordinates. To go from the coordinates \(X\) to distances, we write \[d^2_{i,j} = (x_i^1 - x_j^1)^2 + \dots + (x_i^p - x_j^p)^2.\] We will call the matrix of squared distances DdotD in R and \(D\bullet D\) in the text . We want to find points such that the square of their distances is as close as possible to the \(D\bullet D\) observed.3
3 Here we commit a slight ‘abuse’ by using the longitudes and latitudes of our cities as Cartesian coordinates and ignoring the fact that they are curvilinear coordinates on a sphere-like surface.
2 Here we commit a slight ‘abuse’ by using the longitude and longitude of our cities as Cartesian coordinates and ignoring the curvature of the earth’s surface. Check out the internet for information on the Haversine formula.
X = with(ukraine_coords, cbind(lon, lat * cos(48)))
DdotD = as.matrix(dist(X)^2)The relative distances do not depend on the point of origin of the data. We center the data by using the centering matrix \(H\) defined as \(H=I-\frac{1}{n}{\mathbf{11}}^t\). Let’s check the centering property of \(H\) using:
n = nrow(X)
H = diag(rep(1,n))-(1/n) * matrix(1, nrow = n, ncol = n)
Xc = sweep(X,2,apply(X,2,mean))
Xc[1:2, ] lon
[1,] -1.1722946 -1.184705
[2,] -0.9641429 1.353935HX = H %*% X
HX[1:2, ] lon
[1,] -1.1722946 -1.184705
[2,] -0.9641429 1.353935apply(HX, 2, mean) lon
-1.530923e-15 2.016353e-16 Therefore, given the squared distances between rows (\(D\bullet D\)) and the cross product of the centered matrix \(B=(HX)(HX)^t\), we have shown:
\[ -\frac{1}{2} H(D\bullet D) H=B \tag{9.1}\]
This is always true, and we use it to reverse-engineer an \(X\) which satisfies Equation 9.1 when we are given \(D\bullet D\) to start with.
We can go backwards from a matrix \(D\bullet D\) to \(X\) by taking the eigen-decomposition of \(B\) as defined in Equation 9.1. This also enables us to choose how many coordinates, or columns, we want for the \(X\) matrix. This is very similar to how PCA provides the best rank \(r\) approximation.
Note: As in PCA, we can write this using the singular value decomposition of \(HX\) (or the eigen decomposition of \(HX(HX)^t\)):
\[ HX^{(r)} = US^{(r)}V^t \mbox{ with } S^{(r)} \mbox{ the diagonal matrix of the first } r \mbox{ singular values}, \]
This provides the best approximate representation in an Euclidean space of dimension \(r\). The algorithm gives us the coordinates of points that have approximately the same distances as those provided by the \(D\) matrix.
In summary, given an \(n \times n\) matrix of squared interpoint distances \(D\bullet D\), we can find points and their coordinates \(\tilde{X}\) by the following operations:
Double center the interpoint distance squared and multiply it by \(-\frac{1}{2}\):
\(B = -\frac{1}{2}H D\bullet D H\).
Diagonalize \(B\): \(\quad B = U \Lambda U^t\).
Extract \(\tilde{X}\): \(\quad \tilde{X} = U \Lambda^{1/2}\).
As an example, let’s take objects for which we have similarities (surrogrates for distances) but for which there is no natural underlying Euclidean space.
In a psychology experiment from the 1950s, Ekman (1954) asked 31 subjects to rank the similarities of 14 different colors. His goal was to understand the underlying dimensionality of color perception. The similarity or confusion matrix was scaled to have values between 0 and 1. The colors that were often confused had similarities close to 1. We transform the data into a dissimilarity by subtracting the values from 1:
ekm = read.table("../data/ekman.txt", header=TRUE)
rownames(ekm) = colnames(ekm)
disekm = 1 - ekm - diag(1, ncol(ekm))
disekm[1:5, 1:5] w434 w445 w465 w472 w490
w434 0.00 0.14 0.58 0.58 0.82
w445 0.14 0.00 0.50 0.56 0.78
w465 0.58 0.50 0.00 0.19 0.53
w472 0.58 0.56 0.19 0.00 0.46
w490 0.82 0.78 0.53 0.46 0.00disekm = as.dist(disekm)We compute the MDS coordinates and eigenvalues. We combine the eigenvalues in the screeplot shown in Figure 9.8:
mdsekm = cmdscale(disekm, eig = TRUE)
plotscree(mdsekm)
We plot the different colors using the first two principal coordinates as follows:
dfekm = mdsekm$points[, 1:2] |>
`colnames<-`(paste0("MDS", 1:2)) |>
as_tibble() |>
mutate(
name = rownames(ekm),
rgb = photobiology::w_length2rgb(as.numeric(sub("w", "", name))))
ggplot(dfekm, aes(x = MDS1, y = MDS2)) +
geom_point(col = dfekm$rgb, size = 4) +
geom_text_repel(aes(label = name)) + coord_fixed()
Figure 9.9 shows the Ekman data in the new coordinates. There is a striking pattern that calls for explanation. This horseshoe or arch structure in the points is often an indicator of a sequential latent ordering or gradient in the data (Diaconis et al. 2008). We will revisit this in Section 9.5.
Multidimensional scaling aims to minimize the difference between the squared distances as given by \(D\bullet D\) and the squared distances between the points with their new coordinates. Unfortunately, this objective tends to be sensitive to outliers: one single data point with large distances to everyone else can dominate, and thus skew, the whole analysis. Often, we like to use something that is more robust, and one way to achieve this is to disregard the actual values of the distances and only ask that the relative rankings of the original and the new distances are as similar as possible. Such a rank based approach is robust: its sensitivity to outliers is reduced.
We will use the Ekman data to show how useful robust methods are when we are not quite sure about the ‘scale’ of our measurements. Robust ordination, called non metric multidimensional scaling (NMDS for short) only attempts to embed the points in a new space such that the order of the reconstructed distances in the new map is the same as the ordering of the original distance matrix.
Non metric MDS looks for a transformation \(f\) of the given dissimilarities in the matrix \(d\) and a set of coordinates in a low dimensional space (the map) such that the distance in this new map is \(\tilde{d}\) and \(f(d)\thickapprox \tilde{d}\). The quality of the approximation can be measured by the standardized residual sum of squares (stress) function:
\[ \text{stress}^2=\frac{\sum(f(d)-\tilde{d})^2}{\sum d^2}. \]
NMDS is not sequential in the sense that we have to specify the underlying dimensionality at the outset and the optimization is run to maximize the reconstruction of the distances according to that number. There is no notion of percentage of variation explained by individual axes as provided in PCA. However, we can make a simili-screeplot by running the program for all the successive values of \(k\) (\(k=1, 2, 3, ...\)) and looking at how well the stress drops. Here is an example of looking at these successive approximations and their goodness of fit. As in the case of diagnostics for clustering, we will take the number of axes after the stress has a steep drop.
Because each calculation of a NMDS result requires a new optimization that is both random and dependent on the \(k\) value, we use a similar procedure to what we did for clustering in Chapter 4. We execute the metaMDS function, say, 100 times for each of the four possible values of \(k\) and record the stress values.
library("vegan")
nmds.stress = function(x, sim = 100, kmax = 4) {
sapply(seq_len(kmax), function(k)
replicate(sim, metaMDS(x, k = k, autotransform = FALSE)$stress))
}
stress = nmds.stress(disekm, sim = 100)
dim(stress)Let’s look at the boxplots of the results. This can be a useful diagnostic plot for choosing \(k\) (Figure 9.10).
dfstr = reshape2::melt(stress, varnames = c("replicate","dimensions"))
ggplot(dfstr, aes(y = value, x = dimensions, group = dimensions)) +
geom_boxplot()
We can also compare the distances and their approximations using what is known as a Shepard plot for \(k=2\) for instance, computed with:
nmdsk2 = metaMDS(disekm, k = 2, autotransform = FALSE)
stressplot(nmdsk2, pch = 20)
Both the Shepard’s plot in Figure 9.11 and the screeplot in Figure 9.10 point to a two-dimensional solution for Ekman’s color confusion study. Let us compare the output of the two different MDS programs, the classical metric least squares approximation and the nonmetric rank approximation method. The right panel of Figure 9.12 shows the result from the nonmetric rank approximation, the left panel is the same as Figure 9.9. The projections are almost identical in both cases. For these data, it makes little difference whether we use a Euclidean or nonmetric multidimensional scaling method.
nmdsk2$points[, 1:2] |>
`colnames<-`(paste0("NmMDS", 1:2)) |>
as_tibble() |>
bind_cols(dplyr::select(dfekm, rgb, name)) |>
ggplot(aes(x = NmMDS1, y = NmMDS2)) +
geom_point(col = dfekm$rgb, size = 4) +
geom_text_repel(aes(label = name))
In Chapter 3 we introduced the R data.frame class that enables us to combine heterogeneous data types: categorical factors, text and continuous values. Each row of a dataframe corresponds to an object, or a record, and the columns are the different variables, or features.
Extra information about sample batches, dates of measurement, different protocols are often named metadata; this can be misnomer if it is implied that metadata are somehow less important. Such information is real data that need to be integrated into the analyses. We typically store it in a data.frame or a similar R class and tightly link it to the primary assay data.
Here we show an example of an analysis that was done by Holmes et al. (2011) on bacterial abundance data from Phylochip (Brodie et al. 2006) microarrays. The experiment was designed to detect differences between a group of healthy rats and a group who had Irritable Bowel Disease (Nelson et al. 2010). This example shows a case where the nuisance batch effects become apparent in the analysis of experimental data. It is an illustration of the fact that best practices in data analyses are sequential and that it is better to analyse data as they are collected to adjust for severe problems in the experimental design as they occur, instead of having to deal with them post mortem4.
4 Fisher’s terminology, see Chapter 13.
When data collection started on this project, days 1 and 2 were delivered and we made the plot that appears in Figure 9.14. This showed a definite day effect. When investigating the source of this effect, we found that both the protocol and the array were different in days 1 and 2. This leads to uncertainty in the source of variation, we call this confounding of effects.
day information has been combined with the array data and encoded as a number and could be confused with a continuous variable. We will see in the next section a better practice for storing and manipulating heterogeneous data using a Bioconductor container called SummarizedExperiment.We load the data and the packages we use for this section:
IBDchip = readRDS("../data/vsn28Exprd.rds")
library("ade4")
library("factoextra")
library("sva")The data are normalized abundance measurements of 8634 taxa measured on 28 samples. We use a rank-threshold transformation, giving the top 3000 most abundant taxa scores from 3000 to 1, and letting the remaining (low abundant) ones all have a score of 1. We also separate out the proper assay data from the (awkwardly placed) day variable, which should be considered a factor5:
5 Below, we show how to arrange these data into a Bioconductor SummarizedExperiment, which is a much more sane way of storing such data.
assayIBD = IBDchip[-nrow(IBDchip), ]
day = factor(IBDchip[nrow(IBDchip), ])Instead of using the continuous, somehow normalized data, we use a robust analysis replacing the values by their ranks. The lower values are considered ties encoded as a threshold chosen to reflect the number of expected taxa thought to be present:
rankthreshPCA = function(x, threshold = 3000) {
ranksM = apply(x, 2, rank)
ranksM[ranksM < threshold] = threshold
ranksM = threshold - ranksM
dudi.pca(t(ranksM), scannf = FALSE, nf = 2)
}
pcaDay12 = rankthreshPCA(assayIBD[, day != 3])
fviz_eig(pcaDay12, bar_width = 0.6) + ggtitle("")
day12 = day[ day!=3 ]
rtPCA1 = fviz(pcaDay12, element = "ind", axes = c(1, 2), geom = c("point", "text"),
habillage = day12, repel = TRUE, palette = "Dark2",
addEllipses = TRUE, ellipse.type = "convex") + ggtitle("") +
coord_fixed()
rtPCA1
Figure 9.14 shows that the sample arrange themselves naturally into two different groups according to the day of the samples. After discovering this effect, we delved into the differences that could explain these distinct clusters. There were two different protocols used (protocol 1 on day 1, protocol 2 on day 2) and unfortunately two different provenances for the arrays used on those two days (array 1 on day 1, array 2 on day 2).
A third set of data of four samples had to be collected to deconvolve the confounding effect. Array 2 was used with protocol 2 on Day 3, Figure 9.15 shows the new PCA plot with all the samples created by the following:
pcaDay123 = rankthreshPCA(assayIBD)
fviz(pcaDay123, element = "ind", axes = c(1, 2), geom = c("point", "text"),
habillage = day, repel = TRUE, palette = "Dark2",
addEllipses = TRUE, ellipse.type = "convex") +
ggtitle("") + coord_fixed()
Through this visualization we were able to uncover a flaw in the original experimental design. The first two batches shown in green and brown were both balanced with regards to IBS and healthy rats. They do show very different levels of variability and overall multivariate coordinates. In fact, there are two confounded effects. Both the arrays and protocols were different on those two days. We had to run a third batch of experiments on day 3, represented in purple, this used protocol from day 1 and the arrays from day 2. The third group faithfully overlaps with batch 1, telling us that the change in protocol was responsible for the variability.
Through the combination of the continuous measurements from assayIBD and the supplementary batch number as a factor, the PCA map has provided an invaluable investigation tool. This is a good example of the use of supplementary points6. The mean-barycenter points are created by using the group-means of points in each of the three groups and serve as extra markers on the plot.
6 This is called a supplementary point because the new observation-point is not used in the matrix decomposition.
We can decide to re-align the three groups by subtracting the group means so that all the batches are centered on the origin. A slightly more effective way is to use the ComBat function available in the sva package. This function uses a similar, but slightly more sophisticated method (Empirical Bayes mixture approach (Leek et al. 2010)). We can see its effect on the data by redoing our robust PCA (see the result in Figure 9.17):
model0 = model.matrix(~1, day)
combatIBD = ComBat(dat = assayIBD, batch = day, mod = model0)
pcaDayBatRM = rankthreshPCA(combatIBD)
fviz(pcaDayBatRM, element = "ind", geom = c("point", "text"),
habillage = day, repel=TRUE, palette = "Dark2", addEllipses = TRUE,
ellipse.type = "convex", axes =c(1,2)) + coord_fixed() + ggtitle("")
A more rational way of combining the batch and treatment information into compartments of a composite object is to use the SummarizedExperiment class. It includes special slots for the assay(s) where rows represent features of interest (e.g., genes, transcripts, exons, etc.) and columns represent samples. Supplementary information about the features can be stored in a DataFrame object, accessible using the function rowData. Each row of the DataFrame provides information on the feature in the corresponding row of the SummarizedExperiment object.
DataFrame is not the same as a data.frame.Here we insert the two covariates day and treatment in the colData object and combine it with assay data in a new SummarizedExperiment object.
library("SummarizedExperiment")
treatment = factor(ifelse(grepl("Cntr|^C", colnames(IBDchip)), "CTL", "IBS"))
sampledata = DataFrame(day = day, treatment = treatment)
chipse = SummarizedExperiment(assays = list(abundance = assayIBD),
colData = sampledata)This is the best way to keep all the relevant data together, it will also enable you to quickly filter the data while keeping all the information aligned properly.
chipse, some of the slots are empty.Columns of the DataFrame represent different attributes of the features of interest, e.g., gene or transcript IDs, etc. Here is an example of hybrid data container from single cell experiments (see Bioconductor workflow in Perraudeau et al. (2017) for more details).
corese = readRDS("../data/normse.rds")
norm = assays(corese)$normalizedValuesAfter the pre-processing and normalization steps prescribed in the workflow, we retain the 1000 most variable genes measured on 747 cells.
We can look at a PCA of the normalized values and check graphically that the batch effect has been removed:
respca = dudi.pca(t(norm), nf = 3, scannf = FALSE)
plotscree(respca, 15)
PCS = respca$li[, 1:3]
Since the screeplot in Figure 9.18 shows us that we must not dissociate axes 2 and 3, we will make a three dimensional plot with the rgl package. We use the following interactive code:
library("rgl")In chunk 'screeplotnorm-2': Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display
In chunk 'screeplotnorm-2': Warning: 'rgl.init' failed, will use the null device.See '?rgl.useNULL' for ways to avoid this warning.
batch = colData(corese)$Batch
plot3d(PCS,aspect=sqrt(c(84,24,20)),col=col_clus[batch])
plot3d(PCS,aspect=sqrt(c(84,24,20)),
col = col_clus[as.character(publishedClusters)])
Note: Of course, the book medium is limiting here, as we are showing two static projections that do not do justice to the depth available when looking at the interactive dynamic plots as they appear using the plot3d function. We encourage the reader to experiment extensively with these and other interactive packages and they provide a much more intuitive experience of the data.
Categorical data abound in biological settings: sequence status (CpG/non-CpG), phenotypes, taxa are often coded as factors as we saw in Chapter 2. Cross-tabulation of two such variables gives us a contingency table; the result of counting the co-occurrence of two phenotypes (sex and colorblindness was such an example). We saw that the first step is to look at the independence of the two categorical variables; the standard statistical measure of independence uses the chisquare distance. This quantity will replace the variance we used for continuous measurements.
The columns and rows of the table have the same `status’ and we are not in supervised/regression type setting. We won’t see a sample/variable divide; as a consequence the rows and columns will have the same status and we will ‘center’ both the rows and the columns. This symmetry will also translate in our use of biplots where both dimensions appear on the same plot.
| Patient | Mut1 | Mut2 | Mut3 | ... |
|---|---|---|---|---|
| AHX112 | 0 | 0 | 0 | |
| AHX717 | 1 | 0 | 1 | |
| AHX543 | 1 | 0 | 0 |
If the data are collected as long lists with each subject (or sample) associated to its levels of the categorical variables, we may want to transform them into a contingency table. Here is an example. In Table 9.1 HIV mutations are tabulated as indicator (0/1) binary variables. These data are then transformed into a mutation co-occurrence matrix shown in Table 9.2.
| Patient | Mut1 | Mut2 | Mut3 | ... |
|---|---|---|---|---|
| Mut1 | 853 | 29 | 10 | |
| Mut2 | 29 | 853 | 52 | |
| Mut3 | 10 | 52 | 853 |
Here are some co-occurrence data from the HIV database (Rhee et al. 2003). Some of these mutations have a tendency to co-occur.
Before explaining the details of how correspondence analysis works, let’s look at the output of one of many correspondence analysis functions. We use dudi.coa from the ade4 package to plot the mutations in a lower dimensional projection; the procedure follows what we did for PCA.
cooc = read.delim2("../data/coccurHIV.txt", header = TRUE, sep = ",")
cooc[1:4, 1:11] X4S X6D X6K X11R X20R X21I X35I X35L X35M X35T X39A
4S 0 28 8 0 99 0 22 5 15 3 45
6D 26 0 0 34 131 0 108 4 30 13 84
6K 7 0 0 6 45 0 5 13 38 35 12
11R 0 35 7 0 127 12 60 17 15 6 42HIVca = dudi.coa(cooc, nf = 4, scannf = FALSE)
fviz_eig(HIVca, geom = "bar", bar_width = 0.6) + ggtitle("")
plot3d).
After looking at a screeplot, we see that dimensionality of the underlying variation is definitely three dimensional, we plot these three dimensions. Ideally this would be done with an interactive three-dimensional plotting function such as that provided through the package rgl as shown in Figure 9.21.
fviz_ca_row(HIVca,axes = c(1, 2),geom="text", col.row="purple",
labelsize=3)+ggtitle("") + xlim(-0.55, 1.7) + ylim(-0.53,1.1) +
theme_bw() + coord_fixed()In chunk 'fig-HIVca': Warning: Removed 7 rows containing missing values or values outside the scale range(`geom_text()`).
fviz_ca_row(HIVca,axes = c(1, 3), geom="text",col.row="purple",
labelsize=3)+ggtitle("")+ xlim(-0.55, 1.7)+ylim(-0.5,0.6) +
theme_bw() + coord_fixed()In chunk 'fig-HIVca': Warning: Removed 7 rows containing missing values or values outside the scale range(`geom_text()`).
This first example showed how to map all the different levels of one categorical variable (the mutations) in a similar way to how PCA projects continuous variables. We will now explore how this can be extended to two or more categorical variables.
We will consider a small table, so we can follow the analysis in detail. The data are a contingency table of hair-color and eye-color phenotypic co-occurrence from students as shown in Table 9.3. In Chapter 2, we used a \(\chi^2\) test to detect possible dependencies:
HairColor = HairEyeColor[,,2]
chisq.test(HairColor)
Pearson's Chi-squared test
data: HairColor
X-squared = 106.66, df = 9, p-value < 2.2e-16| Brown | Blue | Hazel | Green | |
|---|---|---|---|---|
| Black | 36 | 9 | 5 | 2 |
| Brown | 66 | 34 | 29 | 14 |
| Red | 16 | 7 | 7 | 7 |
| Blond | 4 | 64 | 5 | 8 |
However, stating non independence between hair and eye color is not enough. We need a more detailed explanation of where the dependencies occur: which hair color occurs more often with green eyes ? Are some of the variable levels independent? In fact we can study the departure from independence using a special weighted version of SVD. This method can be understood as a simple extension of PCA and MDS to contingency tables.
We start by computing the row and column sums; we use these to build the table that would be expected if the two phenotypes were independent. We call this expected table HCexp.
rowsums = as.matrix(apply(HairColor, 1, sum))
rowsums [,1]
Black 52
Brown 143
Red 37
Blond 81colsums = as.matrix(apply(HairColor, 2, sum))
t(colsums) Brown Blue Hazel Green
[1,] 122 114 46 31HCexp = rowsums %*%t (colsums) / sum(colsums)Now we compute the \(\chi^2\) (chi-squared) statistic, which is the sum of the scaled residuals for each of the cells of the table:
sum((HairColor - HCexp)^2/HCexp)[1] 106.6637We can study these residuals from the expected table, first numerically then in Figure 9.23.
round(t(HairColor-HCexp)) Hair
Eye Black Brown Red Blond
Brown 16 10 2 -28
Blue -10 -18 -6 34
Hazel -3 8 2 -7
Green -3 0 3 0library("vcd")
mosaicplot(HairColor, shade=TRUE, las=1, type="pearson", cex.axis=0.7, main="")
Here are the computations we just did in R in a more mathematical form. For a general contingency table \({\mathbf N}\) with \(I\) rows and \(J\) columns and a total sample size of \(n=\sum_{i=1}^I \sum_{j=1}^J n_{ij}= n_{\cdot \cdot}\). If the two categorical variables were independent, each cell frequency would be approximately equal to
\[ n_{ij} = \frac{n_{i \cdot}}{n} \frac{n_{\cdot j}}{n} \times n \]
can also be written:
\[ {\mathbf N} = {\mathbf c r'} \times n, \qquad \mbox{ where } c= \frac{1}{n} {{\mathbf N}} {\mathbb 1}_m \;\mbox{ and }\; r'=\frac{1}{n} {\mathbf N}' {\mathbb 1}_p \]
The departure from independence is measured by the \(\chi^2\) statistic
\[ {\cal X}^2=\sum_{i,j} \frac{\left(n_{ij}-\frac{n_{i\cdot}}{n}\frac{n_{\cdot j}}{n}n\right)^2} {\frac{n_{i\cdot}n_{\cdot j}}{n^2}n} \]
Once we have ascertained that the two variables are not independent, we use a weighted multidimensional scaling using \(\chi^2\) distances to visualize the associations.
CCA in vegan, CA in FactoMineR, ordinate in phyloseq, dudi.coa in ade4.The method is called Correspondence Analysis (CA) or Dual Scaling and there are multiple R packages that implement it.
Here we make a simple biplot of the Hair and Eye colors.
HC = as.data.frame.matrix(HairColor)
coaHC = dudi.coa(HC,scannf=FALSE,nf=2)
round(coaHC$eig[1:3]/sum(coaHC$eig)*100)[1] 89 10 2fviz_ca_biplot(coaHC, repel=TRUE, col.col="brown", col.row="purple") +
ggtitle("") + ylim(c(-0.5,0.5))
CA has a special barycentric property: the biplot scaling is chosen so that the row points are placed at the center of gravity of the column levels with their respective weights. For instance, the Blue eyes column point is at the center gravity of the (Black, Brown, Red, Blond) with weights proportional to (9, 34, 7, 64). The Blond row point is very heavily weighted, this is why Figure 9.24 shows Blond and Blue quite close together.
All the methods we have studied in the last sections are commonly known as ordination methods. In the same way clustering allowed us to detect and interpret a hidden factor/categorical variable, ordination enables us to detect and interpret a hidden ordering, gradient or latent variable in the data.
Ecologists have a long history of interpreting the arches formed by observations points in correspondence analysis and principal components as ecological gradients (Prentice 1977). Let’s illustrate this first with a very simple data set on which we perform a correspondence analysis.
load("../data/lakes.RData")
lakelike[1:3,1:8] plant1 plant2 plant3 plant4 plant5 plant6 plant7 plant8
loc1 6 4 0 3 0 0 0 0
loc2 4 5 5 3 4 2 0 0
loc3 3 4 7 4 5 2 1 1reslake=dudi.coa(lakelike,scannf=FALSE,nf=2)
round(reslake$eig[1:8]/sum(reslake$eig),2)[1] 0.56 0.25 0.09 0.03 0.03 0.02 0.01 0.00We plot both the row-location points (Figure 9.25 (a)) and the biplot of both location and plant species in the lower part of Figure 9.25 (b); this plot was made with:
fviz_ca_row(reslake,repel=TRUE)+ggtitle("")+ylim(c(-0.55,1.7))
fviz_ca_biplot(reslake,repel=TRUE)+ggtitle("")+ylim(c(-0.55,1.7))
We will now analyse a more interesting data set that was published by Moignard et al. (2015). This paper describes the dynamics of blood cell development. The data are single cell gene expression measurements of 3,934 cells with blood and endothelial potential from five populations from between embryonic days E7.0 and E8.25.
Remember from Chapter 4 that several different distances are available for comparing our cells. Here, we start by computing both an \(L_2\) distance and the \(\ell_1\) distance between the 3,934 cells.
Moignard = readRDS("../data/Moignard.rds")
cellt = rowData(Moignard)$celltypes
colsn = c("red", "purple", "orange", "green", "blue")
blom = assay(Moignard)
dist2n.euclid = dist(blom)
dist1n.l1 = dist(blom, "manhattan")The classical multidimensional scaling on these two distances matrices can be carried out using:
ce1Mds = cmdscale(dist1n.l1, k = 20, eig = TRUE)
ce2Mds = cmdscale(dist2n.euclid, k = 20, eig = TRUE)
perc1 = round(100*sum(ce1Mds$eig[1:2])/sum(ce1Mds$eig))
perc2 = round(100*sum(ce2Mds$eig[1:2])/sum(ce2Mds$eig))We look at the underlying dimension and see in Figure 9.27 that two dimensions can provide a substantial fraction of the variance.
plotscree(ce1Mds, m = 4)
plotscree(ce2Mds, m = 4)
The first 2 coordinates account for 78 % of the variability when the \(\ell_1\) distance is used between cells, and 57% when the \(L^2\) distance is used. We see in Figure 9.28 (a) the first plane for the MDS on the \(\ell_1\) distances between cells:
c1mds = ce1Mds$points[, 1:2] |>
`colnames<-`(paste0("L1_PCo", 1:2)) |>
as_tibble()
ggplot(c1mds, aes(x = L1_PCo1, y = L1_PCo2, color = cellt)) +
geom_point(aes(color = cellt), alpha = 0.6) +
scale_colour_manual(values = colsn) + guides(color = "none")
c2mds = ce2Mds$points[, 1:2] |>
`colnames<-`(paste0("L2_PCo", 1:2)) |>
as_tibble()
ggplot(c2mds, aes(x = L2_PCo1, y = L2_PCo2, color = cellt)) +
geom_point(aes(color = cellt), alpha = 0.6) +
scale_colour_manual(values = colsn) + guides(color = "none")
Figure 9.28 (b) is created in the same way and shows the two-dimensional projection created by using MDS on the L2 distances.
Figure 9.28 shows that both distances (L1 and L2) give the same first plane for the MDS with very similar representations of the underlying gradient followed by the cells.
We can see from Figure 9.28 that the cells are not distributed uniformly in the lower dimensions we have been looking at, we see a definite organization of the points. All the cells of type 4SG represented in red form an elongated cluster who are much less mixed with the other cell types.
Multidimensional scaling and non metric multidimensional scaling aims to represent all distances as precisely as possible and the large distances between far away points skew the representations. It can be beneficial when looking for gradients or low dimensional manifolds to restrict ourselves to approximations of points that are close together. This calls for methods that try to represent local (small) distances well and do not try to approximate distances between faraway points with too much accuracy.
There has been substantial progress in such methods in recent years. The use of kernels computed using the calculated interpoint distances allows us to decrease the importance of points that are far apart. A radial basis kernel is of the form
\[ 1-\exp\left(-\frac{d(x,y)^2}{\sigma^2}\right), \quad\mbox{where } \sigma^2 \mbox{ is fixed.} \]
It has the effect of heavily discounting large distances. This can be very useful as the precision of interpoint distances is often better at smaller ranges; several examples of such methods are covered in Exercise 9.6 at the end of this chapter.
This widely used method adds flexibility to the kernel defined above and allows the \(\sigma^2\) parameter to vary locally (there is a normalization step so that it averages to one). The t-SNE method starts out from the positions of the points in the high dimensional space and derives a probability distribution on the set of pairs of points, such that the probabilities are proportional to the points’ proximities or similarities. It then uses this distribution to construct a representation of the dataset in low dimensions. The method is not robust and has the property of separating clusters of points artificially; however, this property can also help clarify a complex situation. One can think of it as a method akin to graph (or network) layout algorithms. They stretch the data to clarify relations between the very close (in the network: connected) points, but the distances between more distal (in the network: unconnected) points cannot be interpreted as being on the same scales in different regions of the plot. In particular, these distances will depend on the local point densities. Here is an example of the output of t-SNE on the cell data:
library("Rtsne")
restsne = Rtsne(blom, dims = 2, perplexity = 30, verbose = FALSE,
max_iter = 900)
dftsne = restsne$Y[, 1:2] |>
`colnames<-`(paste0("axis", 1:2)) |>
as_tibble()
ggplot(dftsne,aes(x = axis1, y = axis2, color = cellt)) +
geom_point(aes(color = cellt), alpha = 0.6) +
scale_color_manual(values = colsn) + guides(color = "none")
restsne3 = Rtsne(blom, dims = 3, perplexity = 30, verbose = FALSE,
max_iter = 900)
dftsne3 = restsne3$Y[, 1:3] |>
`colnames<-`(paste0("axis", 1:3)) |>
as_tibble()
ggplot(dftsne3,aes(x = axis3, y = axis2, group = cellt)) +
geom_point(aes(color = cellt), alpha = 0.6) +
scale_colour_manual(values = colsn) + guides(color = "none")
In this case in order to see the subtle differences between MDS and t-SNE, it is really necessary to use 3d plotting.
Two of these 3d snapshots are shown in Figure 9.30, we see a much stronger grouping of the purple points than in the MDS plots.
Note: A site worth visiting in order to appreciate more about the sensitivity of the t-SNE method to the complexity and \(\sigma\) parameters can be found at http://distill.pub/2016/misread-tsne.
ukraine_tsne = Rtsne(ukraine_dists, is_distance = TRUE, perplexity = 8)
ukraine_tsne_df = tibble(
PCo1 = ukraine_tsne$Y[, 1],
PCo2 = ukraine_tsne$Y[, 2],
labs = attr(ukraine_dists, "Labels")
)
ggplot(ukraine_tsne_df, aes(x = PCo1, y = PCo2, label = labs)) +
geom_point() + geom_text_repel(col = "#0057b7") + coord_fixed() 
There are several other nonlinear methods for estimating nonlinear trajectories followed by points in the relevant state spaces. Here are a few examples.
RDRToolbox Local linear embedding (LLE) and isomap
diffusionMap This package models connections between points as a Markovian kernel.
kernlab Kernel methods
LPCM-package Local principal curves
Current studies often attempt to quantify variation in the microbial, genomic, and metabolic measurements across different experimental conditions. As a result, it is common to perform multiple assays on the same biological samples and ask what features – microbes, genes, or metabolites, for example – are associated with different sample conditions. There are many ways to approach these questions. Which to apply depends on the study’s focus.
As in physics, we define inertia as a sum of distances with ‘weighted’ points. This enables us to compute the inertia of counts in a contingency table as the weighted sum of the squares of distances between observed and expected frequencies (as in the chisquare statistic).
If we want to study two standardized variables measured at the same 10 locations together, we use their covariance. If \(x\) represents the standardized PH, and and \(y\) the standardized humidity, we measure their covariation using the mean
\[ \text{cov}(x,y) = \text{mean}(x1*y1 + x2*y2 + x3*y3 + \cdots + x10*y10). \tag{9.2}\]
If \(x\) and \(y\) co-vary in the same direction, this will be big. We saw how useful the correlation coefficient we defined in Chapter 8 was to our multivariate analyses. Multitable generalizations will be just as useful.
The Mantel coefficient, one of the earliest version of association measures, is probably also the most popular now, especially in ecology (Josse and Holmes 2016). Given two dissimilarity matrices \(D^X\) and \(D^Y\) associated with \(X\) and \(Y\), make these matrices into vectors the way the R dist function does, and compute their linear correlation. A prototypical application is, for instance, to compute \(D^X\) from the soil chemistry at 17 different locations and to use \(D^Y\) to represent dissimilarities in plant occurences as measured by the Jaccard index between the same 17 locations. The Mantel coefficient is defined mathematically as:
\[ r_m(X, Y) = \frac{ \sum_{i\neq j} \left(d_{ij}^X − \bar{d}^X \right) \left(d_{ij}^Y − \bar{d}^Y \right) } { \sqrt{ \left( \sum_{i\neq j} \left(d_{ij}^X − \bar{d}^X \right)^2 \right) \left( \sum_{i\neq j} \left(d_{ij}^Y − \bar{d}^Y \right)^2 \right) } }, \]
with \(\bar{d}^X\) (resp. \(\bar{d}^Y\)) the mean of the upper diagonal elements of the dissimilarity matrix associated to \(d^X\) (resp. to \(d^Y\)). The main difference between the Mantel coefficient and the others such as the RV or the dCov is the absence of double centering. Due to the dependences within distances matrix, the Mantel correlation’s null distribution and its statistical significance cannot be assessed as simply for regular correlation ccoefficients. Instead, it is usually assessed via permutation testing. See Josse and Holmes (2016) for a review with historical background and modern incarnations. The coefficient and associated tests are implemented in several R packages such as ade4 (Chessel et al. 2004), vegan and ecodist (Goslee et al. 2007).
The global measure of similarity of two data tables as opposed to two vectors can be done by a generalization of covariance provided by an inner product between tables that gives the RV coefficient, a number between 0 and 1, like a correlation coefficient, but for tables.
\[ RV(A,B)=\frac{Tr(A'BAB')}{\sqrt{Tr(A'A)}\sqrt{Tr(B'B)}} \]
There are several other measures of matrix correlation available in the package MatrixCorrelation.
If we do ascertain a link between two matrices, we then need to find a way to understand that link. One such method is explained in the next section.
CCA is a method similar to PCA as it was developed by Hotelling in the 1930s to search for associations between two sets of continuous variables \(X\) and \(Y\). Its goal is to find a linear projection of the first set of variables that maximally correlates with a linear projection of the second set of variables.
Finding correlated functions (covariates) of the two views of the same phenomenon by discarding the representation-specific details (noise) is expected to reveal the underlying hidden yet influential factors responsible for the correlation.
Let us consider two matrices:
The \(p\) columns of \(X\) and the \(q\) columns of \(Y\) correspond to variables, and the rows correspond to the same \(n\) experimental units. We denote the \(j\)-th column of the matrix \(X\) by \(X_j\), likewise the \(k\)-th column of \(Y\) by \(Y_k\). Without loss of generality it will be assumed that the columns of \(X\) and \(Y\) are standardized (mean 0 and variance 1).
Classical CCA assumes that \(p \leq n\) and \(q \leq n\), and that the matrices \(X\) and \(Y\) are of full column rank \(p\) and \(q\) respectively. In the following, CCA is presented as a problem solved through an iterative algorithm. The first stage of CCA consists of finding two vectors \(a =(a_1,...,a_p)^t\) and \(b =(b_1,...,b_q)^t\) that maximize the correlation between the linear combinations \(U\) and \(V\) defined as
\[ \begin{aligned} U=Xa&=&a_1 X_1 +a_2 X_2 +\cdots a_p X_p\\ V=Yb&=&b_1 Y_1 +b_2 Y_2 +\cdots a_q Y_q \end{aligned} \]
and assuming that the vectors \(a\) and \(b\) are normalized so that \(\text{var}(U) = \text{var}(V) = 1\). In other words, the problem consists of finding \(a\) and \(b\) such that
\[ \rho_1 = \text{cor}(U, V) = \max_{a,b} \text{cor} (Xa, Yb)\quad \text{subject to}\quad \text{var}(Xa)=\text{var}(Yb) = 1. \tag{9.3}\]
The resulting variables \(U\) and \(V\) are called the first canonical variates and \(\rho_1\) is referred to as the first canonical correlation.
Note: Higher order canonical variates and canonical correlations can be found as a stepwise problem. For \(s = 1,...,p\), we can successively find positive correlations \(\rho_1 \geq \rho_2 \geq ... \geq \rho_p\) with corresponding vectors \((a^1, b^1), ..., (a^p, b^p)\), by maximizing
\[ \rho_s = \text{cor}(U^s,V^s) = \max_{a^s,b^s} \text{cor} (Xa^s,Yb^s)\quad \text{subject to}\quad \text{var}(Xa^s) = \text{var}(Yb^s)=1 \tag{9.4}\]
under the additional restrictions
\[ \text{cor}(U^s,U^t) = \text{cor}(V^s, V^t)=0 \quad\text{for}\quad 1 \leq t < s \leq p. \tag{9.5}\]
We can think of CCA as a generalization of PCA where the variance we maximize is the ‘covariance’ between the two matrices (see Holmes (2006) for more details).
When the number of variables in each table is very large finding two very correlated vectors can be too easy and unstable: we have too many degrees of freedom.
Then it is beneficial to add a penalty maintains the number of non-zero coefficients to a minimum. This approach is called sparse canonical correlation analysis (sparse CCA or sCCA), a method well-suited to both exploratory comparisons between samples and the identification of features with interesting covariation. We will use an implementation from the PMA package.
Here we study a dataset collected by Kashyap et al. (2013) with two tables. One is a contingency table of bacterial abundances and another an abundance table of metabolites. There are 12 samples, so \(n = 12\). The metabolite table has measurements on \(p = 637\) feature and the bacterial abundances had a total of $ q = 20,609$ OTUs, which we will filter down to around 200. We start by loading the data.
library("genefilter")
load("../data/microbe.rda")
metab = read.csv("../data/metabolites.csv", row.names = 1) |> as.matrix()We first filter down to bacteria and metabolites of interest, removing (“by hand”) those that are zero across many samples and giving an upper threshold of 50 to the large values. We transform the data to weaken the heavy tails.
library("phyloseq")
metab = metab[rowSums(metab == 0) <= 3, ]
microbe = prune_taxa(taxa_sums(microbe) > 4, microbe)
microbe = filter_taxa(microbe, filterfun(kOverA(3, 2)), TRUE)
metab = log(1 + metab, base = 10)
X = log(1 + as.matrix(otu_table(microbe)), base = 10)A second step in our preliminary analysis is to look if there is any association between the two matrices using the RV.test from the ade4 package:
colnames(metab) = colnames(X)
pca1 = dudi.pca(t(metab), scale = TRUE, scannf = FALSE)
pca2 = dudi.pca(t(X), scale = TRUE, scannf = FALSE)
rv1 = RV.rtest(pca1$tab, pca2$tab, 999)
rv1Monte-Carlo test
Call: RV.rtest(df1 = pca1$tab, df2 = pca2$tab, nrepet = 999)
Observation: 0.8400429
Based on 999 replicates
Simulated p-value: 0.002
Alternative hypothesis: greater
Std.Obs Expectation Variance
5.911826097 0.317788091 0.007804079 We can now apply sparse CCA. This method compares sets of features across high-dimensional data tables, where there may be more measured features than samples. In the process, it chooses a subset of available features that capture the most covariance – these are the features that reflect signals present across multiple tables. We then apply PCA to this selected subset of features. In this sense, we use sparse CCA as a screening procedure, rather than as an ordination method.
The implementation is below. The parameters penaltyx and penaltyz are sparsity penalties. Smaller values of penaltyx will result in fewer selected microbes, similarly penaltyz modulates the number of selected metabolites. We tune them manually to facilitate subsequent interpretation – we generally prefer more sparsity than the default parameters would provide.
library("PMA")
ccaRes = CCA(t(X), t(metab), penaltyx = 0.15, penaltyz = 0.15,
typex = "standard", typez = "standard")123456789ccaResCall: CCA(x = t(X), z = t(metab), typex = "standard", typez = "standard",
penaltyx = 0.15, penaltyz = 0.15)
Num non-zeros u's: 5
Num non-zeros v's: 16
Type of x: standard
Type of z: standard
Penalty for x: L1 bound is 0.15
Penalty for z: L1 bound is 0.15
Cor(Xu,Zv): 0.9904707With these parameters, 5 bacteria and 16 metabolites were selected based on their ability to explain covariation between tables. Further, these features result in a correlation of 0.99 between the two tables. We interpret this to mean that the microbial and metabolomic data reflect similar underlying signals, and that these signals can be approximated well by the selected features. Be wary of the correlation value, however, since the scores are far from the usual bivariate normal cloud. Further, note that it is possible that other subsets of features could explain the data just as well – sparse CCA has minimized redundancy across features, but makes no guarantee that these are the “true” features in any sense.
Nonetheless, we can still use these 21 features to compress information from the two tables without much loss. To relate the recovered metabolites and OTUs to characteristics of the samples on which they were measured, we use them as input to an ordinary PCA. We have omitted the code we used to generate Figure 9.32, we refer the reader to the online material accompanying the book or the workflow published in Callahan et al. (2016).
Figure 9.32 displays the PCA triplot, where we show different types of samples and the multidomain features (Metabolites and OTUs). This allows comparison across the measured samples – triangles for knockout and circles for wild type –and characterizes the influence the different features – diamonds with text labels. For example, we see that the main variation in the data is across PD and ST samples, which correspond to the different diets. Further, large values of 15 of the features are associated with ST status, while small values for 5 of them indicate PD status.
The advantage of the sparse CCA screening is now clear – we can display most of the variation across samples using a relatively simple plot, and can avoid plotting the hundreds of additional points that would be needed to display all of the features.
cca for their correspondence analyses function.The term constrained correspondence analysis translates the fact that this method is similar to a constrained regression. The method attempts to force the latent variables to be correlated with the environmental variables provided as `explanatory’.
CCpnA creates biplots where the positions of samples are determined by similarity in both species signatures and environmental characteristics. In contrast, principal components analysis or correspondence analysis only look at species signatures. More formally, it ensures that the resulting CCpnA directions lie in the span of the environmental variables. For thorough explanations see Braak (1985; Greenacre 2007).
This method can be run using the function ordinate in phyloseq. In order to use the covariates from the sample data, we provide an extra argument, specifying which of the features to consider.
Here, we take the data we denoised using dada2 in Chapter 4. We will see more details about creating the phyloseq object in Chapter 10. For the time being, we use the otu_table component containing a contingency table of counts for different taxa. We would like to compute the constrained correspondence analyses that explain the taxa abundances by the age and family relationship (both variables are contained in the sample_data slot of the ps1 object).
We would like to make two dimensional plots showing only using the four most abundant taxa (making the biplot easier to read):
ps1=readRDS("../data/ps1.rds")
ps1p=filter_taxa(ps1, function(x) sum(x) > 0, TRUE)
psCCpnA = ordinate(ps1p, "CCA",
formula = ps1p ~ ageBin + family_relationship)To access the positions for the biplot, we can use the scores function in the vegan. Further, to facilitate figure annotation, we also join the site scores with the environmental data in the sample_data slot. Of the 23 total taxonomic orders, we only explicitly annotate the four most abundant – this makes the biplot easier to read.
evalProp = 100 * psCCpnA$CCA$eig[1:2] / sum(psCCpnA$CA$eig)
ggplot() +
geom_point(data = sites,aes(x =CCA2, y =CCA1),shape =2,alpha=0.5) +
geom_point(data = species,aes(x =CCA2,y =CCA1,col = Order),size=1)+
geom_text_repel(data = dplyr::filter(species, CCA2 < (-2)),
aes(x = CCA2, y = CCA1, label = otu_id),
size = 2, segment.size = 0.1) +
facet_grid(. ~ ageBin) +
guides(col = guide_legend(override.aes = list(size = 2))) +
labs(x = sprintf("Axis2 [%s%% variance]", round(evalProp[2])),
y = sprintf("Axis1 [%s%% variance]", round(evalProp[1]))) +
scale_color_brewer(palette = "Set1") + theme(legend.position="bottom")
Figures 9.33 and 9.34 show the plots of these annotated scores, splitting sites by their age bin and litter membership, respectively. Note that to keep the appropriate aspect ratio in the presence of faceting, we have taken the vertical axis as our first canonical component. We have labeled individual bacteria that are outliers along the second CCpnA direction.
Evidently, the first CCpnA direction distinguishes between mice in the two main age bins. Circles on the left and right of the biplot represent bacteria that are characteristic of younger and older mice, respectively. The second CCpnA direction splits off the few mice in the oldest age group; it also partially distinguishes between the two litters. These samples low in the second CCpnA direction have more of the outlier bacteria than the others.
This CCpnA analysis supports the conclusion that the main difference between the microbiome communities of the different mice lies along the age axis. However, in situations where the influence of environmental variables is not so strong, CCA can have more power in detecting such associations. In general, it can be applied whenever it is desirable to incorporate supplemental data, but in a way that (1) is less aggressive than supervised methods, and (2) can use several environmental variables at once.
Heterogeneous data A mixture of many continuous and a few categorical variables can be handled by adding the categorical variables as supplementary information to the PCA. This is done by projecting the mean of all points in a group onto the map.
Using distances Relations between data objects can often be summarized as interpoint distances (whether distances between trees, images, graphs, or other complex objects).
Ordination A useful representation of these distances is available through a method similar to PCA called multidimensional scaling (MDS), otherwise known as PCoA (principal coordinate analysis). It can be helpful to think of the outcome of these analyses as uncovering latent variable. In the case of clustering the latent variables are categorical, in ordination they are latent variables like time or environmental gradients like distance to the water. This is why these methods are often called ordination.
Robust versions can be used when interpoint distances are wildly different. NMDS (nonmetric multidimensional scaling) aims to produce coordinates such that the order of the interpoint distances is respected as closely as possible.
Correspondence analysis: a method for computing low dimensional projections that explain dependencies in categorical data. It decomposes chisquare distance in much the same way that PCA decomposes variance. Correspondence analysis is usually the best way to follow up on a significant chisquare test. Once we have ascertained there are significant dependencies between different levels of categories, we can map them and interpret proximities on this map using plots and biplots.
Permutation test for distances Given two sets of distances between the same points, we can measure whether they are related using the Mantel permutation test.
Generalizations of variance and covariance When dealing with more than one matrix of measurements on the same data, we can generalize the notion of covariance and correlations to vectorial measurements of co-inertia.
Canonical correlation is a method for finding a few linear combinations of variables from each table that are as correlated as possible. When using this method on matrices with large numbers of variables, we use a regularized version with an L1 penalty that reduces the number of non-zero coefficients.
Interpretation of PCoA maps and nonlinear embeddings can also be enhanced the way we did for PCA using generalizations of the supplementary point method, see Trosset and Priebe (2008) or Bengio et al. (2004). We saw in Chapter 7 how we can project one categorical variable onto a PCA. The correspondence analysis framework actually allows us to mix several categorical variables in with any number of continuous variables. This is done through an extension called multiple correspondence analysis (MCA) whereby we can do the same analysis on a large number of binary categorical variables and obtain useful maps. The trick here will be to turn the continuous variables into categorical variables first. For extensive examples using R see for instance the book by Pagès (2016).
A simple extension to PCA that allows for nonlinear principal curve estimates instead of principal directions defined by eigenvectors was proposed in Hastie and Stuetzle (1989) and is available in the package princurve.
Finding curved subspaces containing a high density data for dimensions higher than \(1\) is now called manifold embedding and can be done through Laplacian eigenmaps (Belkin and Niyogi 2003), local linear embedding as in Roweis and Saul (2000) or using the isomap method (Tenenbaum et al. 2000). For textbooks covering nonlinear unsupervised learning methods see Hastie et al. (2008, chap. 14) or Izenman (2008).
A review of many multitable correlation coefficients, and analysis of applications can be found in Josse and Holmes (2016).
7 These values were chosen to give about retain about 3,000 taxa, similar to our previous choice of threshold.
See Figure 9.35.
IBDca = dudi.coa(ibd.pres, scannf = FALSE, nf = 4)
fviz_eig(IBDca, geom = "bar", bar_width = 0.7) +
ylab("Percentage of chisquare") + ggtitle("")
fviz(IBDca, element = "col", axes = c(1, 2), geom = "point",
habillage = day, palette = "Dark2", addEllipses = TRUE, color = day,
ellipse.type = "convex", alpha = 1, col.row.sup = "blue",
select = list(name = NULL, cos2 = NULL, contrib = NULL),
repel = TRUE)
See Figure 9.36. The code is not rendered here, but is shown in the document’s source file.
library("kernlab")
laplacedot1 = laplacedot(sigma = 1/3934)
rbfdot1 = rbfdot(sigma = (1/3934)^2 )
Klaplace_cellsn = kernelMatrix(laplacedot1, blom)
KGauss_cellsn = kernelMatrix(rbfdot1, blom)
Klaplace_rawcells = kernelMatrix(laplacedot1, moignard_raw)
KGauss_rawcells = kernelMatrix(rbfdot1, moignard_raw)Use kernelized distances to protect against outliers and allows discovery of non-linear components.
dist1kr = 1 - Klaplace_rawcells
dist2kr = 1 - KGauss_rawcells
dist1kn = 1 - Klaplace_cellsn
dist2kn = 1 - KGauss_cellsn
cells1.kcmds = cmdscale(dist1kr, k = 20, eig = TRUE)
cells2.kcmds = cmdscale(dist2kr, k = 20, eig = TRUE)
percentage = function(x, n = 4) round(100 * sum(x[seq_len(n)]) / sum(x[x>0]))
kperc1 = percentage(cells1.kcmds$eig)
kperc2 = percentage(cells2.kcmds$eig)
cellsn1.kcmds = cmdscale(dist1kn, k = 20, eig = TRUE)
cellsn2.kcmds = cmdscale(dist2kn, k = 20, eig = TRUE)colc = rowData(Moignard)$cellcol
library("scatterplot3d")
scatterplot3d(cellsn2.kcmds$points[, 1:3], color=colc, pch = 20,
xlab = "Axis k1", ylab = "Axis k2", zlab = "Axis k3", angle=15)
scatterplot3d(cellsn2.kcmds$points[, 1:3], color=colc, pch = 20,
xlab = "Axis k1", ylab = "Axis k2", zlab = "Axis k3", angle = -70)
Page built at 12:54 on 2026-05-18 using R version 4.6.0 (2026-04-24)