2020-10-16

个体发育分析

个体发育，指的是一种生物从其出生到死亡的整个过程中身体出现的形态变化。从幼体到成年，生物的形态变化通常比较明显。在古生物学中研究个体发育，一方面有助于理解不同化石物种的生长变化；另一方面，也有助于在后续研究辅助化石的鉴定。在古生物的个体发育研究中，通常会总结一种生物不同身体部位随年龄增长的变化趋势、这种趋势或是某些器官的增大，或是某些结构数量的增加。由于生物的种类千差万别，因此不同古生物的个体发育过程也存在很大差异，需要根据实际情况选用适宜的衡量标准。

1 数据描述

#先安装包：install.packages("PerformanceAnalytics")
setwd("C:\\Users\\Lenovo\\Desktop\\data_3")
library(PerformanceAnalytics)

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

library(readxl)
data=read_excel("Ontogeny_echinoderms.xlsx")
head(data)

## # A tibble: 6 x 6
##   Specimen_No. thecal_height stem_length `minimun N_brac~ circlets_plates~
##          <dbl>         <dbl>       <dbl>            <dbl>            <dbl>
## 1           10         0.828       0.666                6               10
## 2           15         0.625       1.16                 5                7
## 3           16         0.93        0.871                7                6
## 4           23         0.936       0.892                7                7
## 5           25         1.74        1                   12               10
## 6           27         1.73        1.12                11               11
## # ... with 1 more variable: number_part <dbl>

str(data)

## Classes 'tbl_df', 'tbl' and 'data.frame':    18 obs. of  6 variables:
##  $ Specimen_No.        : num  10 15 16 23 25 27 32 33 39 40 ...
##  $ thecal_height       : num  0.828 0.625 0.93 0.936 1.743 ...
##  $ stem_length         : num  0.666 1.16 0.871 0.892 1 ...
##  $ minimun N_brachioles: num  6 5 7 7 12 11 4 6 7 7 ...
##  $ circlets_plates_No  : num  10 7 6 7 10 11 8 9 7 7 ...
##  $ number_part         : num  7 6 7 6 7 7 6 6 6 6 ...

数据共有18条记录，6个变量，数据类型都为数值型数据。为描述各变量之间的相关性情况，对数据进行形态相关分析如下：

2 个体变量相关性分析

chart.Correlation(data, pch="+")

选择生物的个体发育指标，得到形态相关的相关性值，可以了解到stem_length和thecal_height相关性为0.8，二者相关性强度较高，其次是thecal_height和minimun N_brachioles的相关性最好。

3 个体发育趋势分析

利用贝叶斯信息准则（BIC）作出的频率分布图利用贝叶斯信息准则（BIC）作出的频率分布图，探究各变量之间的正太性和核密度分布，如下图所示：

setwd("C:\\Users\\Lenovo\\Desktop\\data_3")
data1=read_excel("Ontogeny_echinoderms_1.xlsx")
str(data1)

## Classes 'tbl_df', 'tbl' and 'data.frame':    9 obs. of  7 variables:
##  $ X : num  0.2 2.3 3.7 4 5 6 7 8 9.8
##  $ A1: num  -54 -46 -58 -69 -68 -80 -98 -108 -118
##  $ A2: num  -54 -46 -55 -62 -65 -72 -80 -90 -93
##  $ B1: num  -153 -146 NA NA NA NA NA NA NA
##  $ B2: num  -153 -142 -155 -158 -155 -158 -167 -173 -180
##  $ C1: num  -30 -41 -49 NA NA NA NA NA NA
##  $ C2: num  -30 -37 -43 -48 -56 -65 -66 -69 -74

layout(matrix(c(1,2,3,4,5,6),2,3,byrow = TRUE))
plot(data1$X,data1$A1,type="b",yaxt="n",ylab="Precent  surviviing",xlab="Length   (mm)",pch=1,lwd=2,lty=1)
par(new=TRUE)  # 是否叠加新图形，没叠加一次运行一次该命令
plot(data1$X,data1$A2,type="b",ylab="Precent  surviviing",xlab="Length   (mm)",pch=4,lwd=1,lty=6)
text(8, -50,"A")

plot(data1$X,data1$B1,type="b",yaxt="n",ylab="Precent  surviviing",xlab="Length   (mm)",pch=1,lwd=2,lty=1)
par(new=TRUE)  # 是否叠加新图形，没叠加一次运行一次该命令
plot(data1$X,data1$B2,type="b",ylab="Precent  surviviing",xlab="Length   (mm)",pch=4,lwd=1,lty=6)
text(8,-146, "B")

plot(data1$X,data1$C1,type="b",ylab="Precent  surviviing",xlab="Length   (mm)",pch=1,lwd=2,lty=1)
par(new=TRUE)  # 是否叠加新图形，没叠加一次运行一次该命令
plot(data1$X,data1$C2,type="b",ylab="Precent  surviviing",xlab="Length   (mm)",pch=4,lwd=1,lty=6)
text(8, -33, "C")

hist(data$thecal_height, prob=T, xlab='',main='Histogram of maximum pH value')#画直方图
lines(density(data$thecal_height,na.rm=T))#画概率密度曲线
text(1.5, 1.5, "D")

hist(data$circlets_plates_No, prob=T, xlab='',main='Histogram of maximum pH value',ylim=0:1)#画直方图
lines(density(data$circlets_plates_No,na.rm=T))#画概率密度曲线
text(10, 0.9, "E")

hist(data$ stem_length, prob=T, xlab='',main='Histogram of maximum pH value')#画直方图
lines(density(data$ stem_length,na.rm=T))#画概率密度曲线
text(1.5, 1.5, "F")

五组测量变量间的一一对应关系如上所示，说明长度在5到11的呈现右偏分布，0.5到1.5的基本呈现正态分布。

3.1 皮尔逊相关性检验分析

#要先安装basicTrendline包
library(basicTrendline)

## Warning: package 'basicTrendline' was built under R version 3.5.3

trendline(data$thecal_height,data$stem_length, ,model="line3P",Pvalue.corrected = TRUE,
          linecolor = "yellow", lty = 0, lwd = 1, show.equation = FALSE,
          show.Rpvalue = FALSE, Rname = 1, Pname = 0, xname = "x", yname = "y",
          yhat = FALSE, summary = FALSE, ePos.x = NULL, ePos.y = NULL,
          text.col = "black", eDigit = 5, eSize = 1, CI.fill = TRUE,
          CI.level = 0.95, CI.color = "white", CI.alpha = 1, CI.lty = 1,
          CI.lwd = 1, las = 1, xlab = NULL, ylab = NULL,xaxt="n",yaxt="n")

由上图知，皮尔相关系数为0.67609,P值为0.000212.说明变量之间存在很强的显著性差异，变量之间的相关性较好。

3.2 回归系数分析结果

library(basicTrendline)
trendline(log(data$thecal_height),log(data$stem_length), Pvalue.corrected = TRUE,
          linecolor = "black", lty = 1, lwd = 2, show.equation = FALSE,
          show.Rpvalue = TRUE, Rname = 1, Pname = 0, xname = "x", yname = "y",
          yhat = FALSE, summary = FALSE, ePos.x = NULL, ePos.y = NULL,
           eDigit = 5, eSize = 1, CI.fill = TRUE,
          CI.level = 0.95, CI.color = "white", CI.alpha = 1, CI.lty = 1,
          CI.lwd = 1, las = 1, xlab ="", ylab ="",xaxt="n",yaxt="n")
par(new=TRUE)
##95%的结果
mode=lm(log(data$thecal_height)~log(data$stem_length))
trendline(fitted(mode),log(data$stem_length),CI.fill = FALSE,CI.color = "red",CI.lty = 2,linecolor = "red",summary =FALSE,Pvalue.corrected = FALSE, show.equation = FALSE,xlab = "",text.col = "white")

## Warning in summary.lm(fit): essentially perfect fit: summary may be
## unreliable

par(new=TRUE)

trendline(log(data$thecal_height),log(data$stem_length),model="line3P",ePos.x = "topleft",summary =FALSE,eDigit = 5,CI.fill = FALSE,text.col = "white",ylab =NULL,xaxt="n",yaxt="n",Pvalue.corrected = FALSE,show.Rpvalue =FALSE,show.equation = FALSE)

3.3 拟合回归分析

library(basicTrendline)
##OLS回归的结果
mode=lm(log(data$thecal_height)~log(data$stem_length))
trendline(fitted(mode),log(data$stem_length),text.col = "white",CI.fill = FALSE,CI.color = "white",CI.lty = 2,linecolor = "black",xlab="",xaxt="n",yaxt="n",lty=1)

## Warning in summary.lm(fit): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(fit): essentially perfect fit: summary may be
## unreliable

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         Min          1Q      Median          3Q         Max 
## -2.1436e-15  5.5999e-17  1.0999e-16  1.6914e-16  7.3566e-16 
## 
## Coefficients:
##                Estimate  Std. Error     t value  Pr(>|t|)    
## (Intercept) -1.8396e-01  1.4373e-16 -1.2799e+15 < 2.2e-16 ***
## x            1.3685e+00  3.4129e-16  4.0098e+15 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.9521e-16 on 16 degrees of freedom
## Multiple R-squared:       1, Adjusted R-squared:       1 
## F-statistic: 1.6078e+31 on 1 and 16 DF,  p-value: < 2.22e-16
## 
## 
## N: 18 , AIC: -1207.1 , BIC:  -1204.4 
## Residual Sum of Squares:  5.6684e-30

# 95%置信信区间
par(new=TRUE)
trendline(log(data$thecal_height),log(data$stem_length),model="line3P",ePos.x = "topleft",summary =FALSE,eDigit = 5 ,CI.fill = "FALSE",text.col = "white",xaxt="n",linecolor = "white",yaxt="n",CI.color = "green",lty=2)
# 取对数的拟合曲线
par(new=TRUE)
trendline(log(data$thecal_height),log(data$stem_length),CI.fill = FALSE,CI.color = "red",CI.lty = 5,linecolor = "white",Pvalue.corrected = TRUE, show.equation = FALSE,
          show.Rpvalue = FALSE)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.416273 -0.194345 -0.008746  0.136046  0.850577 
## 
## Coefficients:
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -0.214212   0.070382 -3.0435  0.007742 ** 
## x            1.038171   0.146400  7.0914 2.554e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.29314 on 16 degrees of freedom
## Multiple R-squared:  0.75863,    Adjusted R-squared:  0.74354 
## F-statistic: 50.287 on 1 and 16 DF,  p-value: 2.5539e-06
## 
## 
## N: 18 , AIC: 10.786 , BIC:  13.457 
## Residual Sum of Squares:  1.3749

legend("bottomright",c("Observed","OLS fit","95% CL","log() fit"), col=c("black","black","green","red"),
text.col = "green4", pch = c(1, 1,1,1),lty = c(0,1,2,5),bty="n",
merge = F)

通过拟合回归分析得到可决系数值为75.863%，P值远小于0.0001。变量之间的相关系数较好。

4 本章汇总

参数	类别	功能
readxl	函数包	加载readxl包
PerformanceAnalytics	函数包	制作相关性可视化图形
layout	函数	绘制图形分栏布局
text	函数	添加图例的文本标签
lines	函数	绘制线图
par	函数	叠加新图形

5 参考文献

[1]https://blog.csdn.net/tandelin/article/details/99547342

源文件下载在线调试

返回：数据分析教程案例, R 语言