1.1 加载包

library(readxl) #读取xlsx文件

1.2 读取案例(一)数据

re_data <- read_xlsx('C:\\Users\\dell\\Desktop\\北大R语言\\种丰度模型1.xlsx') #读取种丰度模型1
str(re_data)

## Classes 'tbl_df', 'tbl' and 'data.frame':    58 obs. of  9 variables:
##  $ Species : chr  "Lepetidae gen. et sp. indet." "Niveotecturapallida" "Punctrellanobilis" "Lirulariayamadana" ...
##  $ Category: chr  "Undetermined" "Cold-water" "Cosmopolitan" "Cold-water" ...
##  $ W4      : num  0 0 0 5 3 0 4 2 0 0 ...
##  $ W5      : num  0 0 0 0 2 0 2 2 0 0 ...
##  $ W6      : num  5 1 2 0 0 0 0 3 5 0 ...
##  $ W7      : num  7 0 0 0 0 0 0 2 0 1 ...
##  $ W8      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ W9      : num  2 1 1 2 0 1 0 7 0 0 ...
##  $ Total   : num  14 2 3 7 5 1 6 16 5 1 ...

  通过str函数可以发现数据集共有58条数据，9个变量，其中Species和Category为字符串型数据，w2-w9为采样范围，total等于w2-w9的和。

1.3 案例(一)数据预处理

warm_water <- re_data[which(re_data$Category == "Warm-water"),] #筛选Category列类型为Warm-water的数据
cos_cold <- re_data[which(re_data$Category == "Cosmopolitan/\r\nCold-water"),] #筛选Category列类型为Cosmopolitan/Cold-water的数据

  筛选出选Category列类型为Warm-water和Cosmopolitan/Cold-water的数据

total <- apply(re_data[,3:8], 2, sum) #w4-w9列求和

  计算出物种总数

warm_water1 <- apply(warm_water[,3:8], 2, sum) #w4-w9列求和
cos_cold1 <- apply(cos_cold[,3:8], 2, sum) #w4-w9列求和
plot_data <- data.frame()
plot_data[1:6,'Warm-water predator'] <- warm_water1 * 10 / total #计算Warm_water_predator物种丰度
plot_data[1:6,'Cosmopolitan/Cold-water predator'] <- cos_cold1 * 100 / total #计算Cosmopolitan/Cold-water物种丰
rownames(plot_data) <- names(total)
plot_data[4:6, 1] <- 0
head(plot_data)

##    Warm-water predator Cosmopolitan/Cold-water predator
## W4           0.7507862                        0.4716981
## W5           1.0633484                        0.7541478
## W6           0.1312155                        1.1049724
## W7           0.0000000                        1.6997167
## W8           0.0000000                        0.4842615
## W9           0.0000000                        0.4390779

  计算物种相对丰度(计算方法是某个物种的个体数目/ 全部物种的个体数目)

1.4 绘制种丰度条形图

barplot(as.matrix(t(plot_data)),xlab = 'Sampling horizon', ylab = 'Relative abundance(%)', ylim = c(0, 3), col = c('black', 'white'), legend = rownames(t(plot_data)), beside = TRUE)
text(10.5,0.1,'+',cex=1.2);text(13.5,0.1,'-',cex=1.2);text(16.5,0.1,'-',cex=1.2)

  黑色表示的是暖水型捕食者，白色的是世界型/冷水型捕食者，在w4和w5中暖水型捕食者的物种丰度大于世界型/冷水型捕食者(+代表存在，-代表缺失)。

2.1 加载包

library(plotrix) #绘制双坐标图

2.2 读取案例(二)数据

data1 <- read.csv('C:\\Users\\dell\\Desktop\\北大R语言\\种丰度模型2 表1.csv') #读取种丰度模型2 表1.csv数据
str(data1)

## 'data.frame':    5553 obs. of  5 variables:
##  $ SubOrder               : Factor w/ 8 levels "","Adephaga",..: 1 3 1 1 1 1 1 1 1 1 ...
##  $ family_name            : Factor w/ 149 levels "","Acanthocnemidae",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ occurrence.genus_name  : Factor w/ 2139 levels "","Abderina",..: 920 953 972 1000 1215 1420 1472 1472 1472 1488 ...
##  $ occurrence.species_name: Factor w/ 2737 levels "","abacus","abavun",..: 1625 1664 320 2193 2678 2706 796 1238 2678 242 ...
##  $ collection_no          : int  49 49 49 49 49 49 49 49 49 49 ...

data2 <- read.csv('C:\\Users\\dell\\Desktop\\北大R语言\\种丰度模型2 表2.csv') #读取种丰度模型2 表2.csv数据
str(data2)

## 'data.frame':    221 obs. of  8 variables:
##  $ collection_no  : int  2 3 5 7 8 9 10 12 13 14 ...
##  $ collection_name: Factor w/ 221 levels " Unda River",..: 11 2 3 4 29 5 31 6 7 186 ...
##  $ formation      : Factor w/ 157 levels ""," U. Ypresian, Lutetian",..: 96 57 105 1 1 4 1 89 1 61 ...
##  $ epoch          : Factor w/ 9 levels "Cretaceous","Eocene",..: 7 3 5 1 1 2 5 4 3 7 ...
##  $ subepoch       : Factor w/ 94 levels "Early-Middle Eocene",..: 60 19 53 7 12 16 25 50 17 85 ...
##  $ ma_max         : num  259 201.5 28.4 112 136.4 ...
##  $ ma_min         : num  252.3 196.5 23 99.6 130 ...
##  $ lithdescript   : Factor w/ 61 levels "argillite w/ clay lenses",..: 51 8 27 18 52 43 7 12 42 51 ...

  data1是关于甲虫化石信息的数据，共有5553条数据，其中前四个变量分别表示的是昆虫的亚目、科、属、种名，collection_no表示的是化石发现地编号。

  data2是关于化石地点信息的数据，共有221条数据，其中collection_no表示的是化石发现地编号，collection name表示的是发现地名称，formation表示的是岩石地层，epoch和subepoch是关于地层年代信息，ma max表示的是地层年代最大值，ma min表示的是地层年代最小值，lithdescript表示的是采集化石围岩的岩性。

2.3 案例(二)数据预处理

num_col <- data.frame(table(data1$collection_no)) #统计收集地点(collection_no)采集到的昆虫数
names(num_col) <- c('collection_no','occurrences') #变量重新命名

  统计data1中每一个收集地点(collection_no)所采集到的昆虫数。

year <- subset(data2, select = c(collection_no, ma_max, ma_min)) #筛选出data2中的列collection_no, ma_max, ma_min
year['mean_year'] <- (year$ma_max + year$ma_min) / 2 #求取地层年代的中间值
year['time'] <- ceiling(year$mean_year / 25) #以25个百万年间隔划分mean_year(地层年代的中间值)

  计算出data2中每个收集地点的地层年代平均值，并将地层年代平均值以25个百万年间隔划分为12个年代区间，time中的1表示0-25百万年，2表示25-50百万年……12表示的是275-300万年。

new_data <- merge(year, num_col, by = 'collection_no') #合并数据集year，num_col

  合并数据集year，num_col

plot_data2 <- data.frame(time = 1:12)
plot_data2['occurrences'] <- tapply(new_data$occurrences, new_data$time,  sum) #统计在每个年代区间采集到的昆虫
plot_data2['collection'] <- table(new_data$time) #统计在每个年代区间的化石地点数
head(plot_data2)

##   time occurrences collection
## 1    1         889         33
## 2    2        2193         41
## 3    3         268         15
## 4    4          48         10
## 5    5         198         16
## 6    6         761         26

  plot_data2中列occurrences为每个年代区间采集到的昆虫数，列collection为每个年代区间的化石地点数

2.4 绘制种丰度折线图

time <- c(12:1) ; time <- time * 25
twoord.plot(lx = c(12:1), ly = plot_data2$occurrences, rx = c(12:1), ry  = plot_data2$collection, lpch = 1, rpch = 1,xlab = 'time(Ma)', ylab = 'no.occurrences', rylab = 'no.collections', xticklab = time, xtickpos = c(0.65:11.65), lwd = 2, rcol = 'red' , lcol = 'black')
legend('topleft', leg = c('collections','occurrences'), col = c('red', 'black'), lty = 1, bty = 'n', pch = c(1, 1), lwd = 2)

  从图中可以看出在地层年代为25-50百万年昆虫化石数量和化石发现地最多，且昆虫化石数量和化石发现地的波动基本一致。甲虫科的总体取样随着时间的推移而变化，在晚侏罗世到早白垩世和中新生代发现的甲虫的采集和标本出现数量达到高峰。从白垩纪中期到晚白垩世是一个采样率特别低的时期，无论是在采集数量还是出现次数上。大陆脊椎动物化石记录中也记录了类似的模式，至于这是否归因于更大的全球过程，如海洋海侵和海退旋回，仍存在争议。