R在数量生态学分析中的应用（一）
2019-03-22 本文已影响24人
你猜我菜不菜
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")
