LabDSV: Working with Vegetation Tables

Options
Libraries
vegtab
Ordering by Species
Ordering by Environment
Constancy
Importance
Constancy/Cover Tables

Lab 3 — Vegetation Tables and Summaries

Before moving on to formal statistical models it's helpful to review some of the more traditional approaches to vegetation analysis using sorted vegetation tables. Building on the extensive capabilities of R to work with matrices and data.frames, The vegan and LabDSV libraries enable vegetation ecologists to use R for such work (see R for Ecologists if you are unfamiliar with libraries).

In Lab 1 and Lab 2 we explored a dataset from Bryce Canyon National Park. While the Bryce Canyon vegetation data is relatively modest in size (160 plots x 169 species), it's too large to see easily just by listing the whole matrix. For example, if you have loaded the vegetation data file in previous labs, then

veg

will list the entire matrix to your screen, but it's likely to wrap and scroll off. You can control the wrapping aspects with the options() as follows. Suppose you have your window set to a width of 200 characters,

options(width=200)

will produce output that does not wrap unless it's more than 200 characters wide. In addition, to control the periodic insertion of the column headings,

options(height=100)

will cause the column headings to be inserted only every 100 lines. These can of course be combined in a single call.

options(height=100,width=200)

Still, the Bryce vegetation data will scroll uncontrollably. As a first step toward solving this (and other) problem, load the LabDSV and vegan libraries.

library(labdsv)
library(vegan)

This provides several functions we will need. The first is vegtab(), which produces sorted vegetation tables with several specifiers. Specifically,

the vegetation matrix to list
set= a specific subset of plots within that matrix to list (optional)
min= a minimum number of plots a species must occur in to be included (optional)
pltord= a specific order for the plots (optional)
spcord= a specific order for the species (optional)
pltlbl= an alternative plot label for the plots (optional)

The default behavior is to list plots in the order in which they occur in the matrix, but to list species in order from most common to least common. You can specify alternative orders with the optional function arguments. For example, to see the list of most common species that occur on sites with deep soils,

vegtab(veg,site$depth=="deep",min=20)

        chrvis artarb sticom senmul sithys poafen calnut erirac koenit
bcnp__1    0.0    0.0    0.0    0.0    0.0    0.5    0.0    0.0    0.0
bcnp__8    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.0    0.0
bcnp_12    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_13    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_20    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_30    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_31    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.0    0.0
   .        .      .      .      .      .      .      .      .      .
   .        .      .      .      .      .      .      .      .      .
   .        .      .      .      .      .      .      .      .      .
bcnp150    0.0    0.0    0.0    0.0    0.0    0.5    0.0    0.5    0.5
bcnp154    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp158    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp160    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0

where ellipses represent omitted material (to save space). As you can see, there are only 9 species which occur at least 20 times in the set of plots on deep soils (n=56), and the most common is chrvis, followed by artarb, etc.

The plots can be ordered by another variable. For example, to list plots in order of their elevation,

vegtab(veg,site$depth=="deep",min=20,pltord=elev)

         chrvis artarb sticom senmul sithys poafen calnut erirac koenit
bcnp_89     0.0    0.0    0.0    0.0    2.0    0.0    0.0    0.0    0.0
bcnp_90     0.5    0.0    0.0    0.0    5.0    0.0    0.5    0.0    0.0
bcnp_99     0.5    0.0    0.0    0.0    1.0    0.5    0.5    0.0    0.0
bcnp_88     0.5    0.0    0.0    0.0    1.0    0.0    0.5    0.0    0.0
bcnp_98     0.5    0.0    0.0    0.0    1.0    0.0    0.5    0.0    0.0
bcnp_93     0.5    0.0    0.5    0.0    1.0    0.0    0.5    0.0    0.0
      .      .      .      .      .      .      .      .      .      .
      .      .      .      .      .      .      .      .      .      .
      .      .      .      .      .      .      .      .      .      .
bcnp105    2.0    2.0    0.5    0.5    0.5    0.5    0.5    0.5    0.0
bcnp102    1.0    2.0    1.0    0.0    0.5    2.0    0.5    0.0    0.0
bcnp103    0.5    2.0    2.0    0.5    0.5    0.5    0.5    0.5    0.0
bcnp_40    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.5    0.0
bcnp_38    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
     .      .      .      .      .      .      .      .      .      .
     .      .      .      .      .      .      .      .      .      .
     .      .      .      .      .      .      .      .      .      .
bcnp_33    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_31    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.0    0.0
bcnp_30    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_65    0.0    0.0    0.0    0.5    0.5    0.5    0.0    0.0    0.0
bcnp__8    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.0    0.0
bcnp__1    0.0    0.0    0.0    0.0    0.0    0.5    0.0    0.0    0.0

To actually list the elevation, specify pltlbl=elev as well as pltord=elev.

vegtab(veg,site$depth=="deep",pltlbl=elev,pltord=elev,min=20)

         elev chrvis artarb sticom senmul sithys poafen calnut erirac koenit
bcnp_89  6650    0.0    0.0    0.0    0.0    2.0    0.0    0.0    0.0    0.0
bcnp_90  6680    0.5    0.0    0.0    0.0    5.0    0.0    0.5    0.0    0.0
bcnp_99  6720    0.5    0.0    0.0    0.0    1.0    0.5    0.5    0.0    0.0
bcnp_88  6760    0.5    0.0    0.0    0.0    1.0    0.0    0.5    0.0    0.0
bcnp_98  6760    0.5    0.0    0.0    0.0    1.0    0.0    0.5    0.0    0.0
      .    .      .      .      .      .      .      .      .      .      .
      .    .      .      .      .      .      .      .      .      .      .
      .    .      .      .      .      .      .      .      .      .      .
bcnp103  7720    0.5    2.0    2.0    0.5    0.5    0.5    0.5    0.5    0.0
bcnp_40  7800    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.5    0.0
bcnp_38  7840    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp135  7850    1.0    0.5    0.0    0.5    0.0    0.0    0.0    0.0    0.5
      .    .      .      .      .      .      .      .      .      .      .
      .    .      .      .      .      .      .      .      .      .      .
      .    .      .      .      .      .      .      .      .      .      .
bcnp_31  8560    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.0    0.0
bcnp_30  8600    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_65  8660    0.0    0.0    0.0    0.5    0.5    0.5    0.0    0.0    0.0
bcnp__8  8680    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.0    0.0
bcnp__1  8720    0.0    0.0    0.0    0.0    0.0    0.5    0.0    0.0    0.0

Finally, the table can be transposed with the trans=T argument.

vegtab(veg,site$depth=="deep",pltlbl=elev,pltord=elev,min=20,trans=T)

       bcnp_89 bcnp_90 bcnp_99 bcnp_88 bcnp_98 bcnp_93 bcnp158 bcnp160 bcnp148
elev      6650  6680.0  6720.0  6760.0  6760.0  6900.0    6960    7040    7080
chrvis       0     0.5     0.5     0.5     0.5     0.5       0       0       0
artarb       0     0.0     0.0     0.0     0.0     0.0       0       0       0
sticom       0     0.0     0.0     0.0     0.0     0.5       0       0       0
senmul       0     0.0     0.0     0.0     0.0     0.0       0       0       0
sithys       2     5.0     1.0     1.0     1.0     1.0       0       0       0
poafen       0     0.0     0.5     0.0     0.0     0.0       0       0       0
calnut       0     0.5     0.5     0.5     0.5     0.5       0       0       0
erirac       0     0.0     0.0     0.0     0.0     0.0       0       0       0
koenit       0     0.0     0.0     0.0     0.0     0.0       0       0       0
  .          .      .       .       .        .      .        .       .       .  
  .          .      .       .       .        .      .        .       .       .  
  .          .      .       .       .        .      .        .       .       .

It's possible to combine multiple criteria in the set= argument. For example, suppose you wanted to see the vegetation of plots on steep midslopes.

vegtab(veg,site$pos=="mid_slope"&site$slope>20,min=2)

        cermon arcpat ceamar purtri carrss swerad
bcnp_51    0.5      5    0.5    0.5    0.5    0.5
bcnp_75    1.0      3    1.0    0.0    0.5    0.5
bcnp155    1.0      0    0.0    0.5    0.0    0.0

There are only three plots meeting the criteria, so the number of species occurring at least twice is small, but it illustrates the idea. If the criteria get really complicated, you might want to create a set that is the index to the desired plots. For example

demo <- site$elev>7500 & site$slope<10 & site$av>0.5
vegtab(veg,demo,min=20)

        carrss symore senmul berrep arcpat ceamar pacmyr erirac artarb tetcan
bcnp__1    0.5    1.0    0.0    1.0    1.0    0.5    0.0    0.0    0.0    0.0
bcnp__2    0.5    0.5    0.5    0.0    0.5    0.0    0.5    0.0    0.0    0.0
bcnp__3    0.5    0.5    0.5    0.5    1.0    0.5    0.5    0.0    0.0    0.0
bcnp__4    0.5    0.5    0.5    1.0    1.0    0.5    0.0    0.0    0.0    0.0
bcnp__6    0.5    0.5    0.0    1.0    1.0    1.0    0.5    0.0    0.0    0.0

    .       .      .      .      .      .      .      .      .      .      .
    .       .      .      .      .      .      .      .      .      .      .
    .       .      .      .      .      .      .      .      .      .      .
bcnp138    1.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp139    0.0    0.0    0.5    0.0    0.0    0.0    0.0    0.0    0.0    2.0
bcnp140    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    0.5    1.0
bcnp141    2.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0

You can use nested ands and ors as well (remember, & = "and" and | = "or").

demo <- (site$elev>7500 & site$slope<10) | site$pos=="low_slope"
vegtab(veg,demo,pltord=elev,min=35)

        carrss symore senmul arcpat berrep ceamar purtri pacmyr juncom artarb
bcnp_89    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_88    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp157    0.0    0.5    0.0    0.0    0.0    0.0    0.5    0.0    0.0    0.0
bcnp_76    1.0    0.0    0.0    3.0    0.5    0.0    0.0    0.0    0.0    0.0
bcnp_84    1.0    0.5    0.0    0.5    0.5    0.0    0.5    0.0    0.0    0.0
bcnp_77    1.0    0.0    0.0    4.0    0.5    1.0    0.0    0.0    0.0    0.0
     .      .      .      .      .      .      .      .      .      .      . 
     .      .      .      .      .      .      .      .      .      .      .

So far, we have seen how to subset the data by plots and re-order the table by rows. How do we re-order the species?

Ordering Species in the Table

To order the species, we need a vector of length equal to the number of species. For example, back in Lab 1 we created a vector of the number of plots a species occurs in, called spc.pres. We can order the table by that, which is the reverse of the default ordering.

vegtab(veg,spcord=spc.pres,min=40)

          chrvis oryhym sithys pacmyr ceamar purtri berrep arcpat senmul symore
bcnp__1      0.0    0.0    0.0    0.0    0.5    0.0    1.0    1.0    0.0    1.0
bcnp__2      0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.5    0.5    0.5
bcnp__3      0.0    0.0    0.0    0.5    0.5    0.0    0.5    1.0    0.5    0.5
bcnp__4      0.0    0.0    0.0    0.0    0.5    0.0    1.0    1.0    0.5    0.5
bcnp__5      0.0    0.0    0.0    0.5    0.5    0.5    0.5    4.0    0.5    0.5
      .       .      .      .      .      .      .      .      .      .      . 
      .       .      .      .      .      .      .      .      .      .      . 
      .       .      .      .      .      .      .      .      .      .      . 
bcnp156    0.0    0.0    0.0    0.0    0.0    0.5    0.0    0.5    0.5    0.5
bcnp157    0.0    0.0    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.5
bcnp158    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    0.0    0.5
bcnp159    0.0    0.5    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.0
bcnp160    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    0.0    0.0

Much more interesting is to order the species by some ecological factor. To do this, we'll use the wascores() function from the vegan library. wascores stands for "weighted average scores", and is used to compute the weighted average value of some environmental variable for all plots in which a species occurs, weighted by the species abundance. For example,

spc.elev <- wascores(elev,veg)
spc.elev

           [,1]
junost 7335.000
ameuta 8138.462
arcpat 8126.099
arttri 6804.082
atrcan 6792.812
berfre 6928.889
ceamar 8200.122
cerled 7776.667
  .        . 
  .        . 
  .        . 
strcor 7140.000
swerad 8216.400
taroff 7866.000
thafen 8245.000
towmin 7450.000
tradub 7377.692
valacu 8150.000
vicame 8360.000

As you can see, every species has an elevation associated with it, and we can use the rank order of elevation as the ordering variable.

vegtab(veg,spcord=spc.elev,min=40)

    sithys oryhym chrvis purtri carrss arcpat senmul ceamar symore berrep
bcnp__1      0.0    0.0    0.0    0.0    0.5    1.0    0.0    0.5    1.0    1.0
bcnp__2      0.0    0.0    0.0    0.0    0.5    0.5    0.5    0.0    0.5    0.0
bcnp__3      0.0    0.0    0.0    0.0    0.5    1.0    0.5    0.5    0.5    0.5
bcnp__4      0.0    0.0    0.0    0.0    0.5    1.0    0.5    0.5    0.5    1.0
bcnp__5      0.0    0.0    0.0    0.5    1.0    4.0    0.5    0.5    0.5    0.5
      .       .      .      .      .      .      .      .      .      .      . 
      .       .      .      .      .      .      .      .      .      .      . 
bcnp156    0.0    0.0    0.0    0.5    0.0    0.5    0.5    0.0    0.5    0.0
bcnp157    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.0    0.5    0.0
bcnp158    0.0    0.0    0.0    0.0    0.0    0.5    0.0    0.0    0.5    0.0
bcnp159    0.0    0.5    0.0    0.5    0.0    0.0    0.0    0.0    0.0    0.0
bcnp160    0.0    0.0    0.0    0.0    0.0    0.5    0.0    0.0    0.0    0.0

The species are listed in order from low elevation to high for those species which occur at least 40 times. We can order the plots by the same variable and get a sorted listing from low elevation to high in both plot locations and species modes.

vegtab(veg,pltord=elev,spcord=spc.elev,min=40)

         sithys oryhym chrvis purtri carrss arcpat senmul ceamar symore berrep
bcnp_89     2.0    0.5    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_90     5.0    0.5    0.5    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_99     1.0    0.5    0.5    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp100     0.5    0.5    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp101     0.5    0.5    0.5    0.0    0.0    0.0    0.0    0.0    0.0    0.0
      .      .      .      .      .      .      .      .      .      .      . 
      .      .      .      .      .      .      .      .      .      .      . 
      .      .      .      .      .      .      .      .      .      .      . 
bcnp_46     0.0    0.0    0.0    1.0    0.5    0.5    0.5    1.0    0.5    0.5
bcnp_58     0.5    0.5    0.0    0.5    0.5    4.0    0.5    0.0    0.5    0.0
bcnp_61     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    1.0
bcnp_87     0.5    0.5    0.0    0.5    1.0    3.0    0.5    1.0    1.0    1.0
bcnp_60     0.0    0.0    0.0    0.0    0.5    0.5    0.5    0.0    1.0    1.0

Notice how the non-zero values go along the diagonal.

If, for some reason, you want to reverse the order of plots or species, simply put a negative sign in front of the order variable. For example,

vegtab(veg,pltord=-elev,spcord=-spc.elev,min=40)

         pacmyr berrep symore ceamar senmul arcpat carrss purtri chrvis oryhym
bcnp_60     0.5    1.0    1.0    0.0    0.5    0.5    0.5    0.0    0.0    0.0
bcnp_87     0.5    1.0    1.0    1.0    0.5    3.0    1.0    0.5    0.0    0.5
bcnp_61     0.5    1.0    0.5    0.0    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_58     0.5    0.0    0.5    0.0    0.5    4.0    0.5    0.5    0.0    0.5
bcnp_46     0.0    0.5    0.5    1.0    0.5    0.5    0.5    1.0    0.0    0.0
      .      .      .      .      .      .      .      .      .      .      . 
      .      .      .      .      .      .      .      .      .      .      . 
      .      .      .      .      .      .      .      .      .      .      . 
bcnp100     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5
bcnp101     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    0.5
bcnp_99     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    0.5
bcnp_90     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    0.5
bcnp_89     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5

The order here is not a perfect mirror order (e.g. 100,101) due to ties in elevation.

We haven't introduced ordination yet (that comes in later labs). I'll assume, however, that you know what ordination is, and demonstrate ordering plots and species by ordination scores. In Lab 10 we'll study correspondence analysis (CA), which ordinates both plots and species centroids on the same axes. The syntax here will look unfamiliar, but assume that we have stored the CA results for axis I in ca.plt and ca.spc.

vegtab(veg,pltord=ca.plt,spcord=ca.spc,min=40)

         pacmyr berrep arcpat symore ceamar carrss purtri senmul oryhym sithys
bcnp__9     0.5    1.0    0.0    0.5    0.0    0.5    0.0    0.0    0.0    0.0
bcnp_26     0.5    0.5    1.0    1.0    0.5    0.0    0.0    0.0    0.0    0.0
bcnp_22     0.0    0.5    0.5    0.5    0.0    0.0    0.0    0.0    0.0    0.0
bcnp_13     0.5    0.5    0.5    0.5    0.5    0.5    0.0    0.0    0.0    0.0
bcnp_45     0.5    0.5    4.0    0.5    0.5    0.5    0.5    0.0    0.0    0.0
      .      .      .      .      .      .      .      .      .      .      . 
      .      .      .      .      .      .      .      .      .      .      . 
      .      .      .      .      .      .      .      .      .      .      . 
bcnp_92     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    1.0
bcnp_89     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    2.0
bcnp_97     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    1.0
bcnp_96     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    1.0
bcnp_95     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.5    0.5
bcnp_94     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0

Again, as we would expect, the non-zero values go along the diagonal. The species order is not too different from high to low elevation, but the plot order has changed quite a bit. With that teaser, we'll leave vegetation tables for now, and go to vegetation constancy, importance, and discrimination tables.

Constancy and Importance Tables

Even with the ability to subset and order vegetation tables, there is often simply too much information in such a table to visualize and understand. In addition, we often want to compare results from multiple tables, which is difficult with vegtab().

A simple alternative to summarize vegetation tables is a constancy table, which lists for each species the relative frequency of occurrence in a specified sets. Accordingly, the LabDSV library includes a function called const() which behaves similarly to vegtab(). For example

const(veg,min=0.1)

        res
ameuta 0.14
arcpat 0.46
ceamar 0.36
cermon 0.20
chrvis 0.25
juncom 0.24
pacmyr 0.30

   .    .
   .    .
   .    .
opueri 0.10
pedcan 0.21
pencae 0.12
senmul 0.46
swerad 0.15

const(veg,min=0.1) lists the relative frequency of all species which occur with a minimum constancy of 0.1 (10% of the sample plots).

Of much more interest, of course, is to analyze specific variables for effect on the relative frequencies. For example, to look at topographic position,

const(veg,site$pos,min=0.1)

       bottom low_slope mid_slope ridge up_slope
ameuta   0.10      0.12      0.22   .       0.11
arcpat   0.20      0.45      0.54  0.67     0.40
arttri    .         .        0.11   .       0.14
atrcan    .         .        0.11   .       0.14
ceamar   0.10      0.30      0.54  0.56     0.23
  .       .         .         .     .        .
  .       .         .         .     .        .
  .       .         .         .     .        .
sphcoc    .         .        0.15   .        .  
steten    .         .         .    0.11      .  
swerad    .        0.18      0.20  0.22     0.11
taroff   0.15       .         .     .        .  
tradub   0.10       .         .     .        .

In this case, min=0.1 means that a species must achieve a relative frequency of 0.1 in at least one column. The dots replace small values to emphasize more important data. You can specify the threshold with a show=0 or other small value. E.g.

const(veg,site$pos,min=0.1,show=0)

       bottom low_slope mid_slope ridge up_slope
ameuta   0.10      0.12      0.22  0.06     0.11
arcpat   0.20      0.45      0.54  0.67     0.40
arttri   0.00      0.09      0.11  0.06     0.14
atrcan   0.00      0.06      0.11  0.00     0.14
ceamar   0.10      0.30      0.54  0.56     0.23
 .        .         .         .     .        .
 .        .         .         .     .        .
 .        .         .         .     .        .

Constancy is not sensitive to abundance, but rather uses only presence/absence as its criterion. To look at the relative importance of a species, you can use a function called importance() that calculates the average value of a species when it's present (not averaging in the zeroes) in a set. For example,

importance(veg,site$pos,min=0.1)

         bottom low_slope mid_slope ridge up_slope
junost      .         .        0.50   .       0.50
ameuta     0.50      0.50      0.62  0.50     0.50
arcpat     2.38      2.41      2.91  2.00     1.82
arttri      .        1.83      2.17  0.50     1.10
atrcan      .        1.00      1.50   .       1.00
    .       .         .         .     .        .
taroff     0.50      0.50       .     .        .  
thafen      .        0.50      0.50   .        .  
towmin      .        0.50       .     .        .  
tradub     0.50      0.50      1.00   .       0.50
valacu      .        0.50       .     .        .  
vicame      .         .        0.50   .        .

Again, you can include a show=0 argument to see all the values. Because the matrix veg has values between 0 and 6, the "importance" calculated is an average cover code. Way back in Lab 1 were created a matrix of cover midpoints. To get an average cover by species we can use that matrix instead.

importance(cover,site$pos,min=1.0)

       bottom low_slope mid_slope ridge up_slope
ameuta    .         .        1.71   .        .  
arcpat  28.88     32.41     41.26 23.58    19.71
arttri    .       17.67     18.75   .       8.40
atrcan    .        3.00      9.00   .       4.40
ceamar    .        1.50      1.60  2.45     2.62
cerled    .        1.33       .     .        .  
cermon   1.50      1.75      1.57  4.12     3.78
chrpar    .        8.38       .     .        .  
chrvis   2.96      1.25      2.00   .       1.50
  .       .         .         .     .        . 
  .       .         .         .     .        . 
  .       .         .         .     .        .

Here I asked for a minimum of 1.0, which would be an average cover of at least one percent, so few species qualified. Cover class midpoints almost always result in an overestimate of mean cover due to uneven distributions within classes, but that's a more complicated issue than we want to get into here.

Somewhat like vegtab, we can order the species in the table with the spcord= argument.

importance(cover,site$pos,min=1.0,spcord=spc.pres)

       bottom low_slope mid_slope ridge up_slope
carrss  13.67      1.09       .     .        .  
symore    .         .         .    1.46     1.12
arcpat  28.88     32.41     41.26 23.58    19.71
berrep    .        1.12      1.31  1.44     3.63
purtri   1.12       .        3.36  3.96     2.50
ceamar    .        1.50      1.60  2.45     2.62
sithys    .        2.39      4.50   .       2.39
chrvis   2.96      1.25      2.00   .       1.50
  .       .         .         .     .        .
  .       .         .         .     .        .
  .       .         .         .     .        .

Finally, we can combine the constancy and mean abundance estimates in a single table with the concov() command.

concov(cover,site$pos)

          bottom low_slope mid_slope     ridge  up_slope
ameuta 10(  0.5) 12(  0.5) 22(  1.7)  6(  0.5) 11(  0.5)
arcpat 20( 28.9) 45( 31.8) 54( 40.8) 67( 23.6) 40( 19.7)
arttri  0(  0.0)  9( 17.7) 11( 18.8)  6(  0.5) 14(  8.4)
atrcan  0(  0.0)  6(  3.0) 11(  9.0)  0(  0.0) 14(  4.4)
ceamar 10(  0.5) 30(  1.5) 54(  1.6) 56(  2.5) 23(  2.6)
cermon 25(  1.5) 24(  1.8) 13(  1.6) 22(  4.1) 26(  3.8)
  .         .         .         .         .         .
  .         .         .         .         .         .
  .         .         .         .         .         .
steten  0(  0.0)  0(  0.0)  6(  0.5) 11(  0.5)  0(  0.0)
swerad  0(  0.0) 18(  0.5) 20(  0.5) 22(  0.5) 11(  0.5)
taroff 15(  0.5)  6(  0.5)  0(  0.0)  0(  0.0)  0(  0.0)
tradub 10(  0.5)  9(  0.5)  6(  5.3)  0(  0.0)  6(  0.5)

where the number in front of the parentheses is the constancy, and the number inside of the parentheses is the typical cover (not averaging in zeros). If you want the true mean cover (averaging in the zeros) you can say typical=FALSE

concov(cover,site$pos,typical=FALSE)

          bottom low_slope mid_slope     ridge  up_slope
ameuta 10(  0.1) 12(  0.1) 22(  0.4)  6(  0.0) 11(  0.1)
arcpat 20(  5.8) 45( 14.5) 54( 21.9) 67( 15.7) 40(  7.9)
arttri  0(  0.0)  9(  1.6) 11(  2.1)  6(  0.0) 14(  1.2)
atrcan  0(  0.0)  6(  0.2) 11(  1.0)  0(  0.0) 14(  0.6)
ceamar 10(  0.1) 30(  0.5) 54(  0.9) 56(  1.4) 23(  0.6)
cermon 25(  0.4) 24(  0.4) 13(  0.2) 22(  0.9) 26(  1.0)
  .         .         .         .         .         .
  .         .         .         .         .         .
  .         .         .         .         .         .
swerad  0(  0.0) 18(  0.1) 20(  0.1) 22(  0.1) 11(  0.1)
taroff 15(  0.1)  6(  0.0)  0(  0.0)  0(  0.0)  0(  0.0)
tradub 10(  0.1)  9(  0.1)  6(  0.3)  0(  0.0)  6(  0.0)

Species Discrimination Tables

Part of the objective of constancy and importance tables is to look for species which have high constancy or importance in some sets but low in others. You can find these by eye in the tables above (e.g. broano in the importance tables), but it's easier to find them algorithmically using the function spcdisc(). spcdisc calculates the total constancy or importance of a species in a table, and then calculates a statistic based on the Shannon-Weiner index. Species which are uniform across sets have a discrimination of 0.0, and species with all their constancy or importance in a single set have a discrimination of 1.0, and all other species scale in between depending on the concentration of their distribution among types.

To use spcdisc(), you give it the output from const() or importance(). For example

tmp <- const(veg,site$pos,min=0.1,show=0) spcdisc(tmp)

     ameuta      arcpat      arttri      atrcan      ceamar      cermon 
0.054875858 0.040087755 0.165765434 0.351017747 0.089487622 0.015749729 
     chrpar      chrvis      juncom      pacmyr      purtri      quegam 
0.365164787 0.143840862 0.192028393 0.188715185 0.047938940 0.069995689 
     ribcer      sherot      symore      artarb      artfri      berrep 
0.158030413 0.666989397 0.044331310 0.092120032 0.472435047 0.064033580 
       .           .           .           .            .          .
       .           .           .           .            .          .
       .           .           .           .            .          .
     potgra      senmul      sphcoc      steten      swerad      taroff 
0.427443939 0.007709231 0.209069896 0.596599316 0.157004684 0.628274204 
     tradub 
0.155152498

Of the abbreviated list given here, steten and taroff (dandelion) have the highest discrimination. To highlight the best species, use the sort=TRUE argument.

spcdisc(tmp,sort=TRUE)

     irimis      muhric      ericor      junbal      lupkin      astagr 
1.000000000 1.000000000 1.000000000 0.810718950 0.720050581 0.720050581 
     sherot      poapra      broine      taroff      phllon      steten 
0.666989397 0.664352295 0.664352295 0.628274204 0.628274204 0.596599316 
     erifla      poanev      agogla      artfri      potcri      broano 
0.528446070 0.492422067 0.491354331 0.472435047 0.465122958 0.432135526 
       .           .          .            .           .           .
       .           .          .            .           .           .
       .           .          .            .           .           .
     poafen      achmil      carrss      eupfen      cermon      lotuta 
0.027460062 0.026086421 0.023272511 0.020931561 0.015749729 0.008799344 
     senmul 
0.007709231

The first three species are perfect discriminators of topographic position, followed by several other good ones. The bottom of the list are either ubiquitous species (e.g. carrss and senmul) or indiscriminant species.

We haven't covered cluster analysis yet (we'll get to it about Lab 13), but I'm sure you can imagine using such a function to determine the diagnostic species of clusters. For example, assume that clusid contains the cluster membership of each plot in one of ten clusters. Then,

spcdisc(const(veg,clusid,min=0.1,show=0),sort=TRUE)

   taroff    swerad    potgra    potcri    phllon    oenfla    oencor    molpar 
1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 
   lyggra    lupser    lupkin    litmul    leppun    genaff    euplur    eriumb 
1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 
   dessop    compal    cirneo    astmis    astmeg    asthum    astcon    astchi 
1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 
    .          .         .         .         .          .         .        .
    .          .         .         .         .          .         .        .
    .          .         .         .         .          .         .        .
   linlew    koenit    sithys    senmul    berfre    bougra    carrss    eriala 
0.5124234 0.4995058 0.4941893 0.4897064 0.4874319 0.4814601 0.4753900 0.4624004 
   purtri    poafen    chrvis    caslin    elysal    eupfen    hymric    oryhym 
0.4613622 0.4575671 0.4301827 0.4237380 0.4121540 0.4006111 0.2803659 0.2121789

gives us a list of diagnostic species. Notice how I embedded the call to const() inside of spcdisc() rather than saving it as a temporary object.

Using Vegetation Tables

As you can see in this lab, you can use vegetation tables directly to study vegetation, and there is a long history of such work in vegetation ecology. It seems to me, however, that its greatest utility is in determining the details of results from quantitative analyses we will cover in the rest of the labs. Accordingly, we'll leave vegetation tables for a while, but you'll see their use again in subsequent labs.

Functions Used in This Lab

vegtab()

vegtab <- function (comm, set, minval = 1, pltord, 
                  spcord, pltlbl, trans = FALSE) 
{
    if (missing(set)) {
        set <- seq(1:nrow(comm))
    }
    else {
        set <- seq(1:nrow(comm))[set]
        set <- set[!is.na(set)]
    }
    tmp <- comm[set, ]
    spcidx <- apply(tmp > 0, 2, sum)
    tmp <- tmp[, spcidx >= minval]
    if (missing(pltord)) {
        pltord <- seq(1:nrow(tmp))
    }
    else {
        pltord <- pltord[set]
    }
    if (missing(spcord)) {
        x <- apply(tmp > 0, 2, sum)
        y <- apply(tmp, 2, sum)
        spcord <- -1 * (x + y/max(y))
    }
    else {
        spcord <- spcord[spcidx >= minval]
    }
    if (!missing(pltlbl)) {
        if (is.numeric(pltlbl)) {
            tmp <- cbind(pltlbl[set], tmp)
            dimnames(tmp)[[2]][1] <- deparse(substitute(pltlbl))
            spcord <- c(min(spcord) - 1, spcord)
        }
        else {
            dimnames(tmp)[[1]] <- pltlbl
        }
    }
    tmp <- tmp[order(pltord), order(spcord)]
    if (trans == TRUE) {
        tmp <- t(tmp)
    }
    tmp
}

const()

const
function (comm, clustering, minval = 0, show = minval, digits = 2, 
    sort = FALSE, spcord = NULL) 
{
    if (missing(clustering)) {
        const <- apply(comm > 0, 2, sum)/nrow(comm)
        const <- const[const >= minval]
        const <- data.frame(const)
        names(const) <- deparse(substitute(comm))
        return(round(const, digits))
    }
    else if (is.logical(clustering)) {
        comm <- comm[clustering, ]
        const <- apply(comm > 0, 2, sum)/nrow(comm)
        const <- const[const >= minval]
        const <- data.frame(const)
        names(const) <- deparse(substitute(clustering))
        return(round(const, digits))
    }
    clustering <- clustify(clustering)
    if (length(table(clustering)) == 1) {
        const <- apply(comm > 0, 2, sum)/nrow(comm)
        const <- const[const >= minval]
        const <- data.frame(const)
        return(round(const, digits))
    }
    else {
       res <- matrix(0, nrow = ncol(comm), ncol = length(levels(clustering)))
        x <- apply(comm, 2, function(x) {
            tapply(x > 0, clustering, sum)
        })
        y <- as.numeric(table(clustering))
        res <- x/y
        keep <- as.logical(apply(res, 2, max) >= minval)
        res <- res[, keep]
        tmp <- as.data.frame(t(res))
        row.names(tmp) <- names(comm)[keep]
        if (!is.null(spcord)) {
            tmp <- tmp[rev(order(spcord[keep])), ]
        }
        tmp <- format(round(tmp, digits = digits))
        tmp[tmp < show] <- substring(" .  ", 1, digits + 2)
        names(tmp) <- attr(clustering, "levels")
        attr(tmp, "call") <- match.call()
        attr(tmp, "timestamp") <- date()
        if (sort) {
            print(tmp)
            repeat {
                plots <- readline(" enter the species: ")
                if (plots == "") {
                  break
                }
                else {
                  pnt <- readline(" in front of        : ")
                }
                for (i in (strsplit(plots, ",")[[1]])) {
                  ord <- 1:nrow(tmp)
                  x <- match(i, row.names(tmp))
                  if (!is.na(x)) {
                    z <- ord[x]
                    ord <- ord[-x]
                    y <- match(pnt, row.names(tmp))
                    if (!is.na(y)) {
                      if (y > 1) {
                        first <- ord[1:(y - 1)]
                        last <- ord[y:length(ord)]
                        ord <- c(first, z, last)
                      }
                      else {
                        last <- ord[y:length(ord)]
                        ord <- c(z, last)
                      }
                      tmp <- tmp[ord, ]
                      print(tmp)
                    }
                    else {
                      print(paste("species", pnt, "does not exist"))
                    }
                  }
                  else {
                    print(paste("species", i, "does not exist"))
                  }
                }
            }
            return(tmp)
        }
    }
    return(tmp)
}

importance()

importance <- function (comm, clustering, minval = 0, digits = 2, 
    show = minval, sort = FALSE, typical = TRUE, spcord, dots = TRUE) 
{
    if (missing(clustering)) {
        impt <- apply(comm, 2, sum)/nrow(comm)
        impt <- impt[impt >= minval]
        impt <- data.frame(impt)
        names(impt) <- deparse(substitute(comm))
        return(round(impt, digits))
    }
    else if (is.logical(clustering)) {
        comm <- comm[clustering, ]
        impt <- apply(comm, 2, sum)/nrow(comm)
        impt <- impt[impt >= minval]
        impt <- data.frame(impt)
        names(impt) <- deparse(substitute(clustering))
        return(round(impt, digits))
    }
    clustering <- clustify(clustering)
    if (length(table(clustering)) == 1) {
        impt <- apply(comm, 2, sum)/nrow(comm)
        impt <- impt[impt >= minval]
        impt <- data.frame(impt)
        return(round(impt, digits))
    }
    else {
        res <- matrix(0, nrow = ncol(comm), ncol = length(levels(clustering)))
        x <- apply(comm, 2, function(x) {
            tapply(x, clustering, sum)
        })
        if (typical) {
            y <- apply(comm, 2, function(x) {
                tapply(x > 0, clustering, sum)
            })
        }
        else {
            y <- apply(comm, 2, function(x) {
                tapply(x >= 0, clustering, sum)
            })
        }
        y[x == 0] <- 1
        res <- x/y
        keep <- as.logical(apply(res, 2, max) >= minval)
        res <- res[, keep]
        tmp <- as.data.frame(t(res))
        row.names(tmp) <- names(comm)[keep]
        if (!missing(spcord)) {
            tmp <- tmp[rev(order(spcord[keep])), ]
        }
        if (dots) {
            tmpx <- format(round(tmp, digits = digits))
            tmpx[tmp < show] <- substring(" .  ", 1, nchar(tmpx[1, 
                1]))
            print(tmpx)
        }
        if (sort) {
            cat("\nConstancy Table\n\n")
            print(tmp)
            repeat {
                plots <- readline(" enter the species: ")
                if (plots == "") {
                  break
                }
                else {
                  pnt <- readline(" in front of        : ")
                }
                for (i in (strsplit(plots, ",")[[1]])) {
                  ord <- 1:nrow(tmp)
                  x <- match(i, row.names(tmp))
                  if (!is.na(x)) {
                    ord <- ord[-x]
                    y <- match(pnt, row.names(tmp[ord, ]))
                    if (!is.na(y)) {
                      if (y == 1) {
                        ord <- c(x, ord)
                      }
                      else {
                        first <- ord[1:(y - 1)]
                        last <- ord[y:length(ord)]
                        ord <- c(first, x, last)
                      }
                      tmp <- tmp[ord, ]
                      print(tmp)
                    }
                    else {
                      print(paste("species", pnt, "does not exist"))
                    }
                  }
                  else {
                    print(paste("species", i, "does not exist"))
                  }
                }
            }
            attr(tmp, "call") <- match.call()
            attr(tmp, "comm") <- deparse(substitute(comm))
            return(tmp)
        }
    }
}

concov()

concov <- function (comm, clustering, digits = 1, width = 5, 
                 typical = TRUE, thresh = 10) 
{
    if (missing(clustering)) {
        const <- apply(comm > 0, 2, sum)/nrow(comm)
        keep <- const >= thresh/100
        impt <- apply(comm, 2, sum)/nrow(comm)
        a <- formatC(as.numeric(const) * 100, width = 2, format = "d")
        b <- formatC(as.numeric(impt), width = width, digits = digits, 
            format = "f")
        tmp <- NULL
        tmp <- cbind(tmp, paste(a, "(", b, ")", sep = ""))
        tmp <- tmp[keep]
        tmp <- data.frame(tmp)
        row.names(tmp) <- names(comm)[keep]
        attr(tmp, "call") <- match.call()
        attr(tmp, "comm") <- deparse(substitute(comm))
        attr(tmp, "timestamp") <- date()
        return(tmp)
    }
    else if (is.logical(clustering)) {
        comm <- comm[clustering, ]
        comm <- comm[, apply(comm > 0, 2, sum) > 0]
        x <- apply(comm > 0, 2, sum)
        y <- apply(comm, 2, sum)/x
        x <- x/nrow(comm)
        keep <- apply(as.matrix(x), 1, max) >= thresh/100
        a <- formatC(as.numeric(x) * 100, width = 2, format = "d")
        b <- formatC(as.numeric(y), width = width, digits = digits, 
            format = "f")
        tmp <- NULL
        tmp <- cbind(tmp, paste(a, "(", b, ")", sep = ""))
        tmp <- tmp[keep]
        tmp <- data.frame(tmp)
        row.names(tmp) <- names(comm)[keep]
        names(tmp) <- deparse(substitute(clustering))
        attr(tmp, "call") <- match.call()
        attr(tmp, "comm") <- deparse(substitute(comm))
        attr(tmp, "clustering") <- clustering
        attr(tmp, "timestamp") <- date()
        return(tmp)
    }
    clustering <- clustify(clustering)
    if (length(table(clustering)) == 1) {
        const <- apply(comm > 0, 2, sum)/nrow(comm)
        keep <- const >= thresh/100
        impt <- apply(comm, 2, sum)/nrow(comm)
        a <- formatC(as.numeric(const) * 100, width = 2, format = "d")
        b <- formatC(as.numeric(impt), width = width, digits = digits, 
            format = "f")
        tmp <- NULL
        tmp <- cbind(tmp, paste(a, "(", b, ")", sep = ""))
        tmp <- tmp[keep]
        tmp <- data.frame(tmp)
        row.names(tmp) <- names(comm)[keep]
        names(tmp) <- deparse(substitute(clustering))
        attr(tmp, "call") <- match.call()
        attr(tmp, "comm") <- deparse(substitute(comm))
        attr(tmp, "clustering") <- clustering
        attr(tmp, "timestamp") <- date()
        return(tmp)
    }
    else {
        levels <- levels(clustering)
        clustering <- as.integer(clustering)
        counts <- table(clustering)
        const <- t(apply(comm>0,2,function(x){tapply(x,clustering,sum)/counts}))
        cases <- t(apply(comm>0,2,function(x){tapply(x,clustering,sum)}))
        cases[cases==0] <- 1
        if (typical) {
          cov <- t(apply(comm,2,function(x){tapply(x,clustering,sum)}))
          cov <- cov / cases
        } else {
          cov <- t(apply(comm,2,function(x){tapply(x,clustering,mean)}))
        }
        tmp <- NULL
        keep <- apply(as.matrix(const), 1, max) >= thresh/100
        for (i in 1:length(table(clustering))) {
            a <- formatC(as.numeric(const[, i]) * 100, width = 2, 
                format = "d")
            b <- formatC(as.numeric(cov[, i]), width = width, digits = digits, 
                format = "f")
            tmp <- cbind(tmp, paste(a, "(", b, ")", sep = ""))
        }
        tmp <- tmp[keep, ]
        tmp <- data.frame(tmp)
        row.names(tmp) <- names(comm)[keep]
        names(tmp) <- levels
        attr(tmp, "call") <- match.call()
        attr(tmp, "comm") <- deparse(substitute(comm))
        attr(tmp, "clustering") <- clustering
        attr(tmp, "timestamp") <- date()
        return(tmp)
    }
    tmp
}

spcdisc()

spcdisc <- function (x, sort = FALSE) 
{
    shannon <- function(y) {
        y <- as.numeric(y)
        frac <- y/sum(y)
        comp <- sum((-1 * frac * log(frac))[frac > 0])
        comp <- 1 - (comp/log(length(y)))
        comp
    }
    tmp <- apply(x, 1, shannon)
    if (sort) 
        tmp <- tmp[rev(order(tmp))]
    tmp
}