This script reproduces all of the analysis and graphs for betadisper examples in the paper and also includes other analyses not described there. It is set up as an R script that can be “compiled” to HTML, Word, or PDF using knitr::knit(). This is most convenient within R Studio via the File -> Compile Notebook option.

Load packages and the data

NB: This example requires vegan >= 2.4-0

library(vegan)
data(iris)

Calculate distances and run betadisper

dst <- dist(iris[,1:4])
iris.bd <- betadisper(dst, iris$Species)
iris.bd
## 
##  Homogeneity of multivariate dispersions
## 
## Call: betadisper(d = dst, group = iris$Species)
## 
## No. of Positive Eigenvalues: 4
## No. of Negative Eigenvalues: 0
## 
## Average distance to median:
##     setosa versicolor  virginica 
##      0.481      0.706      0.816 
## 
## Eigenvalues for PCoA axes:
##   PCoA1   PCoA2   PCoA3   PCoA4    <NA>    <NA>    <NA>    <NA> 
## 630.008  36.158  11.653   3.551      NA      NA      NA      NA
anova(iris.bd)
## Analysis of Variance Table
## 
## Response: Distances
##            Df Sum Sq Mean Sq F value  Pr(>F)    
## Groups      2   2.91   1.455    10.8 4.4e-05 ***
## Residuals 147  19.89   0.135                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permutest(iris.bd)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##            Df Sum Sq Mean Sq    F N.Perm Pr(>F)    
## Groups      2   2.91   1.455 10.8    999  0.001 ***
## Residuals 147  19.89   0.135                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
labs <- paste("Dimension", 1:4, "(", 
              round(100*iris.bd$eig / sum(iris.bd$eig), 2), "%)")

plot(iris.bd, cex=1, pch=15:17,
     main="Iris data: MDS coordinates", cex.lab=1.25,
     xlab=labs[1], ylab=labs[2],
     hull=FALSE, ellipse=TRUE, conf=0.68, lwd=2)

Boxplot of dispersions

boxplot(iris.bd, xlab="Species", notch=TRUE, col=c("gray", "red", "green"))

Skulls data

data(Skulls, package="heplots")
dst <- dist(Skulls[, -1])
skulls.bd <- betadisper(dst, Skulls$epoch)
skulls.bd
## 
##  Homogeneity of multivariate dispersions
## 
## Call: betadisper(d = dst, group = Skulls$epoch)
## 
## No. of Positive Eigenvalues: 4
## No. of Negative Eigenvalues: 0
## 
## Average distance to median:
## c4000BC c3300BC c1850BC  c200BC  cAD150 
##    8.42    7.80    7.74    7.72    8.95 
## 
## Eigenvalues for PCoA axes:
## PCoA1 PCoA2 PCoA3 PCoA4  <NA>  <NA>  <NA>  <NA> 
##  5301  3536  2825  1380    NA    NA    NA    NA
anova(skulls.bd)
## Analysis of Variance Table
## 
## Response: Distances
##            Df Sum Sq Mean Sq F value Pr(>F)
## Groups      4     35    8.86    0.82   0.51
## Residuals 145   1564   10.79
permutest(skulls.bd)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##            Df Sum Sq Mean Sq    F N.Perm Pr(>F)
## Groups      4     35    8.86 0.82    999   0.55
## Residuals 145   1564   10.79
plot(skulls.bd, main="Skulls data: MDS coordinates", 
     hull=FALSE, ellipse=TRUE, conf=0.68)

boxplot(skulls.bd, notch=TRUE, col="lightblue", xlab="Epoch")

Wine data

data(Wine, package="candisc")
dst <- dist(Wine[, -1])
wine.bd <- betadisper(dst, Wine$Cultivar)
wine.bd
## 
##  Homogeneity of multivariate dispersions
## 
## Call: betadisper(d = dst, group = Wine$Cultivar)
## 
## No. of Positive Eigenvalues: 12
## No. of Negative Eigenvalues: 0
## 
## Average distance to median:
##     barolo grignolino    barbera 
##      176.1      123.1       93.9 
## 
## Eigenvalues for PCoA axes:
##     PCoA1     PCoA2     PCoA3     PCoA4     PCoA5     PCoA6     PCoA7 
## 1.756e+07 3.054e+04 1.671e+03 8.834e+02 2.175e+02 1.489e+02 4.938e+01 
##     PCoA8 
## 2.679e+01
anova(wine.bd)
## Analysis of Variance Table
## 
## Response: Distances
##            Df  Sum Sq Mean Sq F value  Pr(>F)    
## Groups      2  190083   95041    8.33 0.00035 ***
## Residuals 175 1997005   11411                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permutest(wine.bd)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##            Df  Sum Sq Mean Sq    F N.Perm Pr(>F)    
## Groups      2  190083   95041 8.33    999  0.001 ***
## Residuals 175 1997005   11411                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#plot(wine.bd, hull=FALSE, ellipse=TRUE, conf=0.68)
boxplot(wine.bd, xlab="Cultivar", notch=TRUE, col=c("red", "green", "blue"))