library(mvtnorm)
set.seed(123)
covariance <- matrix(c(5,3,0,-3,0, 3,5,0,-3,0, 0,0,5,0,0,
-3,-3,0,6,0, 0,0,0,0,3), nrow=5)
data <- rmvnorm(1000, sigma=covariance)
names(data) <- c("a", "b", "c", "d", "e")8 Unsupervised machine-learning (dimensionality reduction and clustering)
The following material should be classed as extension material (hence being in the appendix), and has been included simply because students have asked for it. This material is normally two days’ worth of lecures and practicals. You have been warned…
Overview
When faced with high-dimensional, continuous data—lots of columns of numbers in a spreadsheet—one of our first thoughts is often how we can focus on the most important dimensions (columns) of data. But if it’s not obvious which column is the most important, perhaps because they are all very similar to one-another, an alternative is to try and compress those many dimensions into fewer dimensions. This compression of data is sort of like finding the fundamental, underlying axes of variation that drive these data: the principal components of variation in the dataset. The three approaches you will learn in this section—PCA, PCoA, and NMDS—are all ways of carrying out such work. A major concern with such approaches is to make sure that you don’t throw the baby out with the bathwater: there is always a loss of information and so you must make sure that the information you’ve retained is still valuable for what you’re trying to understand. You must also be careful, when using these approaches, not to fall in love with them: never mistake the opaque for the profound. Ask yourself constantly and unceasingly what everything you’ve just calculated means. Sure, you’ve compressed 100 columns of data down into 2 columns, but what do those two columns represent? If you don’t know, then you can’t do anything meaningful with them. Check that your results make basic, logical sense in the context of what you know about the biology of the system you’re working with. If it does: congratulations! You might have done it right. Remember that the best way of reducing dimensionality is not with statistics, but with scientific insight. There’s no need to compress hundreds of columns of data if, from first principles, you know which column is the most important one anyway. Beware: these approaches assume that all axes of data are equally important, and so if only one of your columns of data was important in the biology of your system, then there is no guarantee that they will highlight that particular axis.
8.1 Principal Components Analysis—PCA
Principal Components Analysis (PCA) is a way of taking many columns of data and rotating them so as to find their principal axes of variation. It never changes the data, but rather highlights the main ways in which they co-vary. PCA is exactly the same as a major axis regression. PCA loses no information (i.e., it performs no true compression, only highlights correlations), but if you choose to work with only some of its results axes then you explicitly ignore the variation in those axes you do not choose.
8.1.1 An intuitive grasp of PCA
In a standard linear regression, you take a response variable (\(y\)), plot it against an explanatory variable (\(x\)), and find the line through the two that explains the most variation in \(y\) (i.e., that minimises the residual variance in \(y\)). Not simple at all, but something you’re familiar with. From this line, you’re able to get information about the relationship between the two variables (the slope of that line), and you’re able to get residuals—the distance of each \(y\) value from the line and a measure of the error in your prediction. You may recall that you can do the same thing with more than one explanatory variable, in which case the logic is much the same; you’re trying to find a line that minimises the distance between each \(y\) value and the line. If I’m describing this in unfamiliar terms, take a quick look at Figure 8.1.
Imagine now that this line you’re drawing isn’t the line that minimises the error in \(y\) (the distance between the line and each \(y\)), but instead is a line that minimises the total error in both \(y\) and \(x\). If you’re a visual person, you would agree from Figure 8.1 that such a line would minimise the length of the line in red that are at \(90^\circ\) to the line. This sort of regression has lots of names, and you might have heard of it as a ‘major axis regression’, a Deming regression, or … you may not have heard of it. But hopefully it makes a degree of intuitive sense.
If so, congratulations. You now understand a Principal Components Analysis, or PCA, because that is mathematically and precisely exactly what it is. At the end of this lecture you will go and tell your friends this, and someone will tell you that “no, he must be simplifying”: I’m not, this is all there is to it. In the case of two-dimensional data (\(y\) and \(x\) are two variables, so there are two dimensions), the first axis (or first principal component) of a PCA for each data point is how far along each piece of data (\((x,y)\) pair) is along that line. The second axis is how far that data point is from the line itself (a distance that is at \(90^\circ\) to the line—along the red lines in Figure 8.1). A PCA is calculated by trying to find the line that minimises the total variance of that second axis—or, put in a more familiar way, trying to draw the line that minimises the overall error of that line.
In the case that we had more than two variables (dimensions), that line would be still be drawn such that it would minimise the error through the entire dataset. The first axis would be defined in exactly the same way—minimising the error. The second axis would be at right angles—orthogonal—to that first axis (just as in the two dimensional case), and would extend in whatever direction the most variation remained in the data. The third, and final, axis, would be orthogonal to both the first and second axes, and would soak up all the remaining variation. This gives an intuitive reason for two important concepts in PCA: all the axes are orthogonal1 and all of the variation that was in your original data is contained in the new PCA axes.
We’re really just rotating the data in a fancy way so they match onto this line we’ve drawn, and some people prefer to describe PCA as a kind of rotation. In this case, you’re taking the \(x\) and \(y\) variables—the dimensions of the data— and then rotating around so we have a new first axis (the first principal component) and then a second orthogonal to that. This rotation could also happen for as many variables, or dimensions, as you have. Some people find thinking of this as a rotation more natural; if that’s you, then great, but if it’s confusing then focus on the error minimising argument (which is actually more general, and is more useful later in the course).
PCA is really useful, because it means we can take big datasets where everything is correlated with everything else, and collapse that variation down into a number of axes that are all independent—orthogonal—of each other. Because the amount of variance associated with each axis isn’t the same, we can pick fewer numbers of axes that explain “enough” of our data (more on that later), and analyse those. Since we’re just rotating the axes of our real data onto these new axes, we can also figure out how each PCA axis loads onto the original axes. Loading just means “how much of the variation in a given piece of data is associated with a particular PCA axis”. So a variable that loads strongly and positively onto PCA axis 1 is positively correlated with that new axis—the PCA “line” goes through it with a positive slope. A weak negative loading would be a shallow curve with a negative slope. Bear in mind that the sign of the loadings is essentially arbitrary: they only make sense relative to one-another. Thus it’s not important that \(x\) is negatively loaded only PC1, for example, but it would be important if \(x\) were negatively loaded and \(z\) were positively loaded, because it implies that the contributions of \(x\) and \(z\) to the axis are negatively correlated.
8.1.2 Hands on with PCA
OK, let’s simulate some data that we can use to figure out how PCA works. This is probably the scariest part of today’s practical, because we’re going to be drawing data from a multivariate Normal distribution. In practice, what this means is we’re going to draw five variables at random: three of them are going to be strongly correlated with one-another, and the other two are going to vary independently. Take a look at the covariance matrix (covariance): each row and column refers to one of the variables we are going to be randomly drawing from a Normal distribution, and represents the covariance of the two variables. So, for example, row 1 column 2 (3) means “variables 1 and 2 will have a covariance of 3”. Row 1 column 1 (5) means “this variable will have a variance of 5”—covariance with yourself is really just your overall random variation. In cases where the covariance is 0, that means those two variables won’t correlate with each other—they’re independently drawn. The line containing set.seed is just to make sure that the random numbers you generate are the same as the ones I generate (don’t worry about how this works)—it means I can talk about the PCA you’ve generated and be certain that it’s the same as the one I’m working with.2
Notice how I load an R library called mvtnorm using the line library(mvtnorm) at the beginning of this script. If R gives you an error when running this, that means you don’t have that library installed: you can either install it using RStudio’s menus, or you can run the line install.packages("mvtnorm"). Either is fine, but you’ll then have to load the library (installing doesn’t load something in R, just like other programs on your computer). I’m not going to remind you to do that in the future: I will just assume that, if you can’t load a library, you will install it.
Now we’ve got this data, we can run a principal components analysis (PCA) on it. The R code is deceptively simple: the first line below runs the PCA, and the second produces what’s called a biplot. In it we see each point in our data (each data point is labelled according to its row in the data), with arrows that show how each variable (a, b, c, d, and e) loads onto each axis. So you can see that d is negatively correlated with a and b, because they’re both in the same ‘plane’ (i.e., PC1) but on opposite ends. You can also see that c is floating out on its own, orthogonal to everything else, because it was drawn independently from the others. If you were to look at PC axes 2 and 3 (which you can do using the argument, see below) you would see e floating around independently of all the other axes too. It’s very hard to even see e in this particular plot, because it doesn’t load onto either PC1 or PC2.
pca <- prcomp(data)
biplot(pca)
biplot(pca, choices=2:3)
How can we see these ‘loadings’ that I keep talking about? Quite easily: you just print out the PCA as shown below. This gives you the rotations of the original data (rows are original data, columns are PC axes), which shows you how strongly each variable loads onto each principal component. So a, b, and d all load strongly onto PC1 (in opposite directions, hence the sign changes), while c and e are doing their own things on PCs 2 and 3 (notice how big their loadings are—the maximum possible is 1, and they’re dominating the values on those axes). If you want to prove to yourself that, really, PC2 is essentially exactly the same as c then run the second line below, which plots c against PC2. Remember: loadings vary from \(-1\) to \(1\), where \(-1\) means a very strong (perfect) negative correlation between a PC axis and the variable, \(1\) means a perfect positive correlation, and \(0\) means no correlation whatsoever. Think of them as correlations and you’ll do fine.
pca
#> Standard deviations (1, .., p=5):
#> [1] 3.401595 2.212150 1.712135 1.577754 1.423393
#>
#> Rotation (n x k) = (5 x 5):
#> PC1 PC2 PC3 PC4 PC5
#> [1,] 0.574036850 0.014170577 -0.17667112 -0.63925614 -0.48002062
#> [2,] 0.549149941 0.007550034 -0.17999897 -0.06764815 0.81326591
#> [3,] 0.019799941 0.998321195 0.03729037 0.03802989 -0.01122095
#> [4,] -0.606999343 0.052532745 -0.31559566 -0.66980531 0.28381744
#> [5,] 0.008706181 -0.018365568 0.91400027 -0.36971596 0.16583280
plot(pca$x[,2] ~ data[,3], xlab="'c' variable", ylab="PC2")
PCA is meant to be a kind of dimension reduction technique: we’re trying to take a load of correlated data and reduce it down to a smaller set of independent, or orthogonal, principal components. So two critical questions are “how important is each axis” and “how many axes are enough for me to reduce everything down to”? PCA is pretty good at giving us an answer to the first question, but sort of rubbish at the second. In PCA-land, an axis is important if it explains a lot of the variance (variation) in the original data. If we call summary on our PCA object, we can find out how much standard deviation (the square root of variance) is associated with each axis, which sort of gets at that first question. A bigger SD is a more important axis. As for how much is enough… Well, most people plot a scree plot (the second line below) and find where there’s an inflection point in the standard deviations: you find the axes that you can join in a single straight line, and work with those axes plus the next one along. So, in this case, PCs 1 and 2 are on a straight line, PC3 isn’t, so we’d go with PCs 1–3. There is a citation associated with this3, but there’s little more justification to this idea than “it seems a good idea”. I don’t have a definitive answer for you as to how to pick your number of PCA axes, but I will give you two pieces of advice: (1) just because an axis doesn’t have a lot of variance associated with it doesn’t mean it’s not important. The amount of variance is related to all the variables you put into the analysis: if only one of your variables actually mattered to the biology of your system, a PCA isn’t going to help you. The whole point of PCA is that each variable could be as important as the other, and so if you’ve put a load of crap into your PCA then you’re probably going to get a load of crap out and miss the important stuff. (2) As Mick Crawley infamously advises: “never mistake the opaque for the profound”. If you don’t know what a PCA axis represents biologically, it’s essentially useless to you. It’s useful to have a single variable that captures how hot and dry a particular region of environmental space is; it’s not useful to have a single variable that maps onto thirty different climate variables, none of which you understand. User beware…
summary(pca)
#> Importance of components:
#> PC1 PC2 PC3 PC4 PC5
#> Standard deviation 3.4016 2.2121 1.7121 1.5778 1.42339
#> Proportion of Variance 0.4839 0.2047 0.1226 0.1041 0.08473
#> Cumulative Proportion 0.4839 0.6886 0.8112 0.9153 1.00000
plot(pca)
Finally, PCAs are quite sensitive to variance among different axes. This isn’t so much of a problem in this example, but because the whole thing is based around finding axes that explain the most variance, if the variance isn’t equal across your variables you’re going to get skewed results. This is essentially a corollary of piece-of-advice 1 above: if everything’s meant to be equally important, they should all vary equally, or the thing that varies the most will be ‘inflated’ in importance. There’s a simple solution: conduct your PCA on z-transformed variables (scaled to have a standard deviation of 1), and the option for that is simple to use (see below). If you don’t use a scaled version of your variables you will make a mistake: this isn’t an opinion thing, it’s a fact thing. Scale your variables.
pca <- prcomp(data, scale=TRUE)8.2 Principle Co-ordinate Analysis (PCoA) and Non-Metric Multidimensional Scaling (NMDS)
PCoA is a close relative of PCA; whereas PCA works with the underlying data, PCoA works on a distance matrix of the data. This allows PCoA to be applicable to any kind of data you can create a distance matrix for (including categorical data), and in the case that you’re working with the covariance matrix of such data you get exactly the same result as if you’d computed a PCA on that data. PCoA projects a distance matrix onto a Euclidean space of as many dimensions as you specify: it’s up to you to see whether you’re summarised that distance matrix well, and whether it’s even possible to approximate the distance matrix well in Euclidean space. Non-metric Multi-Dimensional Scaling (NMDS) is a popular close-cousin, that differs from PCoA in being able to work with the ranks of distances among data-points. The mathematics of it are, in many ways, simpler, but getting it right is much harder. As always, the choice is yours…
8.2.1 Under pressure: an intuitive feeling for PCoA and NMDS
Principal Co-ordinate Analysis (PCoA) is a method that takes a distance matrix and projects it into Euclidean space. Euclidean space just means “distances among objects work in the way that they do in the real world”: there are lots of ways mathematicians have invented of representing objects in the world that don’t make intuitive sense (e.g., space itself is curved, not straight). So PCoA takes a distance matrix and finds a way of arranging the points that each element of the distance matrix represents so that they can be plotted. If they can be plotted in two dimensions, great, but the data could need as many as three, four, five, six… It proceeds in three steps:
- Centre the distance matrix. This means that, when we map everything into Euclidean space, the centroid of the data will be at the origin (\((0,0)\)) of a graph at the end.
- Find the eigenvalues and eigenvectors of the distance matrix. I will discuss this step in more detail below.
- Scale everything according to the square root of the eigenvalue. This is important, because anything with a negative eigenvalue can’t be represented in Euclidean space (see below), and anything with an eigenvalue of zero doesn’t matter (see below). So the purpose of this step is simply to figure out what what matters.
- Finished! Each row of the resulting matrix is a data-point, and each column the position in space where that point should be plotted.
Bad news: to do PCoA we really do have to understand what eigenvectors and eigenvalues are. Well, they’re the things we were calculating in our PCA: the eigenvectors are regression equations through our data that are orthogonal to each other. The eigenvalues are how much variation in our original data is associated with each of these equations—they’re identical to the standard deviations we worked with during the PCA session. In an ecological context, you can map different environmental variables onto these axes as well, in much the same way that you can look at loadings of variables in a PCA.
The magic of PCoA is that some of these eigenvalues become zero when we run all of this through. That’s really cool, because it means the whole problem of figuring out how many axes are “significant” is done for us: all of them are “significant”4, and so that’s that5. There is, however, a major catch: not all of the variation in the original data can be mapped onto Euclidean space, and there is nothing PCoA can do about that. There are lots of ways of measuring the extent of that missing variation, and they all have names like stress or strain, because they reflect the extent to which your data are resisting being mapped onto this space.
If this loss of information is something that really worries you, it’s possible to use something called Non-Metric Multidimensional Scaling (NMDS; sometimes just called MDS), which sacrifices statistical pleasantry for guaranteeing that it can compress the information in your data into a set number of dimensions. The problem with NMDS is that it requires you to pick a number of dimensions onto which you will project your data, and how do you pick the number of dimensions? Well, there’s no single answer, but you have to check the strain that your data is put under when you map them onto a set number of dimensions, and there’s no clear way to pick that. If this doesn’t sound like a solution, that’s because (in my opinion), frankly it’s not. Either way, you have to decide how much variation you’re happy throwing out. Happily, there is some light at the end of the tunnel: quantile regression lets you work with the original distance matrix itself for statistical analysis, and so all of these problems become purely visualisation issues.
8.2.2 What’s an extra ‘o’ among friends? Hands on with PCoA
OK, let’s get our hands dirty with PCoA to get an intuitive feel for what’s going on. As ever, we’re going to simulate some data, but this time it’s going to be some community ecological data: compositions of species across two environmental gradients. This is the most complex set of data I’m going to simulate for you in this class: we’re going to simulate intercepts and slopes of species’ responses to two environmental gradients, then simulate their abundances across those gradients. If it doesn’t make much sense, just treat it as magic: the details don’t matter, just take the story for granted. There are simpler ways to simulate data such as this, but if you’re interested in the distribution of diversity you might find this a useful starting point for your own investigations.
# Draw species' parameters
intercepts <- rnorm(20, mean=20)
env1 <- rnorm(20)
env2 <- rnorm(20)
# Create environment
env <- expand.grid(env1=seq(-3,3,.5), env2=seq(-3,3,.5))
biomass <- matrix(ncol=20, nrow=nrow(env))
for(i in seq_len(nrow(biomass)))
biomass[i,] <- intercepts + env1*env[i,1] + env2*env[i,2]OK, so we’ve got a lot of data. Notice these are biomasses, thus we’ve got lots of continuous data6. Let’s summarise it using a PCoA:
library(vegan)
#> Loading required package: permute
dist <- dist(biomass)
pcoa <- cmdscale(dist, eig=TRUE)
barplot(pcoa$eig)
The second line creates a distance matrix, as you know. If you’re an ecologist, you’re probably already thinking that it would be a bad idea to work with a Euclidean distance, in which case vegan’s vegdist has got you covered. The third line runs the PCoA, and by asking it to return the eigenvector information we’re guaranteeing that we get all the PCoA axes and lots of additional information. The fourth line plots the eigenvalues of each axis: how much variation they explain in the space. I hope you would agree with me that the first two axes explain the overwhelming majority of information in this space, so we’re fine working with just those. Note that the “scree test” I taught you for PCA does not apply here: that test is based around expectations of eigenvalues under various null models for PCAs, not for PCoAs.
plot(pcoa$points[,1:2], xlab="PCoA1", ylab="PCoA2")
plot(pcoa$points[,1:2], type="n", xlab="PCoA1", ylab="PCoA2")
text(pcoa$points[,1:2]+.25, labels=env[,1], col="red")
text(pcoa$points[,1:2]-.25, labels=env[,2], col="black")
The code above shows how you can plot the various axes (points) of the PCoA. The first line shows how you would plot each site’s position in the two-dimensional space that we have mapped our distance matrix onto; the next three show the values of each site in the first (in red) and second (in black) environmental axes. All I want you to take from this is that we simulated data that varied across two axes, and as a result our PCoA has found two axes across which it can distribute the data.
I’m not going to dwell on two things that are normally mentioned in screeds on PCoA at this point: stress and loadings. Stress, which you’ll recall is how well our distance matrix can be mapped onto a certain number of axes, is manifest in PCoA through negative eigenvalues7. If you want to see an example where there is some stress in a plot, run example(cmdscale) in R. Loadings give you what biplots gave us when we looked at PCAs: they show how variables (e.g., environmental conditions) load onto each of the dimensions the PCoA has estimated. Remember that calculating the PCoA of a distance matrix is the same as running a PCA on the covariances of the same data, and so if you understand biplots in PCA you essentially understand it in PCoA. I describe, in the section on quantile regression, why I think using this sort of thing with a PCoA is a bad idea, however, and so I am not going to spend any more time explaining something I don’t really think you should be doing8.
8.2.3 What’s an extra ‘n’ among friends? Hands on with (N)MDS
You don’t need to have an analytical solution for how to find the Euclidean space onto which a given distance matrix should be mapped. You can, instead, keep moving points around until you find an arrangement where the distance between points in that space is sufficiently close to the distance in the real data. That difference, when scaled in a few different ways, is either called stress or strain depending on the precise method you’re using. The result is quite dependent on the algorithm used to find a good arrangement of points in space, and (you guessed it) there’s no single way to know how many dimensions you should map the data onto or what algorithm to use. With that in mind… Let’s give this method, which is called NMDS, a try.
nmds <- metaMDS(dist)
#> Run 0 stress 1.464206e-16
#> Run 1 stress 3.641244e-05
#> ... Procrustes: rmse 2.710691e-05 max resid 7.139872e-05
#> ... Similar to previous best
#> Run 2 stress 8.417549e-05
#> ... Procrustes: rmse 3.055831e-05 max resid 0.0002167998
#> ... Similar to previous best
#> Run 3 stress 9.234905e-05
#> ... Procrustes: rmse 4.645507e-05 max resid 0.0002139585
#> ... Similar to previous best
#> Run 4 stress 6.740245e-05
#> ... Procrustes: rmse 3.446793e-05 max resid 0.0001473868
#> ... Similar to previous best
#> Run 5 stress 5.119249e-05
#> ... Procrustes: rmse 3.098223e-05 max resid 0.0001082465
#> ... Similar to previous best
#> Run 6 stress 7.36929e-05
#> ... Procrustes: rmse 4.783841e-05 max resid 0.0001722293
#> ... Similar to previous best
#> Run 7 stress 6.053695e-05
#> ... Procrustes: rmse 4.451297e-05 max resid 0.0001591233
#> ... Similar to previous best
#> Run 8 stress 7.950937e-05
#> ... Procrustes: rmse 5.056294e-05 max resid 0.0001607971
#> ... Similar to previous best
#> Run 9 stress 7.856067e-05
#> ... Procrustes: rmse 4.324369e-05 max resid 0.0001565351
#> ... Similar to previous best
#> Run 10 stress 9.298364e-05
#> ... Procrustes: rmse 5.588033e-05 max resid 0.0001933552
#> ... Similar to previous best
#> Run 11 stress 8.298686e-05
#> ... Procrustes: rmse 5.098683e-05 max resid 0.000172562
#> ... Similar to previous best
#> Run 12 stress 9.835348e-05
#> ... Procrustes: rmse 5.855983e-05 max resid 0.0002306781
#> ... Similar to previous best
#> Run 13 stress 9.332286e-05
#> ... Procrustes: rmse 5.565545e-05 max resid 0.000204979
#> ... Similar to previous best
#> Run 14 stress 6.517608e-05
#> ... Procrustes: rmse 4.515782e-05 max resid 0.0002339777
#> ... Similar to previous best
#> Run 15 stress 8.104738e-05
#> ... Procrustes: rmse 3.975073e-05 max resid 0.0001494012
#> ... Similar to previous best
#> Run 16 stress 9.557778e-05
#> ... Procrustes: rmse 5.98405e-05 max resid 0.000199285
#> ... Similar to previous best
#> Run 17 stress 5.080436e-05
#> ... Procrustes: rmse 3.941986e-05 max resid 0.0001529907
#> ... Similar to previous best
#> Run 18 stress 4.887195e-05
#> ... Procrustes: rmse 3.565261e-05 max resid 0.0001328035
#> ... Similar to previous best
#> Run 19 stress 9.518741e-05
#> ... Procrustes: rmse 5.572656e-05 max resid 0.0001927752
#> ... Similar to previous best
#> Run 20 stress 5.038397e-05
#> ... Procrustes: rmse 3.392541e-05 max resid 0.0001308315
#> ... Similar to previous best
#> *** Best solution repeated 20 times
#> Warning in metaMDS(dist): stress is (nearly) zero: you may have insufficient
#> data
plot(nmds)
#> species scores not available
stressplot(nmds)
No great surprise, the answer is essentially the same as before. You will likely get a warning that the stress is very low: this is because these data were simulated to be very well-behaved, and is nothing to worry about. The first two lines are reasonably self-explanatory, but the third line is a Shepard diagram, and can be used to diagnose how well your model has fit. You ideally want a straight 1:1 line, and if you can’t get that then you want the red dots to be close to the blue line. The line is broken up into steps because during NMDS fitting the distances are monotonically transformed to maximise the relationship between the original distances and the new distances. This isn’t a parametric fit—hence the line is broken, as we don’t know what shape it should be so this is the most conservative—and the line can’t decrease—it’s a monotonic fit.
plot(nmds)
#> species scores not available
orditorp(nmds, display="sites")
ordisurf(nmds, env[,1])
#> Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
#> : Fitting terminated with step failure - check results carefully
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
#>
#> Estimated degrees of freedom:
#> 8.7 total = 9.7
#>
#> REML score: -5417.331I’ve already described that I don’t like mapping environmental responses onto these kinds of plots, but the code above will do it for you. All I can say, repeatedly, is that you should be tremendously careful about fitting a smoothed regression to transformed-distance data, because there are so many steps and so many ways for error to creep in (how much variation loads onto each NMDS axis? How well does the smoothed regression fit? etc.). Regardless, above is how you would do it. It’s also possible to fit convex hulls around groups of points that you think correspond to particular treatments: I refuse to show you how to do this because they’re so over-used at conferences, but there are plenty of guides online9
8.3 Coda: quantile regression
NMDS and PCoA are fantastic ways to summarise data, but because you lose information in the process I don’t think they’re good foundations for analysis. I am not as militant about this as I have probably led you to believe over the last few pages, and every analytical tool has its place, but overall I prefer to work with the original distance matrix itself when performing statistics—particularly when the questions related to actual distance matrices themselves. I think this is particular important for analyses of \(\beta\)-diversity, and while NMDS is commonly used to understand how microbial sequencing data varies across time and space, I am not a fan of using those NMDS axes of a distance matrix in analysis. So what do I like doing instead?
mantel(dist, dist(env[,1]))
#>
#> Mantel statistic based on Pearson's product-moment correlation
#>
#> Call:
#> mantel(xdis = dist, ydis = dist(env[, 1]))
#>
#> Mantel statistic r: 0.5224
#> Significance: 0.001
#>
#> Upper quantiles of permutations (null model):
#> 90% 95% 97.5% 99%
#> 0.0238 0.0317 0.0368 0.0420
#> Permutation: free
#> Number of permutations: 999
mantel(dist, dist(env))
#>
#> Mantel statistic based on Pearson's product-moment correlation
#>
#> Call:
#> mantel(xdis = dist, ydis = dist(env))
#>
#> Mantel statistic r: 0.9774
#> Significance: 0.001
#>
#> Upper quantiles of permutations (null model):
#> 90% 95% 97.5% 99%
#> 0.0229 0.0318 0.0372 0.0475
#> Permutation: free
#> Number of permutations: 999A Mantel test is the equivalent of a correlation coefficient for (distance) matrices. There are lots of caveats to its use (read Legendre et al., 2015; Methods in Ecology and Evolution, 6(11), 1239-1247), but if you understand a correlation coefficient then you intuitively understand it. It is rare that you will go wrong with it. In the first line above, I run a Mantel test on a single variable, and in the second I run it on a distance matrix constructed from all the environmental data merged into one matrix. There are lots of version of Mantel tests: I only use one, the basic form of it, and it is perfect for when you want to see if two distance matrices are correlated. I would advise you to ignore things like ‘partial’ Mantel tests, if you come across them, and instead use quantile regressions (see below).
library(quantreg)
#> Loading required package: SparseM
data <- data.frame(dist=as.numeric(dist),
env1=as.numeric(dist(env[,1])), env2=as.numeric(dist(env[,2])))
model <- rq(dist ~ env1 + env2, data=data)
summary(model)
#>
#> Call: rq(formula = dist ~ env1 + env2, data = data)
#>
#> tau: [1] 0.5
#>
#> Coefficients:
#> Value Std. Error t value Pr(>|t|)
#> (Intercept) 0.72399 0.03689 19.62635 0.00000
#> env1 2.40150 0.01017 236.02986 0.00000
#> env2 3.98079 0.01070 371.97971 0.00000If you have more than one variable and you want to know the influence of each on your variable of interest, use a quantile regression. A quantile regression is exactly like a standard regression, only it works on the quantiles of your data—by default, the median. Magically, as long as the pseudo-replication in your data is equal, all its estimates are trustworthy: this is math-speak for saying “it’s fine to use with distance matrices”. If you have more than one explanatory variable, this is the method you should use to understand what is driving pattern in your data.
complex.model <- rq(dist ~ env1 + env2, data=data, tau=c(.2,.5,.8))
summary(complex.model)
#> Warning in summary.rq(xi, U = U, ...): 299 non-positive fis
#>
#> Call: rq(formula = dist ~ env1 + env2, tau = c(0.2, 0.5, 0.8), data = data)
#>
#> tau: [1] 0.2
#>
#> Coefficients:
#> Value Std. Error t value Pr(>|t|)
#> (Intercept) 0.17961 0.00606 29.65255 0.00000
#> env1 2.34742 0.00260 902.97835 0.00000
#> env2 3.98565 0.00437 911.65747 0.00000
#>
#> Call: rq(formula = dist ~ env1 + env2, tau = c(0.2, 0.5, 0.8), data = data)
#>
#> tau: [1] 0.5
#>
#> Coefficients:
#> Value Std. Error t value Pr(>|t|)
#> (Intercept) 0.72399 0.03689 19.62635 0.00000
#> env1 2.40150 0.01017 236.02986 0.00000
#> env2 3.98079 0.01070 371.97971 0.00000
#>
#> Call: rq(formula = dist ~ env1 + env2, tau = c(0.2, 0.5, 0.8), data = data)
#>
#> tau: [1] 0.8
#>
#> Coefficients:
#> Value Std. Error t value Pr(>|t|)
#> (Intercept) 1.72617 0.02797 61.71595 0.00000
#> env1 2.67556 0.01839 145.48183 0.00000
#> env2 4.05886 0.01935 209.71979 0.00000Quantile regression can be run on quantiles other than the median, by varying the \(\tau\) parameter (‘tau’), as I show you above. This makes them particularly cool if you think things are going on in the extremes of your data: I use them all the time, and find them to be very useful in ecology in general, not just distance matrices10.
Let me try to convince you, one last time, to not do what everyone else does and draw circles around groups of species or sites in an NMDS/PCoA plot and write stories about that. PCoA and NMDS are ways of visualising complex data, and as they simplify it they may remove signal from that data. Any relationships or groupings you draw on a PCoA or NMDS plot come from the distance matrix of the data: any story you can tell about those groupings is fundamentally a property of the distances, which you have transformed. Quantile regression works with those original distance matrices, and gives you regression plots based around those original values. Why use something that (1) throws away variation in your data and (2) relies on a transformation of your data into another form that doesn’t represent it? Why not, instead, just plot your data?
simple.eh <- rq(dist ~ env1, data=data)
with(data, plot(dist ~ env1, pch=20,
xlab="Environmental distance", ylab="Site similarity"))
abline(simple.eh, col="red", lwd=3)
8.4 Hierarcihcal cluster analysis
The distinction between dimensionality-reduction and clustering approaches can sometimes be hard to understand, because they both involve simplifying data. Fundamentally, clustering methods focus on detecting groups within data. Whereas something like a PCA tries to summarize lots of continuous variable into fewer variables, cluster analysis tries to detect discrete groupings in that data. Which you end up using depends on what you think is going on with your data: if you’re interested in species delimitation, you’ll probably prefer cluster analysis, for example. Clusters are also easy to present to non-scientist stakeholders (people like to put things into boxes, in my experience), but you must be careful to emphasize the uncertainty associated with clusters. While I’ll give you some examples of some of the techniques that can be used to delimit clusters, the truth is that there is no single way to pick the “correct” numbers of clusters. So be careful, and use your understanding of the biology of the system to drive your approach!
Hierarchical cluster analysis is all about splitting your data into groups. If you have a spreadsheet of data, it will split that data into as many groups as you have rows of data. That doesn’t sound very useful, but the trick is that it makes those groups in a hierarchy—a series of splits—that allow you to stop splitting the data whenever you want, leaving you with as many or as few groups as you like. That leaves, of course, the rather knotty question of where you should stop splitting, but in practice when you look at the output of a hierarchical cluster analysis that’s pretty obvious from eyeballing the data. If it isn’t… then you probably want another method, and we’ll get to that.
8.4.1 Let me break it down for you: an intuitive overview
The aim of hierarchical cluster analysis is to separate your data into a dendrogram—a tree—where each data-point (row) is a tip on the tree. Starting at the root of the tree, each split in the tree takes you closer and closer to data, and as you go further along the tree all the data subtending (falling off from) each node should be increasingly similar. Building such a dendrogram requires the original data to be transformed into a distance matrix: essentially a table where all the rows and columns represent data points in the original data, and all the entries in the table are the distances among data-points. There is an alternative, called model-based cluster analysis, which doesn’t rely on a distance matrix, but we won’t be covering that today.
The devil is in the detail with any cluster analysis, and the most important detail is how you construct the distance matrix. The simplest is to use the Euclidean distance among points: this is just the distance between points if they were plotted according to whatever variables you’ve collected. So, for example, if you had data on how long and how wide a series of boxes were, then the Euclidean distance of these boxes would just be the distances among all the points when plotted in a scatter-plot. The calculation is much the same if you have three variables, or even thirty, but of course you can’t visualize the plot if there are thirty variables! There are lots of different ways to measure distances, and the best one to use boils down to how your data are distributed, but you rarely go too far wrong with Euclidean distance if your data are Normally distributed. If you’re dealing with discrete data, something like the Hamming or binary distance is often useful, which is literally just the count of the number of ways that two things are different (e.g., “Will” vs. “Bill” would have a Hamming distance of 1, because “W” != “B”).
Once you have your distance matrix, you can proceed to form your dendrogram. You start from the ‘bottom’ of the dendrogram, and find the two pieces of data that are most similar to each other. They then form a cluster, and the distance matrix is updated to reflect every other point’s distance from that cluster. You then find the data point or cluster that is closest to another row in the distance matrix, replace that with a cluster, re-calculate, and keep going until everything has been grouped together into one big cluster. By keeping track of all the times you joined things together in the cluster analysis, you’re able to draw a dendrogram—the dendrogram is just a plot of when each data-point joined up with another point, with the lengths of the branches in that tree proportional to the distance between that group/data-point and the group/data-point it joined. It makes a lot more sense once you see it plotted, I promise!
The final complication is how you update the distance matrix when you’re changing everything around. If you replace a group with the average of all its members’ distances, you’re doing an average or Unweighted Pair Group Method with Arithmetic Mean (UPGMA) clustering. This is a popular method if you’re an evolutionary biologist because, under the assumption of a molecular clock, it will produce dendrograms you might reasonably call phylogenies11. The default in R is complete linkage clustering, which replaces a group with its maximum distance to each other data-point/data-point in a group. I’m afraid there’s no “best” method to use; just take a look at what you get and decide for yourself what makes the most sense.
8.4.2 Hands-on with hierarchical cluster analysis
To get started, let’s simulate some data that fall into three fairly obvious groups. As always, if the method of simulation doesn’t make much sense to you, don’t worry about it. Focus on the plot at the end, and make sure you can see that there are roughly three groups.
data <- data.frame(rbind(
cbind(rnorm(50),rnorm(50)),
cbind(rnorm(50,5), rnorm(50,5)),
cbind(rnorm(50,-5), rnorm(50,-5))
))
data$groups <- rep(c("red","blue","black"), each=50)
names(data)[1:2] <- c("x","y")
with(data, plot(y ~ x, pch=20, col=groups))
Great, now we’ve got some data! So, let’s create a distance matrix, using R’s default of a Euclidean distance matrix because we’re dealing with well-behaved Normally-distributed data.
distance <- dist(data[,c("x","y")])
upgma <- hclust(distance, method="average")
plot(upgma)
comp.link <- hclust(distance)
plot(comp.link)
You’d probably agree, looking at this plot, that there are three main groups. Despite being built from the ‘bottom-up’, it is best to read hierarchical dendrograms from the top-down. Starting at the very top, you can see that there are three main splits, and then at the very bottom lots of sub-groups. There are three groups that seem quite distant from each other: look at the height and how there are three splits very high in the tree, representing groups that, when joined together in the cluster, have a very big distance from each other. ‘Lower’ in the tree, data-points within the clusters have lesser distances from each other and so join together at lower heights. You can also see that, for all the fuss that’s made about it, in this case both the UPGMA and the complete linkage methods are forming similar numbers of clusters.
But what if we want to cut our data into groups? Well, luckily R has us covered and it’s quite simple to cut the dendrogram at a certain point or into a certain number of groups.
cut.by.groups <- cutree(upgma, k=2)
cut.by.height <- cutree(upgma, h=8)8.4.3 How many clusters?
The question everyone asks at this point is “but how many groups” and/or “but at what height” should you be doing your cutting? The answer is, surprisingly, quite difficult to get definitively12. The problem can be thought of as a ratio problem: what should the ratio of the average distance among clusters to the average distance within clusters be? Any choice of cut height will determine this ratio, as will a choice in the number of clusters (they are, essentially, the same thing). You might agree that lower intra-cluster distance and higher inter-cluster distance would be the best choice, but in practice it’s difficult to find where the best ‘stopping point’ is (i.e., what ratio of inter- to intra-cluster distance is ‘correct’).
In my experience, one of the best approaches is the DD-weighted gap statistic, as defined by Yan & Ye (2007; Determining the number of clusters using the weighted gap statistic; Biometrics 63: 1031–1037). This metric is appropriate for hierarchical clustering and k-means clustering (more on that next time), and although it tends to work best with spherical clusters so does everything else I’ve seen that was worthwhile. It’s a slightly non-standard install, but I promise you it’s worth the effort:
#install.packages("paran")
#install.packages("ape")
#install.packages("splits", repos="http://R-Forge.R-project.org")
library(splits)
#> Loading required package: ape
#> Loading required package: MASS
#> Loading required package: paran
gap.stat <- ddwtGap(data[,c("x","y")])
#> 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
with(gap.stat, plot(colMeans(DDwGap),pch=15,type='b',
ylim=extendrange(colMeans(DDwGap),f=0.2),
xlab="Number of Clusters", ylab="Weighted Gap Statistic"))
gap.stat$mnGhatWG
#> [1] 3…the take-home from this is the DD-gap-statistic (as you can see from the plot) is at its maximum when we assume there are three clusters in the data, which there were. The DD-gap-statistic (I always think of the “DD” as standing for “double-differenced”) compares the average distances among points within a cluster for a given number of clusters with the distances if there were one fewer, or one more, clusters than a particular value. What this means, in practice, is that when you hit the ‘correct’ number of clusters the DD-gap-statistic peaks, and so to pick your number of clusters you find the peak and that’s your number. This is also statistically testable via a parametric bootstrap approach based vaguely on the distribution of your data (see the help file for details). Or, if you’re feeling lazy, you can just print out that final value (mnGhatWG), which is the “answer” according to this method: 3.
I’m afraid that, fundamentally, there is no universally-accepted definition of what a “good” separation of data looks like, so your choice is either to look at the plot and pick for yourself, or to trust something that someone has told you is quite good and that you like. So your choices are to pick for yourself, read the paper I reference above and see if you agree with me that it’s good, or to just trust me13.
8.5 K-means clustering
K-means is, intuitively, I think much easier to grasp than hierarchical clustering, and it is blindingly fast. However, it also happens to be the method that requires the most assumptions going in. All clustering methods require care in the selection of number of clusters, but this is perhaps most apparently in k-means.
8.5.1 NP-hard? I don’t know the meaning of the word
k-means clustering is so simple that it’s probably the only algorithm in this course where I’m not going to simplify the explanation14. It proceeds in steps:
- Pick as many centers—centroids—at random as you are trying to find clusters in your data. Note that there are many ways of randomly choosing those centroids, and some would argue you shouldn’t do it randomly.
- Assign each data point to the cluster whose centroid it is closest to.
- Calculate the new centroid (center/mean) of all the points in each cluster.
- Go back to step 2 and keep running through and repeating until the assignments of your data to clusters no longer change in step 3.
Amazingly, this often works quite well. If you restart the search lots of times and pick the “best” outcome, where best is defined as the run with the smallest within-group variance, it works even better. Leaving aside the problem that you have to pick the number of clusters to work with, there is one other issue, however: k-means isn’t guaranteed to work as there’s no way to be certain that it will ever return the perfect solution. This is because the problem is what is known as NP-hard in computer science: it means that the only way to be certain you had gotten to the best solution, as far as we know, is to try every single possible assignment of data-points to clusters and check. I mention this solely because, frankly, I find that amazing, and also because you’re going to hear the phrase ‘NP-hard’ quite a lot if you Google around computer science things. There are lots of very, very complicated definitions of NP-hard floating around on the Internet: try and keep a hold of this slightly more intuitive, if not quite as precise, one, because it’ll help you stay afloat15.
8.5.2 Hands-on with k-means
Alright, this time let’s throw something a bit more tricky at our clustering algorithm. It’s not that hierarchical clustering can’t handle it, but rather that now we’re ready for a little more complexity.
data <- data.frame(rbind(
cbind(rnorm(50),rnorm(50)),
cbind(rnorm(50,2.5),rnorm(50,2.5)),
cbind(rnorm(50,-2.5),rnorm(50,-2.5)),
cbind(rnorm(50,5), rnorm(50,5)),
cbind(rnorm(50,-5), rnorm(50,-5))
))
data$groups <- rep(c("red","blue","grey80","grey60","grey20"), each=50)
names(data)[1:2] <- c("x","y")
with(data, plot(y ~ x, pch=20, col=groups))
As you can see, this data is a little bit more messy. Let’s not give our algorithm very many random re-starts, and see how it does. Notice how we’re going to use a contingency table to figure out whether our algorithm correctly assigns points to groups.
k.means <- kmeans(data[,-3], centers=5, nstart=10)
table(k.means$cluster, data$groups)
#>
#> blue grey20 grey60 grey80 red
#> 1 42 0 0 0 3
#> 2 3 0 50 0 0
#> 3 0 50 0 2 0
#> 4 5 0 0 0 45
#> 5 0 0 0 48 2Huh. Even with only ten re-starts (which is vastly too few for a real-world application), we did OK! But the problem is, as I’m sure you’ve recognized, that in order to even get started we had to decide how many groups we wanted. What if we didn’t know? Well, let’s see how our old friend the gap statistic does in this case.
library(splits)
gap.stat <- ddwtGap(data[,c("x","y")])
#> 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
with(gap.stat, plot(colMeans(DDwGap),pch=15,type='b',
ylim=extendrange(colMeans(DDwGap),f=0.2),
xlab="Number of Clusters", ylab="Weighted Gap Statistic"))
gap.stat$mnGhatWG
#> [1] 5You might have noticed that it took a little longer than it did last time; this is because we’ve thrown a lot more data at the method. If you had a much larger dataset, you’d probably have to look into other ways of measuring intra- vs. inter-group variance, but truth be told I don’t think you’d find yourself getting an answer much quicker or much better than the one we have here. Because, as you’ve probably noticed—it worked! Once again, we’ve detected the correct number of clusters!
8.5.3 Fuzzy clustering
If you’re like me, you probably don’t like putting people into boxes and categories. If so, then good news—neither do the people who develop clustering algorithms! Fuzzy clustering algorithms are extensions of standard clustering algorithms where, instead of assigning data-points to categories, we give them membership coefficients (often in a matrix called \(w\) that sum to \(1\) for each data-point) describing how each data-points ‘belongs’ to each category. Using them in R is pretty straightforward, once you’ve installed the relevant packages:
library(fclust)
fuzzy <- FKM(data[,c("x","y")], k=5)
# Get lots of summary statistics
summary(fuzzy)
#>
#> Fuzzy clustering object of class 'fclust'
#>
#> Number of objects:
#> 250
#>
#> Number of clusters:
#> 5
#>
#> Cluster sizes:
#> Clus 1 Clus 2 Clus 3 Clus 4 Clus 5
#> 51 50 52 46 51
#>
#>
#> Clustering index values:
#> SIL.F k=5
#> 0.7605857
#>
#>
#> Closest hard clustering partition:
#> Obj 1 Obj 2 Obj 3 Obj 4 Obj 5 Obj 6 Obj 7 Obj 8 Obj 9 Obj 10
#> 5 5 5 4 5 5 5 5 5 5
#> Obj 11 Obj 12 Obj 13 Obj 14 Obj 15 Obj 16 Obj 17 Obj 18 Obj 19 Obj 20
#> 5 5 5 5 5 5 5 5 5 5
#> Obj 21 Obj 22 Obj 23 Obj 24 Obj 25 Obj 26 Obj 27 Obj 28 Obj 29 Obj 30
#> 5 5 5 5 5 5 5 2 5 2
#> Obj 31 Obj 32 Obj 33 Obj 34 Obj 35 Obj 36 Obj 37 Obj 38 Obj 39 Obj 40
#> 5 5 5 5 4 5 5 5 5 5
#> Obj 41 Obj 42 Obj 43 Obj 44 Obj 45 Obj 46 Obj 47 Obj 48 Obj 49 Obj 50
#> 5 5 5 5 5 5 5 5 5 5
#> Obj 51 Obj 52 Obj 53 Obj 54 Obj 55 Obj 56 Obj 57 Obj 58 Obj 59 Obj 60
#> 4 5 1 4 4 4 4 4 4 4
#> Obj 61 Obj 62 Obj 63 Obj 64 Obj 65 Obj 66 Obj 67 Obj 68 Obj 69 Obj 70
#> 5 4 4 4 4 1 4 4 5 4
#> Obj 71 Obj 72 Obj 73 Obj 74 Obj 75 Obj 76 Obj 77 Obj 78 Obj 79 Obj 80
#> 4 4 4 4 4 4 4 4 5 4
#> Obj 81 Obj 82 Obj 83 Obj 84 Obj 85 Obj 86 Obj 87 Obj 88 Obj 89 Obj 90
#> 4 4 4 5 4 4 4 4 4 4
#> Obj 91 Obj 92 Obj 93 Obj 94 Obj 95 Obj 96 Obj 97 Obj 98 Obj 99 Obj 100
#> 4 4 4 4 4 4 4 4 4 4
#> Obj 101 Obj 102 Obj 103 Obj 104 Obj 105 Obj 106 Obj 107 Obj 108 Obj 109 Obj 110
#> 2 2 2 2 2 2 2 2 2 2
#> Obj 111 Obj 112 Obj 113 Obj 114 Obj 115 Obj 116 Obj 117 Obj 118 Obj 119 Obj 120
#> 2 2 2 2 2 2 2 2 2 2
#> Obj 121 Obj 122 Obj 123 Obj 124 Obj 125 Obj 126 Obj 127 Obj 128 Obj 129 Obj 130
#> 2 2 2 2 2 2 2 2 3 3
#> Obj 131 Obj 132 Obj 133 Obj 134 Obj 135 Obj 136 Obj 137 Obj 138 Obj 139 Obj 140
#> 2 2 2 2 2 2 2 2 2 2
#> Obj 141 Obj 142 Obj 143 Obj 144 Obj 145 Obj 146 Obj 147 Obj 148 Obj 149 Obj 150
#> 2 2 2 2 2 2 2 2 2 2
#> Obj 151 Obj 152 Obj 153 Obj 154 Obj 155 Obj 156 Obj 157 Obj 158 Obj 159 Obj 160
#> 1 1 1 1 1 1 1 1 1 1
#> Obj 161 Obj 162 Obj 163 Obj 164 Obj 165 Obj 166 Obj 167 Obj 168 Obj 169 Obj 170
#> 1 1 1 1 1 1 1 1 1 1
#> Obj 171 Obj 172 Obj 173 Obj 174 Obj 175 Obj 176 Obj 177 Obj 178 Obj 179 Obj 180
#> 1 1 1 1 1 1 1 1 1 1
#> Obj 181 Obj 182 Obj 183 Obj 184 Obj 185 Obj 186 Obj 187 Obj 188 Obj 189 Obj 190
#> 1 1 1 1 1 1 1 1 1 1
#> Obj 191 Obj 192 Obj 193 Obj 194 Obj 195 Obj 196 Obj 197 Obj 198 Obj 199 Obj 200
#> 1 1 1 1 1 4 1 1 1 1
#> Obj 201 Obj 202 Obj 203 Obj 204 Obj 205 Obj 206 Obj 207 Obj 208 Obj 209 Obj 210
#> 3 3 3 3 3 3 3 3 3 3
#> Obj 211 Obj 212 Obj 213 Obj 214 Obj 215 Obj 216 Obj 217 Obj 218 Obj 219 Obj 220
#> 3 3 3 3 3 3 3 3 3 3
#> Obj 221 Obj 222 Obj 223 Obj 224 Obj 225 Obj 226 Obj 227 Obj 228 Obj 229 Obj 230
#> 3 3 3 3 3 3 3 3 3 3
#> Obj 231 Obj 232 Obj 233 Obj 234 Obj 235 Obj 236 Obj 237 Obj 238 Obj 239 Obj 240
#> 3 3 3 3 3 3 3 3 3 3
#> Obj 241 Obj 242 Obj 243 Obj 244 Obj 245 Obj 246 Obj 247 Obj 248 Obj 249 Obj 250
#> 3 3 3 3 3 3 3 3 3 3
#>
#> Cluster memberships:
#> Clus 1 (First 50 objects)
#> [1] "Obj 53" "Obj 66" "Obj 151" "Obj 152" "Obj 153" "Obj 154" "Obj 155"
#> [8] "Obj 156" "Obj 157" "Obj 158" "Obj 159" "Obj 160" "Obj 161" "Obj 162"
#> [15] "Obj 163" "Obj 164" "Obj 165" "Obj 166" "Obj 167" "Obj 168" "Obj 169"
#> [22] "Obj 170" "Obj 171" "Obj 172" "Obj 173" "Obj 174" "Obj 175" "Obj 176"
#> [29] "Obj 177" "Obj 178" "Obj 179" "Obj 180" "Obj 181" "Obj 182" "Obj 183"
#> [36] "Obj 184" "Obj 185" "Obj 186" "Obj 187" "Obj 188" "Obj 189" "Obj 190"
#> [43] "Obj 191" "Obj 192" "Obj 193" "Obj 194" "Obj 195" "Obj 197" "Obj 198"
#> [50] "Obj 199"
#> Clus 2
#> [1] "Obj 28" "Obj 30" "Obj 101" "Obj 102" "Obj 103" "Obj 104" "Obj 105"
#> [8] "Obj 106" "Obj 107" "Obj 108" "Obj 109" "Obj 110" "Obj 111" "Obj 112"
#> [15] "Obj 113" "Obj 114" "Obj 115" "Obj 116" "Obj 117" "Obj 118" "Obj 119"
#> [22] "Obj 120" "Obj 121" "Obj 122" "Obj 123" "Obj 124" "Obj 125" "Obj 126"
#> [29] "Obj 127" "Obj 128" "Obj 131" "Obj 132" "Obj 133" "Obj 134" "Obj 135"
#> [36] "Obj 136" "Obj 137" "Obj 138" "Obj 139" "Obj 140" "Obj 141" "Obj 142"
#> [43] "Obj 143" "Obj 144" "Obj 145" "Obj 146" "Obj 147" "Obj 148" "Obj 149"
#> [50] "Obj 150"
#> Clus 3 (First 50 objects)
#> [1] "Obj 129" "Obj 130" "Obj 201" "Obj 202" "Obj 203" "Obj 204" "Obj 205"
#> [8] "Obj 206" "Obj 207" "Obj 208" "Obj 209" "Obj 210" "Obj 211" "Obj 212"
#> [15] "Obj 213" "Obj 214" "Obj 215" "Obj 216" "Obj 217" "Obj 218" "Obj 219"
#> [22] "Obj 220" "Obj 221" "Obj 222" "Obj 223" "Obj 224" "Obj 225" "Obj 226"
#> [29] "Obj 227" "Obj 228" "Obj 229" "Obj 230" "Obj 231" "Obj 232" "Obj 233"
#> [36] "Obj 234" "Obj 235" "Obj 236" "Obj 237" "Obj 238" "Obj 239" "Obj 240"
#> [43] "Obj 241" "Obj 242" "Obj 243" "Obj 244" "Obj 245" "Obj 246" "Obj 247"
#> [50] "Obj 248"
#> Clus 4
#> [1] "Obj 4" "Obj 35" "Obj 51" "Obj 54" "Obj 55" "Obj 56" "Obj 57"
#> [8] "Obj 58" "Obj 59" "Obj 60" "Obj 62" "Obj 63" "Obj 64" "Obj 65"
#> [15] "Obj 67" "Obj 68" "Obj 70" "Obj 71" "Obj 72" "Obj 73" "Obj 74"
#> [22] "Obj 75" "Obj 76" "Obj 77" "Obj 78" "Obj 80" "Obj 81" "Obj 82"
#> [29] "Obj 83" "Obj 85" "Obj 86" "Obj 87" "Obj 88" "Obj 89" "Obj 90"
#> [36] "Obj 91" "Obj 92" "Obj 93" "Obj 94" "Obj 95" "Obj 96" "Obj 97"
#> [43] "Obj 98" "Obj 99" "Obj 100" "Obj 196"
#> Clus 5 (First 50 objects)
#> [1] "Obj 1" "Obj 2" "Obj 3" "Obj 5" "Obj 6" "Obj 7" "Obj 8" "Obj 9"
#> [9] "Obj 10" "Obj 11" "Obj 12" "Obj 13" "Obj 14" "Obj 15" "Obj 16" "Obj 17"
#> [17] "Obj 18" "Obj 19" "Obj 20" "Obj 21" "Obj 22" "Obj 23" "Obj 24" "Obj 25"
#> [25] "Obj 26" "Obj 27" "Obj 29" "Obj 31" "Obj 32" "Obj 33" "Obj 34" "Obj 36"
#> [33] "Obj 37" "Obj 38" "Obj 39" "Obj 40" "Obj 41" "Obj 42" "Obj 43" "Obj 44"
#> [41] "Obj 45" "Obj 46" "Obj 47" "Obj 48" "Obj 49" "Obj 50" "Obj 52" "Obj 61"
#> [49] "Obj 69" "Obj 79"
#>
#> Number of objects with unclear assignment (maximal membership degree <0.5):
#> 11
#>
#> Objects with unclear assignment:
#> [1] "Obj 4" "Obj 16" "Obj 17" "Obj 49" "Obj 53" "Obj 55" "Obj 104"
#> [8] "Obj 155" "Obj 166" "Obj 188" "Obj 196"
#>
#> Cluster sizes (without unclear assignments):
#> Clus 1 Clus 2 Clus 3 Clus 4 Clus 5 No clus
#> 47 49 52 43 48 11
#>
#> Membership degree matrix (rounded):
#> Clus 1 Clus 2 Clus 3 Clus 4 Clus 5
#> Obj 1 0.03 0.23 0.05 0.09 0.61
#> Obj 2 0.01 0.04 0.01 0.04 0.90
#> Obj 3 0.02 0.11 0.03 0.08 0.75
#> Obj 4 0.06 0.05 0.02 0.49 0.39
#> Obj 5 0.01 0.02 0.01 0.02 0.94
#> Obj 6 0.05 0.09 0.03 0.20 0.63
#> Obj 7 0.04 0.27 0.06 0.11 0.51
#> Obj 8 0.04 0.07 0.02 0.18 0.69
#> Obj 9 0.00 0.01 0.00 0.01 0.98
#> Obj 10 0.01 0.04 0.01 0.03 0.91
#> Obj 11 0.01 0.05 0.01 0.04 0.88
#> Obj 12 0.03 0.22 0.05 0.09 0.61
#> Obj 13 0.06 0.06 0.02 0.34 0.52
#> Obj 14 0.01 0.02 0.01 0.02 0.95
#> Obj 15 0.00 0.02 0.00 0.01 0.96
#> Obj 16 0.08 0.08 0.03 0.41 0.41
#> Obj 17 0.06 0.22 0.07 0.17 0.48
#> Obj 18 0.01 0.03 0.01 0.07 0.88
#> Obj 19 0.07 0.09 0.03 0.30 0.51
#> Obj 20 0.02 0.03 0.01 0.10 0.84
#> Obj 21 0.03 0.17 0.04 0.08 0.69
#> Obj 22 0.07 0.10 0.04 0.24 0.56
#> Obj 23 0.01 0.04 0.01 0.02 0.93
#> Obj 24 0.01 0.02 0.01 0.02 0.95
#> Obj 25 0.01 0.11 0.02 0.04 0.81
#> Obj 26 0.01 0.03 0.01 0.03 0.92
#> Obj 27 0.01 0.11 0.02 0.04 0.82
#> Obj 28 0.02 0.57 0.06 0.05 0.28
#> Obj 29 0.03 0.20 0.04 0.09 0.63
#> Obj 30 0.01 0.85 0.03 0.02 0.09
#> Obj 31 0.02 0.30 0.04 0.06 0.58
#> Obj 32 0.05 0.05 0.02 0.38 0.51
#> Obj 33 0.03 0.09 0.02 0.09 0.77
#> Obj 34 0.03 0.07 0.02 0.11 0.78
#> Obj 35 0.07 0.05 0.02 0.58 0.29
#> Obj 36 0.00 0.00 0.00 0.00 0.99
#> Obj 37 0.04 0.05 0.02 0.20 0.69
#> Obj 38 0.02 0.10 0.02 0.07 0.79
#> Obj 39 0.03 0.05 0.02 0.18 0.73
#> Obj 40 0.02 0.08 0.02 0.05 0.83
#> Obj 41 0.02 0.05 0.01 0.08 0.85
#> Obj 42 0.02 0.09 0.02 0.05 0.82
#> Obj 43 0.01 0.12 0.02 0.05 0.80
#> Obj 44 0.02 0.18 0.03 0.06 0.71
#> Obj 45 0.02 0.03 0.01 0.09 0.85
#> Obj 46 0.02 0.06 0.02 0.09 0.81
#> Obj 47 0.02 0.06 0.01 0.05 0.86
#> Obj 48 0.01 0.03 0.01 0.03 0.93
#> Obj 49 0.05 0.35 0.11 0.12 0.37
#> Obj 50 0.02 0.13 0.03 0.08 0.74
#> Obj 51 0.03 0.01 0.00 0.89 0.06
#> Obj 52 0.04 0.07 0.02 0.19 0.68
#> Obj 53 0.46 0.02 0.01 0.46 0.05
#> Obj 54 0.07 0.02 0.01 0.78 0.11
#> Obj 55 0.09 0.07 0.03 0.47 0.34
#> Obj 56 0.03 0.01 0.01 0.86 0.08
#> Obj 57 0.00 0.00 0.00 0.99 0.00
#> Obj 58 0.04 0.02 0.01 0.81 0.12
#> Obj 59 0.01 0.00 0.00 0.98 0.01
#> Obj 60 0.11 0.06 0.03 0.52 0.28
#> Obj 61 0.03 0.05 0.01 0.19 0.72
#> Obj 62 0.05 0.03 0.01 0.65 0.25
#> Obj 63 0.09 0.01 0.01 0.85 0.04
#> Obj 64 0.06 0.04 0.02 0.57 0.31
#> Obj 65 0.01 0.00 0.00 0.96 0.02
#> Obj 66 0.65 0.02 0.01 0.26 0.06
#> Obj 67 0.03 0.01 0.00 0.90 0.06
#> Obj 68 0.03 0.01 0.00 0.93 0.03
#> Obj 69 0.05 0.05 0.02 0.36 0.51
#> Obj 70 0.03 0.01 0.00 0.94 0.02
#> Obj 71 0.30 0.04 0.02 0.54 0.11
#> Obj 72 0.05 0.01 0.01 0.88 0.05
#> Obj 73 0.13 0.03 0.02 0.70 0.12
#> Obj 74 0.12 0.02 0.01 0.75 0.09
#> Obj 75 0.25 0.02 0.01 0.67 0.06
#> Obj 76 0.05 0.01 0.00 0.90 0.04
#> Obj 77 0.02 0.01 0.00 0.91 0.05
#> Obj 78 0.14 0.01 0.01 0.78 0.05
#> Obj 79 0.04 0.06 0.02 0.17 0.72
#> Obj 80 0.17 0.07 0.03 0.52 0.21
#> Obj 81 0.01 0.00 0.00 0.98 0.01
#> Obj 82 0.06 0.02 0.01 0.85 0.07
#> Obj 83 0.11 0.02 0.01 0.81 0.06
#> Obj 84 0.03 0.04 0.01 0.16 0.75
#> Obj 85 0.10 0.01 0.00 0.85 0.04
#> Obj 86 0.08 0.02 0.01 0.79 0.09
#> Obj 87 0.13 0.01 0.01 0.81 0.04
#> Obj 88 0.09 0.04 0.01 0.71 0.15
#> Obj 89 0.03 0.00 0.00 0.94 0.02
#> Obj 90 0.06 0.01 0.01 0.88 0.05
#> Obj 91 0.05 0.02 0.01 0.82 0.11
#> Obj 92 0.17 0.02 0.01 0.75 0.06
#> Obj 93 0.05 0.04 0.01 0.63 0.27
#> Obj 94 0.15 0.02 0.01 0.75 0.06
#> Obj 95 0.04 0.01 0.00 0.89 0.05
#> Obj 96 0.26 0.02 0.01 0.63 0.07
#> Obj 97 0.38 0.02 0.01 0.53 0.06
#> Obj 98 0.05 0.04 0.01 0.57 0.32
#> Obj 99 0.08 0.03 0.01 0.72 0.15
#> Obj 100 0.09 0.01 0.01 0.85 0.04
#> Obj 101 0.00 0.94 0.03 0.01 0.02
#> Obj 102 0.01 0.85 0.07 0.02 0.06
#> Obj 103 0.01 0.89 0.05 0.01 0.05
#> Obj 104 0.01 0.46 0.45 0.02 0.05
#> Obj 105 0.00 0.95 0.01 0.01 0.02
#> Obj 106 0.00 0.92 0.03 0.01 0.03
#> Obj 107 0.01 0.79 0.05 0.03 0.12
#> Obj 108 0.01 0.84 0.04 0.02 0.09
#> Obj 109 0.00 0.97 0.01 0.00 0.01
#> Obj 110 0.01 0.80 0.07 0.02 0.09
#> Obj 111 0.01 0.73 0.05 0.03 0.17
#> Obj 112 0.00 0.99 0.00 0.00 0.00
#> Obj 113 0.01 0.64 0.27 0.02 0.06
#> Obj 114 0.00 0.95 0.02 0.01 0.03
#> Obj 115 0.01 0.89 0.06 0.01 0.04
#> Obj 116 0.01 0.77 0.16 0.02 0.05
#> Obj 117 0.01 0.72 0.19 0.02 0.06
#> Obj 118 0.02 0.65 0.09 0.05 0.19
#> Obj 119 0.00 0.99 0.01 0.00 0.01
#> Obj 120 0.01 0.86 0.09 0.01 0.03
#> Obj 121 0.00 0.94 0.02 0.01 0.03
#> Obj 122 0.02 0.65 0.08 0.05 0.20
#> Obj 123 0.00 0.92 0.04 0.01 0.03
#> Obj 124 0.01 0.85 0.09 0.01 0.04
#> Obj 125 0.01 0.89 0.05 0.01 0.04
#> Obj 126 0.01 0.74 0.18 0.02 0.05
#> Obj 127 0.00 0.97 0.01 0.00 0.01
#> Obj 128 0.01 0.77 0.07 0.03 0.11
#> Obj 129 0.02 0.27 0.61 0.03 0.07
#> Obj 130 0.01 0.20 0.74 0.01 0.03
#> Obj 131 0.01 0.60 0.32 0.02 0.05
#> Obj 132 0.01 0.70 0.22 0.02 0.05
#> Obj 133 0.00 0.93 0.02 0.01 0.04
#> Obj 134 0.01 0.84 0.07 0.02 0.06
#> Obj 135 0.00 0.99 0.00 0.00 0.01
#> Obj 136 0.01 0.74 0.15 0.02 0.07
#> Obj 137 0.00 0.99 0.01 0.00 0.00
#> Obj 138 0.01 0.88 0.06 0.01 0.04
#> Obj 139 0.02 0.74 0.08 0.03 0.13
#> Obj 140 0.01 0.80 0.04 0.03 0.12
#> Obj 141 0.01 0.82 0.07 0.02 0.08
#> Obj 142 0.00 0.91 0.02 0.01 0.05
#> Obj 143 0.01 0.82 0.11 0.02 0.05
#> Obj 144 0.01 0.55 0.33 0.03 0.08
#> Obj 145 0.00 0.92 0.04 0.01 0.03
#> Obj 146 0.00 1.00 0.00 0.00 0.00
#> Obj 147 0.01 0.75 0.04 0.03 0.16
#> Obj 148 0.02 0.56 0.28 0.04 0.10
#> Obj 149 0.01 0.76 0.06 0.03 0.14
#> Obj 150 0.01 0.87 0.04 0.02 0.07
#> Obj 151 0.82 0.01 0.01 0.14 0.03
#> Obj 152 0.55 0.03 0.01 0.33 0.07
#> Obj 153 0.80 0.01 0.01 0.15 0.03
#> Obj 154 0.65 0.03 0.02 0.22 0.08
#> Obj 155 0.49 0.03 0.01 0.39 0.08
#> Obj 156 0.88 0.01 0.00 0.09 0.02
#> Obj 157 0.85 0.01 0.01 0.11 0.03
#> Obj 158 0.86 0.01 0.00 0.11 0.02
#> Obj 159 0.89 0.01 0.00 0.08 0.02
#> Obj 160 0.93 0.00 0.00 0.05 0.01
#> Obj 161 0.83 0.02 0.01 0.11 0.04
#> Obj 162 0.99 0.00 0.00 0.01 0.00
#> Obj 163 0.96 0.00 0.00 0.03 0.01
#> Obj 164 0.80 0.02 0.01 0.14 0.04
#> Obj 165 0.74 0.02 0.01 0.19 0.05
#> Obj 166 0.49 0.02 0.01 0.44 0.05
#> Obj 167 0.92 0.01 0.00 0.05 0.02
#> Obj 168 0.94 0.01 0.00 0.04 0.01
#> Obj 169 0.86 0.01 0.00 0.11 0.02
#> Obj 170 0.73 0.01 0.01 0.22 0.03
#> Obj 171 0.83 0.02 0.01 0.10 0.04
#> Obj 172 0.53 0.02 0.01 0.38 0.06
#> Obj 173 1.00 0.00 0.00 0.00 0.00
#> Obj 174 0.62 0.03 0.02 0.26 0.07
#> Obj 175 0.84 0.02 0.01 0.10 0.03
#> Obj 176 0.83 0.01 0.00 0.13 0.02
#> Obj 177 0.94 0.01 0.00 0.04 0.01
#> Obj 178 0.70 0.01 0.01 0.24 0.04
#> Obj 179 0.67 0.01 0.01 0.27 0.04
#> Obj 180 0.94 0.00 0.00 0.04 0.01
#> Obj 181 0.85 0.01 0.01 0.09 0.03
#> Obj 182 0.95 0.00 0.00 0.04 0.01
#> Obj 183 0.69 0.02 0.01 0.24 0.04
#> Obj 184 0.76 0.02 0.01 0.16 0.04
#> Obj 185 0.94 0.01 0.00 0.04 0.01
#> Obj 186 0.80 0.02 0.01 0.13 0.04
#> Obj 187 0.79 0.02 0.01 0.13 0.05
#> Obj 188 0.50 0.02 0.01 0.42 0.06
#> Obj 189 0.89 0.01 0.01 0.07 0.02
#> Obj 190 0.85 0.01 0.01 0.09 0.03
#> Obj 191 0.89 0.01 0.01 0.07 0.02
#> Obj 192 0.63 0.02 0.01 0.28 0.06
#> Obj 193 0.94 0.00 0.00 0.05 0.01
#> Obj 194 1.00 0.00 0.00 0.00 0.00
#> Obj 195 0.94 0.01 0.00 0.04 0.01
#> Obj 196 0.42 0.03 0.02 0.44 0.08
#> Obj 197 0.90 0.01 0.01 0.06 0.02
#> Obj 198 0.87 0.01 0.01 0.09 0.03
#> Obj 199 0.93 0.01 0.00 0.05 0.01
#> Obj 200 0.99 0.00 0.00 0.01 0.00
#> Obj 201 0.01 0.16 0.76 0.02 0.05
#> Obj 202 0.00 0.01 0.98 0.00 0.00
#> Obj 203 0.01 0.27 0.65 0.02 0.05
#> Obj 204 0.00 0.03 0.95 0.00 0.01
#> Obj 205 0.00 0.01 0.99 0.00 0.00
#> Obj 206 0.01 0.19 0.76 0.01 0.03
#> Obj 207 0.01 0.17 0.74 0.02 0.05
#> Obj 208 0.01 0.10 0.86 0.01 0.02
#> Obj 209 0.00 0.09 0.89 0.01 0.02
#> Obj 210 0.00 0.00 1.00 0.00 0.00
#> Obj 211 0.00 0.04 0.94 0.01 0.01
#> Obj 212 0.00 0.04 0.94 0.00 0.01
#> Obj 213 0.01 0.06 0.90 0.01 0.02
#> Obj 214 0.01 0.34 0.58 0.02 0.05
#> Obj 215 0.01 0.12 0.82 0.02 0.04
#> Obj 216 0.00 0.06 0.92 0.01 0.01
#> Obj 217 0.00 0.04 0.94 0.01 0.01
#> Obj 218 0.00 0.00 1.00 0.00 0.00
#> Obj 219 0.00 0.00 1.00 0.00 0.00
#> Obj 220 0.00 0.01 0.99 0.00 0.00
#> Obj 221 0.01 0.05 0.92 0.01 0.02
#> Obj 222 0.01 0.09 0.85 0.02 0.03
#> Obj 223 0.01 0.11 0.85 0.01 0.02
#> Obj 224 0.01 0.11 0.81 0.02 0.04
#> Obj 225 0.01 0.23 0.71 0.01 0.04
#> Obj 226 0.00 0.03 0.95 0.00 0.01
#> Obj 227 0.00 0.01 0.99 0.00 0.00
#> Obj 228 0.00 0.04 0.94 0.01 0.01
#> Obj 229 0.00 0.03 0.95 0.00 0.01
#> Obj 230 0.00 0.05 0.93 0.01 0.01
#> Obj 231 0.01 0.14 0.82 0.01 0.03
#> Obj 232 0.00 0.01 0.99 0.00 0.00
#> Obj 233 0.00 0.00 0.99 0.00 0.00
#> Obj 234 0.01 0.11 0.82 0.02 0.04
#> Obj 235 0.00 0.00 1.00 0.00 0.00
#> Obj 236 0.01 0.29 0.61 0.02 0.06
#> Obj 237 0.00 0.04 0.94 0.00 0.01
#> Obj 238 0.01 0.13 0.82 0.01 0.03
#> Obj 239 0.00 0.03 0.96 0.00 0.01
#> Obj 240 0.00 0.02 0.97 0.00 0.01
#> Obj 241 0.00 0.05 0.93 0.00 0.01
#> Obj 242 0.01 0.07 0.89 0.01 0.02
#> Obj 243 0.00 0.04 0.94 0.01 0.01
#> Obj 244 0.01 0.08 0.87 0.01 0.03
#> Obj 245 0.00 0.04 0.94 0.01 0.01
#> Obj 246 0.01 0.11 0.86 0.01 0.02
#> Obj 247 0.01 0.07 0.89 0.01 0.02
#> Obj 248 0.00 0.03 0.96 0.00 0.01
#> Obj 249 0.01 0.07 0.88 0.01 0.02
#> Obj 250 0.02 0.28 0.61 0.03 0.07
#>
#> Cluster summary:
#> Cl.size Min.memb.deg. Max.memb.deg. Av.memb.deg. N.uncl.assignm.
#> Clus 1 51 0.46 1.00 0.81 4
#> Clus 2 50 0.46 1.00 0.82 1
#> Clus 3 52 0.58 1.00 0.88 0
#> Clus 4 46 0.44 0.99 0.76 3
#> Clus 5 51 0.37 0.99 0.75 3
#>
#> Euclidean distance matrix for the prototypes (rounded):
#> Clus 1 Clus 2 Clus 3 Clus 4
#> Clus 2 10.59
#> Clus 3 14.35 3.77
#> Clus 4 3.44 7.16 10.92
#> Clus 5 6.88 3.72 7.48 3.44
#>
#> Available components:
#> [1] "U" "H" "F" "clus" "medoid" "value"
#> [7] "criterion" "iter" "k" "m" "ent" "b"
#> [13] "vp" "delta" "stand" "Xca" "X" "D"
#> [19] "call"
#>
#>
# Plot the predictions out
plot(fuzzy)
This algorithm is essentially the same as k-means, only now we’re measuring the mean intra-cluster distance weighted by all the relative membership coefficients (a weighted mean). Where do those coefficients come from? Well, it depends… But they’re always some kind of weighted distance to the centroid. The precise form of the weighting function is up to you to decide (or go with the default); this is a remarkably active area of research and if you’re interested I recommend the review in the opening of Campello & Hruschka (2006; “A fuzzy extension of the silhouette width criterion for cluster analysis” in Fuzzy Sets and Systems 157: 2858–2875).
There are extensions of the gap statistics I have taught you for fuzzy clusters, but to be honest you’re better off using model-based clustering. Indeed, if you’re going to acknowledge that groupings are fuzzy such that something is almost in several clusters, you might be wondering where you draw the line—what’s the most likely cluster a point is in? Those of you who have taken my other statistics class might be wondering if there is some way to find the most likely cluster (…the maximum likelihood cluster…) for some data. Well, if you are, then you really should use model-based clustering, which is a natural extension of fuzzy methods!
8.6 Model-based clustering
If you’re the sort of person who’s not happy with approximations to things, you’re probably wondering what that “model-based clustering” I was talking about in the last section is all about. Well, it’s actually quite simple to do in R, even if it might not be so easy to understand. You can think of cluster-able data as coming from a series of different statistical distributions, each with their own parameters and possibly functional forms. It should, in principle, be possible to simply fit a load of different potential models to the data, then calculate the fit of those models using something like AIC or BIC and pick the model with the best fit16. Such a thing is possible, and is precisely what the package mclust does for you. Because, for speed, it operates in a frequentist perspective, it uses an expectation maximization algorithm; this basically means it tries to fit all the independent parts of its equation one step at a time, using (1) current parameter estimates to estimate the groupings, then (2) using those groupings to re-estimate the best parameter estimates, and then going back to step (1). This looping repeats until the algorithm has converged on a set of numbers. In well-behaved data, it also works like a charm…
library(mclust)
#> Package 'mclust' version 6.1.3
#> Type 'citation("mclust")' for citing this R package in publications.
#>
#> Attaching package: 'mclust'
#> The following object is masked from 'package:mvtnorm':
#>
#> dmvnorm
model <- Mclust(data[,-3])
summary(model)
#> ----------------------------------------------------
#> Gaussian finite mixture model fitted by EM algorithm
#> ----------------------------------------------------
#>
#> Mclust EII (spherical, equal volume) model with 5 components:
#>
#> log-likelihood n df BIC ICL
#> -1056.169 250 15 -2195.16 -2222.762
#>
#> Clustering table:
#> 1 2 3 4 5
#> 50 45 49 53 53
model
#> 'Mclust' model object: (EII,5)
#>
#> Available components:
#> [1] "call" "data" "modelName" "n"
#> [5] "d" "G" "BIC" "loglik"
#> [9] "df" "bic" "icl" "hypvol"
#> [13] "parameters" "z" "classification" "uncertainty"
plot(model, what="BIC")
If you’re interested in learning more about this approach, I can’t recommend enough reading Ezard et al.’s 2010 paper Algorithmic approaches to aid species’ delimitation in multidimensional morphospace in BMC Evolutionary Biology (10:1, 175). It’s a (remarkably approachable, given the subject matter) introduction to the topic, and it uses these techniques to solve a very interesting species delimitation problem.
8.7 Exercises
Remember, you have two tasks. (1) To read through the material above, make sure you understand it, and that you can interpret and understand everything that is going on in the code examples. (2) Carry out the exercises below. Note also that they involve answering questions and thinking; you could write the code in seconds (as is the case for much of the work in machine learning), and so your task is to make sure you understand what is going on.
- Below is code to download global temperature and precipitation data. Use it to carry out the following exercises.
library(terra)
#> terra 1.9.34
#>
#> Attaching package: 'terra'
#> The following object is masked from 'package:MASS':
#>
#> area
#> The following objects are masked from 'package:ape':
#>
#> rotate, trans, zoom
library(geodata)
r <- worldclim_global(var = "bio", res = 10, path = tempdir())
#> Cached as: /tmp/RtmpMeomFZ/climate/wc2.1_10m//wc2.1_10m_bio.zip
points <- expand.grid(lat = seq(0, 50), long = seq(-50, -150))
v.points <- vect(points, geom = c("long", "lat"), crs = crs(r))
data <- extract(r, v.points, ID = FALSE)
names(data) <- c("temp.mean", "diurnal.range", "isothermality",
"temp.season", "max.temp", "min.temp", "temp.range", "temp.wettest",
"temp.driest", "temp.mean.warmest", "temp.mean.coldest", "precip",
"precip.wettest", "precip.driest", "precip.season", "precip.wettest",
"precip.driest", "precip.warmest", "precip.coldest")
data$lat <- points$lat
data$long <- points$long
data <- na.omit(data)- Perform a scaled PCA on the data.
- Produce a biplot and scree plot of the data.
- How many important axes of variation are there in global temperature and precipitation on the basis of this data? Describe each axis, but if you write more than a sentence per axis you are writing too much.
- Calculate the correlation matrix for these data, and use it to confirm for yourself the conceptual links between a PCA and PCoA.
- Perform an NMDS on these data. Do you learn anything more or less than the PCA or PCoA?
- These questions are based around a classic ecological dataset: the Barro Colorado Island 50-hectare plot. Before the Panama Canal was built, it was a small patch of a much larger forest; now it is a roughly 100-hectare island upon which much of modern ecology has been based. Read into that what you will! You can load the data into
Rby runningdata(BCI)(it’s shipped withvegan).
- Perform a PCoA of the BCI plot dataset using a Euclidean distance matrix.
- Perform an NMDS of the BCI plot dataset.
- Repeat question 1a using a Bray-Curtis distance matrix (use the function
vegdist). Compare your new PCoA with the output from question 1b. Why is it more similar to your answer to 1b?
- These questions are based around a dataset of plant abundances across a gradient in temperature. You will find the code below useful to both load that data, plot it out, and calculate distance matrices from it. Your task is to figure out what is driving differences in community composition within this data: one of your supervisors thinks there are two kinds of communities here (hot vs. cold environments), while the other thinks species are responding to a gradient of temperature. Your task is to help them figure out what’s going on.
# Load the data in
comm <- as.matrix(read.table("hot-sites.txt"))
site.data <- read.table("site-data.txt")
# Build datasets of distances for quantile regression
dist.data <- with(site.data, data.frame(
dist=as.numeric(your.distance.matrix),
temp=as.numeric(dist(temp)),
groups=as.numeric(dist(outer(groups, groups, `==`), method="binary"))
))
# Color plots of data according to a categorical variable
some.plot.function(some.data, col=ifelse(groups=="hot", "red", "blue"))- Perform an NMDS or PCoA of these data.
- Visualise this data, colouring them according to the groupings (hot vs. cold) one of your supervisors favours.
- Add to your visualisation from the above question some contour lines reflecting the environmental gradient.
- Use a quantile regression to determine whether there are two kinds of communities, or whether there is a graded, continuous response to the environment.
- Below is code that will load in over 100 years of electoral data from the US, recording the proportion of the vote in each state that the Republican party received. Notice how we use a logistic transformation on the data to normalize them because percentages are not normally distributed; if you are unfamiliar with this feel free to ask me or not to worry about it and take my class next semester on the analysis of such data. I am asking you to analyze voting data, not to pass judgment on anyone who votes; I am not intending to somehow influence your voting intentions, I’m just showing you some history.
library(cluster)
votes <- na.omit(cluster::votes.repub)
logit <- function(x) log(x / (1-x))
transformed <- logit(votes/100)- Perform a hierarchical cluster analysis of these data.
- Use
ddwtGapwith the argumentgenRndm="uni"to see how many statistically significant clusters there are according to this metric. - Do you think the output of
ddwtGapis correct? Why (not)?17
- Download the file
board-game-mechanics.csvfrom the course website. This is a binary matrix of the top-100 games on https://boardgamegeek.com/, detailing the mechanics of each board game18. Download that data, load it in using the code below, and use it to carry out the following exercises. You do not require any knowledge of board games to answer these questions, and I strongly encourage you not to attempt to develop any in order to answer these questions.
data <- read.csv("board-game-mechanics.csv", row.names=1)- Build a hierarchically-clustered dendrogram of these data. They’re binary data, so create a
binarydistance matrix. - Compare, graphically, your dendrogram from part (a) with a dendrogram built using the
averagemethod. What are the main differences you can see? Which dendrogram do you think would be the best to present to a policy-maker, and why? - It takes too long for us to use
ddwtGapeasily on this dataset19, but try usingcutreeto split your data into groups. Does this change your opinion about which of the methods you would prefer to work with? How do you think the best way to split this tree is?20
- The code below reminds you how to load the
irisdataset intoR; use it to carry out the following exercises.
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- Use
ddwtGapto determine the optimal number of clusters in this dataset. - Run a k-means analysis with at least ten restarts, using the number of clusters you identified above.
- Build a contingency table of your k-means clustering versus the reality of the species definitions. Do you think that k-means has done a better job than PCA did of describing this dataset? Why?
- Repeat this analysis using a model-based clustering approach.
- The following exercise makes use of data from Wright et al. (2004; The worldwide leaf economics spectrum. Nature, 428(6985) 821), and this data should not be shared outside this class in any way21. It concerns data collected on a number of traits related to plant leaves (the details don’t matter); each plant’s average biomass when grown in an experimental garden is also listed in the dataset as a response variable22. You can find the data in the file
glopnet-biomass.csvon your course website.
- Fit a standard linear model (e.g., something like
lm(biomass ~ some.variables, data=data)) to see what plant traits explain biomass in this experiment. What does you model tell you? Why is it not working very well? - You now know at least three different ways of summarizing data for downstream analysis: PCA, hierarchical cluster analysis, and k-means cluster analysis. Explore this dataset using these methods to see what you think summarizes the data best.
- Use whichever method you think performs best to fit a new linear model of the data. Does this perform better or worse than the original model? Why?
I will stop putting this word in italics now, but you must start using it when talking about PCA. Replace it with ‘at right angles to everything else’ in your head if it helps.↩︎
Before doing any serious analysis in
Rafter this practical, you should close downRand start it up again to make sure you get rid of any weirdness my use ofset.seedwould have had. Ask me for more details about how all this works if you’re interested.↩︎Jackson, D. A. (1993). Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches. Ecology, 74(8), 2204-2214. It’s a good paper; if you’re interested, read it.↩︎
Some, of course, are more significant than others…↩︎
As my other footnote indicates, in reality that’s not that but hey ho.↩︎
If you’re interested in extending this to work with abundance data, talk to me. It’s quite fun to do so, but would add extra complexity I don’t want to bother the class as a whole with.↩︎
A negative eigenvalue means there’s no real (in the mathematical sense) solution to the regression equation you’re trying to fit. Thus their presence denotes non-Euclidean variation (variation we can’t find a way to map), and the size of the eigenvalue denotes its importance.↩︎
My class, my rules.↩︎
If you must learn it, Jonathan Lefcheck is sufficiently smart that it might as well be from him: https://jonlefcheck.net/2012/10/24/nmds-tutorial-in-r/↩︎
e.g., Pearse et al. (2013). Ecology, 94(12), 2861-2872.↩︎
Of course, the molecular clock is an absurd over-simplification and isn’t true, but don’t let reality get in the way of a good story…↩︎
If you’re interested, come by my office and borrow my MSc thesis. I assure you that you won’t enjoy the experience of reading it; I certainly didn’t enjoy the experience of writing it.↩︎
Oh, alright, I’ll give you one more piece of advice: Brian Ripley has written reams about this topic, and he’s really quite smart. So take a look at his work.↩︎
There are, however, multiple algorithms that are called the k-means algorithm. This is the most common, and is sometimes called Lloyd’s algorithm (arguably not correctly). As with everything in life, if it’s simple, someone’s got to complicate it.↩︎
There are so many definitions and articles about it because one of the classic unsolved problems in computer science: whether every problem whose answer can quickly be checked (NP problems) can quickly be solved (P problems). If \(P \neq NP\), meaning every quick to check problem can’t quickly be solved, then we’re all fine and dandy. \(P=NP\) then pretty much every intuition we have about computing and statistics is wrong. Gulp.↩︎
I would argue this is no longer a machine learning approach—it’s just a complicated classical statistics problems. Indeed, this is actually precisely how we simulated our example data above↩︎
There is no correct answer to this question. I am trying to get you to develop a healthy suspicion of the output of any machine learning algorithm.↩︎
Those long winter evenings just fly by in the Pearse-Waring household!↩︎
Try it if you’re interested, but be patient and set the
maxClustto be a bit higher than the defaults.↩︎There is no perfect answer to this question. Use your brain, use your intuition, and think about what seems the most intuitively appealing.↩︎
I am required to state this by the conditions of the data release; I would also recommend you not work uncritically with data I have manipulated for a classroom exercise↩︎
Please don’t think about the experiment too much, because it’s not real.↩︎