# Multivariate grouping: clustering, anova, tukey

I have multiple variables (here: weight, horizontal diameter, price and dummy) related to different factors (here: Apple, Orange, Banana and Avocado):

Fruit   Weight      HorDiam     Price       Dummy
Apple   60      60      5       4
Apple   50      70      8       6
Orange  80      75      7       2
Orange  72      70      9       8
Banana  40      30      3       1
Banana  45      35      4       2
Banana  80      50      8       3


I need to test if I can group some species together: are apples and oranges significantly different? ANOVA tells me if weight (or horizontal diameter, or price) is significantly different among species. Tukey test gives me if weight of one species is significantly different from weight of another one (pairwise). Clustering seems only able to group individual observations together, not species. I can't find the appropriate test (or algorithm) to tell me if, for a single variable (weight) or for all of them (weight, horDiam and price), apples can be grouped with oranges and/or with bananas. Any suggestion?

I created a R code for this example:

### CREATE TABLE
Weight<-c(60,50,80,72,40,45,85,90,95)
horDiam<-c(60,70,75,70,30,35,50,60,70)
Price<-c(5,8,7,9,3,4,8,13,14)
Dummy<-c(4,6,2,8,1,2,3,8,6)
myData<-data.frame(Fruit=Fruit, Weight=Weight, horDiam=horDiam, Price=Price, Dummy=Dummy)

### ANOVA
fit.aov<-list()
summaryAOV<-list()
for (i in 1:3){
fit.aov[[i]]<-aov(myData[,i+1]~myData[,1])
summaryAOV[[i]]<-summary(fit.aov[[i]])
}

### TUKEY
par(mfrow=c(1,3))
testTukey<-list()
mainTukey<-c("Weight", "Horiz. Diameter", "Price")
for (i in 1:3){
testTukey[[i]]<-TukeyHSD(fit.aov[[i]], conf.level = 0.95)
plot(testTukey[[i]], main=mainTukey[i])
}

### CLUSTERING
plot( hclust(dist(myData), method="ward") )

### CLUSTERING WITH P-VALUE
fit <- pvclust(t(myData[,-1]), method.hclust="ward", method.dist="euclidean")
plot(fit)
pvrect(fit, alpha=0.95)

-
I don't mean to sound negative, but in your paragraph you've progressed from what sound like interesting and useful questions (do the groups show differences in this or that respect) to a much less interesting one (can they be grouped together--yes or no). What makes you dissatisfied with the nuanced approach and steers you toward that seemingly sterile, either-or question? –  rolando2 Apr 5 '11 at 0:01
The quest for a statistical procedure seems to sweep under the rug what appears to be the key question: how does one compare incommensurate attributes? Specifically, how should a difference of (say) 10 gm in weight be balanced against a difference of 1 in price? If you just throw your data into a PCA/LDA/clustering/MANOVA/whatever grinder, you'll get sophisticated and perhaps interesting answers, but they won't answer this question and each (implicitly and silently) will make its own assumptions about what the answer is--for good or ill. –  whuber Apr 6 '11 at 16:14

You might have a look at the betadisper() function in the vegan package. The function implements the PERMDISP2 procedure (Anderson, 2006) for the analysis of multivariate homogeneity of group dispersions. An example using your data might be the following:

  require(vegan)
distance<-vegdist(myData[,2:5], method="euclidean")
permutest(model, pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions

No. of permutations: 999
Permutation type: free
Permutations are unstratified
Mirrored permutations?: No
Use same permutation within strata?: No

Response: Distances
Df  Sum Sq Mean Sq      F N.Perm Pr(>F)
Groups     3  223.28  74.427 0.3523    999   0.84
Residuals  5 1056.34 211.268

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Apple              8.0000e-03 8.5600e-01  0.004
Banana  6.2093e-01 5.6617e-01             0.797
Orange  1.5060e-31 4.0863e-27 5.6544e-01


Below I have inserted a plot of groups and distances to the group centroid [plot(model)], and a boxplot of the distances to centroid for each group [boxplot(model)].

Hope this helps.

### References

Anderson, M. J. (2006) Distance-based tests for homogeneity of multivariate dispersions. Biometrics 62(1): 245–253.

## Edit

On a second thought I would recommend also a more descriptive approach using linear discriminant analysis (LDA) that could help not only to visualise the spread of objects around their group centroids, but also to find the features that contribute to this configuration. The ade4 package contains the versatile function discrimin() that does this as follows:

require(ade4)
discr <- discrimin(dudi.pca(myData[,2:5], scan = FALSE), myData[,1], scan = FALSE)


Note that the LDA is based on a PCA (function dudi.pca()) of the data so you will need to consider its properties when applying it to your task.

The top left plot represents the coefficients of the linear discriminant functions on the first two axes of the DA. The "Cos(variates, canonical variates)" plot shows the covariances between the object properties projected on the first two axes. Then, on the bottom left is the eigenvalue screeplot demonstrating the contribution of each axis to the variation. The main plot, "Scores and Classes", shows the projections of the individuals on the plane defied by the axes of the DA. Groups are displayed by ellipses where the centers are the means and the ellipses show the variance within the objects. All plots are the result of the plot(discr) command.

A randomisation test (and plot via plot(randtest.discrimin(discr))) of the eigenvalue significance is also available:

randtest.discrimin(discr)
Monte-Carlo test
Call: randtest.discrimin(xtest = discr)

Observation: 0.5074292

Based on 999 replicates
Simulated p-value: 0.052
Alternative hypothesis: greater

Std.Obs Expectation    Variance
1.549410152 0.376074589 0.007187167


-
Thanks, you helped me a lot. I discovered some great tools while discovering betadisper(). Namely, discriminant analysis using package ape4 and principal component analysis with ellipses drawn using function clusplot() from package cluster. –  dernier recours Apr 6 '11 at 1:52
I have just finished a draft for a second answer based on the idea that you might want also a LDA, so I'll post it in a while using your data as the example. –  ils Apr 6 '11 at 11:25
@dernier It will also tell you if there are significant differences between pairs of species: mod <- lm(cbind(Weight, horDiam, Price, Dummy) ~ Fruit, data=myData); summary(mod, test="Pillai"). If your test for differences will be post hoc, you will have to figure out how to implement a multiple comparisons procedure yourself, because neither TukeyHSD nor multcomp seem to work on multivariate models. It is not standard practice, in my field at least, to create groups for further analysis in this way. –  lockedoff Apr 5 '11 at 14:45