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
Avocado 100     60      13      8
Avocado 95      70      14      6

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
Fruit<-c("Apple","Apple","Orange","Orange","Banana","Banana","Banana","Avocado","Avocado")
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)
rownames(myData)<-c("Apple1","Apple2","Orange1","Orange2","Banana1","Banana2","Banana3","Avocado1","Avocado2")


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

 A: 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")
  model<-betadisper(distance, myData[,1])
  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    Avocado     Banana Orange
Apple              8.0000e-03 8.5600e-01  0.004
Avocado 9.2864e-31            7.9800e-01  0.008
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 


A: Are you looking for MANOVA? This is the multivariate generalization of ANOVA where you are testing for differences between mean vectors. In R the manova function will fit this model.
