Movatterモバイル変換

[0]ホーム
URL:

Analyzing FBD Parameters

2025-08-27

This vignette explains how to extract, plot, and statistically testfor differences among BD (Birth-Death) or FBD (fossilized birth-death)parameters (e.g.,net diversification,relative extinction(turnover), andrelative fossilization) across the treewhen using the Skyline BD or Skyline FBD tree models produced by theprogramMr. Bayes orBEAST2 (SA or BDSKY packages).

Load theEvoPhylo package

library(EvoPhylo)

FBD Parameters Statistics and Plots (Mr. Bayes)

Below we demonstrate how to extract evolutionary rate summarystatistics from each node from a Bayesian clock (time-calibrate) summarytree produced by Mr. Bayes, store them in a data frame, produce summarytables, and plots.

1. Import combined log file from all runs.

This is produced by usingcombine_log().The firstargument passed tocombine_log() should be a path to thefolder containing the log files to be imported and combined.

## Import all log (.p) files from all runs and combine them, with burn-in = 25%## and downsampling to 2.5k trees in each log fileposterior3p<-combine_log("LogFiles3p",burnin =0.25,downsample =1000)

Below, we use the posterior datasetposterior3p thataccompaniesEvoPhylo.

data(posterior3p)## Show first 5 lines of combined log filehead(posterior3p,5)

The posterior data must first be transformed from wide to long to beused with the functions described below;FBD_reshape()accomplishes this.

## Reshape imported combined log file from wide to long with FBD_reshapeposterior3p_long<-FBD_reshape(posterior3p,variables =NULL,log.type ="MrBayes")

2. Summarize FBD parameters by time bin

Summary statistics for each FBD parameter by time bin can be quicklysummarized usingFBD_summary():

## Summarize parameters by time bin and analysist3.1<-FBD_summary(posterior3p_long)t3.1
FBD parameters by time bin
parameterTime_binnmeansdminQ1medianQ3max
net_speciation1100000.040.020.000.030.040.060.17
net_speciation2100000.030.020.000.020.030.040.12
net_speciation3100000.020.010.000.010.020.030.12
net_speciation4100000.050.020.000.030.050.060.12
relative_extinction1100000.790.150.080.710.820.901.00
relative_extinction2100000.930.050.550.900.930.961.00
relative_extinction3100000.950.050.180.930.960.981.00
relative_extinction4100000.030.100.000.000.000.010.97
relative_fossilization1100000.040.050.000.010.020.050.72
relative_fossilization2100000.070.040.000.040.060.090.36
relative_fossilization3100000.010.020.000.000.010.020.54
relative_fossilization4100000.040.110.000.000.000.020.99
## Export the tablewrite.csv(t3.1,file ="FBD_summary.csv")

3. Plot the distribution of each FBD parameter

Each of (or all) the FBD parameter distributions can be plotted bytime bin using various plotting alternatives withFBD_dens_plot():

## Plot distribution of the desired FBD parameter by time bin with## kernel density plotFBD_dens_plot(posterior3p_long,parameter ="net_speciation",type ="density",stack =FALSE)

## Plot distribution of the desired FBD parameter by time bin with## stacked kernel density plotFBD_dens_plot(posterior3p_long,parameter ="net_speciation",type ="density",stack =TRUE)

## Plot distribution of the desired FBD parameter by time bin with## a violin plotFBD_dens_plot(posterior3p_long,parameter ="net_speciation",type ="violin",stack =FALSE,color ="red")

## Plot distribution of all FBD parameter by time bin with a violin plotp1<-FBD_dens_plot(posterior3p_long,parameter ="net_speciation",type ="violin",stack =FALSE,color ="red")p2<-FBD_dens_plot(posterior3p_long,parameter ="relative_extinction",type ="violin",stack =FALSE,color ="cyan3")p3<-FBD_dens_plot(posterior3p_long,parameter ="relative_fossilization",type ="violin",stack =FALSE,color ="green3")library(patchwork)p1+ p2+ p3+plot_layout(nrow =1)

## Save your plot to your working directory as a PDFggplot2::ggsave("Plot_regs.pdf",width =12,height =4)

4. Test for assumptions

In this step, users can perform tests for normality andhomoscedasticity in data distribution for each of the FBD parametersunder consideration. The output will determine whether parametric ornonparametric tests will be performed subsequently.

##### Tests for normality and homoscedasticity for each FBD parameter for all time binst3.2<-FBD_tests1(posterior3p_long)
### Export the output table for all testswrite.csv(t3.2,file ="FBD_Tests1_Assum.csv")

The results of the Shapiro-Wilk normality test for each parameter canbe output as seperate tables or as a single combined table.

# Output as separate tablest3.2$shapiro
Shapiro-Wilk normality test
parameterstatisticp-value
Time bin 1net_speciation0.99170
Time bin 2net_speciation0.93850
Time bin 3net_speciation0.92270
Time bin 4net_speciation0.98980
Overallnet_speciation0.95680
Residualsnet_speciation0.98740
parameterstatisticp-value
Time bin 1relative_extinction0.89270
Time bin 2relative_extinction0.92470
Time bin 3relative_extinction0.80440
Time bin 4relative_extinction0.37750
Overallrelative_extinction0.70360
Residualsrelative_extinction0.82380
parameterstatisticp-value
Time bin 1relative_fossilization0.57640
Time bin 2relative_fossilization0.88530
Time bin 3relative_fossilization0.62100
Time bin 4relative_fossilization0.46370
Overallrelative_fossilization0.54730
Residualsrelative_fossilization0.55310
# OR as single merged tablet3.2$shapiro$net_speciation$bin<-row.names(t3.2$shapiro$net_speciation)t3.2$shapiro$relative_extinction$bin<-row.names(t3.2$shapiro$relative_extinction)t3.2$shapiro$relative_fossilization$bin<-row.names(t3.2$shapiro$relative_fossilization)k1all<-rbind(t3.2$shapiro$net_speciation,               t3.2$shapiro$relative_extinction,               t3.2$shapiro$relative_fossilization,make.row.names =FALSE)
k1all
Shapiro-Wilk normality test
parameterstatisticp-valuebin
net_speciation0.99170Time bin 1
net_speciation0.93850Time bin 2
net_speciation0.92270Time bin 3
net_speciation0.98980Time bin 4
net_speciation0.95680Overall
net_speciation0.98740Residuals
relative_extinction0.89270Time bin 1
relative_extinction0.92470Time bin 2
relative_extinction0.80440Time bin 3
relative_extinction0.37750Time bin 4
relative_extinction0.70360Overall
relative_extinction0.82380Residuals
relative_fossilization0.57640Time bin 1
relative_fossilization0.88530Time bin 2
relative_fossilization0.62100Time bin 3
relative_fossilization0.46370Time bin 4
relative_fossilization0.54730Overall
relative_fossilization0.55310Residuals
## Bartlett's test for homogeneity of variancet3.2$bartlett
Bartlett’s test
parameterstatisticp-value
net_speciation3815.4640
relative_extinction18159.2130
relative_fossilization25654.9750
## Fligner-Killeen test for homogeneity of variancet3.2$fligner
Fligner-Killeen test
parameterstatisticp-value
net_speciation3748.1400
relative_extinction12599.8430
relative_fossilization4808.5450

Deviations from normality can be displayed graphically usingFBD_normality_plot():

## Visualize deviations from normality and similarity of variancesFBD_normality_plot(posterior3p_long)

## Save your plot to your working directory as a PDFggplot2::ggsave("Plot_normTests.pdf",width =8,height =6)

5. Test for significant FBD shifts between time bins

Significant shifts in FBD parameters across time bins can be easilytested using parametric (Student’s t-test) and nonparametric(Mann-Whitney test) pairwise comparisons withFBD_tests2().Both are automatically calculated and the preferred pairwise comparisonwill be chosen by the user depending on the results of the assumptiontestsstep #4 (above).

##### Test for significant differences between each time bin for each FBD parametert3.3<-FBD_tests2(posterior3p_long)
### Export the output table for all testswrite.csv(t3.3,file ="FBD_Tests2_Sign.csv")## Pairwise t-tests# Output as separate tablest3.3$t_tests
Significant tests
parameterTime_bin1Time_bin2n1n2p-valuep-value adj
net_speciation12100001000000
net_speciation13100001000000
net_speciation14100001000000
net_speciation23100001000000
net_speciation24100001000000
net_speciation34100001000000
parameterTime_bin1Time_bin2n1n2p-valuep-value adj
relative_extinction12100001000000
relative_extinction13100001000000
relative_extinction14100001000000
relative_extinction23100001000000
relative_extinction24100001000000
relative_extinction34100001000000
parameterTime_bin1Time_bin2n1n2p-valuep-value adj
relative_fossilization12100001000000
relative_fossilization13100001000000
relative_fossilization14100001000000
relative_fossilization23100001000000
relative_fossilization24100001000000
relative_fossilization34100001000000
# OR as single merged tablek3.3a<-rbind(t3.3$t_tests$net_speciation,               t3.3$t_tests$relative_extinction,               t3.3$t_tests$relative_fossilization,make.row.names =FALSE)
k3.3a
Pairwise t-tests
parameterTime_bin1Time_bin2n1n2p-valuep-value adj
net_speciation12100001000000
net_speciation13100001000000
net_speciation14100001000000
net_speciation23100001000000
net_speciation24100001000000
net_speciation34100001000000
relative_extinction12100001000000
relative_extinction13100001000000
relative_extinction14100001000000
relative_extinction23100001000000
relative_extinction24100001000000
relative_extinction34100001000000
relative_fossilization12100001000000
relative_fossilization13100001000000
relative_fossilization14100001000000
relative_fossilization23100001000000
relative_fossilization24100001000000
relative_fossilization34100001000000
## Mann-Whitney tests (use if Tests in step #4 fail assumptions)# Output as separate tablest3.3$mwu_tests
Mann-Whitney tests
parameterTime_bin1Time_bin2n1n2p-valuep-value adj
net_speciation12100001000000
net_speciation13100001000000
net_speciation14100001000000
net_speciation23100001000000
net_speciation24100001000000
net_speciation34100001000000
parameterTime_bin1Time_bin2n1n2p-valuep-value adj
relative_extinction12100001000000
relative_extinction13100001000000
relative_extinction14100001000000
relative_extinction23100001000000
relative_extinction24100001000000
relative_extinction34100001000000
parameterTime_bin1Time_bin2n1n2p-valuep-value adj
relative_fossilization12100001000000
relative_fossilization13100001000000
relative_fossilization14100001000000
relative_fossilization23100001000000
relative_fossilization24100001000000
relative_fossilization34100001000000
# OR as single merged tablek3.3b<-rbind(t3.3$mwu_tests$net_speciation,               t3.3$mwu_tests$relative_extinction,               t3.3$mwu_tests$relative_fossilization,make.row.names =FALSE)
k3.3b
Mann-Whitney tests
parameterTime_bin1Time_bin2n1n2p-valuep-value adj
net_speciation12100001000000
net_speciation13100001000000
net_speciation14100001000000
net_speciation23100001000000
net_speciation24100001000000
net_speciation34100001000000
relative_extinction12100001000000
relative_extinction13100001000000
relative_extinction14100001000000
relative_extinction23100001000000
relative_extinction24100001000000
relative_extinction34100001000000
relative_fossilization12100001000000
relative_fossilization13100001000000
relative_fossilization14100001000000
relative_fossilization23100001000000
relative_fossilization24100001000000
relative_fossilization34100001000000

 
 
 

FBD Parameters Statistics and Plots (BEAST2)

1. Import combined log file from all runs.

The combined posterior log file fromBEAST2is outputted byLogCombiner from their softwarepackage. Our own function to combined log filescombine_logis intended to work with Mr. Bayes posterior files only.

Below, we use the posterior dataset “Penguins_log.log” thataccompaniesEvoPhylo.

posterior<-system.file("extdata","Penguins_log.log",package ="EvoPhylo")posterior<-read.table(posterior,header =TRUE)## Show first 10 lines of combined log filehead(posterior,5)

The posterior data must first be transformed from wide to long to beused with the functions described below;FBD_reshape()accomplishes this.

## Reshape imported combined log file from wide to long with FBD_reshapeposterior_long<-FBD_reshape(posterior,variables =NULL,log.type ="BEAST2")

2. Summarize FBD parameters by time bin

Summary statistics for each FBD parameter by time bin can be quicklysummarized usingFBD_summary():

## Summarize parameters by time bin and analysist3.1<-FBD_summary(posterior_long)t3.1
FBD parameters by time bin
parameterTime_binnmeansdminQ1medianQ3max
diversificationRateFBD15560.040.030.000.020.040.060.17
turnoverFBD15560.860.100.450.810.890.941.00
samplingProportionFBD15560.270.120.060.190.250.340.74
## Export the tablewrite.csv(t3.1,file ="FBD_summary_BEAST2.csv")

3. Plot the distribution of each FBD parameter

Each of (or all) the FBD parameter distributions can be plotted bytime bin using various plotting alternatives withFBD_dens_plot():

## Plot distribution of the desired FBD parameter by time bin with## kernel density plotFBD_dens_plot(posterior_long,parameter ="diversificationRateFBD",type ="density",stack =FALSE)

## Plot distribution of the desired FBD parameter by time bin with## a violin plotFBD_dens_plot(posterior_long,parameter ="diversificationRateFBD",type ="violin",stack =FALSE,color ="red")

## Plot distribution of all FBD parameter by time bin with a violin plotp1<-FBD_dens_plot(posterior_long,parameter ="diversificationRateFBD",type ="violin",stack =FALSE,color ="red")p2<-FBD_dens_plot(posterior_long,parameter ="turnoverFBD",type ="violin",stack =FALSE,color ="cyan3")p3<-FBD_dens_plot(posterior_long,parameter ="samplingProportionFBD",type ="violin",stack =FALSE,color ="green3")library(patchwork)p1+ p2+ p3+plot_layout(nrow =1)

## Save your plot to your working directory as a PDFggplot2::ggsave("Plot_regs.pdf",width =12,height =4)

References

[8]ページ先頭
©2009-2025 Movatter.jp