R在數(shù)量生態(tài)學(xué)分析中的應(yīng)用（一）

Numerical Ecology with R 這本書2018年出了第二版，第一版有中文版的書籍。自己的課題涉及微生物生態(tài)學(xué)，但時常對生態(tài)學(xué)中的一些專業(yè)名詞概念和分析方法不清楚不理解，嚴(yán)重阻礙了對微生物生態(tài)學(xué)問題的解決和分析。R語言是生物統(tǒng)計分析和可視化方面的神器，感覺這本書在R語言和生態(tài)學(xué)的學(xué)習(xí)中可以起到一舉兩得的作用。

1.數(shù)據(jù)的來源
本書的作者提供了兩個數(shù)據(jù)集以及對應(yīng)的代碼，對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)學(xué)研究需要對數(shù)據(jù)進行假設(shè)檢驗和模型化。因此對多維數(shù)據(jù)統(tǒng)計分析的早期是需要使用簡單的統(tǒng)計分析和可視化工具進行數(shù)據(jù)探索性分析的。以期獲得一下信息：

  • 整個數(shù)據(jù)的全貌概況。
  • 對數(shù)據(jù)的變量進行轉(zhuǎn)化和重新編輯
  • 確定進一步分析的方向

2.1數(shù)據(jù)導(dǎo)入和查看

# 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")
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容