bifurcatingr

Introduction

The bifurcating autiegressive (BAR) model is a useful tool for studying binary tree-structured data, often seen in cell-lineage studies. The “bifurcatingr” R package makes it easier to work with these models. It can estimate different orders of BAR models and improve accuracy through various correction methods. These methods include bootstrap-based corrections and linear bias functions. The package also helps generate and visualize data from BAR models and offers eight ways to estimate confidence intervals for the model’s coefficients, with or without bias correction.

Installation

To install the library from CRAN, use:

install.packages("bifurcatingr")

Usage

To use bifurcatingr in your R code, you must first load the library:

library(bifurcatingr)
#> Loading required package: fMultivar

Generating binary tree data

To generate binary tree data based on a the bifurcating autoregressive model (BAR) model, you can use the “bfa.tree.gen” function. Generating a binary tree with 127 observations based on a first-order BAR model with errors that follow a contaminated bivariate normal distribution with 20% contamination. A zero mean vector is used and the standard deviations are 1 and 2, respectively. The intercept is set to 10 and the autoregressive coefficient is 0.7.

dat<-bfa_tree_gen(n=127,p=1,s1=1,s2=2,r1=0.5,r2=0.5,g=0.2,intercept=10,ar_coef=c(0.7))
dat <- round(dat,1)
dat
#>   [1] 30.1 31.7 30.3 31.4 31.6 32.0 31.1 32.4 30.9 32.2 31.2 31.3 32.0 31.8 29.9
#>  [16] 32.9 31.7 30.7 31.7 31.1 30.6 32.9 32.9 31.3 31.9 31.7 31.3 32.2 32.9 29.9
#>  [31] 30.4 32.5 33.2 28.4 31.7 29.8 31.5 31.3 32.3 33.1 33.7 30.9 31.7 32.4 32.3
#>  [46] 31.3 32.5 31.5 31.1 32.2 32.4 33.3 32.0 32.2 29.6 32.9 32.9 33.7 32.9 31.9
#>  [61] 32.1 32.5 30.6 32.8 33.7 32.2 33.1 29.6 30.1 32.9 31.7 30.6 31.3 32.1 32.7
#>  [76] 30.4 31.4 31.0 28.3 36.8 32.7 33.5 34.6 33.3 31.5 32.2 32.2 32.2 32.9 32.2
#>  [91] 32.3 34.0 33.2 31.7 31.9 31.3 32.1 33.8 34.4 32.4 34.4 32.5 32.8 32.4 32.9
#> [106] 32.4 31.5 32.6 32.1 30.0 30.3 32.3 33.7 33.8 32.6 33.4 32.6 32.3 31.0 33.0
#> [121] 32.6 33.8 32.4 34.3 33.8 30.7 31.1

Graphing binary tree data

The function “bfa.tree.plot” is designed to represent the binary tree data usong igraph library. By making use of the cell lineage dataset found in Cowan and Staudte (1986) and stored in the bifurcatingr package as “ecoli”, executing the following command will generate the tree.

bfa_tree_plot(ecoli$lifetime, vertex.size = 10, shape = "circle", text.size = 0.3, vertex.color = "white",vertex.asp = 0.5, arrow.size = 0.2, arrow.width = 0.3, plot.margin = 0)

Graphing scatterplots

The BAR model can be visualized using a scatterplot depicting the relationship between observations at time “t” and their lagged observations. This is done through the “bfa.scatterplot” function. This function requires two inputs: the binary tree data as a numeric vector “z” and the lag order “p”. The following command generates a scatterplot for the genrated binary tree above based on the BAR(1) model:

bfa_scatterplot(dat,1)

The following command generat scatterplot based on the BAR(2) model for the same data above.

bfa_scatterplot(dat,2)

Least squares estimation for BAR models

The “bfa.ls” function is used for obtaining least squares estimates for the parameters of a BAR model of any order. This function yields the LS estimates for the BAR model coefficients along with their corresponding p-values for testing if the coefficient is different from zero, the correlation among the errors, the variance-covariance matrix of the estimates, and both Wald-type asymptotic and standard normal or percentile) bootstrap confidence intervals. The following command illustrates this function on the simulated data above “dat” when fitted by a BAR(1) model.

bfa_ls(dat,p=1,conf=TRUE,cov_matrix=TRUE,conf_level=0.9,p_value=TRUE)
#> $coef
#>      Intercept   X_[t/2]
#> [1,]  12.71179 0.6102418
#> 
#> $p_value
#>        Intercept      X_[t/2]
#> [1,] 0.002054208 2.650564e-06
#> 
#> $error_cor
#> [1] 0.3466633
#> 
#> $cov_matrix
#>            Intercept   X_[t/2]
#> Intercept 2160.09389 -68.02332
#> X_[t/2]    -68.02332   2.14445
#> 
#> $asymptotic_ci
#>                 5%        95%
#> Intercept 5.928166 19.4954124
#> X_[t/2]   0.396503  0.8239806
#> 
#> $bootstarp_ci
#>                  5%       95%
#> Intercept 8.6433373 16.780241
#> X_[t/2]   0.4848966  0.735587
#> 
#> $percentile_ci
#>                  5%        95%
#> Intercept 9.8782880 17.3159189
#> X_[t/2]   0.4703033  0.6983381

Least squares bias correction for the coefficients of BAR models

The “bfa.ls.bc” function provides bias-corrected LS estimates of the BAR model coefficients. Currently, four distinct bias correction techniques are available: single, double, and fast-double bootstrapping, along with the linear-bias-function approach. The “bfa.ls.bc” function to obtain the bias-corrected estimates of the AR coefficient in BAR(1) model based on the single bootstrap bias correction method for the simulated data above “dat” is as follows

bfa_ls_bc(dat,p = 1, method = "boot1",B1 = 999, B2 = 499)
#> [1] 0.6376565

The “bfa.ls.bc” function to obtain the bias-corrected estimates of the AR coefficient in BAR(1) model based on the fast double bootstrap bias correction method for the simulated data above “dat” is as follows

bfa_ls_bc(dat,p = 1, method = "boot2fast",B1 = 999, B2 = 499)
#> [1] 0.6326482

The “bfa.ls.bc” function to obtain the bias-corrected estimates of the AR coefficient in BAR(1) model based on the asymptotic linear bias-correction method for the simulated data above “dat” is as follows

bfa_ls_bc(dat, p=1, method="LBC")
#> [1] 0.6412441

Confidence Intervals

The bifurcatingr package allows for constructing various types of confidence intervals based on the LS estimation of the BAR model coefficients. Both corrected and uncorrected bias confidence intervals are available through the “bfa.ls.bc.ci” function. For instance, the following command can be used to construct a 95% BCa confidence interval for the autoregressive coefficient of the BAR(1) model fitted to the simulated data abouve “dat”.

bfa_ls_bc_ci(dat, p=1, method="BCa", B=4)
#> $Bias_corrected_coef
#>         Estimates
#> X_[t/2] 0.6538157
#> 
#> $BCa_ci
#>              2.5%     97.5%
#> X_[t/2] 0.5156971 0.6655885

To construct a 95% single bootstrap bias-corrected standard normal bootstrap confidence intervals for the autoregressive coefficients of the BAR(1) model for the simulated data above “dat” as follows.

bfa_ls_bc_ci(dat, p=1, method="boot1", B=4)
#> $Bias_corrected_coef
#>           [,1]
#> [1,] 0.5266351
#> 
#> $standard_normal_ci
#>              2.5%     97.5%
#> X_[t/2] 0.4598922 0.5933781

Where to go next

This tutorial is a brief introduction to bifurcatingr in R. We sincerely hope you enjoyed reading it and that it will be useful for your binary tree analysis based on the autoregressive models. For a detailed description of specific functions, see https://cran.r-project.org/package=bifurcatingr. For questions on how to use bifurcatingr, or to report a bug, please email Dr. Tamer Elbayoumi at email: tmelbayoumi@ncat.edu.