Numerical Ecology with R 這本書2018年出了第二版黑毅,第一版有中文版的書籍。自己的課題涉及微生物生態(tài)學灶体,但時常對生態(tài)學中的一些專業(yè)名詞概念和分析方法不清楚不理解阅签,嚴重阻礙了對微生物生態(tài)學問題的解決和分析。R語言是生物統(tǒng)計分析和可視化方面的神器蝎抽,感覺這本書在R語言和生態(tài)學的學習中可以起到一舉兩得的作用政钟。
1.數(shù)據(jù)的來源
本書的作者提供了兩個數(shù)據(jù)集以及對應的代碼,對Doubs River流域30個采樣點的魚類以及其他環(huán)境參數(shù)進行收集樟结,形成了Doubs.RData數(shù)據(jù)集养交,對一個區(qū)域環(huán)境70個采樣點中的螨蟲等微小生物以及其他環(huán)境參數(shù)進行收集,形成了mite.RData數(shù)據(jù)集狭吼。
2.數(shù)據(jù)探索
現(xiàn)如今层坠,大多數(shù)生態(tài)學研究需要對數(shù)據(jù)進行假設檢驗和模型化。因此對多維數(shù)據(jù)統(tǒng)計分析的早期是需要使用簡單的統(tǒng)計分析和可視化工具進行數(shù)據(jù)探索性分析的刁笙。以期獲得一下信息:
- 整個數(shù)據(jù)的全貌概況破花。
- 對數(shù)據(jù)的變量進行轉(zhuǎn)化和重新編輯
- 確定進一步分析的方向
2.1數(shù)據(jù)導入和查看
# 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
地理位置參數(shù)的可視化
> 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展示結(jié)果
> 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")