helixcn/spaa
SPecies Association Analysis

Package index
Search the helixcn/spaa package
Vignettes
Functions
29
Source code
12
Man pages
16
Home
/
GitHub
/
helixcn/spaa
/

In helixcn/spaa: SPecies Association Analysis 

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This tutorial is based on a post (in Chinese) available at: http://blog.sciencenet.cn/blog-255662-1239853.html

Ecological niche refers to resource allocation (Jintun Zhang, 2004), and is important in many aspects, for example, ordination analysis, community assembly etc. At the community scale, niche is usually quantified by the occurrence and co-occurrence of species. Although a great number of studies have focused on this topic in the past three decades ((in particular in the Chinese literature, e.g., C. W. Sun and S. D. Zhao, 1996; C. H. Guo et al. 1997; Z. M. Shi et al. 2001; J. L. Li and J. T. Zhang, 2010; N. J. Wang et al. 2010), a computer programme for conducting the relevant analysis is lacking. Here we describe spaa, which stands for SPecies Association Analysis, an R package for conducting suche analysis.

1. The installation

spaa is an open source package written in R, available on CRAN.

To install, enter the command install.packages("spaa") in the R console, R will ask you to select a CRAN mirror, and the package will be automatically downloaded and installed.

The developing version is available at https://github.com/helixcn/spaa, to install it, type:

library(devtools)
install_github("helixcn/spaa")

2. Functions and datasets in spaa

2.1 Main Functions

2.2 Embedded datasets

3. Calculating niche width

3.1 Converting community data into a matrix

Community data generally needs to be converted to a matrix (Table 1) when calculating diversity indices, species association, niche width, niche overlap, in which each row represents a plot (site, region or track), each column represents a species, and the number indicates the number of individuals.

Table 1 Community species matrix

| | sp1 | sp2 | sp3 | sp4 | sp5 | sp6 | sp7 | |-------|-----|-----|-----|-----|-----|-----|-----| | plot1 | 3 | 6 | 1 | 2 | 1 | 0 | 0 | | plot2 | 8 | 0 | 30 | 0 | 0 | 0 | 0 | | plot3 | 0 | 1 | 0 | 2 | 0 | 1 | 3 |

However, field data are generally recorded in two ways:

  1. Individual-based Field Record Sheet
  2. Plot-based Record Sheet

3.1.2 Individual-based Field Record Sheet

Individual-based Field Record Sheet is often used in forest plots census, in which every tree (usually when DBH >= 1cm) is tagged. Usually, for a woody species, both the main trunk and branches are recorded. The following columns are often included:

  1. plot: indicates the name of the plot (20mX20m) the tree is located at. The first two digits "01" refers to the first 20mX20m from the origin (west to east), the second "01" refers to the first 20mX20m from the origin along the Y axist (south to north)
  2. tag number: the number on the main trunk
  3. species name: the species name of the individual
  4. x-coordinate: the x coordinate of the individual (usually in a local Cartesian System, the unit is often meters easting from the origin)
  5. y-coordinate: the y coordinate of the individual (usually in a local Cartesian System, the unit is often meters northing from the origin)
  6. branch (0 refers to the main trunk, 1 is the second largest branch, 2 is the third largest branch and so on)
  7. diameter at breast height (DBH, usually in cm)
  8. height (usually in meters)
  9. status (alive, dead but standing, collapsed, etc.)
  10. remark, any thing to markdown

Table 1. Field record sheet

| plot | tag | species | x | y | branch | dbh | height | status | remark | |------|------------|---------|-----|------|--------|-----|--------|--------|--------| | 0101 | btm0101001 | sp1 | 1.5 | 5.8 | 0 | 15 | 12 | alive | tree | | 0101 | btm0101002 | sp2 | 3.3 | 12.1 | 0 | 12 | 10 | alive | tree |

... ...

Generating an example community matrix

plot <- c("0101", "0101", "0102")
tag <- c("btm0101001", "btm0101002", "btm0101003" )
species <- c("sp1", "sp2", "sp1")
x <- c(1.5, 3.3, 4.1)
y <- c(5.8, 12.1, 8.9)
branch <- c(0, 0, 0)
dbh  <- c(15, 12, 5)
height   <- c(12, 10, 8)
status  <- c("alive", "alive", "alive")
remark  <- c("", "", "")

btmdata <- data.frame(plot, tag, species, x, y, branch, dbh, height, status, remark)

btmdata

# generate a dataset
btm <- table(btmdata$plot, btmdata$species)
btm

3.1.2 Plot-based Record Sheet

Plot-based Record is much simpler, generally the first column is the name of plot, the second column is species, and the third column is the number of individuals of the species in the sample. Note the value is usually an integer, but could also be a decimal number when representing importance value (Kent M, 2011, page 175).

Table 2 Sample survey form

| plot | species | abundance | |-------|---------|-----------| | plot1 | sp1 | 3 | | plot1 | sp2 | 6 | | plot1 | sp3 | 1 | | plot1 | sp4 | 2 | | plot1 | sp5 | 1 | | plot2 | sp1 | 8 | | plot2 | sp3 | 30 | | plot3 | sp4 | 2 | | plot3 | sp2 | 1 | | plot3 | sp6 | 1 | | plot3 | sp7 | 3 |

The following code converts sample data to a community-species matrix.

library(reshape2)
library(spaa)
data(datasample)

datasample2 <- cbind(plot = row.names(datasample), datasample)

# Convert the modified matrix to the common format of data records
aaa <- melt(data = datasample2,
            id.vars = "plot",
            measure.vars = 2:ncol(datasample2),
            variable.name = "species",
            value.name = "abundance",
            factorsAsStrings = TRUE)

# is 0.00 because the species did not occur in some samples
aaa1 <- aaa[aaa$abundance > 0, ]
head(aaa1)

bbb <- acast(aaa1, 
             formula = plot ~ species, 
             value.var = "abundance", 
             fill = 0)
bbb

When the abundance column only contains integers, it can be converted using data2mat(). Note that the names of the columns should be named "plotname", "species", "abundance" (note R is case sensitive).

library(spaa)
data(testdata)
testdata
data2mat(testdata)

3.2 Calculating species associations

Interspecific associations, including association among all the species and specific association between each pair of species could be computed using spaa (Zhou Xianye et al. 2000; Shi Zuomin et al. 2001; Zhang Siyu and Zheng Shiqun 2002; Zhang Zhiyong et al. 2003; Kang Bing et al. 2005; Wang Wenjin et al. 2007; Wang Naijiang et al. 2010).

3.2.1 Species association of the community as a whole

The function sp.assoc() computes the following indices and the output is a list with the following elements.

An example:

data(testdata)
spmatrix <- data2mat(testdata)
sp.assoc(spmatrix)

1. Variance of species relative frequency (Equation 1) sigmaTsq:

```{R echo=FALSE} knitr::include_graphics("clip_image001.png")

#### 2. Variance of the number of species (Equation 2) `STsq`

```{R echo=FALSE}
knitr::include_graphics("clip_image002.png")

3. Species' relative frequency (Equation 3) pi

```{R echo=FALSE} knitr::include_graphics("clip_image003.png")

#### 4. Variance ratio (Equation 4) `var.ratio`

```{R echo=FALSE}
knitr::include_graphics("clip_image004.png")

5. The W statistic (Equation 5)

The W statistic is used to test the significance of the variance ratio (Qianmei Zhang et al., 2006) with 95% confidence intervals of Chi 0.95,N2 \< W \< Chi 0.05, N2

```{R echo=FALSE} knitr::include_graphics("clip_image005.png")

In Equations 1-5, *N* is the number of plots, *S* is the total number of species, *n* is the number of plots occupied by a species, *T_j* is the number of species in each plot, and *t* is the average of the number of species in all plots.

### 3.2.2 Significance test of inters-pecific associations

The association between each pair of species could be calculated using the following code:

```{R}
data(testdata)
spmatrix <- data2mat(testdata)
result <- sp.pair(spmatrix)

In this test, a 2X2 contingent table is usually used. If there are species X and species Y occurring in the plots, a, b, c, d denote co-occurrence and non-occurrence between them, respectively (Kent M., 2012, page 114):

The Yates-corrected Chi-square is used when determining whether there is a significant association for small samples (Kent M., 2012, page 114).

The output is a list containing the following elements:

The equations are shown below:

1. Yates-corrected chi-square test for significance (Equation 6), see chisq and chisqass

Note: - chisqass is the Chi-square - chisq is the Chi-square with Yates' correction for small samples (Kent M., 2012, page 114)

```{R echo=FALSE} knitr::include_graphics("clip_image006.png")

#### 2. V ratio (Equation 7) `V`

```{R echo=FALSE}
knitr::include_graphics("clip_image007.png")

Indicating whether the species association is positive (>0) or negative (\<0)

3. Jaccard's index (Equation 8) Jaccard:

```{R echo=FALSE} knitr::include_graphics("clip_image008.png")

#### 4. Ochiai index (Equation 9) `Ochiai`:

```{R echo=FALSE, message=FALSE}
knitr::include_graphics("clip_image009.png")

5. Dice Index (Equation 10) Dice:

```{R echo=FALSE} knitr::include_graphics("clip_image010.png")

#### 6. *PCC*: Point Correlation Coefficient `PCC` (Equation 11)

```{R echo=FALSE}
knitr::include_graphics("clip_image011.png")

7. AC Association coefficient (Equation 12) AC

```{R echo=FALSE} knitr::include_graphics("clip_image012.png")

-   If `a*d < b*c` and `d >= a`, then (Equation 13):

```{R echo=FALSE}
knitr::include_graphics("clip_image013.png")

```{R echo=FALSE} knitr::include_graphics("clip_image014.png")

#### 8. Spearman's rank correlation coefficient `Spearman` (Equation 15)

```{R echo=FALSE}
knitr::include_graphics("clip_image015.png")

9. Pearson's correlation coefficient Pearson (Equation 16)

```{R echo=FALSE} knitr::include_graphics("clip_image016.png")

## 3.3 Calculating niche width

Niche width is commonly measured by the Levins or Shannon indices.

### 1. Levins niche width index (Equation 17)

```{R echo=FALSE}
knitr::include_graphics("clip_image017.png")

2. Shannon's niche width index (Equation 18)

```{R echo=FALSE} knitr::include_graphics("clip_image018.png")

where *B_i* is the niche width of the *i*th species, *j* denotes jth plot, and *r* denotes the number of plots.

The function `niche.width()` can calculate the niche width of each species. The arguments are: `niche.width(mat, method = c("shannon", "levins"))`, where `mat` is the sample-species matrix and `method` is eighther `"shannon"` or `"levins"`.

## 3.4 Niche overlap

`niche.overlap()` calculates the niche overlap coefficients between all species pairs, including: `"levins"` (Note this is different from niche width), `"schoener"`, `"petraitis"`, `"pianka"`, `"morisita"`, etc. The output is a distance matrix. These indices are defined as below:

### 1. Levin's niche overlap index (Equation 19)

```{R echo=FALSE}
knitr::include_graphics("clip_image019.png")

2. Schoener's niche overlap index (Equation 20)

```{R echo=FALSE} knitr::include_graphics("clip_image020.png")

### 3. Petraitis's niche overlap index (Equation 21)

```{R echo=FALSE}
knitr::include_graphics("clip_image021.png")

4. Pianka's niche overlap index (Equation 22)

```{R echo=FALSE} knitr::include_graphics("clip_image022.png")

### 5. Czechanowski's niche overlap index (Equation 23)

```{R echo=FALSE}
knitr::include_graphics("clip_image023.png")

6. the simplified Morisita index (Equation 24)

```{R echo=FALSE} knitr::include_graphics("clip_image024.png")

where *O_ik* is the niche overlap of species *i* and species*k*, *P_ij* and *P_kj* are the number of individuals of species *i* and species *k* in the *j*th sample, respectively, *r* is the total number of plots (sites), and *e* is the base of the natural logarithm.

## Example:

```r
data(datasample)
niche.overlap(mat, method = c("levins", "schoener", "petraitis", "pianka", "czech", "morisita"))

where mat is the community matrix and method is either: "levins", "schoener", "petraitis", "pianka", "czech" or "morisita".

niche.overlap.pair() computes the niche overlap between a specific pair. For example:

niche.overlap.pair(vectA, vectB, 
                   method = c("pianka", "schoener", 
                              "petraitis", "czech",  
                              "morisita", "levins"))

Where VectA and vectB are two vectors indicating the number of individuals of species A and species B in the plots, respectively, and method is either "pianka", "schoener", "petraitis", "czech", "morisita", "levins".

3.5 Bootstrap of the confidence intervals for niche width

The confidence interval of niche overlap between each pair of species could be estimated using bootstrap by niche.overlap.boot(), which is show below:

niche.overlap.boot.pair(vectorA, vectorB, method = c("levins",
     "schoener", "petraitis", "pianka", "czech", "morisita"),
      times = 1000, quant = c(0.025, 0.975))

where mat is the input species distribution matrix. method is the niche overlap index, times is the times of bootstrap to be performed, and quant is the quantile of the niche overlap index. By defaults,the 95% confidence interval (0.025 and 0.975) is used.

niche.overlap.boot() is a wrapper of niche.overlap.boot.pair(). In most cases, users do not need to call niche.overlap.boot.pair().

Confidence intervals of species niche overlap analysis

The columns are as follows.

Note: If a species is missing, NaN will appear during in the bootstrap results, which indicates that the denominators is 0.

3.6 Visualisation

3.6.1 Semi-matrix plots

The spaa package also provides some functions for visualising the results of species association.

```{R echo=FALSE} knitr::include_graphics("correlation.jpeg")

The arguments are:

```R
plotlowertri(input, valuename = "r",  
    pchlist = c(19, 17, 15, 1, 5, 2, 7), 
    interval = 6,  cex = 1, ncex = 1, 
    int =1.2, add.number = TRUE, 
    size = FALSE,  add.text = FALSE,  
    show.legend = TRUE, digits = 2)

A correlation coefficient matrix or a distance matrix could easily be visuallized using this function (Figure 1).

3.6.2 Species association network

Networks can be used to show the interconnectedness of species, i.e., in different colors, line widths, or line styles connecting the species to another other, or the similarity of species composition between sites. spaa's plotnetwork() can visualize the relationship (Figure 3). However, network plot should not be used to represent the relationship between a large number of species, for example, more than 10, it could be too crowded.

{R echo=FALSE} knitr::include_graphics("network.jpeg")

A network plot generated using plotnetwork()

3.7 Other functions

3.7.1 Obtain the most frequently occurred species in the community matrix

the sub.sp.matrix(spmatrix, freq = 0.5, common = NULL) can select the species based on relative frequency, for example, freq = 0.5 will only keep the species whose relative frequency greater than 0.5, or common = 10 will select the top 10 species with highest relative frequency.

3.7.2 Other useful functions

Note: These functions were available before the version 0.2.0, but have been removed from the development. I am planning to re-include them in spaa in the near future.

Acknowledgements

We thank Qiong Ding, Jihong Huang, and Dr. Zongshan Li for their help in the development of the package, Jinze Ma, Xueni Zhang, Jun Jiang, Meixiang Gao, Xiaofeng Fang, Yinyou, Wei Li, Alfredo H. Zúñiga Á., ÍsisArantes, Patricia Martínez, Wilson Martins da Silva, Jessica L. Sabo, Maud CHARLERY, Du Zhongyu, Diego Procopio, Zhan Xiaohao, Luis Fernando Gatica Mora, Zhao Wenxi, Xiao Zhengli, Jin Chao, Xing Bingwei, Pan Da, Vicente García-Navas, Jiang Huan, Clara Ruiz. Wang Lulu, Clara Ruiz González, Guo Jiaxing, Du Yuanbao, Zhong Yunjie, Wang Chenhe, Xu Heng, Zhang Wei, Hao Minhui, Suhridam Roy, Wang Yuan, Han Dayong, Yang Fan, Yao Xueqin, Wang Jing, Yang Haitao, Zheng Guiling, Yue Pengpeng, Mary Ann McLean, Bao Zhigui, Russell Bicknell. Joan Giménez Verdugo, Angela Andrea Camargo Sanabria, Xiaodan Zhao, Ramiro Logares Haurié, Simone Cappellari Rabeling, Caitlin Keating-Bitonti. Dr. Kurt Hornik and Prof. Brian Ripley, who provided some comments and suggestions on the compilation of the package, are gratefully acknowledged.

References



helixcn/spaa documentation built on Aug. 13, 2021, 3:21 a.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close