R在数量生态学分析中的应用（一）

Numerical Ecology with R 这本书2018年出了第二版，第一版有中文版的书籍。自己的课题涉及微生物生态学，但时常对生态学中的一些专业名词概念和分析方法不清楚不理解，严重阻碍了对微生物生态学问题的解决和分析。R语言是生物统计分析和可视化方面的神器，感觉这本书在R语言和生态学的学习中可以起到一举两得的作用。

1.数据的来源
本书的作者提供了两个数据集以及对应的代码，对Doubs River流域30个采样点的鱼类以及其他环境参数进行收集，形成了Doubs.RData数据集，对一个区域环境70个采样点中的螨虫等微小生物以及其他环境参数进行收集，形成了mite.RData数据集。

2.数据探索
现如今，大多数生态学研究需要对数据进行假设检验和模型化。因此对多维数据统计分析的早期是需要使用简单的统计分析和可视化工具进行数据探索性分析的。以期获得一下信息：

整个数据的全貌概况。
对数据的变量进行转化和重新编辑
确定进一步分析的方向

2.1数据导入和查看

# Load required packages
> library(vegan)
> library(RgoogleMaps)
> library(googleVis)
> library(labdsv)
# Source additional functions that will be used later in this
# Chapter. Our scripts assume that files to be read are in
# the working directory.
> source("panelutils.R")
# Load the data. File Doubs.Rdata is assumed to be
# in the working directory
> load("Doubs.RData")
# The file Doubs.RData contains the following objects:
# spe: species (community) data frame (fish abundances)
# env: environmental data frame
# spa: spatial data frame – cartesian coordinates
# fishtraits: functional traits of fish species
# latlong: spatial data frame – latitude and longitude

## Exploration of a data frame using basic R functions
> spe # Display the whole data frame in the
> spe
   Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme Gyce
1     0    3    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
2     0    5    4    3    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
3     0    5    5    5    0    0    0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0
4     0    4    5    5    0    0    0    0    0    1    0    0    1    2    2    0    0    0    0    1    0    0    0
5     0    2    3    2    0    0    0    0    5    2    0    0    2    4    4    0    0    2    0    3    0    0    0
6     0    3    4    5    0    0    0    0    1    2    0    0    1    1    1    0    0    0    0    2    0    0    0
7     0    5    4    5    0    0    0    0    1    1    0    0    0    0    0    0    0    0    0    0    0    0    0
8     0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
9     0    0    1    3    0    0    0    0    0    5    0    0    0    0    0    0    0    0    0    1    0    0    0
10    0    1    4    4    0    0    0    0    2    2    0    0    1    0    0    0    0    0    0    0    0    0    0
......
# console
# Not recommended for large datasets!
> spe[1:5, 1:10] # Display only 5 lines and 10 columns
  Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce
1    0    3    0    0    0    0    0    0    0    0
2    0    5    4    3    0    0    0    0    0    0
3    0    5    5    5    0    0    0    0    0    0
4    0    4    5    5    0    0    0    0    0    1
5    0    2    3    2    0    0    0    0    5    2
> head(spe) # Display only the first 6 lines
 Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme Gyce Ruru
1    0    3    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
2    0    5    4    3    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
3    0    5    5    5    0    0    0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0
4    0    4    5    5    0    0    0    0    0    1    0    0    1    2    2    0    0    0    0    1    0    0    0    0
5    0    2    3    2    0    0    0    0    5    2    0    0    2    4    4    0    0    2    0    3    0    0    0    5
6    0    3    4    5    0    0    0    0    1    2    0    0    1    1    1    0    0    0    0    2    0    0    0    1
  Blbj Alal Anan
1    0    0    0
2    0    0    0
3    0    0    0
4    0    0    0
5    0    0    0
6    0    0    0

> tail(spe) # Display only the last 6 rows
   Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme Gyce
25    0    0    0    0    0    0    0    0    1    1    0    0    2    1    0    0    0    1    0    0    0    0    1
26    0    0    0    1    0    0    1    0    1    2    2    1    3    2    1    2    2    1    1    3    2    1    4
27    0    0    0    1    0    0    1    1    2    3    4    1    4    4    1    3    3    1    2    5    3    2    5
28    0    0    0    1    0    0    1    1    2    4    3    1    4    3    2    4    4    2    4    4    3    3    5
29    0    1    1    1    1    1    2    2    3    4    5    3    5    5    4    5    5    2    3    3    4    4    5
30    0    0    0    0    0    0    1    2    3    3    3    5    5    4    5    5    3    5    5    5    5    5    5
   Ruru Blbj Alal Anan
25    1    0    3    0
26    4    2    5    2
27    5    4    5    3
28    5    5    5    4
29    5    4    5    4
30    5    5    5    5

> nrow(spe) # Number of rows (sites)
[1] 30

> ncol(spe) # Number of columns (species)
[1] 27

> dim(spe) # Dimensions of the data frame (rows, columns)
[1] 30 27

> colnames(spe) 
 [1] "Cogo" "Satr" "Phph" "Babl" "Thth" "Teso" "Chna" "Pato" "Lele" "Sqce" "Baba" "Albi" "Gogo" "Eslu" "Pefl" "Rham"
[17] "Legi" "Scer" "Cyca" "Titi" "Abbr" "Icme" "Gyce" "Ruru" "Blbj" "Alal" "Anan"

> rownames(spe) # Row labels (objects = sites)
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15" "16" "17" "18" "19" "20" "21" "22" "23"
[24] "24" "25" "26" "27" "28" "29" "30"

> summary(spe) # Descriptive statistics for columns
      Cogo           Satr           Phph            Babl            Thth           Teso             Chna    
 Min.   :0.00   Min.   :0.00   Min.   :0.000   Min.   :0.000   Min.   :0.00   Min.   :0.0000   Min.   :0.0  
 1st Qu.:0.00   1st Qu.:0.00   1st Qu.:0.000   1st Qu.:1.000   1st Qu.:0.00   1st Qu.:0.0000   1st Qu.:0.0  
 Median :0.00   Median :1.00   Median :3.000   Median :2.000   Median :0.00   Median :0.0000   Median :0.0  
 Mean   :0.50   Mean   :1.90   Mean   :2.267   Mean   :2.433   Mean   :0.50   Mean   :0.6333   Mean   :0.6  
 3rd Qu.:0.75   3rd Qu.:3.75   3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.:0.75   3rd Qu.:0.7500   3rd Qu.:1.0  
 Max.   :3.00   Max.   :5.00   Max.   :5.000   Max.   :5.000   Max.   :4.00   Max.   :5.0000   Max.   :3.0  
      Pato             Lele            Sqce            Baba            Albi          Gogo            Eslu      
 Min.   :0.0000   Min.   :0.000   Min.   :0.000   Min.   :0.000   Min.   :0.0   Min.   :0.000   Min.   :0.000  
 1st Qu.:0.0000   1st Qu.:0.000   1st Qu.:1.000   1st Qu.:0.000   1st Qu.:0.0   1st Qu.:0.000   1st Qu.:0.000  
 Median :0.0000   Median :1.000   Median :2.000   Median :0.000   Median :0.0   Median :1.000   Median :1.000  
 Mean   :0.8667   Mean   :1.433   Mean   :1.867   Mean   :1.433   Mean   :0.9   Mean   :1.833   Mean   :1.333  
 3rd Qu.:2.0000   3rd Qu.:2.000   3rd Qu.:3.000   3rd Qu.:3.000   3rd Qu.:1.0   3rd Qu.:3.750   3rd Qu.:2.000  
 Max.   :4.0000   Max.   :5.000   Max.   :5.000   Max.   :5.000   Max.   :5.0   Max.   :5.000   Max.   :5.000  
      Pefl          Rham          Legi             Scer          Cyca             Titi          Abbr       
 Min.   :0.0   Min.   :0.0   Min.   :0.0000   Min.   :0.0   Min.   :0.0000   Min.   :0.0   Min.   :0.0000  
 1st Qu.:0.0   1st Qu.:0.0   1st Qu.:0.0000   1st Qu.:0.0   1st Qu.:0.0000   1st Qu.:0.0   1st Qu.:0.0000  
 Median :0.5   Median :0.0   Median :0.0000   Median :0.0   Median :0.0000   Median :1.0   Median :0.0000  
 Mean   :1.2   Mean   :1.1   Mean   :0.9667   Mean   :0.7   Mean   :0.8333   Mean   :1.5   Mean   :0.8667  
 3rd Qu.:2.0   3rd Qu.:2.0   3rd Qu.:1.7500   3rd Qu.:1.0   3rd Qu.:1.0000   3rd Qu.:3.0   3rd Qu.:1.0000  
 Max.   :5.0   Max.   :5.0   Max.   :5.0000   Max.   :5.0   Max.   :5.0000   Max.   :5.0   Max.   :5.0000  
      Icme          Gyce            Ruru          Blbj            Alal          Anan     
 Min.   :0.0   Min.   :0.000   Min.   :0.0   Min.   :0.000   Min.   :0.0   Min.   :0.00  
 1st Qu.:0.0   1st Qu.:0.000   1st Qu.:0.0   1st Qu.:0.000   1st Qu.:0.0   1st Qu.:0.00  
 Median :0.0   Median :0.000   Median :1.0   Median :0.000   Median :0.0   Median :0.00  
 Mean   :0.6   Mean   :1.267   Mean   :2.1   Mean   :1.033   Mean   :1.9   Mean   :0.90  
 3rd Qu.:0.0   3rd Qu.:2.000   3rd Qu.:5.0   3rd Qu.:1.750   3rd Qu.:5.0   3rd Qu.:1.75  
 Max.   :5.0   Max.   :5.000   Max.   :5.0   Max.   :5.000   Max.   :5.0   Max.   :5.00  

> range(spe) ## Overall distribution of abundances (dominance codes)
# Minimum and maximum of abundance values in the whole data set
[1] 0 5

> apply(spe, 2, range)
     Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme
[1,]    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
[2,]    3    5    5    5    4    5    3    4    5    5    5    5    5    5    5    5    5    5    5    5    5    5
     Gyce Ruru Blbj Alal Anan
[1,]    0    0    0    0    0
[2,]    5    5    5    5    5

> apply(spe, 2, range)
     Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme
[1,]    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
[2,]    3    5    5    5    4    5    3    4    5    5    5    5    5    5    5    5    5    5    5    5    5    5
     Gyce Ruru Blbj Alal Anan
[1,]    0    0    0    0    0
[2,]    5    5    5    5    5

> (ab <- table(unlist(spe)))  # Count the cases for each abundance class

  0   1   2   3   4   5 
435 108  87  62  54  64 

> barplot(ab, las = 1, xlab = "Abundance class", ylab = "Frequency",col = gray(5 : 0 / 5)) # Barplot of the distribution, all species confounded

> sum(spe == 0) # Number of absences
[1] 435

> sum(spe == 0) / (nrow(spe) * ncol(spe)) # Proportion of zeros in the community data set
[1] 0.537037

地理位置参数的可视化

> plot(spa,
+      asp = 1,
+      type = "n",
+      main = "Site Locations",
+      xlab = "x coordinate (km)",
+      ylab = "y coordinate (km)"
+ )  # Geographic coordinates x and y from the spa data frame
> lines(spa, col = "light blue") # Add a blue line connecting the sites along the Doubs River
> text(spa, row.names(spa), cex = 0.8, col = "red") # Add the site labels
> text(68, 20, "Upstream", cex = 1.2, col = "red") # Add text blocks
> text(15, 35, "Downstream", cex = 1.2, col = "red")

使用Google Maps® map展示结果

> nom <- latlong$Site
> latlong2 <- paste(latlong$LatitudeN, latlong$LongitudeE, sep = ":")
> df <- data.frame(latlong2, nom, stringsAsFactors = FALSE)
> mymap1 <- gvisMap(df,
+                   locationvar = "latlong2",
+                   tipvar = "nom",
+                   options = list(showTip = TRUE)
+ )
> plot(mymap1)

> par(mfrow = c(2,2)) # Divide the plot window into 4 frames, 2 per row
 # Plot four species
> plot(spa,
+      asp = 1,
+      cex.axis = 0.8,
+      col = "brown",
+      cex = spe$Satr,
+      main = "Brown trout",
+      xlab = "x coordinate (km)",
+      ylab = "y coordinate (km)"
+ ) 
> lines(spa, col = "light blue")
> plot(spa,
+      asp = 1,
+      cex.axis = 0.8,
+      col = "brown",
+      cex = spe$Thth,
+      main = "Grayling",
+      xlab = "x coordinate (km)",
+      ylab = "y coordinate (km)"
+ )
> lines(spa, col = "light blue")
> plot(spa,
+      asp = 1,
+      cex.axis = 0.8,
+      col = "brown",
+      cex = spe$Baba,
+      main = "Barbel",
+      xlab = "x coordinate (km)",
+      ylab = "y coordinate (km)"
+ )
> lines(spa, col = "light blue")
> plot(spa,
+      asp = 1,
+      cex.axis = 0.8,
+      col = "brown",
+      cex = spe$Abbr,
+      main = "Common bream",
+      xlab = "x coordinate (km)",
+      ylab = "y coordinate (km)"
+ )
> lines(spa, col = "light blue")

# Compute the number of sites where each species is present
# To sum by columns, the second argument of apply(), MARGIN,
# is set to 2
> spe.pres <- apply(spe > 0, 2, sum)
> spe.pres 
Cogo Satr Phph Babl Thth Teso Chna Pato Lele Sqce Baba Albi Gogo Eslu Pefl Rham Legi Scer Cyca Titi Abbr Icme Gyce 
   8   17   20   24    8    8   12   11   18   25   14   12   20   18   15   11   13   11   12   17    9    7   12 
Ruru Blbj Alal Anan 
  18   10   14   11 
> sort(spe.pres) #排序
Icme Cogo Thth Teso Abbr Blbj Pato Rham Scer Anan Chna Albi Cyca Gyce Legi Baba Alal Pefl Satr Titi Lele Eslu Ruru 
   7    8    8    8    9   10   11   11   11   11   12   12   12   12   13   14   14   15   17   17   18   18   18 
Phph Gogo Babl Sqce 
  20   20   24   25 

> spe.relf <- 100 * spe.pres/nrow(spe) #计算百分比
> spe.relf 
    Cogo     Satr     Phph     Babl     Thth     Teso     Chna     Pato     Lele     Sqce     Baba     Albi     Gogo 
26.66667 56.66667 66.66667 80.00000 26.66667 26.66667 40.00000 36.66667 60.00000 83.33333 46.66667 40.00000 66.66667 
    Eslu     Pefl     Rham     Legi     Scer     Cyca     Titi     Abbr     Icme     Gyce     Ruru     Blbj     Alal 
60.00000 50.00000 36.66667 43.33333 36.66667 40.00000 56.66667 30.00000 23.33333 40.00000 60.00000 33.33333 46.66667 
    Anan 
36.66667 

> par(mfrow = c(1,2))
> hist(spe.pres,
+      main = "Species Occurrences",
+      right = FALSE,
+      las = 1,
+      xlab = "Number of occurrences",
+      ylab = "Number of species",
+      breaks = seq(0, 30, by = 5),
+      col = "bisque"
+ )
> hist(spe.relf,
+      main = "Species Relative Frequencies",
+      right = FALSE,
+      las = 1,
+      xlab = "Frequency of occurrences (%)",
+      ylab = "Number of species",
+      breaks = seq(0, 100, by = 10),
+      col = "bisque"
+ )

> ## Compare sites: species richness
> # Compute the number of species at each site
> # To sum by rows, the second argument of apply(), MARGIN, is
> # set to 1
> sit.pres <- apply(spe > 0, 1, sum)
> sit.pres
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 
 1  3  4  8 11 10  5  0  5  6  6  6  6 10 11 17 22 23 23 22 23 22  3  8  8 21 22 22 26 21 
> # Sort the results in increasing order
> sort(sit.pres)
 8  1  2 23  3  7  9 10 11 12 13  4 24 25  6 14  5 15 16 26 30 17 20 22 27 28 18 19 21 29 
 0  1  3  3  4  5  5  6  6  6  6  8  8  8 10 10 11 11 17 21 21 22 22 22 22 22 23 23 23 26 
> par(mfrow = c(1, 2))
> # Plot species richness vs. position of the sites along the river
> plot(sit.pres,type = "s",
+      las = 1,
+      col = "gray",
+      main = "Species Richness vs. \n Upstream-Downstream Gradient",
+      xlab = "Site numbers",
+      ylab = "Species richness"
+ )
> text(sit.pres, row.names(spe), cex = .8, col = "red")
## Compare sites: species richness
# Compute the number of species at each site
# To sum by rows, the second argument of apply(), MARGIN, is
# set to 1
sit.pres <- apply(spe > 0, 1, sum)
# Sort the results in increasing order
sort(sit.pres)
par(mfrow = c(1, 2))
# Plot species richness vs. position of the sites along the river
plot(sit.pres,type = "s",
las = 1,
col = "gray",
main = "Species Richness vs. \n Upstream-Downstream Gradient",
xlab = "Site numbers",
ylab = "Species richness"
)
text(sit.pres, row.names(spe), cex = .8, col = "red")
# Use geographic coordinates to plot a bubble map
plot(spa,
asp = 1,
main = "Map of Species Richness",
pch = 21,
col = "white",
bg = "brown",
cex = 5 * sit.pres / max(sit.pres),
xlab = "x coordinate (km)",
ylab = "y coordinate (km)"
)
lines(spa, col = "light blue")

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 216,193评论 6赞 498
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 92,306评论 3赞 392
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 162,130评论 0赞 353
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 58,110评论 1赞 292
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 67,118评论 6赞 388
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 51,085评论 1赞 295
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 40,007评论 3赞 417
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 38,844评论 0赞 273
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 45,283评论 1赞 310
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 37,508评论 2赞 332
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 39,667评论 1赞 348
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 35,395评论 5赞 343
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 40,985评论 3赞 325
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 31,630评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,797评论 1赞 268
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 47,653评论 2赞 368
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 44,553评论 2赞 352

R在数量生态学分析中的应用（一）

推荐阅读更多精彩内容