# This code can be used to reproduce the results of the two applications in the manuscript.

library(parallel)
library(dplyr)
library(tidyr)
library(bifurcatingr)

bfa.boot.ls.bc <- function(z, p, method="boot1",burn = 5, B, boot.est=TRUE, boot.data=FALSE){  # burn can be determined by simulation
  n=length(z)
  boot.dat <- matrix(rep(NA, n*B), ncol=n)
  est <- bfa.ls(z,p,resids=TRUE)
  intercept <- est$coef[1,1]
  coef <- t(as.vector(est$coef[1,2:(p+1)]))
  e <- est$resids-mean(est$resids,na.rm=TRUE)
  e.e=matrix(e,ncol=2, byrow=T)
  ee <- e.e[stats::complete.cases(e.e),]
  for (h in 1:B){
    # start-up tree
    n0 = 2^(p + burn) - 1
    term.matrix0 <- matrix(rep(NA,p*n0),ncol=p)
    dat0 <- rep(NA,n0)
    er0 <- matrix(rep(NA,((n0-1)/2)),ncol=2,nrow=((n0-1)/2))
    samp0 <- sample(1:nrow(ee), ((n0-1)/2) ,replace=T)
    for (s in 1: nrow(er0)) {
      l <- samp0[s]
      er0[s,] <- ee[l,]
    }
    error0 <- rep(NA,n0)
    for(i in 1:((n0-1)/2)){
      error0[(i*2)] <- er0[i,1]
      error0[(i*2)+1] <- er0[i,2]
    }
    dat0[1:((2^p)-1)] <- z[1:((2^p)-1)]
    for(j in (2^p) : n0){
      for(k in 1:p){
        term.matrix0[,k] <- coef[1,k]*dat0[floor(j/(2^k))]
      }
      dat0[j] <- intercept+sum(term.matrix0[j,])+error0[j]
    }
    # bootstrap tree
    term.matrix <- matrix(rep(NA,p*n),ncol=p)
    dat <- rep(NA,n)
    er <- matrix(rep(NA,((n-1)/2)),ncol=2,nrow=((n-1)/2))
    samp <- sample(1:nrow(ee), ((n-1)/2) ,replace=T)
    for (s in 1: nrow(er)) {
      l <- samp[s]
      er[s,] <- ee[l,]
    }
    error <- rep(NA,n)
    for(i in 1:((n-1)/2)){
      error[(i*2)] <- er[i,1]
      error[(i*2)+1] <- er[i,2]
    }
    dat[1:((2^p)-1)] <- dat0[bfa.subtree(n0,p)]
    for(j in (2^p) : n){
      for(k in 1:p){
        term.matrix[,k] <- coef[1,k]*dat[floor(j/(2^k))]
      }
      dat[j] <- intercept+sum(term.matrix[j,])+error[j]
    }
    # copying missing values as in the original tree
    for (k in 1:n){
      dat[k] <- ifelse(is.na(z[k])==T, NA, dat[k])
    }
    boot.dat[h,] <- dat
  }
  out <- list()
  if(boot.est)
    out$boot.bcest=matrix(unlist(apply(boot.dat, 1, bfa.ls.bc, p=p, method=method)),nrow=B, byrow=TRUE)
  if(boot.data)
    out$boot.data <- boot.dat
  return(out)
}

##Asymptotic CI
bfa.asy.ci <- function(a1.ls, a1.ls.star, conf.level){
  alpha = 1 - conf.level
  se = sd(a1.ls.star)
  a1.L = a1.ls - qnorm(1-alpha/2)*se
  a1.U = a1.ls + qnorm(1-alpha/2)*se
  return(asy.ci = c(a1.L,a1.U))
}

##Percentile CI
bfa.perc.ci <- function(a1.ls.star, conf.level){
  alpha = 1 - conf.level
  return(perc.ci = quantile(a1.ls.star, prob=c(alpha/2, 1-alpha/2)))
}

##BC Percentile CI
bfa.perc.bc.ci <- function(z, a1.ls, a1.ls.star, conf.level) {
  # generate the bootstrap samples
  n <- length(z)
  B <- length(a1.ls.star)
  # compute z.hat
  z.hat <- qnorm(sum(a1.ls.star < a1.ls)/B)
  # pick out the BC quantiles
  alpha = 1 - conf.level
  p <- c(alpha/2, 1 - alpha/2)
  w <- 2*z.hat + qnorm(p)
  q <- pnorm(w)
  ind <- ceiling(q*B)
  a1.ls.star <- sort(a1.ls.star)
  return(percbc.ci = a1.ls.star[ind])
}

#########################
#BCa Function
# phi1 is the estimated phi
#phis is the vector of bootstrap replicates
##############################
Boot.BCa <- function(z,p,B=99,J=2,conf.level=0.95){
  n=length(z)
  nn=floor((n-1)/J)
  Rs=c(J:((2*J)-1))
  M=length(Rs)
  phi1=bfa.ls(z,1)$coef[1,2]
  dat.matrx = matrix(NA, nrow=J, ncol=nn)
  dat=list()
  # Bootstarp LS estimates  
  phis <- matrix(unlist(apply(matrix(unlist(bfa.boot1.ls(z,1,B=B, boot.est=FALSE,boot.data = TRUE)), ncol=n, nrow = B),1,p=p, bfa.ls)),ncol=3,byrow=TRUE)[,2]
  alpha <- (1 + c(-conf.level, conf.level))/2 
  zalpha <- qnorm(alpha)
  # the bias correction factor
  z0 <- qnorm(sum(as.vector(phis) < phi1) / B)
  # the acceleration factor (jackknife est.) 
  for(R in Rs){
    indx=NULL
    for (i in 0:10){
      ind=((2^i)*R):((2^i)*R+(2^i-1))
      indx=c(indx,if(tail(ind, n=1)<=n){ind})
    }
    dat=append(dat, list(indx))
  }
  for (i in 1:M){
    dat.matrx[i,] = z[dat[[i]]]
  }
  phi.jack <- matrix(unlist(apply(dat.matrx,1,p=p, bfa.ls)),ncol=3,byrow=TRUE)[,2]
  L <- mean(phi.jack) - phi.jack
  a <- sum(L^3)/(6 * sum(L^2)^1.5)
  # BCa conf. limits
  adj.alpha <- pnorm(z0 + (z0+zalpha)/(1-a*(z0+zalpha))) 
  limits <- quantile(phis, adj.alpha, names=F, type=8,na.rm = TRUE)
  return(list("est"=phi1, "BCa"=limits))
  
}
# Cowan & Staudte (1986) data
z<-c(23,17,18,24,27,22,24,27,26,29,30,31,30,26,31,24,29,27,25,28,38,18,23,31,32,27,27,24,24,30,33)

# Staudte et al. (1984) data
z=read.csv("File Path/Tree_means.csv", header = T)
z <- z$Mean.of.each.cell # convert the data to a vector


# Drawing the tree
bfa.tree.plot(z,vertex.size = 7,shape = "circle",text.size = 0.65,vertex.color="white",plot.margin = -0.3)

p=1 # p is the order of the bifurcating autoregressive model
B=999 # This will take too much time. For a quick run, B can be set to 9.
conf.level=0.95 # Confidence interval level 

# The following commands will give results for the Confidence Intervals disscussed in the article including the width of the CIs
ls.fit = bfa.ls(z, p, x.data = TRUE, y.data = TRUE, resids = FALSE, error.cor = TRUE, error.var = FALSE,
                cov.matrix = FALSE, conf = TRUE, conf.level = conf.level, p.value = FALSE)

a1.ls = ls.fit$coef[1,2]
a1.ls.boot1 =bfa.ls.bc(z,p,method="boot1", B1=99, B2=49)
a1.ls.boot2f = bfa.ls.bc(z,p,method="boot2fast", B1=99, B2=49)
a1.ls.star = bfa.boot1.ls(z, p, burn = 5, B, boot.est = TRUE, boot.data = FALSE)
a1.ls.boot1.star = bfa.boot.ls.bc(z, p, method="boot1",burn = 5, B, boot.est=TRUE, boot.data=FALSE)
a1.ls.boot2f.star = bfa.boot.ls.bc(z, p, method="boot2fast",burn = 5, B, boot.est=TRUE, boot.data=FALSE)

asy.fm.ci=ls.fit$conf
asy.fm.ci
asy.ci = bfa.asy.ci(a1.ls, a1.ls.star$boot.est[,2], conf.level)
asy.ci
asy.boot1.ci = bfa.asy.ci(a1.ls.boot1, a1.ls.boot1.star$boot.bcest[,1], conf.level)
asy.boot1.ci
asy.boot2f.ci = bfa.asy.ci(a1.ls.boot2f, a1.ls.boot2f.star$boot.bcest[,1], conf.level)
asy.boot2f.ci
perc.ci = bfa.perc.ci(a1.ls.star$boot.est[,2], conf.level)
perc.ci
perc.boot1.ci = bfa.perc.ci(a1.ls.boot1.star$boot.bcest[,1], conf.level)
perc.boot1.ci
perc.boot2f.ci = bfa.perc.ci(a1.ls.boot2f.star$boot.bcest[,1], conf.level)
perc.boot2f.ci
percbc.ci = bfa.perc.bc.ci(z, a1.ls, a1.ls.star$boot.est[,2], conf.level)
percbc.ci
bca.ci=Boot.BCa(z,p,B,J=4,conf.level)
bca.ci
asy.fm.w=asy.fm.ci[2,2]-asy.fm.ci[2,1]
asy.fm.w
asy.w = asy.ci[2] - asy.ci[1]
asy.w
asy.boot1.w = asy.boot1.ci[2] - asy.boot1.ci[1]
asy.boot1.w
asy.boot2f.w = asy.boot2f.ci[2] - asy.boot2f.ci[1]
asy.boot2f.w
perc.w = perc.ci[2] - perc.ci[1]
perc.w
perc.boot1.w = perc.boot1.ci[2] - perc.boot1.ci[1]
perc.boot1.w
perc.boot2f.w = perc.boot2f.ci[2] - perc.boot2f.ci[1]
perc.boot2f.w
percbc.w = percbc.ci[2] - percbc.ci[1]
percbc.w
bca.w=bca.ci$BCa[2]-bca.ci$BCa[1]
bca.w






