异速生长分析
2020-10-16
2020-10-16
2-1 异速生长分析
1 基本概念
生物在个体发育的不同阶段，其不同器官的生长速率其实是存在差异的。以人类为例，婴儿时期脑部的发育速度最快，这使得人类婴儿的头部相对身体而言显得很大；而随着年龄的增长，四肢的发育速度加快，并在青春期之后使得全身的比例变得协调。
在古生物学中，对于异速生长的研究一方面有助于理解某一类生物的生长发育历史，另一方面也有助于从定量化的角度更准确地认识化石物种在形态上的差异。
2 案例
2.1 案例介绍
性双型现象在古生物中强调得并不多，但是由于胡氏贵州龙的标本数量十分巨大，性双型现象较为明显。性双型现象的本质，其实仍然是异速生长，相当于器官的发育随年龄和性别产生了不等速的发育。但是胡氏贵州龙的性双型从何时开始？是从一出生就有所差异，还是在生长发育过程中出现的？具体在那个阶段产生性双型现象？本研究对此展开了研究。
library(MASS) #加载MASS包,用于fisher判别
library(readxl) #加载包readxl包,用于读取excel数据
library(ggplot2) #加载ggplot2包,用于绘制直方图
library(car) #加载car包,用于绘制置信椭圆2.2 读取数据
data <- read_excel('C:\\Users\\dell\\Desktop\\北大R语言\\2-1 异速生长分析\\Allometry_keichousaurus.xlsx')
str(data)## Classes 'tbl_df', 'tbl' and 'data.frame': 107 obs. of 17 variables:
## $ Number\item: chr "NMNS000933-F0034394" "NMNS-chw-01" "NMNS-chw-02" "NMNS-chw-03" ...
## $ nc.l : chr "28" "11.5" "18.2" "16" ...
## $ ta.l : chr "40" "16" "26" "15" ...
## $ s-v.l : chr "78" "33.5" "48.2" "40" ...
## $ sd.l : chr "7" "2.5" "4" "3" ...
## $ hm.l : chr "6.2" "2" "3.8" "2.8" ...
## $ fm.l : num 6.8 3 4.8 3.5 3 19 12.8 13 21 12 ...
## $ rd.l : chr "3.1" "1" "2" "1.2" ...
## $ ul.l : chr "3.1" "0.9" "2" "1.2" ...
## $ tb.l : chr "2.2999999999999998" "1" "1.8" "1" ...
## $ fb.l : chr "3" "1" "1.6" "1.2" ...
## $ hm.d : chr "2" "0.9" "1.2" "1" ...
## $ hm.m : chr "1" "0.8" "0.8" "0.5" ...
## $ fm.d : chr "1.8" "0.6" "1" "0.9" ...
## $ fm.m : chr "1" "0.5" "0.6" "0.5" ...
## $ Gender 1 : chr "M" "F" "F" "M" ...
## $ Gender 2 : chr "M" "F" "F" "M" ...本案例共有107条数据，17个变量。第2列到第15列变量类型应该为数值型，然而除了变量fm.l，其余的变量都是字符串型，而且还存在缺失值。变量Gender.1的类型应该是因子型。
2.3 数据清洗
for (i in 2:15){
data <- data[grepl('–', unlist(data[, i])) == FALSE, ] #找出第i列数据中含有'-'行,并剔除该行数据
data[, i] <- as.numeric(unlist(data[, i])) #将第i列数据由字符串型改为数值型
}
names(data)[c(1,16:17)] <- c('Number.item', 'Gender.1', 'Gender.2') #变量重命名
data$Gender.1 <- as.factor(data$Gender.1) #将变量Gender.1类型改为因子型
str(data)## Classes 'tbl_df', 'tbl' and 'data.frame': 76 obs. of 17 variables:
## $ Number.item: chr "NMNS000933-F0034394" "NMNS-chw-01" "NMNS-chw-02" "NMNS-chw-03" ...
## $ nc.l : num 28 11.5 18.2 16 76.5 60 57 70 54 81 ...
## $ ta.l : num 40 16 26 15 110 96 82 120 88 127 ...
## $ s-v.l : num 78 33.5 48.2 40 187 ...
## $ sd.l : num 7 2.5 4 3 18.2 13 12.5 15.5 11 18.5 ...
## $ hm.l : num 6.2 2 3.8 2.8 21 13 13 25.5 11.5 22 ...
## $ fm.l : num 6.8 3 4.8 3.5 19 12.8 13 21 12 21 ...
## $ rd.l : num 3.1 1 2 1.2 10 6.5 6.5 12.5 6 10 ...
## $ ul.l : num 3.1 0.9 2 1.2 9 6.2 6 12 6 9.5 ...
## $ tb.l : num 2.3 1 1.8 1 7.5 5 5 8 4.5 7 ...
## $ fb.l : num 3 1 1.6 1.2 8 6 5.5 9 4.5 8 ...
## $ hm.d : num 2 0.9 1.2 1 5 3 3 3 3 5.2 ...
## $ hm.m : num 1 0.8 0.8 0.5 5 2.5 2.5 2 2.2 3.7 ...
## $ fm.d : num 1.8 0.6 1 0.9 4 2.4 2 6 2.2 3 ...
## $ fm.m : num 1 0.5 0.6 0.5 3 1.8 1.8 3 1.3 2.7 ...
## $ Gender.1 : Factor w/ 2 levels "F","M": 2 1 1 2 2 1 1 2 1 1 ...
## $ Gender.2 : chr "M" "F" "F" "M" ...通过for循环，第2列到第15列的变量类型变为数字型，缺失值也被剔除。清洗后的数据共有76条数据，17个变量。
2.4 线性判别分析（LDA）
线性判别分析是一种经典的线性学习方法，在二分类问题上最早由Fisher在1936年提出，亦称Fisher线性判别。
线性判别的思想非常朴素：给定训练样例集，设法将样例投影到一条直线上，使得同类样例的投影点尽可能接近，异样样例的投影点尽可能远离；在对新样本进行分类时，将其投影到同样的直线上，再根据投影点的位置来确定新样本的类别。
ld <- lda(Gender.1~.-Number.item-Gender.1-Gender.2, data = data) #建立fisher模型
pre <- predict(ld, data) #通过fisher模型对数据进行预测建立线性判别分析（LDA）模型，并对原始数据进行预测。
hisdata <- data.frame(sex = ifelse(pre$posterior[,2] > 0.5, 'M', 'F'), pre$x) #建立数据框,pre$posterior[,2]>0.5的性别为M否则性别为F
ggplot(hisdata, aes(x = LD1, fill = sex)) +
labs(title = '', x='Linear discriminant analysis 1', y='Frequency') +
geom_histogram(position = 'identity', colour = 'black')## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
从线性判别分析的结果上看，已有的判别方式可以有效地对胡氏贵州龙的性别进行区分，但是由于亚成年个体的性别比例差异太大，鉴定结果可能不一定准确。
2.5 主成分分析
在用统计分析方法研究多变量的课题时，变量个数太多就会增加课题的复杂性。人们自然希望变量个数较少而得到的信息较多。在很多情形，变量之间是有一定的相关关系的，当两个变量之间有一定相关关系时，可以解释为这两个变量反映此课题的信息有一定的重叠。主成分分析是对于原先提出的所有变量，将重复的变量（关系紧密的变量）删去多余，建立尽可能少的新变量，使得这些新变量是两两不相关的，而且这些新变量在反映课题的信息方面尽可能保持原有的信息。
data.pr <- princomp(data[, 2:15], cor = T) #对2-15列进行主成分分析
summary(data.pr, loadings = T)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4
## Standard deviation 3.5321612 0.74441525 0.54210883 0.46117035
## Proportion of Variance 0.8911545 0.03958243 0.02099157 0.01519129
## Cumulative Proportion 0.8911545 0.93073693 0.95172850 0.96691979
## Comp.5 Comp.6 Comp.7 Comp.8
## Standard deviation 0.41215751 0.295276426 0.249525797 0.204059083
## Proportion of Variance 0.01213384 0.006227726 0.004447366 0.002974294
## Cumulative Proportion 0.97905364 0.985281364 0.989728730 0.992703023
## Comp.9 Comp.10 Comp.11 Comp.12
## Standard deviation 0.192008113 0.14501854 0.129961481 0.1065238524
## Proportion of Variance 0.002633365 0.00150217 0.001206428 0.0008105237
## Cumulative Proportion 0.995336389 0.99683856 0.998044986 0.9988555098
## Comp.13 Comp.14
## Standard deviation 0.094120198 0.0846418997
## Proportion of Variance 0.000632758 0.0005117322
## Cumulative Proportion 0.999488268 1.0000000000
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9
## nc.l -0.268 -0.316 0.129 -0.400 0.210 -0.132 0.424
## ta.l -0.251 -0.297 -0.703 0.577
## s-v.l -0.271 -0.310 0.105 -0.304 0.222 -0.113 0.238
## sd.l -0.272 -0.268 0.101 -0.155 -0.239 0.231 0.113 -0.723
## hm.l -0.277 0.254 -0.194
## fm.l -0.280 -0.108 -0.123 -0.300
## rd.l -0.272 0.311 0.196 -0.500 -0.200
## ul.l -0.273 0.297 0.131 -0.312 0.231
## tb.l -0.275 0.113 -0.261 0.219 -0.443 0.543
## fb.l -0.278 0.110 -0.178 0.107 -0.286 0.237
## hm.d -0.256 0.414 0.332 0.239 0.297 -0.252 0.474 0.411 0.168
## hm.m -0.256 -0.351 0.112 0.557 0.372 -0.278 -0.496 -0.123
## fm.d -0.249 0.163 -0.797 -0.492
## fm.m -0.261 -0.212 -0.425 0.327 0.215 0.657 0.307 0.131
## Comp.10 Comp.11 Comp.12 Comp.13 Comp.14
## nc.l 0.200 0.266 0.539
## ta.l
## s-v.l 0.254 -0.119 -0.210 -0.691
## sd.l -0.368 0.170
## hm.l -0.179 -0.102 0.873
## fm.l 0.697 -0.463 -0.266
## rd.l 0.459 0.504 -0.150
## ul.l -0.259 -0.540 -0.416 0.211 -0.286
## tb.l 0.472 -0.125 -0.131 0.168
## fb.l -0.607 0.452 -0.346 -0.180
## hm.d 0.110
## hm.m
## fm.d
## fm.mProportion of Variance为方差贡献率;Cumulative Proportion为累计方差贡献率。可以看出第1主成分的方差贡献率是89.12%，第2主成分的方差贡献率是3.96%。第1主成分和第2主成分的累计方差贡献率为93.07%。Loadings中的矩阵分别对应主成分与各样本之间的系数关系。
pca_data <- predict(data.pr)
head(pca_data) ## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
## [1,] 4.6745009 0.4678068 -0.518975133 0.0479015 0.07410052 -0.21053246
## [2,] 7.4746625 0.6181172 -0.006923345 0.4037332 0.48414492 0.02608269
## [3,] 6.4245128 0.6580121 -0.111669970 0.1327729 0.25818401 -0.02390274
## [4,] 7.2086854 0.7368848 -0.191516136 0.2187329 0.17812123 -0.07371821
## [5,] -4.2159204 -1.1652225 -0.975272299 1.1968901 0.41194089 -0.34088670
## [6,] 0.2022648 -0.7155536 -0.204886179 -0.2392423 -0.01261650 -0.03078081
## Comp.7 Comp.8 Comp.9 Comp.10 Comp.11
## [1,] 0.26086189 0.16393585 -0.099527885 -0.023550353 -0.02149840
## [2,] -0.14982090 0.06788429 -0.003305722 0.043480192 0.08295511
## [3,] -0.09980026 0.06420871 -0.031795288 0.109972401 0.04901181
## [4,] 0.09247862 0.02541350 0.028266061 0.001031493 0.11155211
## [5,] -0.45793108 0.18436674 -0.151976842 0.179078400 -0.04715084
## [6,] 0.01966840 0.13113553 0.007586874 -0.033939810 -0.08492114
## Comp.12 Comp.13 Comp.14
## [1,] 0.009495784 -0.111905031 -0.08580351
## [2,] -0.086422840 0.008734989 0.05859117
## [3,] -0.118565416 0.059079384 0.01174820
## [4,] -0.053178851 0.032310355 0.06320938
## [5,] -0.007326930 0.073546615 0.03145628
## [6,] 0.174763326 -0.066140591 -0.05782936各样本的主成分的值
dataEllipse(x = pca_data[which(data$Gender.1 == 'F'), 1], y = pca_data[which(data$Gender.1 == 'F'), 2], levels = 0.95,xlab = 'Factor 1', ylab = 'Factor 2', col = 'blue', pch = 1, xlim = c(-10, 10), ylim = c(-2, 2), center.pch = FALSE)
par(new=TRUE)
dataEllipse(x = pca_data[which(data$Gender.1 != 'F'), 1], y = pca_data[which(data$Gender.1 != 'F'), 2], levels = 0.70,xlab = '', ylab = '', col = 'red', pch = 4, xlim = c(-10, 10), ylim = c(-2, 2), center.pch = FALSE)
legend('topleft', legend = c('female', 'male'), col = c('blue', 'red'), pch = c(1, 4), bty = "n") #添加图例
从EFA的结果上看，胡氏贵州龙在体长（Factor 1）和性别（Factor 2）两个维度上的形态差异最大，而体形与年龄呈线性的正相关关系。所以，年龄与性别对胡氏贵州龙的形态影响最大。
2.6 协方差分析（ANCOVA）
协方差分析亦称“共变量（数）分析”。方差分析的引申和扩大。基本原理是将线性回归与方差分析结合起来，调整各组平均数和 F 检验的实验误差项，检验两个或多个调整平均数有无显著差异，以便控制在实验中影响实验效应（因变量）而无法人为控制的协变量（与因变量有密切回归关系的变量）在方差分析中的影响。例如，在研究某种教学方法（实验变量）对学业成绩（实验效应）的影响时，被试的原有知识基础同时影响学业成绩，但往往在实验中难以选取具备相同知识基础的被试参加实验，可用协方差分析从学业成绩的总变异中将归因于被试知识基础差异的部分划分出去，便于确切地分析教学方法对学业成绩的影响，其中被试的知识基础就是协变量。
names(data)[4] <- 's.v.l' #变量's-v.l'名改为's.v.l'
data_f <- data[which(data$Gender.1 == 'F'), ] #筛选出性别为F的数据
data_m <- data[which(data$Gender.1 == 'M'), ] #筛选出性别为M的数据
par(mfrow = c(2, 2), mar = c(4, 4, 1, 1))
#图A
plot(x = data_f$s.v.l, y = data_f$hm.l, xlim = c(0, 250), ylim = c(0, 30), pch = 1,
col = 'blue', xlab = 'Snout-vent length (mm)' , ylab = 'Humerus length (mm)')
abline(lm(data_f$hm.l~data_f$s.v.l), col = 'blue', lwd = 2, h = , v = ) #添加趋势线
par(new = TRUE)
plot(x = data_m$s.v.l, y = data_m$hm.l, xlim = c(0, 250), ylim = c(0, 30), pch = 4,
col = 'red', xlab = '', ylab = '', xaxt = 'n', yaxt = 'n')
abline(lm(data_m$hm.l~data_m$s.v.l), col = 'red', lwd = 2) #添加趋势线
legend(0, 28, legend = c('female', 'male'), col = c('blue', 'red'), pch = c(1, 4), bty = "n") #添加图例
text(15, 29, 'A', cex = 1.5, font = 2) #左上角添加'A'
#图B
plot(x = data_f$s.v.l, y = data_f$fm.l, xlim = c(0, 250), ylim = c(0, 25), pch = 1,
col = 'blue', xlab = 'Snout-vent length (mm)' , ylab = 'Femur length (mm)')
abline(lm(data_f$fm.l~data_f$s.v.l), col = 'blue', lwd = 2, h = , v = ) #添加趋势线
par(new = TRUE)
plot(x = data_m$s.v.l, y = data_m$fm.l, xlim = c(0, 250), ylim = c(0, 25), pch = 4,
col = 'red', xlab = '', ylab = '', xaxt = 'n', yaxt = 'n')
abline(lm(data_m$fm.l~data_m$s.v.l), col = 'red', lwd = 2) #添加趋势线
text(15, 24, 'B', cex = 1.5, font = 2) #左上角添加'B'
#图C
plot(x = data_f$fm.l, y = data_f$hm.l, xlim = c(0, 25), ylim = c(0, 30), pch = 1,
col = 'blue', xlab = 'Femur length (mm)' , ylab = 'Humerus length (mm)')
abline(lm(data_f$hm.l~data_f$fm.l), col = 'blue', lwd = 2, h = , v = ) #添加趋势线
par(new = TRUE)
plot(x = data_m$fm.l, y = data_m$hm.l, xlim = c(0, 25), ylim = c(0, 30), pch = 4,
col = 'red', xlab = '', ylab = '', xaxt = 'n', yaxt = 'n')
abline(lm(data_m$hm.l~data_m$fm.l), col = 'red', lwd = 2) #添加趋势线
text(1, 29, 'C', cex = 1.5, font = 2) #左上角添加'C'
#图D
plot(x = data_f$hm.m, y = data_f$hm.d, xlim = c(0, 5), ylim = c(0, 9), pch = 1,
col = 'blue', xlab = 'Humerus mid-shaft width (mm)' , ylab = 'Humerus distal width (mm)')
abline(lm(data_f$hm.d~data_f$hm.m), col = 'blue', lwd = 2, h = , v = ) #添加趋势线
par(new = TRUE)
plot(x = data_m$hm.m, y = data_m$hm.d, xlim = c(0, 5), ylim = c(0, 9), pch = 4,
col = 'red', xlab = '', ylab = '', xaxt = 'n', yaxt = 'n')
abline(lm(data_m$hm.d~data_m$hm.m), col = 'red', lwd = 2) #添加趋势线
text(0.2, 8.5, 'D', cex = 1.5, font = 2) #左上角添加'D'
协方差分析（ANCOVA）对异速生长的解释。胡氏贵州龙在肱骨/吻臀长度比（图a）、股骨/吻臀长度比（图b）、肱骨/股骨长度比（图c）和肱骨远端宽/肱骨中轴宽度比（图d）四项比列上均体现了雌雄的差异（处在不同的回归线上），但是肱骨远端宽/肱骨中轴宽度的比例从一开始就出现了分离，表明胡氏贵州龙在这一指标的上的差异，可以倒追回胚胎期。换言之，胡氏贵州龙肱骨远端的宽度自出生起就存在着性别的差异。
3 本章汇总
| 参数 | 类别 | 功能 |
|---|---|---|
| readxl | 函数包 | 加载readxl包 |
| MASS | 函数包 | 加载MASS包 |
| ggplot2 | 函数包 | 加载ggplot2包 |
| car | 函数包 | 加载car包 |
| for | 函数 | 循环函数 |
| predict | 函数 | 对数据进行预测 |
| dataEllipse | 函数 | 置信椭圆 |
| abline | 函数 | 添加斜线 |
4 参考文献
[1] Xue, Y., Jiang, D., Motani, R., Rieppel, O., Sun, Y., Sun, Z., Ji, C., and Yang, P. New information on sexual dimorphism and allometric growth in Keichousaurus hui, a pachypleurosaur from the Middle Triassic of Guizhou, South China[J].Acta Palaeontologica Polonica 60 (3): 681–687.2015
[2] https://baike.baidu.com/item/线性判别分析/22657333